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++) {
117: DMLabelMakeValid_Private(label, v);
118: }
119: return 0;
120: }
122: /*
123: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
125: Not collective
127: Input parameter:
128: + label - The DMLabel
129: - v - The stratum value
131: Output parameter:
132: . label - The DMLabel with stratum in hash format
134: Level: developer
136: .seealso: DMLabelCreate()
137: */
138: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
139: {
140: PetscInt p;
141: const PetscInt *points;
143: if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
145: if (label->points[v]) {
146: ISGetIndices(label->points[v], &points);
147: for (p = 0; p < label->stratumSizes[v]; ++p) {
148: PetscHSetIAdd(label->ht[v], points[p]);
149: }
150: ISRestoreIndices(label->points[v],&points);
151: ISDestroy(&(label->points[v]));
152: }
153: label->validIS[v] = PETSC_FALSE;
154: return 0;
155: }
157: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
158: #define DMLABEL_LOOKUP_THRESHOLD 16
159: #endif
161: static inline PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
162: {
163: PetscInt v;
165: *index = -1;
166: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
167: for (v = 0; v < label->numStrata; ++v)
168: if (label->stratumValues[v] == value) {*index = v; break;}
169: } else {
170: PetscHMapIGet(label->hmap, value, index);
171: }
172: if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
173: PetscInt len, loc = -1;
174: PetscHMapIGetSize(label->hmap, &len);
176: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
177: PetscHMapIGet(label->hmap, value, &loc);
178: } else {
179: for (v = 0; v < label->numStrata; ++v)
180: if (label->stratumValues[v] == value) {loc = v; break;}
181: }
183: }
184: return 0;
185: }
187: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
188: {
189: PetscInt v;
190: PetscInt *tmpV;
191: PetscInt *tmpS;
192: PetscHSetI *tmpH, ht;
193: IS *tmpP, is;
194: PetscBool *tmpB;
195: PetscHMapI hmap = label->hmap;
197: v = label->numStrata;
198: tmpV = label->stratumValues;
199: tmpS = label->stratumSizes;
200: tmpH = label->ht;
201: tmpP = label->points;
202: tmpB = label->validIS;
203: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
204: PetscInt *oldV = tmpV;
205: PetscInt *oldS = tmpS;
206: PetscHSetI *oldH = tmpH;
207: IS *oldP = tmpP;
208: PetscBool *oldB = tmpB;
209: PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
210: PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
211: PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
212: PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
213: PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
214: PetscArraycpy(tmpV, oldV, v);
215: PetscArraycpy(tmpS, oldS, v);
216: PetscArraycpy(tmpH, oldH, v);
217: PetscArraycpy(tmpP, oldP, v);
218: PetscArraycpy(tmpB, oldB, v);
219: PetscFree(oldV);
220: PetscFree(oldS);
221: PetscFree(oldH);
222: PetscFree(oldP);
223: PetscFree(oldB);
224: }
225: label->numStrata = v+1;
226: label->stratumValues = tmpV;
227: label->stratumSizes = tmpS;
228: label->ht = tmpH;
229: label->points = tmpP;
230: label->validIS = tmpB;
231: PetscHSetICreate(&ht);
232: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
233: PetscHMapISet(hmap, value, v);
234: tmpV[v] = value;
235: tmpS[v] = 0;
236: tmpH[v] = ht;
237: tmpP[v] = is;
238: tmpB[v] = PETSC_TRUE;
239: PetscObjectStateIncrease((PetscObject) label);
240: *index = v;
241: return 0;
242: }
244: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
245: {
246: DMLabelLookupStratum(label, value, index);
247: if (*index < 0) DMLabelNewStratum(label, value, index);
248: return 0;
249: }
251: static inline PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
252: {
253: *size = 0;
254: if (v < 0) return 0;
255: if (label->validIS[v]) {
256: *size = label->stratumSizes[v];
257: } else {
258: PetscHSetIGetSize(label->ht[v], size);
259: }
260: return 0;
261: }
263: /*@
264: DMLabelAddStratum - Adds a new stratum value in a DMLabel
266: Input Parameters:
267: + label - The DMLabel
268: - value - The stratum value
270: Level: beginner
272: .seealso: DMLabelCreate(), DMLabelDestroy()
273: @*/
274: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
275: {
276: PetscInt v;
279: DMLabelLookupAddStratum(label, value, &v);
280: return 0;
281: }
283: /*@
284: DMLabelAddStrata - Adds new stratum values in a DMLabel
286: Not collective
288: Input Parameters:
289: + label - The DMLabel
290: . numStrata - The number of stratum values
291: - stratumValues - The stratum values
293: Level: beginner
295: .seealso: DMLabelCreate(), DMLabelDestroy()
296: @*/
297: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
298: {
299: PetscInt *values, v;
303: PetscMalloc1(numStrata, &values);
304: PetscArraycpy(values, stratumValues, numStrata);
305: PetscSortRemoveDupsInt(&numStrata, values);
306: if (!label->numStrata) { /* Fast preallocation */
307: PetscInt *tmpV;
308: PetscInt *tmpS;
309: PetscHSetI *tmpH, ht;
310: IS *tmpP, is;
311: PetscBool *tmpB;
312: PetscHMapI hmap = label->hmap;
314: PetscMalloc1(numStrata, &tmpV);
315: PetscMalloc1(numStrata, &tmpS);
316: PetscMalloc1(numStrata, &tmpH);
317: PetscMalloc1(numStrata, &tmpP);
318: PetscMalloc1(numStrata, &tmpB);
319: label->numStrata = numStrata;
320: label->stratumValues = tmpV;
321: label->stratumSizes = tmpS;
322: label->ht = tmpH;
323: label->points = tmpP;
324: label->validIS = tmpB;
325: for (v = 0; v < numStrata; ++v) {
326: PetscHSetICreate(&ht);
327: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
328: PetscHMapISet(hmap, values[v], v);
329: tmpV[v] = values[v];
330: tmpS[v] = 0;
331: tmpH[v] = ht;
332: tmpP[v] = is;
333: tmpB[v] = PETSC_TRUE;
334: }
335: PetscObjectStateIncrease((PetscObject) label);
336: } else {
337: for (v = 0; v < numStrata; ++v) {
338: DMLabelAddStratum(label, values[v]);
339: }
340: }
341: PetscFree(values);
342: return 0;
343: }
345: /*@
346: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
348: Not collective
350: Input Parameters:
351: + label - The DMLabel
352: - valueIS - Index set with stratum values
354: Level: beginner
356: .seealso: DMLabelCreate(), DMLabelDestroy()
357: @*/
358: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
359: {
360: PetscInt numStrata;
361: const PetscInt *stratumValues;
365: ISGetLocalSize(valueIS, &numStrata);
366: ISGetIndices(valueIS, &stratumValues);
367: DMLabelAddStrata(label, numStrata, stratumValues);
368: return 0;
369: }
371: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
372: {
373: PetscInt v;
374: PetscMPIInt rank;
376: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
377: PetscViewerASCIIPushSynchronized(viewer);
378: if (label) {
379: const char *name;
381: PetscObjectGetName((PetscObject) label, &name);
382: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
383: if (label->bt) PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);
384: for (v = 0; v < label->numStrata; ++v) {
385: const PetscInt value = label->stratumValues[v];
386: const PetscInt *points;
387: PetscInt p;
389: ISGetIndices(label->points[v], &points);
390: for (p = 0; p < label->stratumSizes[v]; ++p) {
391: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
392: }
393: ISRestoreIndices(label->points[v],&points);
394: }
395: }
396: PetscViewerFlush(viewer);
397: PetscViewerASCIIPopSynchronized(viewer);
398: return 0;
399: }
401: /*@C
402: DMLabelView - View the label
404: Collective on viewer
406: Input Parameters:
407: + label - The DMLabel
408: - viewer - The PetscViewer
410: Level: intermediate
412: .seealso: DMLabelCreate(), DMLabelDestroy()
413: @*/
414: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
415: {
416: PetscBool iascii;
419: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);
421: if (label) DMLabelMakeAllValid_Private(label);
422: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
423: if (iascii) {
424: DMLabelView_Ascii(label, viewer);
425: }
426: return 0;
427: }
429: /*@
430: DMLabelReset - Destroys internal data structures in a DMLabel
432: Not collective
434: Input Parameter:
435: . label - The DMLabel
437: Level: beginner
439: .seealso: DMLabelDestroy(), DMLabelCreate()
440: @*/
441: PetscErrorCode DMLabelReset(DMLabel label)
442: {
443: PetscInt v;
446: for (v = 0; v < label->numStrata; ++v) {
447: PetscHSetIDestroy(&label->ht[v]);
448: ISDestroy(&label->points[v]);
449: }
450: label->numStrata = 0;
451: PetscFree(label->stratumValues);
452: PetscFree(label->stratumSizes);
453: PetscFree(label->ht);
454: PetscFree(label->points);
455: PetscFree(label->validIS);
456: label->stratumValues = NULL;
457: label->stratumSizes = NULL;
458: label->ht = NULL;
459: label->points = NULL;
460: label->validIS = NULL;
461: PetscHMapIReset(label->hmap);
462: label->pStart = -1;
463: label->pEnd = -1;
464: PetscBTDestroy(&label->bt);
465: return 0;
466: }
468: /*@
469: DMLabelDestroy - Destroys a DMLabel
471: Collective on label
473: Input Parameter:
474: . label - The DMLabel
476: Level: beginner
478: .seealso: DMLabelReset(), DMLabelCreate()
479: @*/
480: PetscErrorCode DMLabelDestroy(DMLabel *label)
481: {
482: if (!*label) return 0;
484: if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return 0;}
485: DMLabelReset(*label);
486: PetscHMapIDestroy(&(*label)->hmap);
487: PetscHeaderDestroy(label);
488: return 0;
489: }
491: /*@
492: DMLabelDuplicate - Duplicates a DMLabel
494: Collective on label
496: Input Parameter:
497: . label - The DMLabel
499: Output Parameter:
500: . labelnew - location to put new vector
502: Level: intermediate
504: .seealso: DMLabelCreate(), DMLabelDestroy()
505: @*/
506: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
507: {
508: const char *name;
509: PetscInt v;
512: DMLabelMakeAllValid_Private(label);
513: PetscObjectGetName((PetscObject) label, &name);
514: DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);
516: (*labelnew)->numStrata = label->numStrata;
517: (*labelnew)->defaultValue = label->defaultValue;
518: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
519: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
520: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
521: PetscMalloc1(label->numStrata, &(*labelnew)->points);
522: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
523: for (v = 0; v < label->numStrata; ++v) {
524: PetscHSetICreate(&(*labelnew)->ht[v]);
525: (*labelnew)->stratumValues[v] = label->stratumValues[v];
526: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
527: PetscObjectReference((PetscObject) (label->points[v]));
528: (*labelnew)->points[v] = label->points[v];
529: (*labelnew)->validIS[v] = PETSC_TRUE;
530: }
531: PetscHMapIDestroy(&(*labelnew)->hmap);
532: PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
533: (*labelnew)->pStart = -1;
534: (*labelnew)->pEnd = -1;
535: (*labelnew)->bt = NULL;
536: return 0;
537: }
539: /*@C
540: DMLabelCompare - Compare two DMLabel objects
542: Collective on comm
544: Input Parameters:
545: + comm - Comm over which to compare labels
546: . l0 - First DMLabel
547: - l1 - Second DMLabel
549: Output Parameters
550: + equal - (Optional) Flag whether the two labels are equal
551: - message - (Optional) Message describing the difference
553: Level: intermediate
555: Notes:
556: The output flag equal is the same on all processes.
557: If it is passed as NULL and difference is found, an error is thrown on all processes.
558: Make sure to pass NULL on all processes.
560: The output message is set independently on each rank.
561: It is set to NULL if no difference was found on the current rank. It must be freed by user.
562: If message is passed as NULL and difference is found, the difference description is printed to stderr in synchronized manner.
563: Make sure to pass NULL on all processes.
565: For the comparison, we ignore the order of stratum values, and strata with no points.
567: The communicator needs to be specified because currently DMLabel can live on PETSC_COMM_SELF even if the underlying DM is parallel.
569: Fortran Notes:
570: This function is currently not available from Fortran.
572: .seealso: DMCompareLabels(), DMLabelGetNumValues(), DMLabelGetDefaultValue(), DMLabelGetNonEmptyStratumValuesIS(), DMLabelGetStratumIS()
573: @*/
574: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
575: {
576: const char *name0, *name1;
577: char msg[PETSC_MAX_PATH_LEN] = "";
578: PetscBool eq;
579: PetscMPIInt rank;
585: MPI_Comm_rank(comm, &rank);
586: PetscObjectGetName((PetscObject)l0, &name0);
587: PetscObjectGetName((PetscObject)l1, &name1);
588: {
589: PetscInt v0, v1;
591: DMLabelGetDefaultValue(l0, &v0);
592: DMLabelGetDefaultValue(l1, &v1);
593: eq = (PetscBool) (v0 == v1);
594: if (!eq) {
595: PetscSNPrintf(msg, sizeof(msg), "Default value of DMLabel l0 \"%s\" = %D != %D = Default value of DMLabel l1 \"%s\"", name0, v0, v1, name1);
596: }
597: MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
598: if (!eq) goto finish;
599: }
600: {
601: IS is0, is1;
603: DMLabelGetNonEmptyStratumValuesIS(l0, &is0);
604: DMLabelGetNonEmptyStratumValuesIS(l1, &is1);
605: ISEqual(is0, is1, &eq);
606: ISDestroy(&is0);
607: ISDestroy(&is1);
608: if (!eq) {
609: PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1);
610: }
611: MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
612: if (!eq) goto finish;
613: }
614: {
615: PetscInt i, nValues;
617: DMLabelGetNumValues(l0, &nValues);
618: for (i=0; i<nValues; i++) {
619: const PetscInt v = l0->stratumValues[i];
620: PetscInt n;
621: IS is0, is1;
623: DMLabelGetStratumSize_Private(l0, i, &n);
624: if (!n) continue;
625: DMLabelGetStratumIS(l0, v, &is0);
626: DMLabelGetStratumIS(l1, v, &is1);
627: ISEqualUnsorted(is0, is1, &eq);
628: ISDestroy(&is0);
629: ISDestroy(&is1);
630: if (!eq) {
631: PetscSNPrintf(msg, sizeof(msg), "Stratum #%D with value %D contains different points in DMLabel l0 \"%s\" and DMLabel l1 \"%s\"", i, v, name0, name1);
632: break;
633: }
634: }
635: MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
636: }
637: finish:
638: /* If message output arg not set, print to stderr */
639: if (message) {
640: *message = NULL;
641: if (msg[0]) {
642: PetscStrallocpy(msg, message);
643: }
644: } else {
645: if (msg[0]) {
646: PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg);
647: }
648: PetscSynchronizedFlush(comm, PETSC_STDERR);
649: }
650: /* If same output arg not ser and labels are not equal, throw error */
651: if (equal) *equal = eq;
653: return 0;
654: }
656: /*@
657: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
659: Not collective
661: Input Parameter:
662: . label - The DMLabel
664: Level: intermediate
666: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
667: @*/
668: PetscErrorCode DMLabelComputeIndex(DMLabel label)
669: {
670: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
673: DMLabelMakeAllValid_Private(label);
674: for (v = 0; v < label->numStrata; ++v) {
675: const PetscInt *points;
676: PetscInt i;
678: ISGetIndices(label->points[v], &points);
679: for (i = 0; i < label->stratumSizes[v]; ++i) {
680: const PetscInt point = points[i];
682: pStart = PetscMin(point, pStart);
683: pEnd = PetscMax(point+1, pEnd);
684: }
685: ISRestoreIndices(label->points[v], &points);
686: }
687: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
688: label->pEnd = pEnd;
689: DMLabelCreateIndex(label, label->pStart, label->pEnd);
690: return 0;
691: }
693: /*@
694: DMLabelCreateIndex - Create an index structure for membership determination
696: Not collective
698: Input Parameters:
699: + label - The DMLabel
700: . pStart - The smallest point
701: - pEnd - The largest point + 1
703: Level: intermediate
705: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
706: @*/
707: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
708: {
709: PetscInt v;
712: DMLabelDestroyIndex(label);
713: DMLabelMakeAllValid_Private(label);
714: label->pStart = pStart;
715: label->pEnd = pEnd;
716: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
717: PetscBTCreate(pEnd - pStart, &label->bt);
718: for (v = 0; v < label->numStrata; ++v) {
719: const PetscInt *points;
720: PetscInt i;
722: ISGetIndices(label->points[v], &points);
723: for (i = 0; i < label->stratumSizes[v]; ++i) {
724: const PetscInt point = points[i];
727: PetscBTSet(label->bt, point - pStart);
728: }
729: ISRestoreIndices(label->points[v], &points);
730: }
731: return 0;
732: }
734: /*@
735: DMLabelDestroyIndex - Destroy the index structure
737: Not collective
739: Input Parameter:
740: . label - the DMLabel
742: Level: intermediate
744: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
745: @*/
746: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
747: {
749: label->pStart = -1;
750: label->pEnd = -1;
751: PetscBTDestroy(&label->bt);
752: return 0;
753: }
755: /*@
756: DMLabelGetBounds - Return the smallest and largest point in the label
758: Not collective
760: Input Parameter:
761: . label - the DMLabel
763: Output Parameters:
764: + pStart - The smallest point
765: - pEnd - The largest point + 1
767: Note: This will compute an index for the label if one does not exist.
769: Level: intermediate
771: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
772: @*/
773: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
774: {
776: if ((label->pStart == -1) && (label->pEnd == -1)) DMLabelComputeIndex(label);
777: if (pStart) {
779: *pStart = label->pStart;
780: }
781: if (pEnd) {
783: *pEnd = label->pEnd;
784: }
785: return 0;
786: }
788: /*@
789: DMLabelHasValue - Determine whether a label assigns the value to any point
791: Not collective
793: Input Parameters:
794: + label - the DMLabel
795: - value - the value
797: Output Parameter:
798: . contains - Flag indicating whether the label maps this value to any point
800: Level: developer
802: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
803: @*/
804: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
805: {
806: PetscInt v;
810: DMLabelLookupStratum(label, value, &v);
811: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
812: return 0;
813: }
815: /*@
816: DMLabelHasPoint - Determine whether a label assigns a value to a point
818: Not collective
820: Input Parameters:
821: + label - the DMLabel
822: - point - the point
824: Output Parameter:
825: . contains - Flag indicating whether the label maps this point to a value
827: Note: The user must call DMLabelCreateIndex() before this function.
829: Level: developer
831: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
832: @*/
833: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
834: {
838: DMLabelMakeAllValid_Private(label);
839: if (PetscDefined(USE_DEBUG)) {
842: }
843: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
844: return 0;
845: }
847: /*@
848: DMLabelStratumHasPoint - Return true if the stratum contains a point
850: Not collective
852: Input Parameters:
853: + label - the DMLabel
854: . value - the stratum value
855: - point - the point
857: Output Parameter:
858: . contains - true if the stratum contains the point
860: Level: intermediate
862: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
863: @*/
864: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
865: {
866: PetscInt v;
871: *contains = PETSC_FALSE;
872: DMLabelLookupStratum(label, value, &v);
873: if (v < 0) return 0;
875: if (label->validIS[v]) {
876: PetscInt i;
878: ISLocate(label->points[v], point, &i);
879: if (i >= 0) *contains = PETSC_TRUE;
880: } else {
881: PetscBool has;
883: PetscHSetIHas(label->ht[v], point, &has);
884: if (has) *contains = PETSC_TRUE;
885: }
886: return 0;
887: }
889: /*@
890: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
891: When a label is created, it is initialized to -1.
893: Not collective
895: Input parameter:
896: . label - a DMLabel object
898: Output parameter:
899: . defaultValue - the default value
901: Level: beginner
903: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
904: @*/
905: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
906: {
908: *defaultValue = label->defaultValue;
909: return 0;
910: }
912: /*@
913: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
914: When a label is created, it is initialized to -1.
916: Not collective
918: Input parameter:
919: . label - a DMLabel object
921: Output parameter:
922: . defaultValue - the default value
924: Level: beginner
926: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
927: @*/
928: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
929: {
931: label->defaultValue = defaultValue;
932: return 0;
933: }
935: /*@
936: 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())
938: Not collective
940: Input Parameters:
941: + label - the DMLabel
942: - point - the point
944: Output Parameter:
945: . value - The point value, or the default value (-1 by default)
947: Level: intermediate
949: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
950: @*/
951: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
952: {
953: PetscInt v;
958: *value = label->defaultValue;
959: for (v = 0; v < label->numStrata; ++v) {
960: if (label->validIS[v]) {
961: PetscInt i;
963: ISLocate(label->points[v], point, &i);
964: if (i >= 0) {
965: *value = label->stratumValues[v];
966: break;
967: }
968: } else {
969: PetscBool has;
971: PetscHSetIHas(label->ht[v], point, &has);
972: if (has) {
973: *value = label->stratumValues[v];
974: break;
975: }
976: }
977: }
978: return 0;
979: }
981: /*@
982: 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.
984: Not collective
986: Input Parameters:
987: + label - the DMLabel
988: . point - the point
989: - value - The point value
991: Level: intermediate
993: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
994: @*/
995: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
996: {
997: PetscInt v;
1000: /* Find label value, add new entry if needed */
1001: if (value == label->defaultValue) return 0;
1002: DMLabelLookupAddStratum(label, value, &v);
1003: /* Set key */
1004: DMLabelMakeInvalid_Private(label, v);
1005: PetscHSetIAdd(label->ht[v], point);
1006: return 0;
1007: }
1009: /*@
1010: DMLabelClearValue - Clear the value a label assigns to a point
1012: Not collective
1014: Input Parameters:
1015: + label - the DMLabel
1016: . point - the point
1017: - value - The point value
1019: Level: intermediate
1021: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
1022: @*/
1023: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1024: {
1025: PetscInt v;
1028: /* Find label value */
1029: DMLabelLookupStratum(label, value, &v);
1030: if (v < 0) return 0;
1032: if (label->bt) {
1034: PetscBTClear(label->bt, point - label->pStart);
1035: }
1037: /* Delete key */
1038: DMLabelMakeInvalid_Private(label, v);
1039: PetscHSetIDel(label->ht[v], point);
1040: return 0;
1041: }
1043: /*@
1044: DMLabelInsertIS - Set all points in the IS to a value
1046: Not collective
1048: Input Parameters:
1049: + label - the DMLabel
1050: . is - the point IS
1051: - value - The point value
1053: Level: intermediate
1055: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1056: @*/
1057: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1058: {
1059: PetscInt v, n, p;
1060: const PetscInt *points;
1064: /* Find label value, add new entry if needed */
1065: if (value == label->defaultValue) return 0;
1066: DMLabelLookupAddStratum(label, value, &v);
1067: /* Set keys */
1068: DMLabelMakeInvalid_Private(label, v);
1069: ISGetLocalSize(is, &n);
1070: ISGetIndices(is, &points);
1071: for (p = 0; p < n; ++p) PetscHSetIAdd(label->ht[v], points[p]);
1072: ISRestoreIndices(is, &points);
1073: return 0;
1074: }
1076: /*@
1077: DMLabelGetNumValues - Get the number of values that the DMLabel takes
1079: Not collective
1081: Input Parameter:
1082: . label - the DMLabel
1084: Output Parameter:
1085: . numValues - the number of values
1087: Level: intermediate
1089: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1090: @*/
1091: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1092: {
1095: *numValues = label->numStrata;
1096: return 0;
1097: }
1099: /*@
1100: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1102: Not collective
1104: Input Parameter:
1105: . label - the DMLabel
1107: Output Parameter:
1108: . is - the value IS
1110: Level: intermediate
1112: Notes:
1113: The output IS should be destroyed when no longer needed.
1114: Strata which are allocated but empty [DMLabelGetStratumSize() yields 0] are counted.
1115: If you need to count only nonempty strata, use DMLabelGetNonEmptyStratumValuesIS().
1117: .seealso: DMLabelGetNonEmptyStratumValuesIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1118: @*/
1119: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1120: {
1123: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1124: return 0;
1125: }
1127: /*@
1128: DMLabelGetNonEmptyStratumValuesIS - Get an IS of all values that the DMlabel takes
1130: Not collective
1132: Input Parameter:
1133: . label - the DMLabel
1135: Output Paramater:
1136: . is - the value IS
1138: Level: intermediate
1140: Notes:
1141: The output IS should be destroyed when no longer needed.
1142: This is similar to DMLabelGetValueIS() but counts only nonempty strata.
1144: .seealso: DMLabelGetValueIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1145: @*/
1146: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1147: {
1148: PetscInt i, j;
1149: PetscInt *valuesArr;
1153: PetscMalloc1(label->numStrata, &valuesArr);
1154: for (i = 0, j = 0; i < label->numStrata; i++) {
1155: PetscInt n;
1157: DMLabelGetStratumSize_Private(label, i, &n);
1158: if (n) valuesArr[j++] = label->stratumValues[i];
1159: }
1160: if (j == label->numStrata) {
1161: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1162: } else {
1163: ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values);
1164: }
1165: PetscFree(valuesArr);
1166: return 0;
1167: }
1169: /*@
1170: DMLabelGetValueIndex - Get the index of a given value in the list of values for the DMlabel, or -1 if it is not present
1172: Not collective
1174: Input Parameters:
1175: + label - the DMLabel
1176: - value - the value
1178: Output Parameter:
1179: . index - the index of value in the list of values
1181: Level: intermediate
1183: .seealso: DMLabelGetValueIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1184: @*/
1185: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1186: {
1187: PetscInt v;
1191: /* Do not assume they are sorted */
1192: for (v = 0; v < label->numStrata; ++v) if (label->stratumValues[v] == value) break;
1193: if (v >= label->numStrata) *index = -1;
1194: else *index = v;
1195: return 0;
1196: }
1198: /*@
1199: DMLabelHasStratum - Determine whether points exist with the given value
1201: Not collective
1203: Input Parameters:
1204: + label - the DMLabel
1205: - value - the stratum value
1207: Output Parameter:
1208: . exists - Flag saying whether points exist
1210: Level: intermediate
1212: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1213: @*/
1214: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1215: {
1216: PetscInt v;
1220: DMLabelLookupStratum(label, value, &v);
1221: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1222: return 0;
1223: }
1225: /*@
1226: DMLabelGetStratumSize - Get the size of a stratum
1228: Not collective
1230: Input Parameters:
1231: + label - the DMLabel
1232: - value - the stratum value
1234: Output Parameter:
1235: . size - The number of points in the stratum
1237: Level: intermediate
1239: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1240: @*/
1241: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1242: {
1243: PetscInt v;
1247: DMLabelLookupStratum(label, value, &v);
1248: DMLabelGetStratumSize_Private(label, v, size);
1249: return 0;
1250: }
1252: /*@
1253: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1255: Not collective
1257: Input Parameters:
1258: + label - the DMLabel
1259: - value - the stratum value
1261: Output Parameters:
1262: + start - the smallest point in the stratum
1263: - end - the largest point in the stratum
1265: Level: intermediate
1267: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1268: @*/
1269: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1270: {
1271: PetscInt v, min, max;
1276: DMLabelLookupStratum(label, value, &v);
1277: if (v < 0) return 0;
1278: DMLabelMakeValid_Private(label, v);
1279: if (label->stratumSizes[v] <= 0) return 0;
1280: ISGetMinMax(label->points[v], &min, &max);
1281: if (start) *start = min;
1282: if (end) *end = max+1;
1283: return 0;
1284: }
1286: /*@
1287: DMLabelGetStratumIS - Get an IS with the stratum points
1289: Not collective
1291: Input Parameters:
1292: + label - the DMLabel
1293: - value - the stratum value
1295: Output Parameter:
1296: . points - The stratum points
1298: Level: intermediate
1300: Notes:
1301: The output IS should be destroyed when no longer needed.
1302: Returns NULL if the stratum is empty.
1304: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1305: @*/
1306: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1307: {
1308: PetscInt v;
1312: *points = NULL;
1313: DMLabelLookupStratum(label, value, &v);
1314: if (v < 0) return 0;
1315: DMLabelMakeValid_Private(label, v);
1316: PetscObjectReference((PetscObject) label->points[v]);
1317: *points = label->points[v];
1318: return 0;
1319: }
1321: /*@
1322: DMLabelSetStratumIS - Set the stratum points using an IS
1324: Not collective
1326: Input Parameters:
1327: + label - the DMLabel
1328: . value - the stratum value
1329: - points - The stratum points
1331: Level: intermediate
1333: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1334: @*/
1335: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1336: {
1337: PetscInt v;
1341: DMLabelLookupAddStratum(label, value, &v);
1342: if (is == label->points[v]) return 0;
1343: DMLabelClearStratum(label, value);
1344: ISGetLocalSize(is, &(label->stratumSizes[v]));
1345: PetscObjectReference((PetscObject)is);
1346: ISDestroy(&(label->points[v]));
1347: label->points[v] = is;
1348: label->validIS[v] = PETSC_TRUE;
1349: PetscObjectStateIncrease((PetscObject) label);
1350: if (label->bt) {
1351: const PetscInt *points;
1352: PetscInt p;
1354: ISGetIndices(is,&points);
1355: for (p = 0; p < label->stratumSizes[v]; ++p) {
1356: const PetscInt point = points[p];
1359: PetscBTSet(label->bt, point - label->pStart);
1360: }
1361: }
1362: return 0;
1363: }
1365: /*@
1366: DMLabelClearStratum - Remove a stratum
1368: Not collective
1370: Input Parameters:
1371: + label - the DMLabel
1372: - value - the stratum value
1374: Level: intermediate
1376: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1377: @*/
1378: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1379: {
1380: PetscInt v;
1383: DMLabelLookupStratum(label, value, &v);
1384: if (v < 0) return 0;
1385: if (label->validIS[v]) {
1386: if (label->bt) {
1387: PetscInt i;
1388: const PetscInt *points;
1390: ISGetIndices(label->points[v], &points);
1391: for (i = 0; i < label->stratumSizes[v]; ++i) {
1392: const PetscInt point = points[i];
1395: PetscBTClear(label->bt, point - label->pStart);
1396: }
1397: ISRestoreIndices(label->points[v], &points);
1398: }
1399: label->stratumSizes[v] = 0;
1400: ISDestroy(&label->points[v]);
1401: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1402: PetscObjectSetName((PetscObject) label->points[v], "indices");
1403: PetscObjectStateIncrease((PetscObject) label);
1404: } else {
1405: PetscHSetIClear(label->ht[v]);
1406: }
1407: return 0;
1408: }
1410: /*@
1411: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1413: Not collective
1415: Input Parameters:
1416: + label - The DMLabel
1417: . value - The label value for all points
1418: . pStart - The first point
1419: - pEnd - A point beyond all marked points
1421: Note: The marks points are [pStart, pEnd), and only the bounds are stored.
1423: Level: intermediate
1425: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1426: @*/
1427: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1428: {
1429: IS pIS;
1431: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1432: DMLabelSetStratumIS(label, value, pIS);
1433: ISDestroy(&pIS);
1434: return 0;
1435: }
1437: /*@
1438: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1440: Not collective
1442: Input Parameters:
1443: + label - The DMLabel
1444: . value - The label value
1445: - p - A point with this value
1447: Output Parameter:
1448: . 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
1450: Level: intermediate
1452: .seealso: DMLabelGetValueIndex(), DMLabelGetStratumIS(), DMLabelCreate()
1453: @*/
1454: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1455: {
1456: const PetscInt *indices;
1457: PetscInt v;
1461: *index = -1;
1462: DMLabelLookupStratum(label, value, &v);
1463: if (v < 0) return 0;
1464: DMLabelMakeValid_Private(label, v);
1465: ISGetIndices(label->points[v], &indices);
1466: PetscFindInt(p, label->stratumSizes[v], indices, index);
1467: ISRestoreIndices(label->points[v], &indices);
1468: return 0;
1469: }
1471: /*@
1472: DMLabelFilter - Remove all points outside of [start, end)
1474: Not collective
1476: Input Parameters:
1477: + label - the DMLabel
1478: . start - the first point kept
1479: - end - one more than the last point kept
1481: Level: intermediate
1483: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1484: @*/
1485: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1486: {
1487: PetscInt v;
1490: DMLabelDestroyIndex(label);
1491: DMLabelMakeAllValid_Private(label);
1492: for (v = 0; v < label->numStrata; ++v) {
1493: ISGeneralFilter(label->points[v], start, end);
1494: }
1495: DMLabelCreateIndex(label, start, end);
1496: return 0;
1497: }
1499: /*@
1500: DMLabelPermute - Create a new label with permuted points
1502: Not collective
1504: Input Parameters:
1505: + label - the DMLabel
1506: - permutation - the point permutation
1508: Output Parameter:
1509: . labelnew - the new label containing the permuted points
1511: Level: intermediate
1513: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1514: @*/
1515: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1516: {
1517: const PetscInt *perm;
1518: PetscInt numValues, numPoints, v, q;
1522: DMLabelMakeAllValid_Private(label);
1523: DMLabelDuplicate(label, labelNew);
1524: DMLabelGetNumValues(*labelNew, &numValues);
1525: ISGetLocalSize(permutation, &numPoints);
1526: ISGetIndices(permutation, &perm);
1527: for (v = 0; v < numValues; ++v) {
1528: const PetscInt size = (*labelNew)->stratumSizes[v];
1529: const PetscInt *points;
1530: PetscInt *pointsNew;
1532: ISGetIndices((*labelNew)->points[v],&points);
1533: PetscMalloc1(size,&pointsNew);
1534: for (q = 0; q < size; ++q) {
1535: const PetscInt point = points[q];
1538: pointsNew[q] = perm[point];
1539: }
1540: ISRestoreIndices((*labelNew)->points[v],&points);
1541: PetscSortInt(size, pointsNew);
1542: ISDestroy(&((*labelNew)->points[v]));
1543: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1544: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1545: PetscFree(pointsNew);
1546: } else {
1547: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1548: }
1549: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1550: }
1551: ISRestoreIndices(permutation, &perm);
1552: if (label->bt) {
1553: PetscBTDestroy(&label->bt);
1554: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1555: }
1556: return 0;
1557: }
1559: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1560: {
1561: MPI_Comm comm;
1562: PetscInt s, l, nroots, nleaves, offset, size;
1563: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1564: PetscSection rootSection;
1565: PetscSF labelSF;
1567: if (label) DMLabelMakeAllValid_Private(label);
1568: PetscObjectGetComm((PetscObject)sf, &comm);
1569: /* Build a section of stratum values per point, generate the according SF
1570: and distribute point-wise stratum values to leaves. */
1571: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1572: PetscSectionCreate(comm, &rootSection);
1573: PetscSectionSetChart(rootSection, 0, nroots);
1574: if (label) {
1575: for (s = 0; s < label->numStrata; ++s) {
1576: const PetscInt *points;
1578: ISGetIndices(label->points[s], &points);
1579: for (l = 0; l < label->stratumSizes[s]; l++) {
1580: PetscSectionAddDof(rootSection, points[l], 1);
1581: }
1582: ISRestoreIndices(label->points[s], &points);
1583: }
1584: }
1585: PetscSectionSetUp(rootSection);
1586: /* Create a point-wise array of stratum values */
1587: PetscSectionGetStorageSize(rootSection, &size);
1588: PetscMalloc1(size, &rootStrata);
1589: PetscCalloc1(nroots, &rootIdx);
1590: if (label) {
1591: for (s = 0; s < label->numStrata; ++s) {
1592: const PetscInt *points;
1594: ISGetIndices(label->points[s], &points);
1595: for (l = 0; l < label->stratumSizes[s]; l++) {
1596: const PetscInt p = points[l];
1597: PetscSectionGetOffset(rootSection, p, &offset);
1598: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1599: }
1600: ISRestoreIndices(label->points[s], &points);
1601: }
1602: }
1603: /* Build SF that maps label points to remote processes */
1604: PetscSectionCreate(comm, leafSection);
1605: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1606: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1607: PetscFree(remoteOffsets);
1608: /* Send the strata for each point over the derived SF */
1609: PetscSectionGetStorageSize(*leafSection, &size);
1610: PetscMalloc1(size, leafStrata);
1611: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1612: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1613: /* Clean up */
1614: PetscFree(rootStrata);
1615: PetscFree(rootIdx);
1616: PetscSectionDestroy(&rootSection);
1617: PetscSFDestroy(&labelSF);
1618: return 0;
1619: }
1621: /*@
1622: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1624: Collective on sf
1626: Input Parameters:
1627: + label - the DMLabel
1628: - sf - the map from old to new distribution
1630: Output Parameter:
1631: . labelnew - the new redistributed label
1633: Level: intermediate
1635: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1636: @*/
1637: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1638: {
1639: MPI_Comm comm;
1640: PetscSection leafSection;
1641: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1642: PetscInt *leafStrata, *strataIdx;
1643: PetscInt **points;
1644: const char *lname = NULL;
1645: char *name;
1646: PetscInt nameSize;
1647: PetscHSetI stratumHash;
1648: size_t len = 0;
1649: PetscMPIInt rank;
1652: if (label) {
1654: DMLabelMakeAllValid_Private(label);
1655: }
1656: PetscObjectGetComm((PetscObject)sf, &comm);
1657: MPI_Comm_rank(comm, &rank);
1658: /* Bcast name */
1659: if (rank == 0) {
1660: PetscObjectGetName((PetscObject) label, &lname);
1661: PetscStrlen(lname, &len);
1662: }
1663: nameSize = len;
1664: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1665: PetscMalloc1(nameSize+1, &name);
1666: if (rank == 0) PetscArraycpy(name, lname, nameSize+1);
1667: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1668: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1669: PetscFree(name);
1670: /* Bcast defaultValue */
1671: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1672: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1673: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1674: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1675: /* Determine received stratum values and initialise new label*/
1676: PetscHSetICreate(&stratumHash);
1677: PetscSectionGetStorageSize(leafSection, &size);
1678: for (p = 0; p < size; ++p) PetscHSetIAdd(stratumHash, leafStrata[p]);
1679: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1680: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1681: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1682: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1683: /* Turn leafStrata into indices rather than stratum values */
1684: offset = 0;
1685: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1686: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1687: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1688: PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1689: }
1690: for (p = 0; p < size; ++p) {
1691: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1692: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1693: }
1694: }
1695: /* Rebuild the point strata on the receiver */
1696: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1697: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1698: for (p=pStart; p<pEnd; p++) {
1699: PetscSectionGetDof(leafSection, p, &dof);
1700: PetscSectionGetOffset(leafSection, p, &offset);
1701: for (s=0; s<dof; s++) {
1702: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1703: }
1704: }
1705: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1706: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1707: PetscMalloc1((*labelNew)->numStrata,&points);
1708: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1709: PetscHSetICreate(&(*labelNew)->ht[s]);
1710: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1711: }
1712: /* Insert points into new strata */
1713: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1714: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1715: for (p=pStart; p<pEnd; p++) {
1716: PetscSectionGetDof(leafSection, p, &dof);
1717: PetscSectionGetOffset(leafSection, p, &offset);
1718: for (s=0; s<dof; s++) {
1719: stratum = leafStrata[offset+s];
1720: points[stratum][strataIdx[stratum]++] = p;
1721: }
1722: }
1723: for (s = 0; s < (*labelNew)->numStrata; s++) {
1724: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1725: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1726: }
1727: PetscFree(points);
1728: PetscHSetIDestroy(&stratumHash);
1729: PetscFree(leafStrata);
1730: PetscFree(strataIdx);
1731: PetscSectionDestroy(&leafSection);
1732: return 0;
1733: }
1735: /*@
1736: DMLabelGather - Gather all label values from leafs into roots
1738: Collective on sf
1740: Input Parameters:
1741: + label - the DMLabel
1742: - sf - the Star Forest point communication map
1744: Output Parameters:
1745: . labelNew - the new DMLabel with localised leaf values
1747: Level: developer
1749: Note: This is the inverse operation to DMLabelDistribute.
1751: .seealso: DMLabelDistribute()
1752: @*/
1753: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1754: {
1755: MPI_Comm comm;
1756: PetscSection rootSection;
1757: PetscSF sfLabel;
1758: PetscSFNode *rootPoints, *leafPoints;
1759: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1760: const PetscInt *rootDegree, *ilocal;
1761: PetscInt *rootStrata;
1762: const char *lname;
1763: char *name;
1764: PetscInt nameSize;
1765: size_t len = 0;
1766: PetscMPIInt rank, size;
1770: PetscObjectGetComm((PetscObject)sf, &comm);
1771: MPI_Comm_rank(comm, &rank);
1772: MPI_Comm_size(comm, &size);
1773: /* Bcast name */
1774: if (rank == 0) {
1775: PetscObjectGetName((PetscObject) label, &lname);
1776: PetscStrlen(lname, &len);
1777: }
1778: nameSize = len;
1779: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1780: PetscMalloc1(nameSize+1, &name);
1781: if (rank == 0) PetscArraycpy(name, lname, nameSize+1);
1782: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1783: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1784: PetscFree(name);
1785: /* Gather rank/index pairs of leaves into local roots to build
1786: an inverse, multi-rooted SF. Note that this ignores local leaf
1787: indexing due to the use of the multiSF in PetscSFGather. */
1788: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1789: PetscMalloc1(nroots, &leafPoints);
1790: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1791: for (p = 0; p < nleaves; p++) {
1792: PetscInt ilp = ilocal ? ilocal[p] : p;
1794: leafPoints[ilp].index = ilp;
1795: leafPoints[ilp].rank = rank;
1796: }
1797: PetscSFComputeDegreeBegin(sf, &rootDegree);
1798: PetscSFComputeDegreeEnd(sf, &rootDegree);
1799: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1800: PetscMalloc1(nmultiroots, &rootPoints);
1801: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1802: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1803: PetscSFCreate(comm,& sfLabel);
1804: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1805: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1806: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1807: /* Rebuild the point strata on the receiver */
1808: for (p = 0, idx = 0; p < nroots; p++) {
1809: for (d = 0; d < rootDegree[p]; d++) {
1810: PetscSectionGetDof(rootSection, idx+d, &dof);
1811: PetscSectionGetOffset(rootSection, idx+d, &offset);
1812: for (s = 0; s < dof; s++) DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);
1813: }
1814: idx += rootDegree[p];
1815: }
1816: PetscFree(leafPoints);
1817: PetscFree(rootStrata);
1818: PetscSectionDestroy(&rootSection);
1819: PetscSFDestroy(&sfLabel);
1820: return 0;
1821: }
1823: /*@
1824: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1826: Not collective
1828: Input Parameter:
1829: . label - the DMLabel
1831: Output Parameters:
1832: + section - the section giving offsets for each stratum
1833: - is - An IS containing all the label points
1835: Level: developer
1837: .seealso: DMLabelDistribute()
1838: @*/
1839: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1840: {
1841: IS vIS;
1842: const PetscInt *values;
1843: PetscInt *points;
1844: PetscInt nV, vS = 0, vE = 0, v, N;
1847: DMLabelGetNumValues(label, &nV);
1848: DMLabelGetValueIS(label, &vIS);
1849: ISGetIndices(vIS, &values);
1850: if (nV) {vS = values[0]; vE = values[0]+1;}
1851: for (v = 1; v < nV; ++v) {
1852: vS = PetscMin(vS, values[v]);
1853: vE = PetscMax(vE, values[v]+1);
1854: }
1855: PetscSectionCreate(PETSC_COMM_SELF, section);
1856: PetscSectionSetChart(*section, vS, vE);
1857: for (v = 0; v < nV; ++v) {
1858: PetscInt n;
1860: DMLabelGetStratumSize(label, values[v], &n);
1861: PetscSectionSetDof(*section, values[v], n);
1862: }
1863: PetscSectionSetUp(*section);
1864: PetscSectionGetStorageSize(*section, &N);
1865: PetscMalloc1(N, &points);
1866: for (v = 0; v < nV; ++v) {
1867: IS is;
1868: const PetscInt *spoints;
1869: PetscInt dof, off, p;
1871: PetscSectionGetDof(*section, values[v], &dof);
1872: PetscSectionGetOffset(*section, values[v], &off);
1873: DMLabelGetStratumIS(label, values[v], &is);
1874: ISGetIndices(is, &spoints);
1875: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1876: ISRestoreIndices(is, &spoints);
1877: ISDestroy(&is);
1878: }
1879: ISRestoreIndices(vIS, &values);
1880: ISDestroy(&vIS);
1881: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1882: return 0;
1883: }
1885: /*@
1886: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1887: the local section and an SF describing the section point overlap.
1889: Collective on sf
1891: Input Parameters:
1892: + s - The PetscSection for the local field layout
1893: . sf - The SF describing parallel layout of the section points
1894: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1895: . label - The label specifying the points
1896: - labelValue - The label stratum specifying the points
1898: Output Parameter:
1899: . gsection - The PetscSection for the global field layout
1901: Note: This gives negative sizes and offsets to points not owned by this process
1903: Level: developer
1905: .seealso: PetscSectionCreate()
1906: @*/
1907: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1908: {
1909: PetscInt *neg = NULL, *tmpOff = NULL;
1910: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1915: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1916: PetscSectionGetChart(s, &pStart, &pEnd);
1917: PetscSectionSetChart(*gsection, pStart, pEnd);
1918: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1919: if (nroots >= 0) {
1921: PetscCalloc1(nroots, &neg);
1922: if (nroots > pEnd-pStart) {
1923: PetscCalloc1(nroots, &tmpOff);
1924: } else {
1925: tmpOff = &(*gsection)->atlasDof[-pStart];
1926: }
1927: }
1928: /* Mark ghost points with negative dof */
1929: for (p = pStart; p < pEnd; ++p) {
1930: PetscInt value;
1932: DMLabelGetValue(label, p, &value);
1933: if (value != labelValue) continue;
1934: PetscSectionGetDof(s, p, &dof);
1935: PetscSectionSetDof(*gsection, p, dof);
1936: PetscSectionGetConstraintDof(s, p, &cdof);
1937: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
1938: if (neg) neg[p] = -(dof+1);
1939: }
1940: PetscSectionSetUpBC(*gsection);
1941: if (nroots >= 0) {
1942: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1943: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1944: if (nroots > pEnd-pStart) {
1945: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1946: }
1947: }
1948: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1949: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1950: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1951: (*gsection)->atlasOff[p] = off;
1952: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1953: }
1954: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1955: globalOff -= off;
1956: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1957: (*gsection)->atlasOff[p] += globalOff;
1958: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1959: }
1960: /* Put in negative offsets for ghost points */
1961: if (nroots >= 0) {
1962: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1963: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1964: if (nroots > pEnd-pStart) {
1965: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1966: }
1967: }
1968: if (nroots >= 0 && nroots > pEnd-pStart) PetscFree(tmpOff);
1969: PetscFree(neg);
1970: return 0;
1971: }
1973: typedef struct _n_PetscSectionSym_Label
1974: {
1975: DMLabel label;
1976: PetscCopyMode *modes;
1977: PetscInt *sizes;
1978: const PetscInt ***perms;
1979: const PetscScalar ***rots;
1980: PetscInt (*minMaxOrients)[2];
1981: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1982: } PetscSectionSym_Label;
1984: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1985: {
1986: PetscInt i, j;
1987: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1989: for (i = 0; i <= sl->numStrata; i++) {
1990: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1991: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1992: if (sl->perms[i]) PetscFree(sl->perms[i][j]);
1993: if (sl->rots[i]) PetscFree(sl->rots[i][j]);
1994: }
1995: if (sl->perms[i]) {
1996: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1998: PetscFree(perms);
1999: }
2000: if (sl->rots[i]) {
2001: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
2003: PetscFree(rots);
2004: }
2005: }
2006: }
2007: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
2008: DMLabelDestroy(&sl->label);
2009: sl->numStrata = 0;
2010: return 0;
2011: }
2013: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2014: {
2015: PetscSectionSymLabelReset(sym);
2016: PetscFree(sym->data);
2017: return 0;
2018: }
2020: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2021: {
2022: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
2023: PetscBool isAscii;
2024: DMLabel label = sl->label;
2025: const char *name;
2027: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
2028: if (isAscii) {
2029: PetscInt i, j, k;
2030: PetscViewerFormat format;
2032: PetscViewerGetFormat(viewer,&format);
2033: if (label) {
2034: PetscViewerGetFormat(viewer,&format);
2035: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2036: PetscViewerASCIIPushTab(viewer);
2037: DMLabelView(label, viewer);
2038: PetscViewerASCIIPopTab(viewer);
2039: } else {
2040: PetscObjectGetName((PetscObject) sl->label, &name);
2041: PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);
2042: }
2043: } else {
2044: PetscViewerASCIIPrintf(viewer, "No label given\n");
2045: }
2046: PetscViewerASCIIPushTab(viewer);
2047: for (i = 0; i <= sl->numStrata; i++) {
2048: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
2050: if (!(sl->perms[i] || sl->rots[i])) {
2051: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
2052: } else {
2053: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
2054: PetscViewerASCIIPushTab(viewer);
2055: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
2056: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2057: PetscViewerASCIIPushTab(viewer);
2058: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2059: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2060: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
2061: } else {
2062: PetscInt tab;
2064: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
2065: PetscViewerASCIIPushTab(viewer);
2066: PetscViewerASCIIGetTab(viewer,&tab);
2067: if (sl->perms[i] && sl->perms[i][j]) {
2068: PetscViewerASCIIPrintf(viewer,"Permutation:");
2069: PetscViewerASCIISetTab(viewer,0);
2070: for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);
2071: PetscViewerASCIIPrintf(viewer,"\n");
2072: PetscViewerASCIISetTab(viewer,tab);
2073: }
2074: if (sl->rots[i] && sl->rots[i][j]) {
2075: PetscViewerASCIIPrintf(viewer,"Rotations: ");
2076: PetscViewerASCIISetTab(viewer,0);
2077: #if defined(PETSC_USE_COMPLEX)
2078: for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer," %+f+i*%+f",PetscRealPart(sl->rots[i][j][k]),PetscImaginaryPart(sl->rots[i][j][k]));
2079: #else
2080: for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);
2081: #endif
2082: PetscViewerASCIIPrintf(viewer,"\n");
2083: PetscViewerASCIISetTab(viewer,tab);
2084: }
2085: PetscViewerASCIIPopTab(viewer);
2086: }
2087: }
2088: PetscViewerASCIIPopTab(viewer);
2089: }
2090: PetscViewerASCIIPopTab(viewer);
2091: }
2092: }
2093: PetscViewerASCIIPopTab(viewer);
2094: }
2095: return 0;
2096: }
2098: /*@
2099: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2101: Logically collective on sym
2103: Input parameters:
2104: + sym - the section symmetries
2105: - label - the DMLabel describing the types of points
2107: Level: developer:
2109: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
2110: @*/
2111: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2112: {
2113: PetscSectionSym_Label *sl;
2116: sl = (PetscSectionSym_Label *) sym->data;
2117: if (sl->label && sl->label != label) PetscSectionSymLabelReset(sym);
2118: if (label) {
2119: sl->label = label;
2120: PetscObjectReference((PetscObject) label);
2121: DMLabelGetNumValues(label,&sl->numStrata);
2122: 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);
2123: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
2124: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
2125: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
2126: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
2127: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
2128: }
2129: return 0;
2130: }
2132: /*@C
2133: PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum
2135: Logically collective on sym
2137: Input Parameters:
2138: + sym - the section symmetries
2139: - stratum - the stratum value in the label that we are assigning symmetries for
2141: Output Parameters:
2142: + size - the number of dofs for points in the stratum of the label
2143: . minOrient - the smallest orientation for a point in this stratum
2144: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2145: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2146: - 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
2148: Level: developer
2150: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2151: @*/
2152: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2153: {
2154: PetscSectionSym_Label *sl;
2155: const char *name;
2156: PetscInt i;
2159: sl = (PetscSectionSym_Label *) sym->data;
2161: for (i = 0; i <= sl->numStrata; i++) {
2162: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2164: if (stratum == value) break;
2165: }
2166: PetscObjectGetName((PetscObject) sl->label, &name);
2173: return 0;
2174: }
2176: /*@C
2177: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2179: Logically collective on sym
2181: InputParameters:
2182: + sym - the section symmetries
2183: . stratum - the stratum value in the label that we are assigning symmetries for
2184: . size - the number of dofs for points in the stratum of the label
2185: . minOrient - the smallest orientation for a point in this stratum
2186: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2187: . mode - how sym should copy the perms and rots arrays
2188: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2189: - 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
2191: Level: developer
2193: .seealso: PetscSectionSymLabelGetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2194: @*/
2195: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2196: {
2197: PetscSectionSym_Label *sl;
2198: const char *name;
2199: PetscInt i, j, k;
2202: sl = (PetscSectionSym_Label *) sym->data;
2204: for (i = 0; i <= sl->numStrata; i++) {
2205: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2207: if (stratum == value) break;
2208: }
2209: PetscObjectGetName((PetscObject) sl->label, &name);
2211: sl->sizes[i] = size;
2212: sl->modes[i] = mode;
2213: sl->minMaxOrients[i][0] = minOrient;
2214: sl->minMaxOrients[i][1] = maxOrient;
2215: if (mode == PETSC_COPY_VALUES) {
2216: if (perms) {
2217: PetscInt **ownPerms;
2219: PetscCalloc1(maxOrient - minOrient,&ownPerms);
2220: for (j = 0; j < maxOrient-minOrient; j++) {
2221: if (perms[j]) {
2222: PetscMalloc1(size,&ownPerms[j]);
2223: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2224: }
2225: }
2226: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2227: }
2228: if (rots) {
2229: PetscScalar **ownRots;
2231: PetscCalloc1(maxOrient - minOrient,&ownRots);
2232: for (j = 0; j < maxOrient-minOrient; j++) {
2233: if (rots[j]) {
2234: PetscMalloc1(size,&ownRots[j]);
2235: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2236: }
2237: }
2238: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2239: }
2240: } else {
2241: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2242: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2243: }
2244: return 0;
2245: }
2247: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2248: {
2249: PetscInt i, j, numStrata;
2250: PetscSectionSym_Label *sl;
2251: DMLabel label;
2253: sl = (PetscSectionSym_Label *) sym->data;
2254: numStrata = sl->numStrata;
2255: label = sl->label;
2256: for (i = 0; i < numPoints; i++) {
2257: PetscInt point = points[2*i];
2258: PetscInt ornt = points[2*i+1];
2260: for (j = 0; j < numStrata; j++) {
2261: if (label->validIS[j]) {
2262: PetscInt k;
2264: ISLocate(label->points[j],point,&k);
2265: if (k >= 0) break;
2266: } else {
2267: PetscBool has;
2269: PetscHSetIHas(label->ht[j], point, &has);
2270: if (has) break;
2271: }
2272: }
2274: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2275: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2276: }
2277: return 0;
2278: }
2280: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2281: {
2282: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) nsym->data;
2283: IS valIS;
2284: const PetscInt *values;
2285: PetscInt Nv, v;
2287: DMLabelGetNumValues(sl->label, &Nv);
2288: DMLabelGetValueIS(sl->label, &valIS);
2289: ISGetIndices(valIS, &values);
2290: for (v = 0; v < Nv; ++v) {
2291: const PetscInt val = values[v];
2292: PetscInt size, minOrient, maxOrient;
2293: const PetscInt **perms;
2294: const PetscScalar **rots;
2296: PetscSectionSymLabelGetStratum(sym, val, &size, &minOrient, &maxOrient, &perms, &rots);
2297: PetscSectionSymLabelSetStratum(nsym, val, size, minOrient, maxOrient, PETSC_COPY_VALUES, perms, rots);
2298: }
2299: ISDestroy(&valIS);
2300: return 0;
2301: }
2303: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2304: {
2305: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
2306: DMLabel dlabel;
2308: DMLabelDistribute(sl->label, migrationSF, &dlabel);
2309: PetscSectionSymCreateLabel(PetscObjectComm((PetscObject) sym), dlabel, dsym);
2310: DMLabelDestroy(&dlabel);
2311: PetscSectionSymCopy(sym, *dsym);
2312: return 0;
2313: }
2315: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2316: {
2317: PetscSectionSym_Label *sl;
2319: PetscNewLog(sym,&sl);
2320: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2321: sym->ops->distribute = PetscSectionSymDistribute_Label;
2322: sym->ops->copy = PetscSectionSymCopy_Label;
2323: sym->ops->view = PetscSectionSymView_Label;
2324: sym->ops->destroy = PetscSectionSymDestroy_Label;
2325: sym->data = (void *) sl;
2326: return 0;
2327: }
2329: /*@
2330: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2332: Collective
2334: Input Parameters:
2335: + comm - the MPI communicator for the new symmetry
2336: - label - the label defining the strata
2338: Output Parameters:
2339: . sym - the section symmetries
2341: Level: developer
2343: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2344: @*/
2345: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2346: {
2347: DMInitializePackage();
2348: PetscSectionSymCreate(comm,sym);
2349: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2350: PetscSectionSymLabelSetLabel(*sym,label);
2351: return 0;
2352: }