Actual source code: dmlabel.c
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/sectionimpl.h>
4: #include <petscsf.h>
5: #include <petscsection.h>
7: PetscFunctionList DMLabelList = NULL;
8: PetscBool DMLabelRegisterAllCalled = PETSC_FALSE;
10: /*@
11: DMLabelCreate - Create a `DMLabel` object, which is a multimap
13: Collective
15: Input Parameters:
16: + comm - The communicator, usually `PETSC_COMM_SELF`
17: - name - The label name
19: Output Parameter:
20: . label - The `DMLabel`
22: Level: beginner
24: Notes:
25: The label name is actually usually the `PetscObject` name.
26: One can get/set it with `PetscObjectGetName()`/`PetscObjectSetName()`.
28: .seealso: `DMLabel`, `DM`, `DMLabelDestroy()`
29: @*/
30: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
31: {
32: PetscFunctionBegin;
33: PetscAssertPointer(label, 3);
34: PetscCall(DMInitializePackage());
36: PetscCall(PetscHeaderCreate(*label, DMLABEL_CLASSID, "DMLabel", "DMLabel", "DM", comm, DMLabelDestroy, DMLabelView));
37: (*label)->numStrata = 0;
38: (*label)->defaultValue = -1;
39: (*label)->stratumValues = NULL;
40: (*label)->validIS = NULL;
41: (*label)->stratumSizes = NULL;
42: (*label)->points = NULL;
43: (*label)->ht = NULL;
44: (*label)->pStart = -1;
45: (*label)->pEnd = -1;
46: (*label)->bt = NULL;
47: PetscCall(PetscHMapICreate(&(*label)->hmap));
48: PetscCall(PetscObjectSetName((PetscObject)*label, name));
49: PetscCall(DMLabelSetType(*label, DMLABELCONCRETE));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*@
54: DMLabelSetUp - SetUp a `DMLabel` object
56: Collective
58: Input Parameters:
59: . label - The `DMLabel`
61: Level: intermediate
63: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
64: @*/
65: PetscErrorCode DMLabelSetUp(DMLabel label)
66: {
67: PetscFunctionBegin;
69: PetscTryTypeMethod(label, setup);
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*
74: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
76: Not collective
78: Input parameter:
79: + label - The `DMLabel`
80: - v - The stratum value
82: Output parameter:
83: . label - The `DMLabel` with stratum in sorted list format
85: Level: developer
87: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
88: */
89: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
90: {
91: IS is;
92: PetscInt off = 0, *pointArray, p;
94: PetscFunctionBegin;
95: if ((PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
96: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeValid_Private", v);
97: PetscCall(PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]));
98: PetscCall(PetscMalloc1(label->stratumSizes[v], &pointArray));
99: PetscCall(PetscHSetIGetElems(label->ht[v], &off, pointArray));
100: PetscCall(PetscHSetIClear(label->ht[v]));
101: PetscCall(PetscSortInt(label->stratumSizes[v], pointArray));
102: if (label->bt) {
103: for (p = 0; p < label->stratumSizes[v]; ++p) {
104: const PetscInt point = pointArray[p];
105: PetscCheck(!(point < label->pStart) && !(point >= label->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %" PetscInt_FMT " is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, label->pStart, label->pEnd);
106: PetscCall(PetscBTSet(label->bt, point - label->pStart));
107: }
108: }
109: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
110: PetscCall(ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is));
111: PetscCall(PetscFree(pointArray));
112: } else {
113: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is));
114: }
115: PetscCall(ISSetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, PETSC_TRUE));
116: PetscCall(PetscObjectSetName((PetscObject)is, "indices"));
117: label->points[v] = is;
118: label->validIS[v] = PETSC_TRUE;
119: PetscCall(PetscObjectStateIncrease((PetscObject)label));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*
124: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
126: Not Collective
128: Input parameter:
129: . label - The `DMLabel`
131: Output parameter:
132: . label - The `DMLabel` with all strata in sorted list format
134: Level: developer
136: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
137: */
138: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
139: {
140: PetscInt v;
142: PetscFunctionBegin;
143: for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeValid_Private(label, v));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*
148: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
150: Not Collective
152: Input parameter:
153: + label - The `DMLabel`
154: - v - The stratum value
156: Output parameter:
157: . label - The `DMLabel` with stratum in hash format
159: Level: developer
161: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`
162: */
163: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
164: {
165: PetscInt p;
166: const PetscInt *points;
168: PetscFunctionBegin;
169: if ((PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) || label->readonly) PetscFunctionReturn(PETSC_SUCCESS);
170: PetscCheck(v >= 0 && v < label->numStrata, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %" PetscInt_FMT " in DMLabelMakeInvalid_Private", v);
171: if (label->points[v]) {
172: PetscCall(ISGetIndices(label->points[v], &points));
173: for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
174: PetscCall(ISRestoreIndices(label->points[v], &points));
175: PetscCall(ISDestroy(&label->points[v]));
176: }
177: label->validIS[v] = PETSC_FALSE;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label)
182: {
183: PetscInt v;
185: PetscFunctionBegin;
186: for (v = 0; v < label->numStrata; v++) PetscCall(DMLabelMakeInvalid_Private(label, v));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
191: #define DMLABEL_LOOKUP_THRESHOLD 16
192: #endif
194: PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
195: {
196: PetscInt v;
198: PetscFunctionBegin;
199: *index = -1;
200: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD || label->readonly) {
201: for (v = 0; v < label->numStrata; ++v)
202: if (label->stratumValues[v] == value) {
203: *index = v;
204: break;
205: }
206: } else {
207: PetscCall(PetscHMapIGet(label->hmap, value, index));
208: }
209: if (PetscDefined(USE_DEBUG) && !label->readonly) { /* Check strata hash map consistency */
210: PetscInt len, loc = -1;
211: PetscCall(PetscHMapIGetSize(label->hmap, &len));
212: PetscCheck(len == label->numStrata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
213: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
214: PetscCall(PetscHMapIGet(label->hmap, value, &loc));
215: } else {
216: for (v = 0; v < label->numStrata; ++v)
217: if (label->stratumValues[v] == value) {
218: loc = v;
219: break;
220: }
221: }
222: PetscCheck(loc == *index, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
223: }
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
228: {
229: PetscInt v;
230: PetscInt *tmpV;
231: PetscInt *tmpS;
232: PetscHSetI *tmpH, ht;
233: IS *tmpP, is;
234: PetscBool *tmpB;
235: PetscHMapI hmap = label->hmap;
237: PetscFunctionBegin;
238: v = label->numStrata;
239: tmpV = label->stratumValues;
240: tmpS = label->stratumSizes;
241: tmpH = label->ht;
242: tmpP = label->points;
243: tmpB = label->validIS;
244: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
245: PetscInt *oldV = tmpV;
246: PetscInt *oldS = tmpS;
247: PetscHSetI *oldH = tmpH;
248: IS *oldP = tmpP;
249: PetscBool *oldB = tmpB;
250: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpV), &tmpV));
251: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpS), &tmpS));
252: PetscCall(PetscCalloc((v + 1) * sizeof(*tmpH), &tmpH));
253: PetscCall(PetscCalloc((v + 1) * sizeof(*tmpP), &tmpP));
254: PetscCall(PetscMalloc((v + 1) * sizeof(*tmpB), &tmpB));
255: PetscCall(PetscArraycpy(tmpV, oldV, v));
256: PetscCall(PetscArraycpy(tmpS, oldS, v));
257: PetscCall(PetscArraycpy(tmpH, oldH, v));
258: PetscCall(PetscArraycpy(tmpP, oldP, v));
259: PetscCall(PetscArraycpy(tmpB, oldB, v));
260: PetscCall(PetscFree(oldV));
261: PetscCall(PetscFree(oldS));
262: PetscCall(PetscFree(oldH));
263: PetscCall(PetscFree(oldP));
264: PetscCall(PetscFree(oldB));
265: }
266: label->numStrata = v + 1;
267: label->stratumValues = tmpV;
268: label->stratumSizes = tmpS;
269: label->ht = tmpH;
270: label->points = tmpP;
271: label->validIS = tmpB;
272: PetscCall(PetscHSetICreate(&ht));
273: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
274: PetscCall(PetscHMapISet(hmap, value, v));
275: tmpV[v] = value;
276: tmpS[v] = 0;
277: tmpH[v] = ht;
278: tmpP[v] = is;
279: tmpB[v] = PETSC_TRUE;
280: PetscCall(PetscObjectStateIncrease((PetscObject)label));
281: *index = v;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
286: {
287: PetscFunctionBegin;
288: PetscCall(DMLabelLookupStratum(label, value, index));
289: if (*index < 0) PetscCall(DMLabelNewStratum(label, value, index));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
294: {
295: PetscFunctionBegin;
296: *size = 0;
297: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
298: if (label->readonly || label->validIS[v]) {
299: *size = label->stratumSizes[v];
300: } else {
301: PetscCall(PetscHSetIGetSize(label->ht[v], size));
302: }
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@
307: DMLabelAddStratum - Adds a new stratum value in a `DMLabel`
309: Input Parameters:
310: + label - The `DMLabel`
311: - value - The stratum value
313: Level: beginner
315: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
316: @*/
317: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
318: {
319: PetscInt v;
321: PetscFunctionBegin;
323: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
324: PetscCall(DMLabelLookupAddStratum(label, value, &v));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: DMLabelAddStrata - Adds new stratum values in a `DMLabel`
331: Not Collective
333: Input Parameters:
334: + label - The `DMLabel`
335: . numStrata - The number of stratum values
336: - stratumValues - The stratum values
338: Level: beginner
340: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
341: @*/
342: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
343: {
344: PetscInt *values, v;
346: PetscFunctionBegin;
348: if (numStrata) PetscAssertPointer(stratumValues, 3);
349: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
350: PetscCall(PetscMalloc1(numStrata, &values));
351: PetscCall(PetscArraycpy(values, stratumValues, numStrata));
352: PetscCall(PetscSortRemoveDupsInt(&numStrata, values));
353: if (!label->numStrata) { /* Fast preallocation */
354: PetscInt *tmpV;
355: PetscInt *tmpS;
356: PetscHSetI *tmpH, ht;
357: IS *tmpP, is;
358: PetscBool *tmpB;
359: PetscHMapI hmap = label->hmap;
361: PetscCall(PetscMalloc1(numStrata, &tmpV));
362: PetscCall(PetscMalloc1(numStrata, &tmpS));
363: PetscCall(PetscCalloc1(numStrata, &tmpH));
364: PetscCall(PetscCalloc1(numStrata, &tmpP));
365: PetscCall(PetscMalloc1(numStrata, &tmpB));
366: label->numStrata = numStrata;
367: label->stratumValues = tmpV;
368: label->stratumSizes = tmpS;
369: label->ht = tmpH;
370: label->points = tmpP;
371: label->validIS = tmpB;
372: for (v = 0; v < numStrata; ++v) {
373: PetscCall(PetscHSetICreate(&ht));
374: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is));
375: PetscCall(PetscHMapISet(hmap, values[v], v));
376: tmpV[v] = values[v];
377: tmpS[v] = 0;
378: tmpH[v] = ht;
379: tmpP[v] = is;
380: tmpB[v] = PETSC_TRUE;
381: }
382: PetscCall(PetscObjectStateIncrease((PetscObject)label));
383: } else {
384: for (v = 0; v < numStrata; ++v) PetscCall(DMLabelAddStratum(label, values[v]));
385: }
386: PetscCall(PetscFree(values));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@
391: DMLabelAddStrataIS - Adds new stratum values in a `DMLabel`
393: Not Collective
395: Input Parameters:
396: + label - The `DMLabel`
397: - valueIS - Index set with stratum values
399: Level: beginner
401: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
402: @*/
403: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
404: {
405: PetscInt numStrata;
406: const PetscInt *stratumValues;
408: PetscFunctionBegin;
411: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
412: PetscCall(ISGetLocalSize(valueIS, &numStrata));
413: PetscCall(ISGetIndices(valueIS, &stratumValues));
414: PetscCall(DMLabelAddStrata(label, numStrata, stratumValues));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: static PetscErrorCode DMLabelView_Concrete_Ascii(DMLabel label, PetscViewer viewer)
419: {
420: PetscInt v;
421: PetscMPIInt rank;
423: PetscFunctionBegin;
424: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
425: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
426: if (label) {
427: const char *name;
429: PetscCall(PetscObjectGetName((PetscObject)label, &name));
430: PetscCall(PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name));
431: if (label->bt) PetscCall(PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", label->pStart, label->pEnd));
432: for (v = 0; v < label->numStrata; ++v) {
433: const PetscInt value = label->stratumValues[v];
434: const PetscInt *points;
435: PetscInt p;
437: PetscCall(ISGetIndices(label->points[v], &points));
438: for (p = 0; p < label->stratumSizes[v]; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value));
439: PetscCall(ISRestoreIndices(label->points[v], &points));
440: }
441: }
442: PetscCall(PetscViewerFlush(viewer));
443: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode DMLabelView_Concrete(DMLabel label, PetscViewer viewer)
448: {
449: PetscBool isascii;
451: PetscFunctionBegin;
452: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
453: if (isascii) PetscCall(DMLabelView_Concrete_Ascii(label, viewer));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: DMLabelView - View the label
460: Collective
462: Input Parameters:
463: + label - The `DMLabel`
464: - viewer - The `PetscViewer`
466: Level: intermediate
468: .seealso: `DMLabel`, `PetscViewer`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
469: @*/
470: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
471: {
472: PetscFunctionBegin;
474: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer));
476: PetscCall(DMLabelMakeAllValid_Private(label));
477: PetscUseTypeMethod(label, view, viewer);
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: 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: PetscCall(PetscHMapIReset(label->hmap));
510: label->pStart = -1;
511: label->pEnd = -1;
512: PetscCall(PetscBTDestroy(&label->bt));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: /*@
517: DMLabelDestroy - Destroys a `DMLabel`
519: Collective
521: Input Parameter:
522: . label - The `DMLabel`
524: Level: beginner
526: .seealso: `DMLabel`, `DM`, `DMLabelReset()`, `DMLabelCreate()`
527: @*/
528: PetscErrorCode DMLabelDestroy(DMLabel *label)
529: {
530: PetscFunctionBegin;
531: if (!*label) PetscFunctionReturn(PETSC_SUCCESS);
533: if (--((PetscObject)*label)->refct > 0) {
534: *label = NULL;
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
537: PetscCall(DMLabelReset(*label));
538: PetscCall(PetscHMapIDestroy(&(*label)->hmap));
539: PetscCall(PetscHeaderDestroy(label));
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: static PetscErrorCode DMLabelDuplicate_Concrete(DMLabel label, DMLabel *labelnew)
544: {
545: PetscFunctionBegin;
546: for (PetscInt v = 0; v < label->numStrata; ++v) {
547: PetscCall(PetscHSetICreate(&(*labelnew)->ht[v]));
548: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
549: (*labelnew)->points[v] = label->points[v];
550: }
551: PetscCall(PetscHMapIDestroy(&(*labelnew)->hmap));
552: PetscCall(PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: /*@
557: DMLabelDuplicate - Duplicates a `DMLabel`
559: Collective
561: Input Parameter:
562: . label - The `DMLabel`
564: Output Parameter:
565: . labelnew - new label
567: Level: intermediate
569: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelDestroy()`
570: @*/
571: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
572: {
573: const char *name;
575: PetscFunctionBegin;
577: PetscCall(DMLabelMakeAllValid_Private(label));
578: PetscCall(PetscObjectGetName((PetscObject)label, &name));
579: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew));
581: (*labelnew)->numStrata = label->numStrata;
582: (*labelnew)->defaultValue = label->defaultValue;
583: (*labelnew)->readonly = label->readonly;
584: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues));
585: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes));
586: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->ht));
587: PetscCall(PetscCalloc1(label->numStrata, &(*labelnew)->points));
588: PetscCall(PetscMalloc1(label->numStrata, &(*labelnew)->validIS));
589: for (PetscInt v = 0; v < label->numStrata; ++v) {
590: (*labelnew)->stratumValues[v] = label->stratumValues[v];
591: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
592: (*labelnew)->validIS[v] = PETSC_TRUE;
593: }
594: (*labelnew)->pStart = -1;
595: (*labelnew)->pEnd = -1;
596: (*labelnew)->bt = NULL;
597: PetscUseTypeMethod(label, duplicate, labelnew);
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*@C
602: DMLabelCompare - Compare two `DMLabel` objects
604: Collective; No Fortran Support
606: Input Parameters:
607: + comm - Comm over which to compare labels
608: . l0 - First `DMLabel`
609: - l1 - Second `DMLabel`
611: Output Parameters:
612: + equal - (Optional) Flag whether the two labels are equal
613: - message - (Optional) Message describing the difference
615: Level: intermediate
617: Notes:
618: The output flag equal is the same on all processes.
619: If it is passed as `NULL` and difference is found, an error is thrown on all processes.
620: Make sure to pass `NULL` on all processes.
622: The output message is set independently on each rank.
623: It is set to `NULL` if no difference was found on the current rank. It must be freed by user.
624: If message is passed as `NULL` and difference is found, the difference description is printed to stderr in synchronized manner.
625: Make sure to pass `NULL` on all processes.
627: For the comparison, we ignore the order of stratum values, and strata with no points.
629: The communicator needs to be specified because currently `DMLabel` can live on `PETSC_COMM_SELF` even if the underlying `DM` is parallel.
631: Developer Note:
632: Fortran stub cannot be generated automatically because `message` must be freed with `PetscFree()`
634: .seealso: `DMLabel`, `DM`, `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
635: @*/
636: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
637: {
638: const char *name0, *name1;
639: char msg[PETSC_MAX_PATH_LEN] = "";
640: PetscBool eq;
641: PetscMPIInt rank;
643: PetscFunctionBegin;
646: if (equal) PetscAssertPointer(equal, 4);
647: if (message) PetscAssertPointer(message, 5);
648: PetscCallMPI(MPI_Comm_rank(comm, &rank));
649: PetscCall(PetscObjectGetName((PetscObject)l0, &name0));
650: PetscCall(PetscObjectGetName((PetscObject)l1, &name1));
651: {
652: PetscInt v0, v1;
654: PetscCall(DMLabelGetDefaultValue(l0, &v0));
655: PetscCall(DMLabelGetDefaultValue(l1, &v1));
656: eq = (PetscBool)(v0 == v1);
657: 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));
658: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
659: if (!eq) goto finish;
660: }
661: {
662: IS is0, is1;
664: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l0, &is0));
665: PetscCall(DMLabelGetNonEmptyStratumValuesIS(l1, &is1));
666: PetscCall(ISEqual(is0, is1, &eq));
667: PetscCall(ISDestroy(&is0));
668: PetscCall(ISDestroy(&is1));
669: if (!eq) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1));
670: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
671: if (!eq) goto finish;
672: }
673: {
674: PetscInt i, nValues;
676: PetscCall(DMLabelGetNumValues(l0, &nValues));
677: for (i = 0; i < nValues; i++) {
678: const PetscInt v = l0->stratumValues[i];
679: PetscInt n;
680: IS is0, is1;
682: PetscCall(DMLabelGetStratumSize_Private(l0, i, &n));
683: if (!n) continue;
684: PetscCall(DMLabelGetStratumIS(l0, v, &is0));
685: PetscCall(DMLabelGetStratumIS(l1, v, &is1));
686: PetscCall(ISEqualUnsorted(is0, is1, &eq));
687: PetscCall(ISDestroy(&is0));
688: PetscCall(ISDestroy(&is1));
689: if (!eq) {
690: 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));
691: break;
692: }
693: }
694: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm));
695: }
696: finish:
697: /* If message output arg not set, print to stderr */
698: if (message) {
699: *message = NULL;
700: if (msg[0]) PetscCall(PetscStrallocpy(msg, message));
701: } else {
702: if (msg[0]) PetscCall(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg));
703: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDERR));
704: }
705: /* If same output arg not ser and labels are not equal, throw error */
706: if (equal) *equal = eq;
707: else PetscCheck(eq, comm, PETSC_ERR_ARG_INCOMP, "DMLabels l0 \"%s\" and l1 \"%s\" are not equal", name0, name1);
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: /*@
712: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
714: Not Collective
716: Input Parameter:
717: . label - The `DMLabel`
719: Level: intermediate
721: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
722: @*/
723: PetscErrorCode DMLabelComputeIndex(DMLabel label)
724: {
725: PetscInt pStart = PETSC_INT_MAX, pEnd = -1, v;
727: PetscFunctionBegin;
729: PetscCall(DMLabelMakeAllValid_Private(label));
730: for (v = 0; v < label->numStrata; ++v) {
731: const PetscInt *points;
732: PetscInt i;
734: PetscCall(ISGetIndices(label->points[v], &points));
735: for (i = 0; i < label->stratumSizes[v]; ++i) {
736: const PetscInt point = points[i];
738: pStart = PetscMin(point, pStart);
739: pEnd = PetscMax(point + 1, pEnd);
740: }
741: PetscCall(ISRestoreIndices(label->points[v], &points));
742: }
743: label->pStart = pStart == PETSC_INT_MAX ? -1 : pStart;
744: label->pEnd = pEnd;
745: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /*@
750: DMLabelCreateIndex - Create an index structure for membership determination
752: Not Collective
754: Input Parameters:
755: + label - The `DMLabel`
756: . pStart - The smallest point
757: - pEnd - The largest point + 1
759: Level: intermediate
761: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
762: @*/
763: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
764: {
765: PetscInt v;
767: PetscFunctionBegin;
769: PetscCall(DMLabelDestroyIndex(label));
770: PetscCall(DMLabelMakeAllValid_Private(label));
771: label->pStart = pStart;
772: label->pEnd = pEnd;
773: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
774: PetscCall(PetscBTCreate(pEnd - pStart, &label->bt));
775: for (v = 0; v < label->numStrata; ++v) {
776: IS pointIS;
777: const PetscInt *points;
778: PetscInt i;
780: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
781: PetscCall(ISGetIndices(pointIS, &points));
782: for (i = 0; i < label->stratumSizes[v]; ++i) {
783: const PetscInt point = points[i];
785: 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);
786: PetscCall(PetscBTSet(label->bt, point - pStart));
787: }
788: PetscCall(ISRestoreIndices(label->points[v], &points));
789: PetscCall(ISDestroy(&pointIS));
790: }
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@
795: DMLabelDestroyIndex - Destroy the index structure
797: Not Collective
799: Input Parameter:
800: . label - the `DMLabel`
802: Level: intermediate
804: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
805: @*/
806: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
807: {
808: PetscFunctionBegin;
810: label->pStart = -1;
811: label->pEnd = -1;
812: PetscCall(PetscBTDestroy(&label->bt));
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: /*@
817: DMLabelGetBounds - Return the smallest and largest point in the label
819: Not Collective
821: Input Parameter:
822: . label - the `DMLabel`
824: Output Parameters:
825: + pStart - The smallest point
826: - pEnd - The largest point + 1
828: Level: intermediate
830: Note:
831: This will compute an index for the label if one does not exist.
833: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
834: @*/
835: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
836: {
837: PetscFunctionBegin;
839: if ((label->pStart == -1) && (label->pEnd == -1)) PetscCall(DMLabelComputeIndex(label));
840: if (pStart) {
841: PetscAssertPointer(pStart, 2);
842: *pStart = label->pStart;
843: }
844: if (pEnd) {
845: PetscAssertPointer(pEnd, 3);
846: *pEnd = label->pEnd;
847: }
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: /*@
852: DMLabelHasValue - Determine whether a label assigns the value to any point
854: Not Collective
856: Input Parameters:
857: + label - the `DMLabel`
858: - value - the value
860: Output Parameter:
861: . contains - Flag indicating whether the label maps this value to any point
863: Level: developer
865: .seealso: `DMLabel`, `DM`, `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
866: @*/
867: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
868: {
869: PetscInt v;
871: PetscFunctionBegin;
873: PetscAssertPointer(contains, 3);
874: PetscCall(DMLabelLookupStratum(label, value, &v));
875: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: /*@
880: DMLabelHasPoint - Determine whether a label assigns a value to a point
882: Not Collective
884: Input Parameters:
885: + label - the `DMLabel`
886: - point - the point
888: Output Parameter:
889: . contains - Flag indicating whether the label maps this point to a value
891: Level: developer
893: Note:
894: The user must call `DMLabelCreateIndex()` before this function.
896: .seealso: `DMLabel`, `DM`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
897: @*/
898: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
899: {
900: PetscFunctionBeginHot;
902: PetscAssertPointer(contains, 3);
903: PetscCall(DMLabelMakeAllValid_Private(label));
904: if (PetscDefined(USE_DEBUG)) {
905: PetscCheck(label->bt, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
906: 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);
907: }
908: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: /*@
913: DMLabelStratumHasPoint - Return true if the stratum contains a point
915: Not Collective
917: Input Parameters:
918: + label - the `DMLabel`
919: . value - the stratum value
920: - point - the point
922: Output Parameter:
923: . contains - true if the stratum contains the point
925: Level: intermediate
927: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
928: @*/
929: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
930: {
931: PetscFunctionBeginHot;
933: PetscAssertPointer(contains, 4);
934: if (value == label->defaultValue) {
935: PetscInt pointVal;
937: PetscCall(DMLabelGetValue(label, point, &pointVal));
938: *contains = (PetscBool)(pointVal == value);
939: } else {
940: PetscInt v;
942: PetscCall(DMLabelLookupStratum(label, value, &v));
943: if (v >= 0) {
944: if (label->validIS[v] || label->readonly) {
945: IS is;
946: PetscInt i;
948: PetscUseTypeMethod(label, getstratumis, v, &is);
949: PetscCall(ISLocate(is, point, &i));
950: PetscCall(ISDestroy(&is));
951: *contains = (PetscBool)(i >= 0);
952: } else {
953: PetscCall(PetscHSetIHas(label->ht[v], point, contains));
954: }
955: } else { // value is not present
956: *contains = PETSC_FALSE;
957: }
958: }
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: /*@
963: DMLabelGetDefaultValue - Get the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
964: When a label is created, it is initialized to -1.
966: Not Collective
968: Input Parameter:
969: . label - a `DMLabel` object
971: Output Parameter:
972: . defaultValue - the default value
974: Level: beginner
976: .seealso: `DMLabel`, `DM`, `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
977: @*/
978: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
979: {
980: PetscFunctionBegin;
982: *defaultValue = label->defaultValue;
983: PetscFunctionReturn(PETSC_SUCCESS);
984: }
986: /*@
987: DMLabelSetDefaultValue - Set the default value returned by `DMLabelGetValue()` if a point has not been explicitly given a value.
988: When a label is created, it is initialized to -1.
990: Not Collective
992: Input Parameter:
993: . label - a `DMLabel` object
995: Output Parameter:
996: . defaultValue - the default value
998: Level: beginner
1000: .seealso: `DMLabel`, `DM`, `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1001: @*/
1002: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
1003: {
1004: PetscFunctionBegin;
1006: label->defaultValue = defaultValue;
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }
1010: /*@
1011: 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
1012: `DMLabelSetDefaultValue()`)
1014: Not Collective
1016: Input Parameters:
1017: + label - the `DMLabel`
1018: - point - the point
1020: Output Parameter:
1021: . value - The point value, or the default value (-1 by default)
1023: Level: intermediate
1025: Note:
1026: A label may assign multiple values to a point. No guarantees are made about which value is returned in that case.
1027: Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
1029: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1030: @*/
1031: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
1032: {
1033: PetscInt v;
1035: PetscFunctionBeginHot;
1037: PetscAssertPointer(value, 3);
1038: *value = label->defaultValue;
1039: for (v = 0; v < label->numStrata; ++v) {
1040: if (label->validIS[v] || label->readonly) {
1041: IS is;
1042: PetscInt i;
1044: PetscUseTypeMethod(label, getstratumis, v, &is);
1045: PetscCall(ISLocate(label->points[v], point, &i));
1046: PetscCall(ISDestroy(&is));
1047: if (i >= 0) {
1048: *value = label->stratumValues[v];
1049: break;
1050: }
1051: } else {
1052: PetscBool has;
1054: PetscCall(PetscHSetIHas(label->ht[v], point, &has));
1055: if (has) {
1056: *value = label->stratumValues[v];
1057: break;
1058: }
1059: }
1060: }
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: /*@
1065: 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
1066: be changed with `DMLabelSetDefaultValue()` to something different), then this function will do nothing.
1068: Not Collective
1070: Input Parameters:
1071: + label - the `DMLabel`
1072: . point - the point
1073: - value - The point value
1075: Level: intermediate
1077: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1078: @*/
1079: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1080: {
1081: PetscInt v;
1083: PetscFunctionBegin;
1085: /* Find label value, add new entry if needed */
1086: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1087: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1088: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1089: /* Set key */
1090: PetscCall(DMLabelMakeInvalid_Private(label, v));
1091: PetscCall(PetscHSetIAdd(label->ht[v], point));
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /*@
1096: DMLabelClearValue - Clear the value a label assigns to a point
1098: Not Collective
1100: Input Parameters:
1101: + label - the `DMLabel`
1102: . point - the point
1103: - value - The point value
1105: Level: intermediate
1107: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1108: @*/
1109: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1110: {
1111: PetscInt v;
1113: PetscFunctionBegin;
1115: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1116: /* Find label value */
1117: PetscCall(DMLabelLookupStratum(label, value, &v));
1118: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1120: if (label->bt) {
1121: 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);
1122: PetscCall(PetscBTClear(label->bt, point - label->pStart));
1123: }
1125: /* Delete key */
1126: PetscCall(DMLabelMakeInvalid_Private(label, v));
1127: PetscCall(PetscHSetIDel(label->ht[v], point));
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: /*@
1132: DMLabelInsertIS - Set all points in the `IS` to a value
1134: Not Collective
1136: Input Parameters:
1137: + label - the `DMLabel`
1138: . is - the point `IS`
1139: - value - The point value
1141: Level: intermediate
1143: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1144: @*/
1145: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1146: {
1147: PetscInt v, n, p;
1148: const PetscInt *points;
1150: PetscFunctionBegin;
1153: /* Find label value, add new entry if needed */
1154: if (value == label->defaultValue) PetscFunctionReturn(PETSC_SUCCESS);
1155: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1156: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1157: /* Set keys */
1158: PetscCall(DMLabelMakeInvalid_Private(label, v));
1159: PetscCall(ISGetLocalSize(is, &n));
1160: PetscCall(ISGetIndices(is, &points));
1161: for (p = 0; p < n; ++p) PetscCall(PetscHSetIAdd(label->ht[v], points[p]));
1162: PetscCall(ISRestoreIndices(is, &points));
1163: PetscFunctionReturn(PETSC_SUCCESS);
1164: }
1166: /*@
1167: DMLabelGetNumValues - Get the number of values that the `DMLabel` takes
1169: Not Collective
1171: Input Parameter:
1172: . label - the `DMLabel`
1174: Output Parameter:
1175: . numValues - the number of values
1177: Level: intermediate
1179: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1180: @*/
1181: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1182: {
1183: PetscFunctionBegin;
1185: PetscAssertPointer(numValues, 2);
1186: *numValues = label->numStrata;
1187: PetscFunctionReturn(PETSC_SUCCESS);
1188: }
1190: /*@
1191: DMLabelGetValueIS - Get an `IS` of all values that the `DMlabel` takes
1193: Not Collective
1195: Input Parameter:
1196: . label - the `DMLabel`
1198: Output Parameter:
1199: . values - the value `IS`
1201: Level: intermediate
1203: Notes:
1204: The output `IS` should be destroyed when no longer needed.
1205: Strata which are allocated but empty [`DMLabelGetStratumSize()` yields 0] are counted.
1206: If you need to count only nonempty strata, use `DMLabelGetNonEmptyStratumValuesIS()`.
1208: .seealso: `DMLabel`, `DM`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1209: @*/
1210: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1211: {
1212: PetscFunctionBegin;
1214: PetscAssertPointer(values, 2);
1215: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1216: PetscFunctionReturn(PETSC_SUCCESS);
1217: }
1219: /*@
1220: DMLabelGetValueBounds - Return the smallest and largest value in the label
1222: Not Collective
1224: Input Parameter:
1225: . label - the `DMLabel`
1227: Output Parameters:
1228: + minValue - The smallest value
1229: - maxValue - The largest value
1231: Level: intermediate
1233: .seealso: `DMLabel`, `DM`, `DMLabelGetBounds()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1234: @*/
1235: PetscErrorCode DMLabelGetValueBounds(DMLabel label, PetscInt *minValue, PetscInt *maxValue)
1236: {
1237: PetscInt min = PETSC_INT_MAX, max = PETSC_INT_MIN;
1239: PetscFunctionBegin;
1241: for (PetscInt v = 0; v < label->numStrata; ++v) {
1242: min = PetscMin(min, label->stratumValues[v]);
1243: max = PetscMax(max, label->stratumValues[v]);
1244: }
1245: if (minValue) {
1246: PetscAssertPointer(minValue, 2);
1247: *minValue = min;
1248: }
1249: if (maxValue) {
1250: PetscAssertPointer(maxValue, 3);
1251: *maxValue = max;
1252: }
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@
1257: DMLabelGetNonEmptyStratumValuesIS - Get an `IS` of all values that the `DMlabel` takes
1259: Not Collective
1261: Input Parameter:
1262: . label - the `DMLabel`
1264: Output Parameter:
1265: . values - the value `IS`
1267: Level: intermediate
1269: Notes:
1270: The output `IS` should be destroyed when no longer needed.
1271: This is similar to `DMLabelGetValueIS()` but counts only nonempty strata.
1273: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1274: @*/
1275: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1276: {
1277: PetscInt i, j;
1278: PetscInt *valuesArr;
1280: PetscFunctionBegin;
1282: PetscAssertPointer(values, 2);
1283: PetscCall(PetscMalloc1(label->numStrata, &valuesArr));
1284: for (i = 0, j = 0; i < label->numStrata; i++) {
1285: PetscInt n;
1287: PetscCall(DMLabelGetStratumSize_Private(label, i, &n));
1288: if (n) valuesArr[j++] = label->stratumValues[i];
1289: }
1290: if (j == label->numStrata) {
1291: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values));
1292: } else {
1293: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values));
1294: }
1295: PetscCall(PetscFree(valuesArr));
1296: PetscFunctionReturn(PETSC_SUCCESS);
1297: }
1299: /*@
1300: DMLabelGetValueIndex - Get the index of a given value in the list of values for the `DMlabel`, or -1 if it is not present
1302: Not Collective
1304: Input Parameters:
1305: + label - the `DMLabel`
1306: - value - the value
1308: Output Parameter:
1309: . index - the index of value in the list of values
1311: Level: intermediate
1313: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1314: @*/
1315: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1316: {
1317: PetscInt v;
1319: PetscFunctionBegin;
1321: PetscAssertPointer(index, 3);
1322: /* Do not assume they are sorted */
1323: for (v = 0; v < label->numStrata; ++v)
1324: if (label->stratumValues[v] == value) break;
1325: if (v >= label->numStrata) *index = -1;
1326: else *index = v;
1327: PetscFunctionReturn(PETSC_SUCCESS);
1328: }
1330: /*@
1331: DMLabelHasStratum - Determine whether points exist with the given value
1333: Not Collective
1335: Input Parameters:
1336: + label - the `DMLabel`
1337: - value - the stratum value
1339: Output Parameter:
1340: . exists - Flag saying whether points exist
1342: Level: intermediate
1344: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1345: @*/
1346: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1347: {
1348: PetscInt v;
1350: PetscFunctionBegin;
1352: PetscAssertPointer(exists, 3);
1353: PetscCall(DMLabelLookupStratum(label, value, &v));
1354: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: /*@
1359: DMLabelGetStratumSize - Get the size of a stratum
1361: Not Collective
1363: Input Parameters:
1364: + label - the `DMLabel`
1365: - value - the stratum value
1367: Output Parameter:
1368: . size - The number of points in the stratum
1370: Level: intermediate
1372: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1373: @*/
1374: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1375: {
1376: PetscInt v;
1378: PetscFunctionBegin;
1380: PetscAssertPointer(size, 3);
1381: PetscCall(DMLabelLookupStratum(label, value, &v));
1382: PetscCall(DMLabelGetStratumSize_Private(label, v, size));
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: /*@
1387: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1389: Not Collective
1391: Input Parameters:
1392: + label - the `DMLabel`
1393: - value - the stratum value
1395: Output Parameters:
1396: + start - the smallest point in the stratum
1397: - end - the largest point in the stratum
1399: Level: intermediate
1401: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1402: @*/
1403: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1404: {
1405: IS is;
1406: PetscInt v, min, max;
1408: PetscFunctionBegin;
1410: if (start) {
1411: PetscAssertPointer(start, 3);
1412: *start = -1;
1413: }
1414: if (end) {
1415: PetscAssertPointer(end, 4);
1416: *end = -1;
1417: }
1418: PetscCall(DMLabelLookupStratum(label, value, &v));
1419: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1420: PetscCall(DMLabelMakeValid_Private(label, v));
1421: if (label->stratumSizes[v] <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1422: PetscUseTypeMethod(label, getstratumis, v, &is);
1423: PetscCall(ISGetMinMax(is, &min, &max));
1424: PetscCall(ISDestroy(&is));
1425: if (start) *start = min;
1426: if (end) *end = max + 1;
1427: PetscFunctionReturn(PETSC_SUCCESS);
1428: }
1430: static PetscErrorCode DMLabelGetStratumIS_Concrete(DMLabel label, PetscInt v, IS *pointIS)
1431: {
1432: PetscFunctionBegin;
1433: PetscCall(PetscObjectReference((PetscObject)label->points[v]));
1434: *pointIS = label->points[v];
1435: PetscFunctionReturn(PETSC_SUCCESS);
1436: }
1438: /*@
1439: DMLabelGetStratumIS - Get an `IS` with the stratum points
1441: Not Collective
1443: Input Parameters:
1444: + label - the `DMLabel`
1445: - value - the stratum value
1447: Output Parameter:
1448: . points - The stratum points
1450: Level: intermediate
1452: Notes:
1453: The output `IS` should be destroyed when no longer needed.
1454: Returns `NULL` if the stratum is empty.
1456: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1457: @*/
1458: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1459: {
1460: PetscInt v;
1462: PetscFunctionBegin;
1464: PetscAssertPointer(points, 3);
1465: *points = NULL;
1466: PetscCall(DMLabelLookupStratum(label, value, &v));
1467: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1468: PetscCall(DMLabelMakeValid_Private(label, v));
1469: PetscUseTypeMethod(label, getstratumis, v, points);
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: /*@
1474: DMLabelSetStratumIS - Set the stratum points using an `IS`
1476: Not Collective
1478: Input Parameters:
1479: + label - the `DMLabel`
1480: . value - the stratum value
1481: - is - The stratum points
1483: Level: intermediate
1485: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1486: @*/
1487: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1488: {
1489: PetscInt v;
1491: PetscFunctionBegin;
1494: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1495: PetscCall(DMLabelLookupAddStratum(label, value, &v));
1496: if (is == label->points[v]) PetscFunctionReturn(PETSC_SUCCESS);
1497: PetscCall(DMLabelClearStratum(label, value));
1498: PetscCall(ISGetLocalSize(is, &label->stratumSizes[v]));
1499: PetscCall(PetscObjectReference((PetscObject)is));
1500: PetscCall(ISDestroy(&label->points[v]));
1501: label->points[v] = is;
1502: label->validIS[v] = PETSC_TRUE;
1503: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1504: if (label->bt) {
1505: const PetscInt *points;
1506: PetscInt p;
1508: PetscCall(ISGetIndices(is, &points));
1509: for (p = 0; p < label->stratumSizes[v]; ++p) {
1510: const PetscInt point = points[p];
1512: 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);
1513: PetscCall(PetscBTSet(label->bt, point - label->pStart));
1514: }
1515: }
1516: PetscFunctionReturn(PETSC_SUCCESS);
1517: }
1519: /*@
1520: DMLabelClearStratum - Remove a stratum
1522: Not Collective
1524: Input Parameters:
1525: + label - the `DMLabel`
1526: - value - the stratum value
1528: Level: intermediate
1530: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1531: @*/
1532: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1533: {
1534: PetscInt v;
1536: PetscFunctionBegin;
1538: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1539: PetscCall(DMLabelLookupStratum(label, value, &v));
1540: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1541: if (label->validIS[v]) {
1542: if (label->bt) {
1543: PetscInt i;
1544: const PetscInt *points;
1546: PetscCall(ISGetIndices(label->points[v], &points));
1547: for (i = 0; i < label->stratumSizes[v]; ++i) {
1548: const PetscInt point = points[i];
1550: 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);
1551: PetscCall(PetscBTClear(label->bt, point - label->pStart));
1552: }
1553: PetscCall(ISRestoreIndices(label->points[v], &points));
1554: }
1555: label->stratumSizes[v] = 0;
1556: PetscCall(ISDestroy(&label->points[v]));
1557: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]));
1558: PetscCall(PetscObjectSetName((PetscObject)label->points[v], "indices"));
1559: PetscCall(PetscObjectStateIncrease((PetscObject)label));
1560: } else {
1561: PetscCall(PetscHSetIClear(label->ht[v]));
1562: }
1563: PetscFunctionReturn(PETSC_SUCCESS);
1564: }
1566: /*@
1567: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1569: Not Collective
1571: Input Parameters:
1572: + label - The `DMLabel`
1573: . value - The label value for all points
1574: . pStart - The first point
1575: - pEnd - A point beyond all marked points
1577: Level: intermediate
1579: Note:
1580: The marks points are [`pStart`, `pEnd`), and only the bounds are stored.
1582: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1583: @*/
1584: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1585: {
1586: IS pIS;
1588: PetscFunctionBegin;
1589: PetscCall(ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS));
1590: PetscCall(DMLabelSetStratumIS(label, value, pIS));
1591: PetscCall(ISDestroy(&pIS));
1592: PetscFunctionReturn(PETSC_SUCCESS);
1593: }
1595: /*@
1596: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1598: Not Collective
1600: Input Parameters:
1601: + label - The `DMLabel`
1602: . value - The label value
1603: - p - A point with this value
1605: Output Parameter:
1606: . 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
1608: Level: intermediate
1610: .seealso: `DMLabel`, `DM`, `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1611: @*/
1612: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1613: {
1614: IS pointIS;
1615: PetscInt v;
1617: PetscFunctionBegin;
1619: PetscAssertPointer(index, 4);
1620: *index = -1;
1621: PetscCall(DMLabelLookupStratum(label, value, &v));
1622: if (v < 0) PetscFunctionReturn(PETSC_SUCCESS);
1623: PetscCall(DMLabelMakeValid_Private(label, v));
1624: PetscUseTypeMethod(label, getstratumis, v, &pointIS);
1625: PetscCall(ISLocate(pointIS, p, index));
1626: PetscCall(ISDestroy(&pointIS));
1627: PetscFunctionReturn(PETSC_SUCCESS);
1628: }
1630: /*@
1631: DMLabelFilter - Remove all points outside of [`start`, `end`)
1633: Not Collective
1635: Input Parameters:
1636: + label - the `DMLabel`
1637: . start - the first point kept
1638: - end - one more than the last point kept
1640: Level: intermediate
1642: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1643: @*/
1644: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1645: {
1646: PetscInt v;
1648: PetscFunctionBegin;
1650: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1651: PetscCall(DMLabelDestroyIndex(label));
1652: PetscCall(DMLabelMakeAllValid_Private(label));
1653: for (v = 0; v < label->numStrata; ++v) {
1654: PetscCall(ISGeneralFilter(label->points[v], start, end));
1655: PetscCall(ISGetLocalSize(label->points[v], &label->stratumSizes[v]));
1656: }
1657: PetscCall(DMLabelCreateIndex(label, start, end));
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: /*@
1662: DMLabelPermute - Create a new label with permuted points
1664: Not Collective
1666: Input Parameters:
1667: + label - the `DMLabel`
1668: - permutation - the point permutation
1670: Output Parameter:
1671: . labelNew - the new label containing the permuted points
1673: Level: intermediate
1675: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1676: @*/
1677: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1678: {
1679: const PetscInt *perm;
1680: PetscInt numValues, numPoints, v, q;
1682: PetscFunctionBegin;
1685: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1686: PetscCall(DMLabelMakeAllValid_Private(label));
1687: PetscCall(DMLabelDuplicate(label, labelNew));
1688: PetscCall(DMLabelGetNumValues(*labelNew, &numValues));
1689: PetscCall(ISGetLocalSize(permutation, &numPoints));
1690: PetscCall(ISGetIndices(permutation, &perm));
1691: for (v = 0; v < numValues; ++v) {
1692: const PetscInt size = (*labelNew)->stratumSizes[v];
1693: const PetscInt *points;
1694: PetscInt *pointsNew;
1696: PetscCall(ISGetIndices((*labelNew)->points[v], &points));
1697: PetscCall(PetscCalloc1(size, &pointsNew));
1698: for (q = 0; q < size; ++q) {
1699: const PetscInt point = points[q];
1701: 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);
1702: pointsNew[q] = perm[point];
1703: }
1704: PetscCall(ISRestoreIndices((*labelNew)->points[v], &points));
1705: PetscCall(PetscSortInt(size, pointsNew));
1706: PetscCall(ISDestroy(&(*labelNew)->points[v]));
1707: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1708: PetscCall(ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v])));
1709: PetscCall(PetscFree(pointsNew));
1710: } else {
1711: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v])));
1712: }
1713: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices"));
1714: }
1715: PetscCall(ISRestoreIndices(permutation, &perm));
1716: if (label->bt) {
1717: PetscCall(PetscBTDestroy(&label->bt));
1718: PetscCall(DMLabelCreateIndex(label, label->pStart, label->pEnd));
1719: }
1720: PetscFunctionReturn(PETSC_SUCCESS);
1721: }
1723: /*@
1724: DMLabelPermuteValues - Permute the values in a label
1726: Not collective
1728: Input Parameters:
1729: + label - the `DMLabel`
1730: - permutation - the value permutation, permutation[old value] = new value
1732: Output Parameter:
1733: . label - the `DMLabel` now with permuted values
1735: Note:
1736: The modification is done in-place
1738: Level: intermediate
1740: .seealso: `DMLabelRewriteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1741: @*/
1742: PetscErrorCode DMLabelPermuteValues(DMLabel label, IS permutation)
1743: {
1744: PetscInt Nv, Np;
1746: PetscFunctionBegin;
1749: PetscCall(DMLabelGetNumValues(label, &Nv));
1750: PetscCall(ISGetLocalSize(permutation, &Np));
1751: PetscCheck(Np == Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " != %" PetscInt_FMT " number of label values", Np, Nv);
1752: if (PetscDefined(USE_DEBUG)) {
1753: PetscBool flg;
1754: PetscCall(ISGetInfo(permutation, IS_PERMUTATION, IS_LOCAL, PETSC_TRUE, &flg));
1755: PetscCheck(flg, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "IS is not a permutation");
1756: }
1757: PetscCall(DMLabelRewriteValues(label, permutation));
1758: PetscFunctionReturn(PETSC_SUCCESS);
1759: }
1761: /*@
1762: DMLabelRewriteValues - Permute the values in a label, but some may be omitted
1764: Not collective
1766: Input Parameters:
1767: + label - the `DMLabel`
1768: - permutation - the value permutation, permutation[old value] = new value, but some maybe omitted
1770: Output Parameter:
1771: . label - the `DMLabel` now with permuted values
1773: Note:
1774: The modification is done in-place
1776: Level: intermediate
1778: .seealso: `DMLabelPermuteValues()`, `DMLabel`, `DM`, `DMLabelPermute()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1779: @*/
1780: PetscErrorCode DMLabelRewriteValues(DMLabel label, IS permutation)
1781: {
1782: const PetscInt *perm;
1783: PetscInt Nv, Np;
1785: PetscFunctionBegin;
1788: PetscCall(DMLabelMakeAllValid_Private(label));
1789: PetscCall(DMLabelGetNumValues(label, &Nv));
1790: PetscCall(ISGetLocalSize(permutation, &Np));
1791: PetscCheck(Np >= Nv, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_SIZ, "Permutation has size %" PetscInt_FMT " < %" PetscInt_FMT " number of label values", Np, Nv);
1792: PetscCall(ISGetIndices(permutation, &perm));
1793: for (PetscInt v = 0; v < Nv; ++v) label->stratumValues[v] = perm[label->stratumValues[v]];
1794: PetscCall(ISRestoreIndices(permutation, &perm));
1795: PetscFunctionReturn(PETSC_SUCCESS);
1796: }
1798: static PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1799: {
1800: MPI_Comm comm;
1801: PetscInt s, l, nroots, nleaves, offset, size;
1802: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1803: PetscSection rootSection;
1804: PetscSF labelSF;
1806: PetscFunctionBegin;
1807: if (label) PetscCall(DMLabelMakeAllValid_Private(label));
1808: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1809: /* Build a section of stratum values per point, generate the according SF
1810: and distribute point-wise stratum values to leaves. */
1811: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
1812: PetscCall(PetscSectionCreate(comm, &rootSection));
1813: PetscCall(PetscSectionSetChart(rootSection, 0, nroots));
1814: if (label) {
1815: for (s = 0; s < label->numStrata; ++s) {
1816: const PetscInt *points;
1818: PetscCall(ISGetIndices(label->points[s], &points));
1819: for (l = 0; l < label->stratumSizes[s]; l++) PetscCall(PetscSectionAddDof(rootSection, points[l], 1));
1820: PetscCall(ISRestoreIndices(label->points[s], &points));
1821: }
1822: }
1823: PetscCall(PetscSectionSetUp(rootSection));
1824: /* Create a point-wise array of stratum values */
1825: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
1826: PetscCall(PetscMalloc1(size, &rootStrata));
1827: PetscCall(PetscCalloc1(nroots, &rootIdx));
1828: if (label) {
1829: for (s = 0; s < label->numStrata; ++s) {
1830: const PetscInt *points;
1832: PetscCall(ISGetIndices(label->points[s], &points));
1833: for (l = 0; l < label->stratumSizes[s]; l++) {
1834: const PetscInt p = points[l];
1835: PetscCall(PetscSectionGetOffset(rootSection, p, &offset));
1836: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1837: }
1838: PetscCall(ISRestoreIndices(label->points[s], &points));
1839: }
1840: }
1841: /* Build SF that maps label points to remote processes */
1842: PetscCall(PetscSectionCreate(comm, leafSection));
1843: PetscCall(PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection));
1844: PetscCall(PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF));
1845: PetscCall(PetscFree(remoteOffsets));
1846: /* Send the strata for each point over the derived SF */
1847: PetscCall(PetscSectionGetStorageSize(*leafSection, &size));
1848: PetscCall(PetscMalloc1(size, leafStrata));
1849: PetscCall(PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1850: PetscCall(PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE));
1851: /* Clean up */
1852: PetscCall(PetscFree(rootStrata));
1853: PetscCall(PetscFree(rootIdx));
1854: PetscCall(PetscSectionDestroy(&rootSection));
1855: PetscCall(PetscSFDestroy(&labelSF));
1856: PetscFunctionReturn(PETSC_SUCCESS);
1857: }
1859: /*@
1860: DMLabelDistribute - Create a new label pushed forward over the `PetscSF`
1862: Collective
1864: Input Parameters:
1865: + label - the `DMLabel`
1866: - sf - the map from old to new distribution
1868: Output Parameter:
1869: . labelNew - the new redistributed label
1871: Level: intermediate
1873: .seealso: `DMLabel`, `DM`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1874: @*/
1875: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1876: {
1877: MPI_Comm comm;
1878: PetscSection leafSection;
1879: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1880: PetscInt *leafStrata, *strataIdx;
1881: PetscInt **points;
1882: const char *lname = NULL;
1883: char *name;
1884: PetscMPIInt nameSize;
1885: PetscHSetI stratumHash;
1886: size_t len = 0;
1887: PetscMPIInt rank;
1889: PetscFunctionBegin;
1891: if (label) {
1893: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
1894: PetscCall(DMLabelMakeAllValid_Private(label));
1895: }
1896: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
1897: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1898: /* Bcast name */
1899: if (rank == 0) {
1900: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
1901: PetscCall(PetscStrlen(lname, &len));
1902: }
1903: PetscCall(PetscMPIIntCast(len, &nameSize));
1904: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
1905: PetscCall(PetscMalloc1(nameSize + 1, &name));
1906: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
1907: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
1908: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
1909: PetscCall(PetscFree(name));
1910: /* Bcast defaultValue */
1911: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1912: PetscCallMPI(MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm));
1913: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1914: PetscCall(DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata));
1915: /* Determine received stratum values and initialise new label*/
1916: PetscCall(PetscHSetICreate(&stratumHash));
1917: PetscCall(PetscSectionGetStorageSize(leafSection, &size));
1918: for (p = 0; p < size; ++p) PetscCall(PetscHSetIAdd(stratumHash, leafStrata[p]));
1919: PetscCall(PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata));
1920: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS));
1921: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1922: PetscCall(PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues));
1923: /* Turn leafStrata into indices rather than stratum values */
1924: offset = 0;
1925: PetscCall(PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues));
1926: PetscCall(PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues));
1927: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscCall(PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s));
1928: for (p = 0; p < size; ++p) {
1929: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1930: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
1931: leafStrata[p] = s;
1932: break;
1933: }
1934: }
1935: }
1936: /* Rebuild the point strata on the receiver */
1937: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes));
1938: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1939: for (p = pStart; p < pEnd; p++) {
1940: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1941: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1942: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
1943: }
1944: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht));
1945: PetscCall(PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->points));
1946: PetscCall(PetscCalloc1((*labelNew)->numStrata, &points));
1947: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1948: PetscCall(PetscHSetICreate(&(*labelNew)->ht[s]));
1949: PetscCall(PetscMalloc1((*labelNew)->stratumSizes[s], &points[s]));
1950: }
1951: /* Insert points into new strata */
1952: PetscCall(PetscCalloc1((*labelNew)->numStrata, &strataIdx));
1953: PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
1954: for (p = pStart; p < pEnd; p++) {
1955: PetscCall(PetscSectionGetDof(leafSection, p, &dof));
1956: PetscCall(PetscSectionGetOffset(leafSection, p, &offset));
1957: for (s = 0; s < dof; s++) {
1958: stratum = leafStrata[offset + s];
1959: points[stratum][strataIdx[stratum]++] = p;
1960: }
1961: }
1962: for (s = 0; s < (*labelNew)->numStrata; s++) {
1963: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &points[s][0], PETSC_OWN_POINTER, &((*labelNew)->points[s])));
1964: PetscCall(PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices"));
1965: }
1966: PetscCall(PetscFree(points));
1967: PetscCall(PetscHSetIDestroy(&stratumHash));
1968: PetscCall(PetscFree(leafStrata));
1969: PetscCall(PetscFree(strataIdx));
1970: PetscCall(PetscSectionDestroy(&leafSection));
1971: PetscFunctionReturn(PETSC_SUCCESS);
1972: }
1974: /*@
1975: DMLabelGather - Gather all label values from leafs into roots
1977: Collective
1979: Input Parameters:
1980: + label - the `DMLabel`
1981: - sf - the `PetscSF` communication map
1983: Output Parameter:
1984: . labelNew - the new `DMLabel` with localised leaf values
1986: Level: developer
1988: Note:
1989: This is the inverse operation to `DMLabelDistribute()`.
1991: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
1992: @*/
1993: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1994: {
1995: MPI_Comm comm;
1996: PetscSection rootSection;
1997: PetscSF sfLabel;
1998: PetscSFNode *rootPoints, *leafPoints;
1999: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
2000: const PetscInt *rootDegree, *ilocal;
2001: PetscInt *rootStrata;
2002: const char *lname;
2003: char *name;
2004: PetscMPIInt nameSize;
2005: size_t len = 0;
2006: PetscMPIInt rank, size;
2008: PetscFunctionBegin;
2011: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2012: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
2013: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2014: PetscCallMPI(MPI_Comm_size(comm, &size));
2015: /* Bcast name */
2016: if (rank == 0) {
2017: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
2018: PetscCall(PetscStrlen(lname, &len));
2019: }
2020: PetscCall(PetscMPIIntCast(len, &nameSize));
2021: PetscCallMPI(MPI_Bcast(&nameSize, 1, MPI_INT, 0, comm));
2022: PetscCall(PetscMalloc1(nameSize + 1, &name));
2023: if (rank == 0) PetscCall(PetscArraycpy(name, lname, nameSize + 1));
2024: PetscCallMPI(MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm));
2025: PetscCall(DMLabelCreate(PETSC_COMM_SELF, name, labelNew));
2026: PetscCall(PetscFree(name));
2027: /* Gather rank/index pairs of leaves into local roots to build
2028: an inverse, multi-rooted SF. Note that this ignores local leaf
2029: indexing due to the use of the multiSF in PetscSFGather. */
2030: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL));
2031: PetscCall(PetscMalloc1(nroots, &leafPoints));
2032: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
2033: for (p = 0; p < nleaves; p++) {
2034: PetscInt ilp = ilocal ? ilocal[p] : p;
2036: leafPoints[ilp].index = ilp;
2037: leafPoints[ilp].rank = rank;
2038: }
2039: PetscCall(PetscSFComputeDegreeBegin(sf, &rootDegree));
2040: PetscCall(PetscSFComputeDegreeEnd(sf, &rootDegree));
2041: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
2042: PetscCall(PetscMalloc1(nmultiroots, &rootPoints));
2043: PetscCall(PetscSFGatherBegin(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2044: PetscCall(PetscSFGatherEnd(sf, MPIU_SF_NODE, leafPoints, rootPoints));
2045: PetscCall(PetscSFCreate(comm, &sfLabel));
2046: PetscCall(PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER));
2047: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
2048: PetscCall(DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata));
2049: /* Rebuild the point strata on the receiver */
2050: for (p = 0, idx = 0; p < nroots; p++) {
2051: for (d = 0; d < rootDegree[p]; d++) {
2052: PetscCall(PetscSectionGetDof(rootSection, idx + d, &dof));
2053: PetscCall(PetscSectionGetOffset(rootSection, idx + d, &offset));
2054: for (s = 0; s < dof; s++) PetscCall(DMLabelSetValue(*labelNew, p, rootStrata[offset + s]));
2055: }
2056: idx += rootDegree[p];
2057: }
2058: PetscCall(PetscFree(leafPoints));
2059: PetscCall(PetscFree(rootStrata));
2060: PetscCall(PetscSectionDestroy(&rootSection));
2061: PetscCall(PetscSFDestroy(&sfLabel));
2062: PetscFunctionReturn(PETSC_SUCCESS);
2063: }
2065: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
2066: {
2067: const PetscInt *degree;
2068: const PetscInt *points;
2069: PetscInt Nr, r, Nl, l, val, defVal;
2071: PetscFunctionBegin;
2072: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2073: /* Add in leaves */
2074: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2075: for (l = 0; l < Nl; ++l) {
2076: PetscCall(DMLabelGetValue(label, points[l], &val));
2077: if (val != defVal) valArray[points[l]] = val;
2078: }
2079: /* Add in shared roots */
2080: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2081: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2082: for (r = 0; r < Nr; ++r) {
2083: if (degree[r]) {
2084: PetscCall(DMLabelGetValue(label, r, &val));
2085: if (val != defVal) valArray[r] = val;
2086: }
2087: }
2088: PetscFunctionReturn(PETSC_SUCCESS);
2089: }
2091: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
2092: {
2093: const PetscInt *degree;
2094: const PetscInt *points;
2095: PetscInt Nr, r, Nl, l, val, defVal;
2097: PetscFunctionBegin;
2098: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2099: /* Read out leaves */
2100: PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL));
2101: for (l = 0; l < Nl; ++l) {
2102: const PetscInt p = points[l];
2103: const PetscInt cval = valArray[p];
2105: if (cval != defVal) {
2106: PetscCall(DMLabelGetValue(label, p, &val));
2107: if (val == defVal) {
2108: PetscCall(DMLabelSetValue(label, p, cval));
2109: if (markPoint) PetscCall((*markPoint)(label, p, cval, ctx));
2110: }
2111: }
2112: }
2113: /* Read out shared roots */
2114: PetscCall(PetscSFComputeDegreeBegin(pointSF, °ree));
2115: PetscCall(PetscSFComputeDegreeEnd(pointSF, °ree));
2116: for (r = 0; r < Nr; ++r) {
2117: if (degree[r]) {
2118: const PetscInt cval = valArray[r];
2120: if (cval != defVal) {
2121: PetscCall(DMLabelGetValue(label, r, &val));
2122: if (val == defVal) {
2123: PetscCall(DMLabelSetValue(label, r, cval));
2124: if (markPoint) PetscCall((*markPoint)(label, r, cval, ctx));
2125: }
2126: }
2127: }
2128: }
2129: PetscFunctionReturn(PETSC_SUCCESS);
2130: }
2132: /*@
2133: DMLabelPropagateBegin - Setup a cycle of label propagation
2135: Collective
2137: Input Parameters:
2138: + label - The `DMLabel` to propagate across processes
2139: - sf - The `PetscSF` describing parallel layout of the label points
2141: Level: intermediate
2143: .seealso: `DMLabel`, `DM`, `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
2144: @*/
2145: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
2146: {
2147: PetscInt Nr, r, defVal;
2148: PetscMPIInt size;
2150: PetscFunctionBegin;
2151: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2152: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
2153: if (size > 1) {
2154: PetscCall(DMLabelGetDefaultValue(label, &defVal));
2155: PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL));
2156: if (Nr >= 0) PetscCall(PetscMalloc1(Nr, &label->propArray));
2157: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
2158: }
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: /*@
2163: DMLabelPropagateEnd - Tear down a cycle of label propagation
2165: Collective
2167: Input Parameters:
2168: + label - The `DMLabel` to propagate across processes
2169: - pointSF - The `PetscSF` describing parallel layout of the label points
2171: Level: intermediate
2173: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
2174: @*/
2175: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
2176: {
2177: PetscFunctionBegin;
2178: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2179: PetscCall(PetscFree(label->propArray));
2180: label->propArray = NULL;
2181: PetscFunctionReturn(PETSC_SUCCESS);
2182: }
2184: /*@C
2185: DMLabelPropagatePush - Tear down a cycle of label propagation
2187: Collective
2189: Input Parameters:
2190: + label - The `DMLabel` to propagate across processes
2191: . pointSF - The `PetscSF` describing parallel layout of the label points
2192: . markPoint - An optional callback that is called when a point is marked, or `NULL`
2193: - ctx - An optional user context for the callback, or `NULL`
2195: Calling sequence of `markPoint`:
2196: + label - The `DMLabel`
2197: . p - The point being marked
2198: . val - The label value for `p`
2199: - ctx - An optional user context
2201: Level: intermediate
2203: .seealso: `DMLabel`, `DM`, `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
2204: @*/
2205: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel label, PetscInt p, PetscInt val, void *ctx), void *ctx)
2206: {
2207: PetscInt *valArray = label->propArray, Nr;
2208: PetscMPIInt size;
2210: PetscFunctionBegin;
2211: PetscCheck(!label->readonly, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_WRONG, "Read-only labels cannot be altered");
2212: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size));
2213: PetscCall(PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL));
2214: if (size > 1 && Nr >= 0) {
2215: /* Communicate marked edges
2216: The current implementation allocates an array the size of the number of root. We put the label values into the
2217: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
2219: TODO: We could use in-place communication with a different SF
2220: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
2221: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
2223: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
2224: 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
2225: edge to the queue.
2226: */
2227: PetscCall(DMLabelPropagateInit_Internal(label, pointSF, valArray));
2228: PetscCall(PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2229: PetscCall(PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX));
2230: PetscCall(PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2231: PetscCall(PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE));
2232: PetscCall(DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx));
2233: }
2234: PetscFunctionReturn(PETSC_SUCCESS);
2235: }
2237: /*@
2238: DMLabelConvertToSection - Make a `PetscSection`/`IS` pair that encodes the label
2240: Not Collective
2242: Input Parameter:
2243: . label - the `DMLabel`
2245: Output Parameters:
2246: + section - the section giving offsets for each stratum
2247: - is - An `IS` containing all the label points
2249: Level: developer
2251: .seealso: `DMLabel`, `DM`, `DMLabelDistribute()`
2252: @*/
2253: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2254: {
2255: IS vIS;
2256: const PetscInt *values;
2257: PetscInt *points;
2258: PetscInt nV, vS = 0, vE = 0, v, N;
2260: PetscFunctionBegin;
2262: PetscCall(DMLabelGetNumValues(label, &nV));
2263: PetscCall(DMLabelGetValueIS(label, &vIS));
2264: PetscCall(ISGetIndices(vIS, &values));
2265: if (nV) {
2266: vS = values[0];
2267: vE = values[0] + 1;
2268: }
2269: for (v = 1; v < nV; ++v) {
2270: vS = PetscMin(vS, values[v]);
2271: vE = PetscMax(vE, values[v] + 1);
2272: }
2273: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, section));
2274: PetscCall(PetscSectionSetChart(*section, vS, vE));
2275: for (v = 0; v < nV; ++v) {
2276: PetscInt n;
2278: PetscCall(DMLabelGetStratumSize(label, values[v], &n));
2279: PetscCall(PetscSectionSetDof(*section, values[v], n));
2280: }
2281: PetscCall(PetscSectionSetUp(*section));
2282: PetscCall(PetscSectionGetStorageSize(*section, &N));
2283: PetscCall(PetscMalloc1(N, &points));
2284: for (v = 0; v < nV; ++v) {
2285: IS is;
2286: const PetscInt *spoints;
2287: PetscInt dof, off, p;
2289: PetscCall(PetscSectionGetDof(*section, values[v], &dof));
2290: PetscCall(PetscSectionGetOffset(*section, values[v], &off));
2291: PetscCall(DMLabelGetStratumIS(label, values[v], &is));
2292: PetscCall(ISGetIndices(is, &spoints));
2293: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2294: PetscCall(ISRestoreIndices(is, &spoints));
2295: PetscCall(ISDestroy(&is));
2296: }
2297: PetscCall(ISRestoreIndices(vIS, &values));
2298: PetscCall(ISDestroy(&vIS));
2299: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is));
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2303: /*@C
2304: DMLabelRegister - Adds a new label component implementation
2306: Not Collective
2308: Input Parameters:
2309: + name - The name of a new user-defined creation routine
2310: - create_func - The creation routine itself
2312: Notes:
2313: `DMLabelRegister()` may be called multiple times to add several user-defined labels
2315: Example Usage:
2316: .vb
2317: DMLabelRegister("my_label", MyLabelCreate);
2318: .ve
2320: Then, your label type can be chosen with the procedural interface via
2321: .vb
2322: DMLabelCreate(MPI_Comm, DMLabel *);
2323: DMLabelSetType(DMLabel, "my_label");
2324: .ve
2325: or at runtime via the option
2326: .vb
2327: -dm_label_type my_label
2328: .ve
2330: Level: advanced
2332: .seealso: `DMLabel`, `DM`, `DMLabelType`, `DMLabelRegisterAll()`, `DMLabelRegisterDestroy()`
2333: @*/
2334: PetscErrorCode DMLabelRegister(const char name[], PetscErrorCode (*create_func)(DMLabel))
2335: {
2336: PetscFunctionBegin;
2337: PetscCall(DMInitializePackage());
2338: PetscCall(PetscFunctionListAdd(&DMLabelList, name, create_func));
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel);
2343: PETSC_EXTERN PetscErrorCode DMLabelCreate_Ephemeral(DMLabel);
2345: /*@C
2346: DMLabelRegisterAll - Registers all of the `DMLabel` implementations in the `DM` package.
2348: Not Collective
2350: Level: advanced
2352: .seealso: `DMLabel`, `DM`, `DMRegisterAll()`, `DMLabelRegisterDestroy()`
2353: @*/
2354: PetscErrorCode DMLabelRegisterAll(void)
2355: {
2356: PetscFunctionBegin;
2357: if (DMLabelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
2358: DMLabelRegisterAllCalled = PETSC_TRUE;
2360: PetscCall(DMLabelRegister(DMLABELCONCRETE, DMLabelCreate_Concrete));
2361: PetscCall(DMLabelRegister(DMLABELEPHEMERAL, DMLabelCreate_Ephemeral));
2362: PetscFunctionReturn(PETSC_SUCCESS);
2363: }
2365: /*@C
2366: DMLabelRegisterDestroy - This function destroys the `DMLabel` registry. It is called from `PetscFinalize()`.
2368: Level: developer
2370: .seealso: `DMLabel`, `DM`, `PetscInitialize()`
2371: @*/
2372: PetscErrorCode DMLabelRegisterDestroy(void)
2373: {
2374: PetscFunctionBegin;
2375: PetscCall(PetscFunctionListDestroy(&DMLabelList));
2376: DMLabelRegisterAllCalled = PETSC_FALSE;
2377: PetscFunctionReturn(PETSC_SUCCESS);
2378: }
2380: /*@
2381: DMLabelSetType - Sets the particular implementation for a label.
2383: Collective
2385: Input Parameters:
2386: + label - The label
2387: - method - The name of the label type
2389: Options Database Key:
2390: . -dm_label_type <type> - Sets the label type; use -help for a list of available types or see `DMLabelType`
2392: Level: intermediate
2394: .seealso: `DMLabel`, `DM`, `DMLabelGetType()`, `DMLabelCreate()`
2395: @*/
2396: PetscErrorCode DMLabelSetType(DMLabel label, DMLabelType method)
2397: {
2398: PetscErrorCode (*r)(DMLabel);
2399: PetscBool match;
2401: PetscFunctionBegin;
2403: PetscCall(PetscObjectTypeCompare((PetscObject)label, method, &match));
2404: if (match) PetscFunctionReturn(PETSC_SUCCESS);
2406: PetscCall(DMLabelRegisterAll());
2407: PetscCall(PetscFunctionListFind(DMLabelList, method, &r));
2408: PetscCheck(r, PetscObjectComm((PetscObject)label), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMLabel type: %s", method);
2410: PetscTryTypeMethod(label, destroy);
2411: PetscCall(PetscMemzero(label->ops, sizeof(*label->ops)));
2412: PetscCall(PetscObjectChangeTypeName((PetscObject)label, method));
2413: PetscCall((*r)(label));
2414: PetscFunctionReturn(PETSC_SUCCESS);
2415: }
2417: /*@
2418: DMLabelGetType - Gets the type name (as a string) from the label.
2420: Not Collective
2422: Input Parameter:
2423: . label - The `DMLabel`
2425: Output Parameter:
2426: . type - The `DMLabel` type name
2428: Level: intermediate
2430: .seealso: `DMLabel`, `DM`, `DMLabelSetType()`, `DMLabelCreate()`
2431: @*/
2432: PetscErrorCode DMLabelGetType(DMLabel label, DMLabelType *type)
2433: {
2434: PetscFunctionBegin;
2436: PetscAssertPointer(type, 2);
2437: PetscCall(DMLabelRegisterAll());
2438: *type = ((PetscObject)label)->type_name;
2439: PetscFunctionReturn(PETSC_SUCCESS);
2440: }
2442: static PetscErrorCode DMLabelInitialize_Concrete(DMLabel label)
2443: {
2444: PetscFunctionBegin;
2445: label->ops->view = DMLabelView_Concrete;
2446: label->ops->setup = NULL;
2447: label->ops->duplicate = DMLabelDuplicate_Concrete;
2448: label->ops->getstratumis = DMLabelGetStratumIS_Concrete;
2449: PetscFunctionReturn(PETSC_SUCCESS);
2450: }
2452: PETSC_EXTERN PetscErrorCode DMLabelCreate_Concrete(DMLabel label)
2453: {
2454: PetscFunctionBegin;
2456: PetscCall(DMLabelInitialize_Concrete(label));
2457: PetscFunctionReturn(PETSC_SUCCESS);
2458: }
2460: /*@
2461: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2462: the local section and an `PetscSF` describing the section point overlap.
2464: Collective
2466: Input Parameters:
2467: + s - The `PetscSection` for the local field layout
2468: . sf - The `PetscSF` describing parallel layout of the section points
2469: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
2470: . label - The label specifying the points
2471: - labelValue - The label stratum specifying the points
2473: Output Parameter:
2474: . gsection - The `PetscSection` for the global field layout
2476: Level: developer
2478: Note:
2479: This gives negative sizes and offsets to points not owned by this process
2481: .seealso: `DMLabel`, `DM`, `PetscSectionCreate()`
2482: @*/
2483: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2484: {
2485: PetscInt *neg = NULL, *tmpOff = NULL;
2486: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2488: PetscFunctionBegin;
2492: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
2493: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2494: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
2495: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
2496: if (nroots >= 0) {
2497: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
2498: PetscCall(PetscCalloc1(nroots, &neg));
2499: if (nroots > pEnd - pStart) {
2500: PetscCall(PetscCalloc1(nroots, &tmpOff));
2501: } else {
2502: tmpOff = &(*gsection)->atlasDof[-pStart];
2503: }
2504: }
2505: /* Mark ghost points with negative dof */
2506: for (p = pStart; p < pEnd; ++p) {
2507: PetscInt value;
2509: PetscCall(DMLabelGetValue(label, p, &value));
2510: if (value != labelValue) continue;
2511: PetscCall(PetscSectionGetDof(s, p, &dof));
2512: PetscCall(PetscSectionSetDof(*gsection, p, dof));
2513: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2514: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
2515: if (neg) neg[p] = -(dof + 1);
2516: }
2517: PetscCall(PetscSectionSetUpBC(*gsection));
2518: if (nroots >= 0) {
2519: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2520: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2521: if (nroots > pEnd - pStart) {
2522: for (p = pStart; p < pEnd; ++p) {
2523: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2524: }
2525: }
2526: }
2527: /* Calculate new sizes, get process offset, and calculate point offsets */
2528: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2529: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2530: (*gsection)->atlasOff[p] = off;
2531: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2532: }
2533: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
2534: globalOff -= off;
2535: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2536: (*gsection)->atlasOff[p] += globalOff;
2537: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2538: }
2539: /* Put in negative offsets for ghost points */
2540: if (nroots >= 0) {
2541: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2542: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
2543: if (nroots > pEnd - pStart) {
2544: for (p = pStart; p < pEnd; ++p) {
2545: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2546: }
2547: }
2548: }
2549: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
2550: PetscCall(PetscFree(neg));
2551: PetscFunctionReturn(PETSC_SUCCESS);
2552: }
2554: typedef struct _n_PetscSectionSym_Label {
2555: DMLabel label;
2556: PetscCopyMode *modes;
2557: PetscInt *sizes;
2558: const PetscInt ***perms;
2559: const PetscScalar ***rots;
2560: PetscInt (*minMaxOrients)[2];
2561: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2562: } PetscSectionSym_Label;
2564: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2565: {
2566: PetscInt i, j;
2567: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2569: PetscFunctionBegin;
2570: for (i = 0; i <= sl->numStrata; i++) {
2571: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2572: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2573: if (sl->perms[i]) PetscCall(PetscFree(sl->perms[i][j]));
2574: if (sl->rots[i]) PetscCall(PetscFree(sl->rots[i][j]));
2575: }
2576: if (sl->perms[i]) {
2577: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2579: PetscCall(PetscFree(perms));
2580: }
2581: if (sl->rots[i]) {
2582: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2584: PetscCall(PetscFree(rots));
2585: }
2586: }
2587: }
2588: PetscCall(PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients));
2589: PetscCall(DMLabelDestroy(&sl->label));
2590: sl->numStrata = 0;
2591: PetscFunctionReturn(PETSC_SUCCESS);
2592: }
2594: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2595: {
2596: PetscFunctionBegin;
2597: PetscCall(PetscSectionSymLabelReset(sym));
2598: PetscCall(PetscFree(sym->data));
2599: PetscFunctionReturn(PETSC_SUCCESS);
2600: }
2602: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2603: {
2604: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2605: PetscBool isAscii;
2606: DMLabel label = sl->label;
2607: const char *name;
2609: PetscFunctionBegin;
2610: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii));
2611: if (isAscii) {
2612: PetscInt i, j, k;
2613: PetscViewerFormat format;
2615: PetscCall(PetscViewerGetFormat(viewer, &format));
2616: if (label) {
2617: PetscCall(PetscViewerGetFormat(viewer, &format));
2618: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2619: PetscCall(PetscViewerASCIIPushTab(viewer));
2620: PetscCall(DMLabelView(label, viewer));
2621: PetscCall(PetscViewerASCIIPopTab(viewer));
2622: } else {
2623: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2624: PetscCall(PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name));
2625: }
2626: } else {
2627: PetscCall(PetscViewerASCIIPrintf(viewer, "No label given\n"));
2628: }
2629: PetscCall(PetscViewerASCIIPushTab(viewer));
2630: for (i = 0; i <= sl->numStrata; i++) {
2631: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2633: if (!(sl->perms[i] || sl->rots[i])) {
2634: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]));
2635: } else {
2636: PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]));
2637: PetscCall(PetscViewerASCIIPushTab(viewer));
2638: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]));
2639: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2640: PetscCall(PetscViewerASCIIPushTab(viewer));
2641: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2642: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2643: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j));
2644: } else {
2645: PetscInt tab;
2647: PetscCall(PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j));
2648: PetscCall(PetscViewerASCIIPushTab(viewer));
2649: PetscCall(PetscViewerASCIIGetTab(viewer, &tab));
2650: if (sl->perms[i] && sl->perms[i][j]) {
2651: PetscCall(PetscViewerASCIIPrintf(viewer, "Permutation:"));
2652: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2653: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]));
2654: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2655: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2656: }
2657: if (sl->rots[i] && sl->rots[i][j]) {
2658: PetscCall(PetscViewerASCIIPrintf(viewer, "Rotations: "));
2659: PetscCall(PetscViewerASCIISetTab(viewer, 0));
2660: #if defined(PETSC_USE_COMPLEX)
2661: 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])));
2662: #else
2663: for (k = 0; k < sl->sizes[i]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]));
2664: #endif
2665: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2666: PetscCall(PetscViewerASCIISetTab(viewer, tab));
2667: }
2668: PetscCall(PetscViewerASCIIPopTab(viewer));
2669: }
2670: }
2671: PetscCall(PetscViewerASCIIPopTab(viewer));
2672: }
2673: PetscCall(PetscViewerASCIIPopTab(viewer));
2674: }
2675: }
2676: PetscCall(PetscViewerASCIIPopTab(viewer));
2677: }
2678: PetscFunctionReturn(PETSC_SUCCESS);
2679: }
2681: /*@
2682: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2684: Logically
2686: Input Parameters:
2687: + sym - the section symmetries
2688: - label - the `DMLabel` describing the types of points
2690: Level: developer:
2692: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2693: @*/
2694: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2695: {
2696: PetscSectionSym_Label *sl;
2698: PetscFunctionBegin;
2700: sl = (PetscSectionSym_Label *)sym->data;
2701: if (sl->label && sl->label != label) PetscCall(PetscSectionSymLabelReset(sym));
2702: if (label) {
2703: sl->label = label;
2704: PetscCall(PetscObjectReference((PetscObject)label));
2705: PetscCall(DMLabelGetNumValues(label, &sl->numStrata));
2706: 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));
2707: PetscCall(PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode)));
2708: PetscCall(PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt)));
2709: PetscCall(PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **)));
2710: PetscCall(PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **)));
2711: PetscCall(PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2])));
2712: }
2713: PetscFunctionReturn(PETSC_SUCCESS);
2714: }
2716: /*@C
2717: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2719: Logically Collective
2721: Input Parameters:
2722: + sym - the section symmetries
2723: - stratum - the stratum value in the label that we are assigning symmetries for
2725: Output Parameters:
2726: + size - the number of dofs for points in the `stratum` of the label
2727: . minOrient - the smallest orientation for a point in this `stratum`
2728: . maxOrient - one greater than the largest orientation for a ppoint in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2729: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2730: - 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
2732: Level: developer
2734: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2735: @*/
2736: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2737: {
2738: PetscSectionSym_Label *sl;
2739: const char *name;
2740: PetscInt i;
2742: PetscFunctionBegin;
2744: sl = (PetscSectionSym_Label *)sym->data;
2745: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2746: for (i = 0; i <= sl->numStrata; i++) {
2747: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2749: if (stratum == value) break;
2750: }
2751: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2752: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2753: if (size) {
2754: PetscAssertPointer(size, 3);
2755: *size = sl->sizes[i];
2756: }
2757: if (minOrient) {
2758: PetscAssertPointer(minOrient, 4);
2759: *minOrient = sl->minMaxOrients[i][0];
2760: }
2761: if (maxOrient) {
2762: PetscAssertPointer(maxOrient, 5);
2763: *maxOrient = sl->minMaxOrients[i][1];
2764: }
2765: if (perms) {
2766: PetscAssertPointer(perms, 6);
2767: *perms = PetscSafePointerPlusOffset(sl->perms[i], sl->minMaxOrients[i][0]);
2768: }
2769: if (rots) {
2770: PetscAssertPointer(rots, 7);
2771: *rots = PetscSafePointerPlusOffset(sl->rots[i], sl->minMaxOrients[i][0]);
2772: }
2773: PetscFunctionReturn(PETSC_SUCCESS);
2774: }
2776: /*@C
2777: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2779: Logically
2781: Input Parameters:
2782: + sym - the section symmetries
2783: . stratum - the stratum value in the label that we are assigning symmetries for
2784: . size - the number of dofs for points in the `stratum` of the label
2785: . minOrient - the smallest orientation for a point in this `stratum`
2786: . maxOrient - one greater than the largest orientation for a point in this `stratum` (i.e., orientations are in the range [`minOrient`, `maxOrient`))
2787: . mode - how `sym` should copy the `perms` and `rots` arrays
2788: . perms - `NULL` if there are no permutations, or (`maxOrient` - `minOrient`) permutations, one for each orientation. A `NULL` permutation is the identity
2789: - 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
2791: Level: developer
2793: .seealso: `DMLabel`, `DM`, `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2794: @*/
2795: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2796: {
2797: PetscSectionSym_Label *sl;
2798: const char *name;
2799: PetscInt i, j, k;
2801: PetscFunctionBegin;
2803: sl = (PetscSectionSym_Label *)sym->data;
2804: PetscCheck(sl->label, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_WRONGSTATE, "No label set yet");
2805: for (i = 0; i <= sl->numStrata; i++) {
2806: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2808: if (stratum == value) break;
2809: }
2810: PetscCall(PetscObjectGetName((PetscObject)sl->label, &name));
2811: PetscCheck(i <= sl->numStrata, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_OUTOFRANGE, "Stratum %" PetscInt_FMT " not found in label %s", stratum, name);
2812: sl->sizes[i] = size;
2813: sl->modes[i] = mode;
2814: sl->minMaxOrients[i][0] = minOrient;
2815: sl->minMaxOrients[i][1] = maxOrient;
2816: if (mode == PETSC_COPY_VALUES) {
2817: if (perms) {
2818: PetscInt **ownPerms;
2820: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownPerms));
2821: for (j = 0; j < maxOrient - minOrient; j++) {
2822: if (perms[j]) {
2823: PetscCall(PetscMalloc1(size, &ownPerms[j]));
2824: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2825: }
2826: }
2827: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2828: }
2829: if (rots) {
2830: PetscScalar **ownRots;
2832: PetscCall(PetscCalloc1(maxOrient - minOrient, &ownRots));
2833: for (j = 0; j < maxOrient - minOrient; j++) {
2834: if (rots[j]) {
2835: PetscCall(PetscMalloc1(size, &ownRots[j]));
2836: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2837: }
2838: }
2839: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2840: }
2841: } else {
2842: sl->perms[i] = PetscSafePointerPlusOffset(perms, -minOrient);
2843: sl->rots[i] = PetscSafePointerPlusOffset(rots, -minOrient);
2844: }
2845: PetscFunctionReturn(PETSC_SUCCESS);
2846: }
2848: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2849: {
2850: PetscInt i, j, numStrata;
2851: PetscSectionSym_Label *sl;
2852: DMLabel label;
2854: PetscFunctionBegin;
2855: sl = (PetscSectionSym_Label *)sym->data;
2856: numStrata = sl->numStrata;
2857: label = sl->label;
2858: for (i = 0; i < numPoints; i++) {
2859: PetscInt point = points[2 * i];
2860: PetscInt ornt = points[2 * i + 1];
2862: for (j = 0; j < numStrata; j++) {
2863: if (label->validIS[j]) {
2864: PetscInt k;
2866: PetscCall(ISLocate(label->points[j], point, &k));
2867: if (k >= 0) break;
2868: } else {
2869: PetscBool has;
2871: PetscCall(PetscHSetIHas(label->ht[j], point, &has));
2872: if (has) break;
2873: }
2874: }
2875: 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],
2876: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2877: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2878: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2879: }
2880: PetscFunctionReturn(PETSC_SUCCESS);
2881: }
2883: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2884: {
2885: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2886: IS valIS;
2887: const PetscInt *values;
2888: PetscInt Nv, v;
2890: PetscFunctionBegin;
2891: PetscCall(DMLabelGetNumValues(sl->label, &Nv));
2892: PetscCall(DMLabelGetValueIS(sl->label, &valIS));
2893: PetscCall(ISGetIndices(valIS, &values));
2894: for (v = 0; v < Nv; ++v) {
2895: const PetscInt val = values[v];
2896: PetscInt size, minOrient, maxOrient;
2897: const PetscInt **perms;
2898: const PetscScalar **rots;
2900: PetscCall(PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots));
2901: PetscCall(PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots));
2902: }
2903: PetscCall(ISDestroy(&valIS));
2904: PetscFunctionReturn(PETSC_SUCCESS);
2905: }
2907: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2908: {
2909: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2910: DMLabel dlabel;
2912: PetscFunctionBegin;
2913: PetscCall(DMLabelDistribute(sl->label, migrationSF, &dlabel));
2914: PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym));
2915: PetscCall(DMLabelDestroy(&dlabel));
2916: PetscCall(PetscSectionSymCopy(sym, *dsym));
2917: PetscFunctionReturn(PETSC_SUCCESS);
2918: }
2920: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2921: {
2922: PetscSectionSym_Label *sl;
2924: PetscFunctionBegin;
2925: PetscCall(PetscNew(&sl));
2926: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2927: sym->ops->distribute = PetscSectionSymDistribute_Label;
2928: sym->ops->copy = PetscSectionSymCopy_Label;
2929: sym->ops->view = PetscSectionSymView_Label;
2930: sym->ops->destroy = PetscSectionSymDestroy_Label;
2931: sym->data = (void *)sl;
2932: PetscFunctionReturn(PETSC_SUCCESS);
2933: }
2935: /*@
2936: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2938: Collective
2940: Input Parameters:
2941: + comm - the MPI communicator for the new symmetry
2942: - label - the label defining the strata
2944: Output Parameter:
2945: . sym - the section symmetries
2947: Level: developer
2949: .seealso: `DMLabel`, `DM`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
2950: @*/
2951: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2952: {
2953: PetscFunctionBegin;
2954: PetscCall(DMInitializePackage());
2955: PetscCall(PetscSectionSymCreate(comm, sym));
2956: PetscCall(PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL));
2957: PetscCall(PetscSectionSymLabelSetLabel(*sym, label));
2958: PetscFunctionReturn(PETSC_SUCCESS);
2959: }