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: /*@C
8: DMLabelCreate - Create a DMLabel object, which is a multimap
10: Collective
12: Input parameters:
13: + comm - The communicator, usually PETSC_COMM_SELF
14: - name - The label name
16: Output parameter:
17: . label - The DMLabel
19: Level: beginner
21: Notes:
22: The label name is actually usual PetscObject name.
23: One can get/set it with PetscObjectGetName()/PetscObjectSetName().
25: .seealso: `DMLabelDestroy()`
26: @*/
27: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
28: {
30: DMInitializePackage();
32: PetscHeaderCreate(*label, DMLABEL_CLASSID, "DMLabel", "DMLabel", "DM", comm, DMLabelDestroy, DMLabelView);
34: (*label)->numStrata = 0;
35: (*label)->defaultValue = -1;
36: (*label)->stratumValues = NULL;
37: (*label)->validIS = NULL;
38: (*label)->stratumSizes = NULL;
39: (*label)->points = NULL;
40: (*label)->ht = NULL;
41: (*label)->pStart = -1;
42: (*label)->pEnd = -1;
43: (*label)->bt = NULL;
44: PetscHMapICreate(&(*label)->hmap);
45: PetscObjectSetName((PetscObject)*label, name);
46: return 0;
47: }
49: /*
50: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
52: Not collective
54: Input parameter:
55: + label - The DMLabel
56: - v - The stratum value
58: Output parameter:
59: . label - The DMLabel with stratum in sorted list format
61: Level: developer
63: .seealso: `DMLabelCreate()`
64: */
65: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
66: {
67: IS is;
68: PetscInt off = 0, *pointArray, p;
70: if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
72: PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
73: PetscMalloc1(label->stratumSizes[v], &pointArray);
74: PetscHSetIGetElems(label->ht[v], &off, pointArray);
75: PetscHSetIClear(label->ht[v]);
76: PetscSortInt(label->stratumSizes[v], pointArray);
77: if (label->bt) {
78: for (p = 0; p < label->stratumSizes[v]; ++p) {
79: const PetscInt point = pointArray[p];
81: PetscBTSet(label->bt, point - label->pStart);
82: }
83: }
84: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v] - 1] == pointArray[0] + label->stratumSizes[v] - 1) {
85: ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
86: PetscFree(pointArray);
87: } else {
88: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
89: }
90: PetscObjectSetName((PetscObject)is, "indices");
91: label->points[v] = is;
92: label->validIS[v] = PETSC_TRUE;
93: PetscObjectStateIncrease((PetscObject)label);
94: return 0;
95: }
97: /*
98: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
100: Not collective
102: Input parameter:
103: . label - The DMLabel
105: Output parameter:
106: . label - The DMLabel with all strata in sorted list format
108: Level: developer
110: .seealso: `DMLabelCreate()`
111: */
112: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
113: {
114: PetscInt v;
116: for (v = 0; v < label->numStrata; v++) DMLabelMakeValid_Private(label, v);
117: return 0;
118: }
120: /*
121: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
123: Not collective
125: Input parameter:
126: + label - The DMLabel
127: - v - The stratum value
129: Output parameter:
130: . label - The DMLabel with stratum in hash format
132: Level: developer
134: .seealso: `DMLabelCreate()`
135: */
136: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
137: {
138: PetscInt p;
139: const PetscInt *points;
141: if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
143: if (label->points[v]) {
144: ISGetIndices(label->points[v], &points);
145: for (p = 0; p < label->stratumSizes[v]; ++p) PetscHSetIAdd(label->ht[v], points[p]);
146: ISRestoreIndices(label->points[v], &points);
147: ISDestroy(&(label->points[v]));
148: }
149: label->validIS[v] = PETSC_FALSE;
150: return 0;
151: }
153: PetscErrorCode DMLabelMakeAllInvalid_Internal(DMLabel label)
154: {
155: PetscInt v;
157: for (v = 0; v < label->numStrata; v++) DMLabelMakeInvalid_Private(label, v);
158: return 0;
159: }
161: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
162: #define DMLABEL_LOOKUP_THRESHOLD 16
163: #endif
165: static inline PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
166: {
167: PetscInt v;
169: *index = -1;
170: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
171: for (v = 0; v < label->numStrata; ++v)
172: if (label->stratumValues[v] == value) {
173: *index = v;
174: break;
175: }
176: } else {
177: PetscHMapIGet(label->hmap, value, index);
178: }
179: if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
180: PetscInt len, loc = -1;
181: PetscHMapIGetSize(label->hmap, &len);
183: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
184: PetscHMapIGet(label->hmap, value, &loc);
185: } else {
186: for (v = 0; v < label->numStrata; ++v)
187: if (label->stratumValues[v] == value) {
188: loc = v;
189: break;
190: }
191: }
193: }
194: return 0;
195: }
197: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
198: {
199: PetscInt v;
200: PetscInt *tmpV;
201: PetscInt *tmpS;
202: PetscHSetI *tmpH, ht;
203: IS *tmpP, is;
204: PetscBool *tmpB;
205: PetscHMapI hmap = label->hmap;
207: v = label->numStrata;
208: tmpV = label->stratumValues;
209: tmpS = label->stratumSizes;
210: tmpH = label->ht;
211: tmpP = label->points;
212: tmpB = label->validIS;
213: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
214: PetscInt *oldV = tmpV;
215: PetscInt *oldS = tmpS;
216: PetscHSetI *oldH = tmpH;
217: IS *oldP = tmpP;
218: PetscBool *oldB = tmpB;
219: PetscMalloc((v + 1) * sizeof(*tmpV), &tmpV);
220: PetscMalloc((v + 1) * sizeof(*tmpS), &tmpS);
221: PetscMalloc((v + 1) * sizeof(*tmpH), &tmpH);
222: PetscMalloc((v + 1) * sizeof(*tmpP), &tmpP);
223: PetscMalloc((v + 1) * sizeof(*tmpB), &tmpB);
224: PetscArraycpy(tmpV, oldV, v);
225: PetscArraycpy(tmpS, oldS, v);
226: PetscArraycpy(tmpH, oldH, v);
227: PetscArraycpy(tmpP, oldP, v);
228: PetscArraycpy(tmpB, oldB, v);
229: PetscFree(oldV);
230: PetscFree(oldS);
231: PetscFree(oldH);
232: PetscFree(oldP);
233: PetscFree(oldB);
234: }
235: label->numStrata = v + 1;
236: label->stratumValues = tmpV;
237: label->stratumSizes = tmpS;
238: label->ht = tmpH;
239: label->points = tmpP;
240: label->validIS = tmpB;
241: PetscHSetICreate(&ht);
242: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is);
243: PetscHMapISet(hmap, value, v);
244: tmpV[v] = value;
245: tmpS[v] = 0;
246: tmpH[v] = ht;
247: tmpP[v] = is;
248: tmpB[v] = PETSC_TRUE;
249: PetscObjectStateIncrease((PetscObject)label);
250: *index = v;
251: return 0;
252: }
254: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
255: {
256: DMLabelLookupStratum(label, value, index);
257: if (*index < 0) DMLabelNewStratum(label, value, index);
258: return 0;
259: }
261: static inline PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
262: {
263: *size = 0;
264: if (v < 0) return 0;
265: if (label->validIS[v]) {
266: *size = label->stratumSizes[v];
267: } else {
268: PetscHSetIGetSize(label->ht[v], size);
269: }
270: return 0;
271: }
273: /*@
274: DMLabelAddStratum - Adds a new stratum value in a DMLabel
276: Input Parameters:
277: + label - The DMLabel
278: - value - The stratum value
280: Level: beginner
282: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
283: @*/
284: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
285: {
286: PetscInt v;
289: DMLabelLookupAddStratum(label, value, &v);
290: return 0;
291: }
293: /*@
294: DMLabelAddStrata - Adds new stratum values in a DMLabel
296: Not collective
298: Input Parameters:
299: + label - The DMLabel
300: . numStrata - The number of stratum values
301: - stratumValues - The stratum values
303: Level: beginner
305: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
306: @*/
307: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
308: {
309: PetscInt *values, v;
313: PetscMalloc1(numStrata, &values);
314: PetscArraycpy(values, stratumValues, numStrata);
315: PetscSortRemoveDupsInt(&numStrata, values);
316: if (!label->numStrata) { /* Fast preallocation */
317: PetscInt *tmpV;
318: PetscInt *tmpS;
319: PetscHSetI *tmpH, ht;
320: IS *tmpP, is;
321: PetscBool *tmpB;
322: PetscHMapI hmap = label->hmap;
324: PetscMalloc1(numStrata, &tmpV);
325: PetscMalloc1(numStrata, &tmpS);
326: PetscMalloc1(numStrata, &tmpH);
327: PetscMalloc1(numStrata, &tmpP);
328: PetscMalloc1(numStrata, &tmpB);
329: label->numStrata = numStrata;
330: label->stratumValues = tmpV;
331: label->stratumSizes = tmpS;
332: label->ht = tmpH;
333: label->points = tmpP;
334: label->validIS = tmpB;
335: for (v = 0; v < numStrata; ++v) {
336: PetscHSetICreate(&ht);
337: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &is);
338: PetscHMapISet(hmap, values[v], v);
339: tmpV[v] = values[v];
340: tmpS[v] = 0;
341: tmpH[v] = ht;
342: tmpP[v] = is;
343: tmpB[v] = PETSC_TRUE;
344: }
345: PetscObjectStateIncrease((PetscObject)label);
346: } else {
347: for (v = 0; v < numStrata; ++v) DMLabelAddStratum(label, values[v]);
348: }
349: PetscFree(values);
350: return 0;
351: }
353: /*@
354: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
356: Not collective
358: Input Parameters:
359: + label - The DMLabel
360: - valueIS - Index set with stratum values
362: Level: beginner
364: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
365: @*/
366: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
367: {
368: PetscInt numStrata;
369: const PetscInt *stratumValues;
373: ISGetLocalSize(valueIS, &numStrata);
374: ISGetIndices(valueIS, &stratumValues);
375: DMLabelAddStrata(label, numStrata, stratumValues);
376: return 0;
377: }
379: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
380: {
381: PetscInt v;
382: PetscMPIInt rank;
384: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
385: PetscViewerASCIIPushSynchronized(viewer);
386: if (label) {
387: const char *name;
389: PetscObjectGetName((PetscObject)label, &name);
390: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
391: if (label->bt) PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", label->pStart, label->pEnd);
392: for (v = 0; v < label->numStrata; ++v) {
393: const PetscInt value = label->stratumValues[v];
394: const PetscInt *points;
395: PetscInt p;
397: ISGetIndices(label->points[v], &points);
398: for (p = 0; p < label->stratumSizes[v]; ++p) PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %" PetscInt_FMT " (%" PetscInt_FMT ")\n", rank, points[p], value);
399: ISRestoreIndices(label->points[v], &points);
400: }
401: }
402: PetscViewerFlush(viewer);
403: PetscViewerASCIIPopSynchronized(viewer);
404: return 0;
405: }
407: /*@C
408: DMLabelView - View the label
410: Collective on viewer
412: Input Parameters:
413: + label - The DMLabel
414: - viewer - The PetscViewer
416: Level: intermediate
418: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
419: @*/
420: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
421: {
422: PetscBool iascii;
425: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);
427: if (label) DMLabelMakeAllValid_Private(label);
428: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
429: if (iascii) DMLabelView_Ascii(label, viewer);
430: return 0;
431: }
433: /*@
434: DMLabelReset - Destroys internal data structures in a DMLabel
436: Not collective
438: Input Parameter:
439: . label - The DMLabel
441: Level: beginner
443: .seealso: `DMLabelDestroy()`, `DMLabelCreate()`
444: @*/
445: PetscErrorCode DMLabelReset(DMLabel label)
446: {
447: PetscInt v;
450: for (v = 0; v < label->numStrata; ++v) {
451: PetscHSetIDestroy(&label->ht[v]);
452: ISDestroy(&label->points[v]);
453: }
454: label->numStrata = 0;
455: PetscFree(label->stratumValues);
456: PetscFree(label->stratumSizes);
457: PetscFree(label->ht);
458: PetscFree(label->points);
459: PetscFree(label->validIS);
460: label->stratumValues = NULL;
461: label->stratumSizes = NULL;
462: label->ht = NULL;
463: label->points = NULL;
464: label->validIS = NULL;
465: PetscHMapIReset(label->hmap);
466: label->pStart = -1;
467: label->pEnd = -1;
468: PetscBTDestroy(&label->bt);
469: return 0;
470: }
472: /*@
473: DMLabelDestroy - Destroys a DMLabel
475: Collective on label
477: Input Parameter:
478: . label - The DMLabel
480: Level: beginner
482: .seealso: `DMLabelReset()`, `DMLabelCreate()`
483: @*/
484: PetscErrorCode DMLabelDestroy(DMLabel *label)
485: {
486: if (!*label) return 0;
488: if (--((PetscObject)(*label))->refct > 0) {
489: *label = NULL;
490: return 0;
491: }
492: DMLabelReset(*label);
493: PetscHMapIDestroy(&(*label)->hmap);
494: PetscHeaderDestroy(label);
495: return 0;
496: }
498: /*@
499: DMLabelDuplicate - Duplicates a DMLabel
501: Collective on label
503: Input Parameter:
504: . label - The DMLabel
506: Output Parameter:
507: . labelnew - location to put new vector
509: Level: intermediate
511: .seealso: `DMLabelCreate()`, `DMLabelDestroy()`
512: @*/
513: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
514: {
515: const char *name;
516: PetscInt v;
519: DMLabelMakeAllValid_Private(label);
520: PetscObjectGetName((PetscObject)label, &name);
521: DMLabelCreate(PetscObjectComm((PetscObject)label), name, labelnew);
523: (*labelnew)->numStrata = label->numStrata;
524: (*labelnew)->defaultValue = label->defaultValue;
525: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
526: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
527: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
528: PetscMalloc1(label->numStrata, &(*labelnew)->points);
529: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
530: for (v = 0; v < label->numStrata; ++v) {
531: PetscHSetICreate(&(*labelnew)->ht[v]);
532: (*labelnew)->stratumValues[v] = label->stratumValues[v];
533: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
534: PetscObjectReference((PetscObject)(label->points[v]));
535: (*labelnew)->points[v] = label->points[v];
536: (*labelnew)->validIS[v] = PETSC_TRUE;
537: }
538: PetscHMapIDestroy(&(*labelnew)->hmap);
539: PetscHMapIDuplicate(label->hmap, &(*labelnew)->hmap);
540: (*labelnew)->pStart = -1;
541: (*labelnew)->pEnd = -1;
542: (*labelnew)->bt = NULL;
543: return 0;
544: }
546: /*@C
547: DMLabelCompare - Compare two DMLabel objects
549: Collective on comm
551: Input Parameters:
552: + comm - Comm over which to compare labels
553: . l0 - First DMLabel
554: - l1 - Second DMLabel
556: Output Parameters
557: + equal - (Optional) Flag whether the two labels are equal
558: - message - (Optional) Message describing the difference
560: Level: intermediate
562: Notes:
563: The output flag equal is the same on all processes.
564: If it is passed as NULL and difference is found, an error is thrown on all processes.
565: Make sure to pass NULL on all processes.
567: The output message is set independently on each rank.
568: It is set to NULL if no difference was found on the current rank. It must be freed by user.
569: If message is passed as NULL and difference is found, the difference description is printed to stderr in synchronized manner.
570: Make sure to pass NULL on all processes.
572: For the comparison, we ignore the order of stratum values, and strata with no points.
574: The communicator needs to be specified because currently DMLabel can live on PETSC_COMM_SELF even if the underlying DM is parallel.
576: Fortran Notes:
577: This function is currently not available from Fortran.
579: .seealso: `DMCompareLabels()`, `DMLabelGetNumValues()`, `DMLabelGetDefaultValue()`, `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelGetStratumIS()`
580: @*/
581: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
582: {
583: const char *name0, *name1;
584: char msg[PETSC_MAX_PATH_LEN] = "";
585: PetscBool eq;
586: PetscMPIInt rank;
592: MPI_Comm_rank(comm, &rank);
593: PetscObjectGetName((PetscObject)l0, &name0);
594: PetscObjectGetName((PetscObject)l1, &name1);
595: {
596: PetscInt v0, v1;
598: DMLabelGetDefaultValue(l0, &v0);
599: DMLabelGetDefaultValue(l1, &v1);
600: eq = (PetscBool)(v0 == v1);
601: if (!eq) PetscSNPrintf(msg, sizeof(msg), "Default value of DMLabel l0 \"%s\" = %" PetscInt_FMT " != %" PetscInt_FMT " = Default value of DMLabel l1 \"%s\"", name0, v0, v1, name1);
602: MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
603: if (!eq) goto finish;
604: }
605: {
606: IS is0, is1;
608: DMLabelGetNonEmptyStratumValuesIS(l0, &is0);
609: DMLabelGetNonEmptyStratumValuesIS(l1, &is1);
610: ISEqual(is0, is1, &eq);
611: ISDestroy(&is0);
612: ISDestroy(&is1);
613: if (!eq) PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1);
614: MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
615: if (!eq) goto finish;
616: }
617: {
618: PetscInt i, nValues;
620: DMLabelGetNumValues(l0, &nValues);
621: for (i = 0; i < nValues; i++) {
622: const PetscInt v = l0->stratumValues[i];
623: PetscInt n;
624: IS is0, is1;
626: DMLabelGetStratumSize_Private(l0, i, &n);
627: if (!n) continue;
628: DMLabelGetStratumIS(l0, v, &is0);
629: DMLabelGetStratumIS(l1, v, &is1);
630: ISEqualUnsorted(is0, is1, &eq);
631: ISDestroy(&is0);
632: ISDestroy(&is1);
633: if (!eq) {
634: 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);
635: break;
636: }
637: }
638: MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
639: }
640: finish:
641: /* If message output arg not set, print to stderr */
642: if (message) {
643: *message = NULL;
644: if (msg[0]) PetscStrallocpy(msg, message);
645: } else {
646: if (msg[0]) PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg);
647: PetscSynchronizedFlush(comm, PETSC_STDERR);
648: }
649: /* If same output arg not ser and labels are not equal, throw error */
650: if (equal) *equal = eq;
652: return 0;
653: }
655: /*@
656: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
658: Not collective
660: Input Parameter:
661: . label - The DMLabel
663: Level: intermediate
665: .seealso: `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
666: @*/
667: PetscErrorCode DMLabelComputeIndex(DMLabel label)
668: {
669: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
672: DMLabelMakeAllValid_Private(label);
673: for (v = 0; v < label->numStrata; ++v) {
674: const PetscInt *points;
675: PetscInt i;
677: ISGetIndices(label->points[v], &points);
678: for (i = 0; i < label->stratumSizes[v]; ++i) {
679: const PetscInt point = points[i];
681: pStart = PetscMin(point, pStart);
682: pEnd = PetscMax(point + 1, pEnd);
683: }
684: ISRestoreIndices(label->points[v], &points);
685: }
686: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
687: label->pEnd = pEnd;
688: DMLabelCreateIndex(label, label->pStart, label->pEnd);
689: return 0;
690: }
692: /*@
693: DMLabelCreateIndex - Create an index structure for membership determination
695: Not collective
697: Input Parameters:
698: + label - The DMLabel
699: . pStart - The smallest point
700: - pEnd - The largest point + 1
702: Level: intermediate
704: .seealso: `DMLabelHasPoint()`, `DMLabelComputeIndex()`, `DMLabelDestroyIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
705: @*/
706: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
707: {
708: PetscInt v;
711: DMLabelDestroyIndex(label);
712: DMLabelMakeAllValid_Private(label);
713: label->pStart = pStart;
714: label->pEnd = pEnd;
715: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
716: PetscBTCreate(pEnd - pStart, &label->bt);
717: for (v = 0; v < label->numStrata; ++v) {
718: const PetscInt *points;
719: PetscInt i;
721: ISGetIndices(label->points[v], &points);
722: for (i = 0; i < label->stratumSizes[v]; ++i) {
723: const PetscInt point = points[i];
726: PetscBTSet(label->bt, point - pStart);
727: }
728: ISRestoreIndices(label->points[v], &points);
729: }
730: return 0;
731: }
733: /*@
734: DMLabelDestroyIndex - Destroy the index structure
736: Not collective
738: Input Parameter:
739: . label - the DMLabel
741: Level: intermediate
743: .seealso: `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
744: @*/
745: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
746: {
748: label->pStart = -1;
749: label->pEnd = -1;
750: PetscBTDestroy(&label->bt);
751: return 0;
752: }
754: /*@
755: DMLabelGetBounds - Return the smallest and largest point in the label
757: Not collective
759: Input Parameter:
760: . label - the DMLabel
762: Output Parameters:
763: + pStart - The smallest point
764: - pEnd - The largest point + 1
766: Note: This will compute an index for the label if one does not exist.
768: Level: intermediate
770: .seealso: `DMLabelHasPoint()`, `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
771: @*/
772: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
773: {
775: if ((label->pStart == -1) && (label->pEnd == -1)) DMLabelComputeIndex(label);
776: if (pStart) {
778: *pStart = label->pStart;
779: }
780: if (pEnd) {
782: *pEnd = label->pEnd;
783: }
784: return 0;
785: }
787: /*@
788: DMLabelHasValue - Determine whether a label assigns the value to any point
790: Not collective
792: Input Parameters:
793: + label - the DMLabel
794: - value - the value
796: Output Parameter:
797: . contains - Flag indicating whether the label maps this value to any point
799: Level: developer
801: .seealso: `DMLabelHasPoint()`, `DMLabelGetValue()`, `DMLabelSetValue()`
802: @*/
803: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
804: {
805: PetscInt v;
809: DMLabelLookupStratum(label, value, &v);
810: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
811: return 0;
812: }
814: /*@
815: DMLabelHasPoint - Determine whether a label assigns a value to a point
817: Not collective
819: Input Parameters:
820: + label - the DMLabel
821: - point - the point
823: Output Parameter:
824: . contains - Flag indicating whether the label maps this point to a value
826: Note: The user must call DMLabelCreateIndex() before this function.
828: Level: developer
830: .seealso: `DMLabelCreateIndex()`, `DMLabelGetValue()`, `DMLabelSetValue()`
831: @*/
832: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
833: {
837: DMLabelMakeAllValid_Private(label);
838: if (PetscDefined(USE_DEBUG)) {
841: }
842: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
843: return 0;
844: }
846: /*@
847: DMLabelStratumHasPoint - Return true if the stratum contains a point
849: Not collective
851: Input Parameters:
852: + label - the DMLabel
853: . value - the stratum value
854: - point - the point
856: Output Parameter:
857: . contains - true if the stratum contains the point
859: Level: intermediate
861: .seealso: `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`
862: @*/
863: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
864: {
868: if (value == label->defaultValue) {
869: PetscInt pointVal;
871: DMLabelGetValue(label, point, &pointVal);
872: *contains = (PetscBool)(pointVal == value);
873: } else {
874: PetscInt v;
876: DMLabelLookupStratum(label, value, &v);
877: if (v >= 0) {
878: if (label->validIS[v]) {
879: PetscInt i;
881: ISLocate(label->points[v], point, &i);
882: *contains = (PetscBool)(i >= 0);
883: } else {
884: PetscHSetIHas(label->ht[v], point, contains);
885: }
886: } else { // value is not present
887: *contains = PETSC_FALSE;
888: }
889: }
890: return 0;
891: }
893: /*@
894: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
895: When a label is created, it is initialized to -1.
897: Not collective
899: Input parameter:
900: . label - a DMLabel object
902: Output parameter:
903: . defaultValue - the default value
905: Level: beginner
907: .seealso: `DMLabelSetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
908: @*/
909: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
910: {
912: *defaultValue = label->defaultValue;
913: return 0;
914: }
916: /*@
917: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
918: When a label is created, it is initialized to -1.
920: Not collective
922: Input parameter:
923: . label - a DMLabel object
925: Output parameter:
926: . defaultValue - the default value
928: Level: beginner
930: .seealso: `DMLabelGetDefaultValue()`, `DMLabelGetValue()`, `DMLabelSetValue()`
931: @*/
932: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
933: {
935: label->defaultValue = defaultValue;
936: return 0;
937: }
939: /*@
940: 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 DMLabelSetDefaultValue())
942: Not collective
944: Input Parameters:
945: + label - the DMLabel
946: - point - the point
948: Output Parameter:
949: . value - The point value, or the default value (-1 by default)
951: Note: a label may assign multiple values to a point. No guarantees are made about which value is returned in that case. Use `DMLabelStratumHasPoint()` to check for inclusion in a specific value stratum.
953: Level: intermediate
955: .seealso: `DMLabelCreate()`, `DMLabelSetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
956: @*/
957: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
958: {
959: PetscInt v;
964: *value = label->defaultValue;
965: for (v = 0; v < label->numStrata; ++v) {
966: if (label->validIS[v]) {
967: PetscInt i;
969: ISLocate(label->points[v], point, &i);
970: if (i >= 0) {
971: *value = label->stratumValues[v];
972: break;
973: }
974: } else {
975: PetscBool has;
977: PetscHSetIHas(label->ht[v], point, &has);
978: if (has) {
979: *value = label->stratumValues[v];
980: break;
981: }
982: }
983: }
984: return 0;
985: }
987: /*@
988: 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 be changed with DMLabelSetDefaultValue() to something different), then this function will do nothing.
990: Not collective
992: Input Parameters:
993: + label - the DMLabel
994: . point - the point
995: - value - The point value
997: Level: intermediate
999: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelClearValue()`, `DMLabelGetDefaultValue()`, `DMLabelSetDefaultValue()`
1000: @*/
1001: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
1002: {
1003: PetscInt v;
1006: /* Find label value, add new entry if needed */
1007: if (value == label->defaultValue) return 0;
1008: DMLabelLookupAddStratum(label, value, &v);
1009: /* Set key */
1010: DMLabelMakeInvalid_Private(label, v);
1011: PetscHSetIAdd(label->ht[v], point);
1012: return 0;
1013: }
1015: /*@
1016: DMLabelClearValue - Clear the value a label assigns to a point
1018: Not collective
1020: Input Parameters:
1021: + label - the DMLabel
1022: . point - the point
1023: - value - The point value
1025: Level: intermediate
1027: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`
1028: @*/
1029: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1030: {
1031: PetscInt v;
1034: /* Find label value */
1035: DMLabelLookupStratum(label, value, &v);
1036: if (v < 0) return 0;
1038: if (label->bt) {
1040: PetscBTClear(label->bt, point - label->pStart);
1041: }
1043: /* Delete key */
1044: DMLabelMakeInvalid_Private(label, v);
1045: PetscHSetIDel(label->ht[v], point);
1046: return 0;
1047: }
1049: /*@
1050: DMLabelInsertIS - Set all points in the IS to a value
1052: Not collective
1054: Input Parameters:
1055: + label - the DMLabel
1056: . is - the point IS
1057: - value - The point value
1059: Level: intermediate
1061: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1062: @*/
1063: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1064: {
1065: PetscInt v, n, p;
1066: const PetscInt *points;
1070: /* Find label value, add new entry if needed */
1071: if (value == label->defaultValue) return 0;
1072: DMLabelLookupAddStratum(label, value, &v);
1073: /* Set keys */
1074: DMLabelMakeInvalid_Private(label, v);
1075: ISGetLocalSize(is, &n);
1076: ISGetIndices(is, &points);
1077: for (p = 0; p < n; ++p) PetscHSetIAdd(label->ht[v], points[p]);
1078: ISRestoreIndices(is, &points);
1079: return 0;
1080: }
1082: /*@
1083: DMLabelGetNumValues - Get the number of values that the DMLabel takes
1085: Not collective
1087: Input Parameter:
1088: . label - the DMLabel
1090: Output Parameter:
1091: . numValues - the number of values
1093: Level: intermediate
1095: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1096: @*/
1097: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1098: {
1101: *numValues = label->numStrata;
1102: return 0;
1103: }
1105: /*@
1106: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1108: Not collective
1110: Input Parameter:
1111: . label - the DMLabel
1113: Output Parameter:
1114: . is - the value IS
1116: Level: intermediate
1118: Notes:
1119: The output IS should be destroyed when no longer needed.
1120: Strata which are allocated but empty [DMLabelGetStratumSize() yields 0] are counted.
1121: If you need to count only nonempty strata, use DMLabelGetNonEmptyStratumValuesIS().
1123: .seealso: `DMLabelGetNonEmptyStratumValuesIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1124: @*/
1125: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1126: {
1129: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1130: return 0;
1131: }
1133: /*@
1134: DMLabelGetNonEmptyStratumValuesIS - Get an IS of all values that the DMlabel takes
1136: Not collective
1138: Input Parameter:
1139: . label - the DMLabel
1141: Output Parameter:
1142: . is - the value IS
1144: Level: intermediate
1146: Notes:
1147: The output IS should be destroyed when no longer needed.
1148: This is similar to DMLabelGetValueIS() but counts only nonempty strata.
1150: .seealso: `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1151: @*/
1152: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1153: {
1154: PetscInt i, j;
1155: PetscInt *valuesArr;
1159: PetscMalloc1(label->numStrata, &valuesArr);
1160: for (i = 0, j = 0; i < label->numStrata; i++) {
1161: PetscInt n;
1163: DMLabelGetStratumSize_Private(label, i, &n);
1164: if (n) valuesArr[j++] = label->stratumValues[i];
1165: }
1166: if (j == label->numStrata) {
1167: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1168: } else {
1169: ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values);
1170: }
1171: PetscFree(valuesArr);
1172: return 0;
1173: }
1175: /*@
1176: DMLabelGetValueIndex - Get the index of a given value in the list of values for the DMlabel, or -1 if it is not present
1178: Not collective
1180: Input Parameters:
1181: + label - the DMLabel
1182: - value - the value
1184: Output Parameter:
1185: . index - the index of value in the list of values
1187: Level: intermediate
1189: .seealso: `DMLabelGetValueIS()`, `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1190: @*/
1191: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1192: {
1193: PetscInt v;
1197: /* Do not assume they are sorted */
1198: for (v = 0; v < label->numStrata; ++v)
1199: if (label->stratumValues[v] == value) break;
1200: if (v >= label->numStrata) *index = -1;
1201: else *index = v;
1202: return 0;
1203: }
1205: /*@
1206: DMLabelHasStratum - Determine whether points exist with the given value
1208: Not collective
1210: Input Parameters:
1211: + label - the DMLabel
1212: - value - the stratum value
1214: Output Parameter:
1215: . exists - Flag saying whether points exist
1217: Level: intermediate
1219: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1220: @*/
1221: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1222: {
1223: PetscInt v;
1227: DMLabelLookupStratum(label, value, &v);
1228: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1229: return 0;
1230: }
1232: /*@
1233: DMLabelGetStratumSize - Get the size of a stratum
1235: Not collective
1237: Input Parameters:
1238: + label - the DMLabel
1239: - value - the stratum value
1241: Output Parameter:
1242: . size - The number of points in the stratum
1244: Level: intermediate
1246: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1247: @*/
1248: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1249: {
1250: PetscInt v;
1254: DMLabelLookupStratum(label, value, &v);
1255: DMLabelGetStratumSize_Private(label, v, size);
1256: return 0;
1257: }
1259: /*@
1260: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1262: Not collective
1264: Input Parameters:
1265: + label - the DMLabel
1266: - value - the stratum value
1268: Output Parameters:
1269: + start - the smallest point in the stratum
1270: - end - the largest point in the stratum
1272: Level: intermediate
1274: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1275: @*/
1276: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1277: {
1278: PetscInt v, min, max;
1281: if (start) {
1283: *start = -1;
1284: }
1285: if (end) {
1287: *end = -1;
1288: }
1289: DMLabelLookupStratum(label, value, &v);
1290: if (v < 0) return 0;
1291: DMLabelMakeValid_Private(label, v);
1292: if (label->stratumSizes[v] <= 0) return 0;
1293: ISGetMinMax(label->points[v], &min, &max);
1294: if (start) *start = min;
1295: if (end) *end = max + 1;
1296: return 0;
1297: }
1299: /*@
1300: DMLabelGetStratumIS - Get an IS with the stratum points
1302: Not collective
1304: Input Parameters:
1305: + label - the DMLabel
1306: - value - the stratum value
1308: Output Parameter:
1309: . points - The stratum points
1311: Level: intermediate
1313: Notes:
1314: The output IS should be destroyed when no longer needed.
1315: Returns NULL if the stratum is empty.
1317: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1318: @*/
1319: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1320: {
1321: PetscInt v;
1325: *points = NULL;
1326: DMLabelLookupStratum(label, value, &v);
1327: if (v < 0) return 0;
1328: DMLabelMakeValid_Private(label, v);
1329: PetscObjectReference((PetscObject)label->points[v]);
1330: *points = label->points[v];
1331: return 0;
1332: }
1334: /*@
1335: DMLabelSetStratumIS - Set the stratum points using an IS
1337: Not collective
1339: Input Parameters:
1340: + label - the DMLabel
1341: . value - the stratum value
1342: - points - The stratum points
1344: Level: intermediate
1346: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1347: @*/
1348: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1349: {
1350: PetscInt v;
1354: DMLabelLookupAddStratum(label, value, &v);
1355: if (is == label->points[v]) return 0;
1356: DMLabelClearStratum(label, value);
1357: ISGetLocalSize(is, &(label->stratumSizes[v]));
1358: PetscObjectReference((PetscObject)is);
1359: ISDestroy(&(label->points[v]));
1360: label->points[v] = is;
1361: label->validIS[v] = PETSC_TRUE;
1362: PetscObjectStateIncrease((PetscObject)label);
1363: if (label->bt) {
1364: const PetscInt *points;
1365: PetscInt p;
1367: ISGetIndices(is, &points);
1368: for (p = 0; p < label->stratumSizes[v]; ++p) {
1369: const PetscInt point = points[p];
1372: PetscBTSet(label->bt, point - label->pStart);
1373: }
1374: }
1375: return 0;
1376: }
1378: /*@
1379: DMLabelClearStratum - Remove a stratum
1381: Not collective
1383: Input Parameters:
1384: + label - the DMLabel
1385: - value - the stratum value
1387: Level: intermediate
1389: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1390: @*/
1391: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1392: {
1393: PetscInt v;
1396: DMLabelLookupStratum(label, value, &v);
1397: if (v < 0) return 0;
1398: if (label->validIS[v]) {
1399: if (label->bt) {
1400: PetscInt i;
1401: const PetscInt *points;
1403: ISGetIndices(label->points[v], &points);
1404: for (i = 0; i < label->stratumSizes[v]; ++i) {
1405: const PetscInt point = points[i];
1408: PetscBTClear(label->bt, point - label->pStart);
1409: }
1410: ISRestoreIndices(label->points[v], &points);
1411: }
1412: label->stratumSizes[v] = 0;
1413: ISDestroy(&label->points[v]);
1414: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1415: PetscObjectSetName((PetscObject)label->points[v], "indices");
1416: PetscObjectStateIncrease((PetscObject)label);
1417: } else {
1418: PetscHSetIClear(label->ht[v]);
1419: }
1420: return 0;
1421: }
1423: /*@
1424: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1426: Not collective
1428: Input Parameters:
1429: + label - The DMLabel
1430: . value - The label value for all points
1431: . pStart - The first point
1432: - pEnd - A point beyond all marked points
1434: Note: The marks points are [pStart, pEnd), and only the bounds are stored.
1436: Level: intermediate
1438: .seealso: `DMLabelCreate()`, `DMLabelSetStratumIS()`, `DMLabelGetStratumIS()`
1439: @*/
1440: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1441: {
1442: IS pIS;
1444: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1445: DMLabelSetStratumIS(label, value, pIS);
1446: ISDestroy(&pIS);
1447: return 0;
1448: }
1450: /*@
1451: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1453: Not collective
1455: Input Parameters:
1456: + label - The DMLabel
1457: . value - The label value
1458: - p - A point with this value
1460: Output Parameter:
1461: . 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
1463: Level: intermediate
1465: .seealso: `DMLabelGetValueIndex()`, `DMLabelGetStratumIS()`, `DMLabelCreate()`
1466: @*/
1467: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1468: {
1469: const PetscInt *indices;
1470: PetscInt v;
1474: *index = -1;
1475: DMLabelLookupStratum(label, value, &v);
1476: if (v < 0) return 0;
1477: DMLabelMakeValid_Private(label, v);
1478: ISGetIndices(label->points[v], &indices);
1479: PetscFindInt(p, label->stratumSizes[v], indices, index);
1480: ISRestoreIndices(label->points[v], &indices);
1481: return 0;
1482: }
1484: /*@
1485: DMLabelFilter - Remove all points outside of [start, end)
1487: Not collective
1489: Input Parameters:
1490: + label - the DMLabel
1491: . start - the first point kept
1492: - end - one more than the last point kept
1494: Level: intermediate
1496: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1497: @*/
1498: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1499: {
1500: PetscInt v;
1503: DMLabelDestroyIndex(label);
1504: DMLabelMakeAllValid_Private(label);
1505: for (v = 0; v < label->numStrata; ++v) ISGeneralFilter(label->points[v], start, end);
1506: DMLabelCreateIndex(label, start, end);
1507: return 0;
1508: }
1510: /*@
1511: DMLabelPermute - Create a new label with permuted points
1513: Not collective
1515: Input Parameters:
1516: + label - the DMLabel
1517: - permutation - the point permutation
1519: Output Parameter:
1520: . labelnew - the new label containing the permuted points
1522: Level: intermediate
1524: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1525: @*/
1526: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1527: {
1528: const PetscInt *perm;
1529: PetscInt numValues, numPoints, v, q;
1533: DMLabelMakeAllValid_Private(label);
1534: DMLabelDuplicate(label, labelNew);
1535: DMLabelGetNumValues(*labelNew, &numValues);
1536: ISGetLocalSize(permutation, &numPoints);
1537: ISGetIndices(permutation, &perm);
1538: for (v = 0; v < numValues; ++v) {
1539: const PetscInt size = (*labelNew)->stratumSizes[v];
1540: const PetscInt *points;
1541: PetscInt *pointsNew;
1543: ISGetIndices((*labelNew)->points[v], &points);
1544: PetscMalloc1(size, &pointsNew);
1545: for (q = 0; q < size; ++q) {
1546: const PetscInt point = points[q];
1549: pointsNew[q] = perm[point];
1550: }
1551: ISRestoreIndices((*labelNew)->points[v], &points);
1552: PetscSortInt(size, pointsNew);
1553: ISDestroy(&((*labelNew)->points[v]));
1554: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1555: ISCreateStride(PETSC_COMM_SELF, size, pointsNew[0], 1, &((*labelNew)->points[v]));
1556: PetscFree(pointsNew);
1557: } else {
1558: ISCreateGeneral(PETSC_COMM_SELF, size, pointsNew, PETSC_OWN_POINTER, &((*labelNew)->points[v]));
1559: }
1560: PetscObjectSetName((PetscObject)((*labelNew)->points[v]), "indices");
1561: }
1562: ISRestoreIndices(permutation, &perm);
1563: if (label->bt) {
1564: PetscBTDestroy(&label->bt);
1565: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1566: }
1567: return 0;
1568: }
1570: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1571: {
1572: MPI_Comm comm;
1573: PetscInt s, l, nroots, nleaves, offset, size;
1574: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1575: PetscSection rootSection;
1576: PetscSF labelSF;
1578: if (label) DMLabelMakeAllValid_Private(label);
1579: PetscObjectGetComm((PetscObject)sf, &comm);
1580: /* Build a section of stratum values per point, generate the according SF
1581: and distribute point-wise stratum values to leaves. */
1582: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1583: PetscSectionCreate(comm, &rootSection);
1584: PetscSectionSetChart(rootSection, 0, nroots);
1585: if (label) {
1586: for (s = 0; s < label->numStrata; ++s) {
1587: const PetscInt *points;
1589: ISGetIndices(label->points[s], &points);
1590: for (l = 0; l < label->stratumSizes[s]; l++) PetscSectionAddDof(rootSection, points[l], 1);
1591: ISRestoreIndices(label->points[s], &points);
1592: }
1593: }
1594: PetscSectionSetUp(rootSection);
1595: /* Create a point-wise array of stratum values */
1596: PetscSectionGetStorageSize(rootSection, &size);
1597: PetscMalloc1(size, &rootStrata);
1598: PetscCalloc1(nroots, &rootIdx);
1599: if (label) {
1600: for (s = 0; s < label->numStrata; ++s) {
1601: const PetscInt *points;
1603: ISGetIndices(label->points[s], &points);
1604: for (l = 0; l < label->stratumSizes[s]; l++) {
1605: const PetscInt p = points[l];
1606: PetscSectionGetOffset(rootSection, p, &offset);
1607: rootStrata[offset + rootIdx[p]++] = label->stratumValues[s];
1608: }
1609: ISRestoreIndices(label->points[s], &points);
1610: }
1611: }
1612: /* Build SF that maps label points to remote processes */
1613: PetscSectionCreate(comm, leafSection);
1614: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1615: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1616: PetscFree(remoteOffsets);
1617: /* Send the strata for each point over the derived SF */
1618: PetscSectionGetStorageSize(*leafSection, &size);
1619: PetscMalloc1(size, leafStrata);
1620: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE);
1621: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata, MPI_REPLACE);
1622: /* Clean up */
1623: PetscFree(rootStrata);
1624: PetscFree(rootIdx);
1625: PetscSectionDestroy(&rootSection);
1626: PetscSFDestroy(&labelSF);
1627: return 0;
1628: }
1630: /*@
1631: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1633: Collective on sf
1635: Input Parameters:
1636: + label - the DMLabel
1637: - sf - the map from old to new distribution
1639: Output Parameter:
1640: . labelnew - the new redistributed label
1642: Level: intermediate
1644: .seealso: `DMLabelCreate()`, `DMLabelGetValue()`, `DMLabelSetValue()`, `DMLabelClearValue()`
1645: @*/
1646: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1647: {
1648: MPI_Comm comm;
1649: PetscSection leafSection;
1650: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1651: PetscInt *leafStrata, *strataIdx;
1652: PetscInt **points;
1653: const char *lname = NULL;
1654: char *name;
1655: PetscInt nameSize;
1656: PetscHSetI stratumHash;
1657: size_t len = 0;
1658: PetscMPIInt rank;
1661: if (label) {
1663: DMLabelMakeAllValid_Private(label);
1664: }
1665: PetscObjectGetComm((PetscObject)sf, &comm);
1666: MPI_Comm_rank(comm, &rank);
1667: /* Bcast name */
1668: if (rank == 0) {
1669: PetscObjectGetName((PetscObject)label, &lname);
1670: PetscStrlen(lname, &len);
1671: }
1672: nameSize = len;
1673: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1674: PetscMalloc1(nameSize + 1, &name);
1675: if (rank == 0) PetscArraycpy(name, lname, nameSize + 1);
1676: MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm);
1677: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1678: PetscFree(name);
1679: /* Bcast defaultValue */
1680: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1681: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1682: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1683: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1684: /* Determine received stratum values and initialise new label*/
1685: PetscHSetICreate(&stratumHash);
1686: PetscSectionGetStorageSize(leafSection, &size);
1687: for (p = 0; p < size; ++p) PetscHSetIAdd(stratumHash, leafStrata[p]);
1688: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1689: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1690: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1691: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1692: /* Turn leafStrata into indices rather than stratum values */
1693: offset = 0;
1694: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1695: PetscSortInt((*labelNew)->numStrata, (*labelNew)->stratumValues);
1696: for (s = 0; s < (*labelNew)->numStrata; ++s) PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1697: for (p = 0; p < size; ++p) {
1698: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1699: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {
1700: leafStrata[p] = s;
1701: break;
1702: }
1703: }
1704: }
1705: /* Rebuild the point strata on the receiver */
1706: PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->stratumSizes);
1707: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1708: for (p = pStart; p < pEnd; p++) {
1709: PetscSectionGetDof(leafSection, p, &dof);
1710: PetscSectionGetOffset(leafSection, p, &offset);
1711: for (s = 0; s < dof; s++) (*labelNew)->stratumSizes[leafStrata[offset + s]]++;
1712: }
1713: PetscCalloc1((*labelNew)->numStrata, &(*labelNew)->ht);
1714: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->points);
1715: PetscMalloc1((*labelNew)->numStrata, &points);
1716: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1717: PetscHSetICreate(&(*labelNew)->ht[s]);
1718: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1719: }
1720: /* Insert points into new strata */
1721: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1722: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1723: for (p = pStart; p < pEnd; p++) {
1724: PetscSectionGetDof(leafSection, p, &dof);
1725: PetscSectionGetOffset(leafSection, p, &offset);
1726: for (s = 0; s < dof; s++) {
1727: stratum = leafStrata[offset + s];
1728: points[stratum][strataIdx[stratum]++] = p;
1729: }
1730: }
1731: for (s = 0; s < (*labelNew)->numStrata; s++) {
1732: ISCreateGeneral(PETSC_COMM_SELF, (*labelNew)->stratumSizes[s], &(points[s][0]), PETSC_OWN_POINTER, &((*labelNew)->points[s]));
1733: PetscObjectSetName((PetscObject)((*labelNew)->points[s]), "indices");
1734: }
1735: PetscFree(points);
1736: PetscHSetIDestroy(&stratumHash);
1737: PetscFree(leafStrata);
1738: PetscFree(strataIdx);
1739: PetscSectionDestroy(&leafSection);
1740: return 0;
1741: }
1743: /*@
1744: DMLabelGather - Gather all label values from leafs into roots
1746: Collective on sf
1748: Input Parameters:
1749: + label - the DMLabel
1750: - sf - the Star Forest point communication map
1752: Output Parameters:
1753: . labelNew - the new DMLabel with localised leaf values
1755: Level: developer
1757: Note: This is the inverse operation to DMLabelDistribute.
1759: .seealso: `DMLabelDistribute()`
1760: @*/
1761: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1762: {
1763: MPI_Comm comm;
1764: PetscSection rootSection;
1765: PetscSF sfLabel;
1766: PetscSFNode *rootPoints, *leafPoints;
1767: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1768: const PetscInt *rootDegree, *ilocal;
1769: PetscInt *rootStrata;
1770: const char *lname;
1771: char *name;
1772: PetscInt nameSize;
1773: size_t len = 0;
1774: PetscMPIInt rank, size;
1778: PetscObjectGetComm((PetscObject)sf, &comm);
1779: MPI_Comm_rank(comm, &rank);
1780: MPI_Comm_size(comm, &size);
1781: /* Bcast name */
1782: if (rank == 0) {
1783: PetscObjectGetName((PetscObject)label, &lname);
1784: PetscStrlen(lname, &len);
1785: }
1786: nameSize = len;
1787: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1788: PetscMalloc1(nameSize + 1, &name);
1789: if (rank == 0) PetscArraycpy(name, lname, nameSize + 1);
1790: MPI_Bcast(name, nameSize + 1, MPI_CHAR, 0, comm);
1791: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1792: PetscFree(name);
1793: /* Gather rank/index pairs of leaves into local roots to build
1794: an inverse, multi-rooted SF. Note that this ignores local leaf
1795: indexing due to the use of the multiSF in PetscSFGather. */
1796: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1797: PetscMalloc1(nroots, &leafPoints);
1798: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1799: for (p = 0; p < nleaves; p++) {
1800: PetscInt ilp = ilocal ? ilocal[p] : p;
1802: leafPoints[ilp].index = ilp;
1803: leafPoints[ilp].rank = rank;
1804: }
1805: PetscSFComputeDegreeBegin(sf, &rootDegree);
1806: PetscSFComputeDegreeEnd(sf, &rootDegree);
1807: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1808: PetscMalloc1(nmultiroots, &rootPoints);
1809: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1810: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1811: PetscSFCreate(comm, &sfLabel);
1812: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1813: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1814: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1815: /* Rebuild the point strata on the receiver */
1816: for (p = 0, idx = 0; p < nroots; p++) {
1817: for (d = 0; d < rootDegree[p]; d++) {
1818: PetscSectionGetDof(rootSection, idx + d, &dof);
1819: PetscSectionGetOffset(rootSection, idx + d, &offset);
1820: for (s = 0; s < dof; s++) DMLabelSetValue(*labelNew, p, rootStrata[offset + s]);
1821: }
1822: idx += rootDegree[p];
1823: }
1824: PetscFree(leafPoints);
1825: PetscFree(rootStrata);
1826: PetscSectionDestroy(&rootSection);
1827: PetscSFDestroy(&sfLabel);
1828: return 0;
1829: }
1831: static PetscErrorCode DMLabelPropagateInit_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[])
1832: {
1833: const PetscInt *degree;
1834: const PetscInt *points;
1835: PetscInt Nr, r, Nl, l, val, defVal;
1837: DMLabelGetDefaultValue(label, &defVal);
1838: /* Add in leaves */
1839: PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL);
1840: for (l = 0; l < Nl; ++l) {
1841: DMLabelGetValue(label, points[l], &val);
1842: if (val != defVal) valArray[points[l]] = val;
1843: }
1844: /* Add in shared roots */
1845: PetscSFComputeDegreeBegin(pointSF, °ree);
1846: PetscSFComputeDegreeEnd(pointSF, °ree);
1847: for (r = 0; r < Nr; ++r) {
1848: if (degree[r]) {
1849: DMLabelGetValue(label, r, &val);
1850: if (val != defVal) valArray[r] = val;
1851: }
1852: }
1853: return 0;
1854: }
1856: static PetscErrorCode DMLabelPropagateFini_Internal(DMLabel label, PetscSF pointSF, PetscInt valArray[], PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1857: {
1858: const PetscInt *degree;
1859: const PetscInt *points;
1860: PetscInt Nr, r, Nl, l, val, defVal;
1862: DMLabelGetDefaultValue(label, &defVal);
1863: /* Read out leaves */
1864: PetscSFGetGraph(pointSF, &Nr, &Nl, &points, NULL);
1865: for (l = 0; l < Nl; ++l) {
1866: const PetscInt p = points[l];
1867: const PetscInt cval = valArray[p];
1869: if (cval != defVal) {
1870: DMLabelGetValue(label, p, &val);
1871: if (val == defVal) {
1872: DMLabelSetValue(label, p, cval);
1873: if (markPoint) (*markPoint)(label, p, cval, ctx);
1874: }
1875: }
1876: }
1877: /* Read out shared roots */
1878: PetscSFComputeDegreeBegin(pointSF, °ree);
1879: PetscSFComputeDegreeEnd(pointSF, °ree);
1880: for (r = 0; r < Nr; ++r) {
1881: if (degree[r]) {
1882: const PetscInt cval = valArray[r];
1884: if (cval != defVal) {
1885: DMLabelGetValue(label, r, &val);
1886: if (val == defVal) {
1887: DMLabelSetValue(label, r, cval);
1888: if (markPoint) (*markPoint)(label, r, cval, ctx);
1889: }
1890: }
1891: }
1892: }
1893: return 0;
1894: }
1896: /*@
1897: DMLabelPropagateBegin - Setup a cycle of label propagation
1899: Collective on sf
1901: Input Parameters:
1902: + label - The DMLabel to propagate across processes
1903: - sf - The SF describing parallel layout of the label points
1905: Level: intermediate
1907: .seealso: `DMLabelPropagateEnd()`, `DMLabelPropagatePush()`
1908: @*/
1909: PetscErrorCode DMLabelPropagateBegin(DMLabel label, PetscSF sf)
1910: {
1911: PetscInt Nr, r, defVal;
1912: PetscMPIInt size;
1914: MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size);
1915: if (size > 1) {
1916: DMLabelGetDefaultValue(label, &defVal);
1917: PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL);
1918: if (Nr >= 0) PetscMalloc1(Nr, &label->propArray);
1919: for (r = 0; r < Nr; ++r) label->propArray[r] = defVal;
1920: }
1921: return 0;
1922: }
1924: /*@
1925: DMLabelPropagateEnd - Tear down a cycle of label propagation
1927: Collective on sf
1929: Input Parameters:
1930: + label - The DMLabel to propagate across processes
1931: - sf - The SF describing parallel layout of the label points
1933: Level: intermediate
1935: .seealso: `DMLabelPropagateBegin()`, `DMLabelPropagatePush()`
1936: @*/
1937: PetscErrorCode DMLabelPropagateEnd(DMLabel label, PetscSF pointSF)
1938: {
1939: PetscFree(label->propArray);
1940: label->propArray = NULL;
1941: return 0;
1942: }
1944: /*@C
1945: DMLabelPropagatePush - Tear down a cycle of label propagation
1947: Collective on sf
1949: Input Parameters:
1950: + label - The DMLabel to propagate across processes
1951: . sf - The SF describing parallel layout of the label points
1952: . markPoint - An optional callback that is called when a point is marked, or NULL
1953: - ctx - An optional user context for the callback, or NULL
1955: Calling sequence of markPoint:
1956: $ markPoint(DMLabel label, PetscInt p, PetscInt val, void *ctx);
1958: + label - The DMLabel
1959: . p - The point being marked
1960: . val - The label value for p
1961: - ctx - An optional user context
1963: Level: intermediate
1965: .seealso: `DMLabelPropagateBegin()`, `DMLabelPropagateEnd()`
1966: @*/
1967: PetscErrorCode DMLabelPropagatePush(DMLabel label, PetscSF pointSF, PetscErrorCode (*markPoint)(DMLabel, PetscInt, PetscInt, void *), void *ctx)
1968: {
1969: PetscInt *valArray = label->propArray, Nr;
1970: PetscMPIInt size;
1972: MPI_Comm_size(PetscObjectComm((PetscObject)pointSF), &size);
1973: PetscSFGetGraph(pointSF, &Nr, NULL, NULL, NULL);
1974: if (size > 1 && Nr >= 0) {
1975: /* Communicate marked edges
1976: The current implementation allocates an array the size of the number of root. We put the label values into the
1977: array, and then call PetscSFReduce()+PetscSFBcast() to make the marks consistent.
1979: TODO: We could use in-place communication with a different SF
1980: We use MPI_SUM for the Reduce, and check the result against the rootdegree. If sum >= rootdegree+1, then the edge has
1981: already been marked. If not, it might have been handled on the process in this round, but we add it anyway.
1983: In order to update the queue with the new edges from the label communication, we use BcastAnOp(MPI_SUM), so that new
1984: 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
1985: edge to the queue.
1986: */
1987: DMLabelPropagateInit_Internal(label, pointSF, valArray);
1988: PetscSFReduceBegin(pointSF, MPIU_INT, valArray, valArray, MPI_MAX);
1989: PetscSFReduceEnd(pointSF, MPIU_INT, valArray, valArray, MPI_MAX);
1990: PetscSFBcastBegin(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE);
1991: PetscSFBcastEnd(pointSF, MPIU_INT, valArray, valArray, MPI_REPLACE);
1992: DMLabelPropagateFini_Internal(label, pointSF, valArray, markPoint, ctx);
1993: }
1994: return 0;
1995: }
1997: /*@
1998: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
2000: Not collective
2002: Input Parameter:
2003: . label - the DMLabel
2005: Output Parameters:
2006: + section - the section giving offsets for each stratum
2007: - is - An IS containing all the label points
2009: Level: developer
2011: .seealso: `DMLabelDistribute()`
2012: @*/
2013: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
2014: {
2015: IS vIS;
2016: const PetscInt *values;
2017: PetscInt *points;
2018: PetscInt nV, vS = 0, vE = 0, v, N;
2021: DMLabelGetNumValues(label, &nV);
2022: DMLabelGetValueIS(label, &vIS);
2023: ISGetIndices(vIS, &values);
2024: if (nV) {
2025: vS = values[0];
2026: vE = values[0] + 1;
2027: }
2028: for (v = 1; v < nV; ++v) {
2029: vS = PetscMin(vS, values[v]);
2030: vE = PetscMax(vE, values[v] + 1);
2031: }
2032: PetscSectionCreate(PETSC_COMM_SELF, section);
2033: PetscSectionSetChart(*section, vS, vE);
2034: for (v = 0; v < nV; ++v) {
2035: PetscInt n;
2037: DMLabelGetStratumSize(label, values[v], &n);
2038: PetscSectionSetDof(*section, values[v], n);
2039: }
2040: PetscSectionSetUp(*section);
2041: PetscSectionGetStorageSize(*section, &N);
2042: PetscMalloc1(N, &points);
2043: for (v = 0; v < nV; ++v) {
2044: IS is;
2045: const PetscInt *spoints;
2046: PetscInt dof, off, p;
2048: PetscSectionGetDof(*section, values[v], &dof);
2049: PetscSectionGetOffset(*section, values[v], &off);
2050: DMLabelGetStratumIS(label, values[v], &is);
2051: ISGetIndices(is, &spoints);
2052: for (p = 0; p < dof; ++p) points[off + p] = spoints[p];
2053: ISRestoreIndices(is, &spoints);
2054: ISDestroy(&is);
2055: }
2056: ISRestoreIndices(vIS, &values);
2057: ISDestroy(&vIS);
2058: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
2059: return 0;
2060: }
2062: /*@
2063: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
2064: the local section and an SF describing the section point overlap.
2066: Collective on sf
2068: Input Parameters:
2069: + s - The PetscSection for the local field layout
2070: . sf - The SF describing parallel layout of the section points
2071: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
2072: . label - The label specifying the points
2073: - labelValue - The label stratum specifying the points
2075: Output Parameter:
2076: . gsection - The PetscSection for the global field layout
2078: Note: This gives negative sizes and offsets to points not owned by this process
2080: Level: developer
2082: .seealso: `PetscSectionCreate()`
2083: @*/
2084: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
2085: {
2086: PetscInt *neg = NULL, *tmpOff = NULL;
2087: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
2092: PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection);
2093: PetscSectionGetChart(s, &pStart, &pEnd);
2094: PetscSectionSetChart(*gsection, pStart, pEnd);
2095: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
2096: if (nroots >= 0) {
2098: PetscCalloc1(nroots, &neg);
2099: if (nroots > pEnd - pStart) {
2100: PetscCalloc1(nroots, &tmpOff);
2101: } else {
2102: tmpOff = &(*gsection)->atlasDof[-pStart];
2103: }
2104: }
2105: /* Mark ghost points with negative dof */
2106: for (p = pStart; p < pEnd; ++p) {
2107: PetscInt value;
2109: DMLabelGetValue(label, p, &value);
2110: if (value != labelValue) continue;
2111: PetscSectionGetDof(s, p, &dof);
2112: PetscSectionSetDof(*gsection, p, dof);
2113: PetscSectionGetConstraintDof(s, p, &cdof);
2114: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
2115: if (neg) neg[p] = -(dof + 1);
2116: }
2117: PetscSectionSetUpBC(*gsection);
2118: if (nroots >= 0) {
2119: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2120: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2121: if (nroots > pEnd - pStart) {
2122: for (p = pStart; p < pEnd; ++p) {
2123: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
2124: }
2125: }
2126: }
2127: /* Calculate new sizes, get process offset, and calculate point offsets */
2128: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2129: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
2130: (*gsection)->atlasOff[p] = off;
2131: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p] - cdof : 0;
2132: }
2133: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
2134: globalOff -= off;
2135: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
2136: (*gsection)->atlasOff[p] += globalOff;
2137: if (neg) neg[p] = -((*gsection)->atlasOff[p] + 1);
2138: }
2139: /* Put in negative offsets for ghost points */
2140: if (nroots >= 0) {
2141: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2142: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
2143: if (nroots > pEnd - pStart) {
2144: for (p = pStart; p < pEnd; ++p) {
2145: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
2146: }
2147: }
2148: }
2149: if (nroots >= 0 && nroots > pEnd - pStart) PetscFree(tmpOff);
2150: PetscFree(neg);
2151: return 0;
2152: }
2154: typedef struct _n_PetscSectionSym_Label {
2155: DMLabel label;
2156: PetscCopyMode *modes;
2157: PetscInt *sizes;
2158: const PetscInt ***perms;
2159: const PetscScalar ***rots;
2160: PetscInt (*minMaxOrients)[2];
2161: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
2162: } PetscSectionSym_Label;
2164: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
2165: {
2166: PetscInt i, j;
2167: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2169: for (i = 0; i <= sl->numStrata; i++) {
2170: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
2171: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2172: if (sl->perms[i]) PetscFree(sl->perms[i][j]);
2173: if (sl->rots[i]) PetscFree(sl->rots[i][j]);
2174: }
2175: if (sl->perms[i]) {
2176: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
2178: PetscFree(perms);
2179: }
2180: if (sl->rots[i]) {
2181: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2183: PetscFree(rots);
2184: }
2185: }
2186: }
2187: PetscFree5(sl->modes, sl->sizes, sl->perms, sl->rots, sl->minMaxOrients);
2188: DMLabelDestroy(&sl->label);
2189: sl->numStrata = 0;
2190: return 0;
2191: }
2193: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2194: {
2195: PetscSectionSymLabelReset(sym);
2196: PetscFree(sym->data);
2197: return 0;
2198: }
2200: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2201: {
2202: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2203: PetscBool isAscii;
2204: DMLabel label = sl->label;
2205: const char *name;
2207: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii);
2208: if (isAscii) {
2209: PetscInt i, j, k;
2210: PetscViewerFormat format;
2212: PetscViewerGetFormat(viewer, &format);
2213: if (label) {
2214: PetscViewerGetFormat(viewer, &format);
2215: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2216: PetscViewerASCIIPushTab(viewer);
2217: DMLabelView(label, viewer);
2218: PetscViewerASCIIPopTab(viewer);
2219: } else {
2220: PetscObjectGetName((PetscObject)sl->label, &name);
2221: PetscViewerASCIIPrintf(viewer, " Label '%s'\n", name);
2222: }
2223: } else {
2224: PetscViewerASCIIPrintf(viewer, "No label given\n");
2225: }
2226: PetscViewerASCIIPushTab(viewer);
2227: for (i = 0; i <= sl->numStrata; i++) {
2228: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2230: if (!(sl->perms[i] || sl->rots[i])) {
2231: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point): no symmetries\n", value, sl->sizes[i]);
2232: } else {
2233: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %" PetscInt_FMT " (%" PetscInt_FMT " dofs per point):\n", value, sl->sizes[i]);
2234: PetscViewerASCIIPushTab(viewer);
2235: PetscViewerASCIIPrintf(viewer, "Orientation range: [%" PetscInt_FMT ", %" PetscInt_FMT ")\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
2236: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2237: PetscViewerASCIIPushTab(viewer);
2238: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2239: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2240: PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ": identity\n", j);
2241: } else {
2242: PetscInt tab;
2244: PetscViewerASCIIPrintf(viewer, "Orientation %" PetscInt_FMT ":\n", j);
2245: PetscViewerASCIIPushTab(viewer);
2246: PetscViewerASCIIGetTab(viewer, &tab);
2247: if (sl->perms[i] && sl->perms[i][j]) {
2248: PetscViewerASCIIPrintf(viewer, "Permutation:");
2249: PetscViewerASCIISetTab(viewer, 0);
2250: for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sl->perms[i][j][k]);
2251: PetscViewerASCIIPrintf(viewer, "\n");
2252: PetscViewerASCIISetTab(viewer, tab);
2253: }
2254: if (sl->rots[i] && sl->rots[i][j]) {
2255: PetscViewerASCIIPrintf(viewer, "Rotations: ");
2256: PetscViewerASCIISetTab(viewer, 0);
2257: #if defined(PETSC_USE_COMPLEX)
2258: for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer, " %+g+i*%+g", (double)PetscRealPart(sl->rots[i][j][k]), (double)PetscImaginaryPart(sl->rots[i][j][k]));
2259: #else
2260: for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer, " %+g", (double)sl->rots[i][j][k]);
2261: #endif
2262: PetscViewerASCIIPrintf(viewer, "\n");
2263: PetscViewerASCIISetTab(viewer, tab);
2264: }
2265: PetscViewerASCIIPopTab(viewer);
2266: }
2267: }
2268: PetscViewerASCIIPopTab(viewer);
2269: }
2270: PetscViewerASCIIPopTab(viewer);
2271: }
2272: }
2273: PetscViewerASCIIPopTab(viewer);
2274: }
2275: return 0;
2276: }
2278: /*@
2279: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2281: Logically collective on sym
2283: Input parameters:
2284: + sym - the section symmetries
2285: - label - the DMLabel describing the types of points
2287: Level: developer:
2289: .seealso: `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreateLabel()`, `PetscSectionGetPointSyms()`
2290: @*/
2291: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2292: {
2293: PetscSectionSym_Label *sl;
2296: sl = (PetscSectionSym_Label *)sym->data;
2297: if (sl->label && sl->label != label) PetscSectionSymLabelReset(sym);
2298: if (label) {
2299: sl->label = label;
2300: PetscObjectReference((PetscObject)label);
2301: DMLabelGetNumValues(label, &sl->numStrata);
2302: 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);
2303: PetscMemzero((void *)sl->modes, (sl->numStrata + 1) * sizeof(PetscCopyMode));
2304: PetscMemzero((void *)sl->sizes, (sl->numStrata + 1) * sizeof(PetscInt));
2305: PetscMemzero((void *)sl->perms, (sl->numStrata + 1) * sizeof(const PetscInt **));
2306: PetscMemzero((void *)sl->rots, (sl->numStrata + 1) * sizeof(const PetscScalar **));
2307: PetscMemzero((void *)sl->minMaxOrients, (sl->numStrata + 1) * sizeof(PetscInt[2]));
2308: }
2309: return 0;
2310: }
2312: /*@C
2313: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2315: Logically collective on sym
2317: Input Parameters:
2318: + sym - the section symmetries
2319: - stratum - the stratum value in the label that we are assigning symmetries for
2321: Output Parameters:
2322: + size - the number of dofs for points in the stratum of the label
2323: . minOrient - the smallest orientation for a point in this stratum
2324: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2325: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2326: - 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
2328: Level: developer
2330: .seealso: `PetscSectionSymLabelSetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2331: @*/
2332: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2333: {
2334: PetscSectionSym_Label *sl;
2335: const char *name;
2336: PetscInt i;
2339: sl = (PetscSectionSym_Label *)sym->data;
2341: for (i = 0; i <= sl->numStrata; i++) {
2342: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2344: if (stratum == value) break;
2345: }
2346: PetscObjectGetName((PetscObject)sl->label, &name);
2348: if (size) {
2350: *size = sl->sizes[i];
2351: }
2352: if (minOrient) {
2354: *minOrient = sl->minMaxOrients[i][0];
2355: }
2356: if (maxOrient) {
2358: *maxOrient = sl->minMaxOrients[i][1];
2359: }
2360: if (perms) {
2362: *perms = sl->perms[i] ? &sl->perms[i][sl->minMaxOrients[i][0]] : NULL;
2363: }
2364: if (rots) {
2366: *rots = sl->rots[i] ? &sl->rots[i][sl->minMaxOrients[i][0]] : NULL;
2367: }
2368: return 0;
2369: }
2371: /*@C
2372: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2374: Logically collective on sym
2376: InputParameters:
2377: + sym - the section symmetries
2378: . stratum - the stratum value in the label that we are assigning symmetries for
2379: . size - the number of dofs for points in the stratum of the label
2380: . minOrient - the smallest orientation for a point in this stratum
2381: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2382: . mode - how sym should copy the perms and rots arrays
2383: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2384: - 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
2386: Level: developer
2388: .seealso: `PetscSectionSymLabelGetStratum()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreateLabel()`
2389: @*/
2390: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2391: {
2392: PetscSectionSym_Label *sl;
2393: const char *name;
2394: PetscInt i, j, k;
2397: sl = (PetscSectionSym_Label *)sym->data;
2399: for (i = 0; i <= sl->numStrata; i++) {
2400: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2402: if (stratum == value) break;
2403: }
2404: PetscObjectGetName((PetscObject)sl->label, &name);
2406: sl->sizes[i] = size;
2407: sl->modes[i] = mode;
2408: sl->minMaxOrients[i][0] = minOrient;
2409: sl->minMaxOrients[i][1] = maxOrient;
2410: if (mode == PETSC_COPY_VALUES) {
2411: if (perms) {
2412: PetscInt **ownPerms;
2414: PetscCalloc1(maxOrient - minOrient, &ownPerms);
2415: for (j = 0; j < maxOrient - minOrient; j++) {
2416: if (perms[j]) {
2417: PetscMalloc1(size, &ownPerms[j]);
2418: for (k = 0; k < size; k++) ownPerms[j][k] = perms[j][k];
2419: }
2420: }
2421: sl->perms[i] = (const PetscInt **)&ownPerms[-minOrient];
2422: }
2423: if (rots) {
2424: PetscScalar **ownRots;
2426: PetscCalloc1(maxOrient - minOrient, &ownRots);
2427: for (j = 0; j < maxOrient - minOrient; j++) {
2428: if (rots[j]) {
2429: PetscMalloc1(size, &ownRots[j]);
2430: for (k = 0; k < size; k++) ownRots[j][k] = rots[j][k];
2431: }
2432: }
2433: sl->rots[i] = (const PetscScalar **)&ownRots[-minOrient];
2434: }
2435: } else {
2436: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2437: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2438: }
2439: return 0;
2440: }
2442: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2443: {
2444: PetscInt i, j, numStrata;
2445: PetscSectionSym_Label *sl;
2446: DMLabel label;
2448: sl = (PetscSectionSym_Label *)sym->data;
2449: numStrata = sl->numStrata;
2450: label = sl->label;
2451: for (i = 0; i < numPoints; i++) {
2452: PetscInt point = points[2 * i];
2453: PetscInt ornt = points[2 * i + 1];
2455: for (j = 0; j < numStrata; j++) {
2456: if (label->validIS[j]) {
2457: PetscInt k;
2459: ISLocate(label->points[j], point, &k);
2460: if (k >= 0) break;
2461: } else {
2462: PetscBool has;
2464: PetscHSetIHas(label->ht[j], point, &has);
2465: if (has) break;
2466: }
2467: }
2469: j < numStrata ? label->stratumValues[j] : label->defaultValue);
2470: if (perms) perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;
2471: if (rots) rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;
2472: }
2473: return 0;
2474: }
2476: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2477: {
2478: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)nsym->data;
2479: IS valIS;
2480: const PetscInt *values;
2481: PetscInt Nv, v;
2483: DMLabelGetNumValues(sl->label, &Nv);
2484: DMLabelGetValueIS(sl->label, &valIS);
2485: ISGetIndices(valIS, &values);
2486: for (v = 0; v < Nv; ++v) {
2487: const PetscInt val = values[v];
2488: PetscInt size, minOrient, maxOrient;
2489: const PetscInt **perms;
2490: const PetscScalar **rots;
2492: PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots);
2493: PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots);
2494: }
2495: ISDestroy(&valIS);
2496: return 0;
2497: }
2499: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2500: {
2501: PetscSectionSym_Label *sl = (PetscSectionSym_Label *)sym->data;
2502: DMLabel dlabel;
2504: DMLabelDistribute(sl->label, migrationSF, &dlabel);
2505: PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)sym), dlabel, dsym);
2506: DMLabelDestroy(&dlabel);
2507: PetscSectionSymCopy(sym, *dsym);
2508: return 0;
2509: }
2511: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2512: {
2513: PetscSectionSym_Label *sl;
2515: PetscNew(&sl);
2516: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2517: sym->ops->distribute = PetscSectionSymDistribute_Label;
2518: sym->ops->copy = PetscSectionSymCopy_Label;
2519: sym->ops->view = PetscSectionSymView_Label;
2520: sym->ops->destroy = PetscSectionSymDestroy_Label;
2521: sym->data = (void *)sl;
2522: return 0;
2523: }
2525: /*@
2526: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2528: Collective
2530: Input Parameters:
2531: + comm - the MPI communicator for the new symmetry
2532: - label - the label defining the strata
2534: Output Parameters:
2535: . sym - the section symmetries
2537: Level: developer
2539: .seealso: `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
2540: @*/
2541: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2542: {
2543: DMInitializePackage();
2544: PetscSectionSymCreate(comm, sym);
2545: PetscSectionSymSetType(*sym, PETSCSECTIONSYMLABEL);
2546: PetscSectionSymLabelSetLabel(*sym, label);
2547: return 0;
2548: }