Actual source code: dmlabel.c
petsc-3.12.5 2020-03-29
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: {
33: DMInitializePackage();
35: PetscHeaderCreate(*label,DMLABEL_CLASSID,"DMLabel","DMLabel","DM",comm,DMLabelDestroy,DMLabelView);
37: (*label)->numStrata = 0;
38: (*label)->defaultValue = -1;
39: (*label)->stratumValues = NULL;
40: (*label)->validIS = NULL;
41: (*label)->stratumSizes = NULL;
42: (*label)->points = NULL;
43: (*label)->ht = NULL;
44: (*label)->pStart = -1;
45: (*label)->pEnd = -1;
46: (*label)->bt = NULL;
47: PetscHMapICreate(&(*label)->hmap);
48: PetscObjectSetName((PetscObject) *label, name);
49: return(0);
50: }
52: /*
53: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
55: Not collective
57: Input parameter:
58: + label - The DMLabel
59: - v - The stratum value
61: Output parameter:
62: . label - The DMLabel with stratum in sorted list format
64: Level: developer
66: .seealso: DMLabelCreate()
67: */
68: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
69: {
70: IS is;
71: PetscInt off = 0, *pointArray, p;
74: if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
76: if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
77: PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
78: PetscMalloc1(label->stratumSizes[v], &pointArray);
79: PetscHSetIGetElems(label->ht[v], &off, pointArray);
80: PetscHSetIClear(label->ht[v]);
81: PetscSortInt(label->stratumSizes[v], pointArray);
82: if (label->bt) {
83: for (p = 0; p < label->stratumSizes[v]; ++p) {
84: const PetscInt point = pointArray[p];
85: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
86: PetscBTSet(label->bt, point - label->pStart);
87: }
88: }
89: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
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;
118: for (v = 0; v < label->numStrata; v++){
119: DMLabelMakeValid_Private(label, v);
120: }
121: return(0);
122: }
124: /*
125: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
127: Not collective
129: Input parameter:
130: + label - The DMLabel
131: - v - The stratum value
133: Output parameter:
134: . label - The DMLabel with stratum in hash format
136: Level: developer
138: .seealso: DMLabelCreate()
139: */
140: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
141: {
142: PetscInt p;
143: const PetscInt *points;
146: if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
148: if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeInvalid_Private\n", v);
149: if (label->points[v]) {
150: ISGetIndices(label->points[v], &points);
151: for (p = 0; p < label->stratumSizes[v]; ++p) {
152: PetscHSetIAdd(label->ht[v], points[p]);
153: }
154: ISRestoreIndices(label->points[v],&points);
155: ISDestroy(&(label->points[v]));
156: }
157: label->validIS[v] = PETSC_FALSE;
158: return(0);
159: }
161: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
162: #define DMLABEL_LOOKUP_THRESHOLD 16
163: #endif
165: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
166: {
167: PetscInt v;
171: *index = -1;
172: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
173: for (v = 0; v < label->numStrata; ++v)
174: if (label->stratumValues[v] == value) {*index = v; break;}
175: } else {
176: PetscHMapIGet(label->hmap, value, index);
177: }
178: #if defined(PETSC_USE_DEBUG)
179: { /* Check strata hash map consistency */
180: PetscInt len, loc = -1;
181: PetscHMapIGetSize(label->hmap, &len);
182: if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
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) {loc = v; break;}
188: }
189: if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
190: }
191: #endif
192: return(0);
193: }
195: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
196: {
197: PetscInt v;
198: PetscInt *tmpV;
199: PetscInt *tmpS;
200: PetscHSetI *tmpH, ht;
201: IS *tmpP, is;
202: PetscBool *tmpB;
203: 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: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
255: {
258: DMLabelLookupStratum(label, value, index);
259: if (*index < 0) {DMLabelNewStratum(label, value, index);}
260: return(0);
261: }
263: /*@
264: DMLabelAddStratum - Adds a new stratum value in a DMLabel
266: Input Parameter:
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;
281: DMLabelLookupAddStratum(label, value, &v);
282: return(0);
283: }
285: /*@
286: DMLabelAddStrata - Adds new stratum values in a DMLabel
288: Not collective
290: Input Parameter:
291: + label - The DMLabel
292: . numStrata - The number of stratum values
293: - stratumValues - The stratum values
295: Level: beginner
297: .seealso: DMLabelCreate(), DMLabelDestroy()
298: @*/
299: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
300: {
301: PetscInt *values, v;
307: PetscMalloc1(numStrata, &values);
308: PetscArraycpy(values, stratumValues, numStrata);
309: PetscSortRemoveDupsInt(&numStrata, values);
310: if (!label->numStrata) { /* Fast preallocation */
311: PetscInt *tmpV;
312: PetscInt *tmpS;
313: PetscHSetI *tmpH, ht;
314: IS *tmpP, is;
315: PetscBool *tmpB;
316: PetscHMapI hmap = label->hmap;
318: PetscMalloc1(numStrata, &tmpV);
319: PetscMalloc1(numStrata, &tmpS);
320: PetscMalloc1(numStrata, &tmpH);
321: PetscMalloc1(numStrata, &tmpP);
322: PetscMalloc1(numStrata, &tmpB);
323: label->numStrata = numStrata;
324: label->stratumValues = tmpV;
325: label->stratumSizes = tmpS;
326: label->ht = tmpH;
327: label->points = tmpP;
328: label->validIS = tmpB;
329: for (v = 0; v < numStrata; ++v) {
330: PetscHSetICreate(&ht);
331: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
332: PetscHMapISet(hmap, values[v], v);
333: tmpV[v] = values[v];
334: tmpS[v] = 0;
335: tmpH[v] = ht;
336: tmpP[v] = is;
337: tmpB[v] = PETSC_TRUE;
338: }
339: PetscObjectStateIncrease((PetscObject) label);
340: } else {
341: for (v = 0; v < numStrata; ++v) {
342: DMLabelAddStratum(label, values[v]);
343: }
344: }
345: PetscFree(values);
346: return(0);
347: }
349: /*@
350: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
352: Not collective
354: Input Parameter:
355: + label - The DMLabel
356: - valueIS - Index set with stratum values
358: Level: beginner
360: .seealso: DMLabelCreate(), DMLabelDestroy()
361: @*/
362: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
363: {
364: PetscInt numStrata;
365: const PetscInt *stratumValues;
371: ISGetLocalSize(valueIS, &numStrata);
372: ISGetIndices(valueIS, &stratumValues);
373: DMLabelAddStrata(label, numStrata, stratumValues);
374: return(0);
375: }
377: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
378: {
379: PetscInt v;
380: 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 [%D, %D)\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) {
399: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
400: }
401: ISRestoreIndices(label->points[v],&points);
402: }
403: }
404: PetscViewerFlush(viewer);
405: PetscViewerASCIIPopSynchronized(viewer);
406: return(0);
407: }
409: /*@C
410: DMLabelView - View the label
412: Collective on viewer
414: Input Parameters:
415: + label - The DMLabel
416: - viewer - The PetscViewer
418: Level: intermediate
420: .seealso: DMLabelCreate(), DMLabelDestroy()
421: @*/
422: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
423: {
424: PetscBool iascii;
429: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
431: if (label) {DMLabelMakeAllValid_Private(label);}
432: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
433: if (iascii) {
434: DMLabelView_Ascii(label, viewer);
435: }
436: return(0);
437: }
439: /*@
440: DMLabelReset - Destroys internal data structures in a DMLabel
442: Not collective
444: Input Parameter:
445: . label - The DMLabel
447: Level: beginner
449: .seealso: DMLabelDestroy(), DMLabelCreate()
450: @*/
451: PetscErrorCode DMLabelReset(DMLabel label)
452: {
453: PetscInt v;
458: for (v = 0; v < label->numStrata; ++v) {
459: PetscHSetIDestroy(&label->ht[v]);
460: ISDestroy(&label->points[v]);
461: }
462: label->numStrata = 0;
463: PetscFree(label->stratumValues);
464: PetscFree(label->stratumSizes);
465: PetscFree(label->ht);
466: PetscFree(label->points);
467: PetscFree(label->validIS);
468: PetscHMapIReset(label->hmap);
469: label->pStart = -1;
470: label->pEnd = -1;
471: PetscBTDestroy(&label->bt);
472: return(0);
473: }
475: /*@
476: DMLabelDestroy - Destroys a DMLabel
478: Collective on label
480: Input Parameter:
481: . label - The DMLabel
483: Level: beginner
485: .seealso: DMLabelReset(), DMLabelCreate()
486: @*/
487: PetscErrorCode DMLabelDestroy(DMLabel *label)
488: {
492: if (!*label) return(0);
494: if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
495: DMLabelReset(*label);
496: PetscHMapIDestroy(&(*label)->hmap);
497: PetscHeaderDestroy(label);
498: return(0);
499: }
501: /*@
502: DMLabelDuplicate - Duplicates a DMLabel
504: Collective on label
506: Input Parameter:
507: . label - The DMLabel
509: Output Parameter:
510: . labelnew - location to put new vector
512: Level: intermediate
514: .seealso: DMLabelCreate(), DMLabelDestroy()
515: @*/
516: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
517: {
518: const char *name;
519: PetscInt v;
524: DMLabelMakeAllValid_Private(label);
525: PetscObjectGetName((PetscObject) label, &name);
526: DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);
528: (*labelnew)->numStrata = label->numStrata;
529: (*labelnew)->defaultValue = label->defaultValue;
530: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
531: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
532: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
533: PetscMalloc1(label->numStrata, &(*labelnew)->points);
534: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
535: for (v = 0; v < label->numStrata; ++v) {
536: PetscHSetICreate(&(*labelnew)->ht[v]);
537: (*labelnew)->stratumValues[v] = label->stratumValues[v];
538: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
539: PetscObjectReference((PetscObject) (label->points[v]));
540: (*labelnew)->points[v] = label->points[v];
541: (*labelnew)->validIS[v] = PETSC_TRUE;
542: }
543: PetscHMapIDestroy(&(*labelnew)->hmap);
544: PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
545: (*labelnew)->pStart = -1;
546: (*labelnew)->pEnd = -1;
547: (*labelnew)->bt = NULL;
548: return(0);
549: }
551: /*@
552: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
554: Not collective
556: Input Parameter:
557: . label - The DMLabel
559: Level: intermediate
561: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
562: @*/
563: PetscErrorCode DMLabelComputeIndex(DMLabel label)
564: {
565: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
570: DMLabelMakeAllValid_Private(label);
571: for (v = 0; v < label->numStrata; ++v) {
572: const PetscInt *points;
573: PetscInt i;
575: ISGetIndices(label->points[v], &points);
576: for (i = 0; i < label->stratumSizes[v]; ++i) {
577: const PetscInt point = points[i];
579: pStart = PetscMin(point, pStart);
580: pEnd = PetscMax(point+1, pEnd);
581: }
582: ISRestoreIndices(label->points[v], &points);
583: }
584: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
585: label->pEnd = pEnd;
586: DMLabelCreateIndex(label, label->pStart, label->pEnd);
587: return(0);
588: }
590: /*@
591: DMLabelCreateIndex - Create an index structure for membership determination
593: Not collective
595: Input Parameters:
596: + label - The DMLabel
597: . pStart - The smallest point
598: - pEnd - The largest point + 1
600: Level: intermediate
602: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
603: @*/
604: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
605: {
606: PetscInt v;
611: DMLabelDestroyIndex(label);
612: DMLabelMakeAllValid_Private(label);
613: label->pStart = pStart;
614: label->pEnd = pEnd;
615: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
616: PetscBTCreate(pEnd - pStart, &label->bt);
617: for (v = 0; v < label->numStrata; ++v) {
618: const PetscInt *points;
619: PetscInt i;
621: ISGetIndices(label->points[v], &points);
622: for (i = 0; i < label->stratumSizes[v]; ++i) {
623: const PetscInt point = points[i];
625: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
626: PetscBTSet(label->bt, point - pStart);
627: }
628: ISRestoreIndices(label->points[v], &points);
629: }
630: return(0);
631: }
633: /*@
634: DMLabelDestroyIndex - Destroy the index structure
636: Not collective
638: Input Parameter:
639: . label - the DMLabel
641: Level: intermediate
643: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
644: @*/
645: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
646: {
651: label->pStart = -1;
652: label->pEnd = -1;
653: PetscBTDestroy(&label->bt);
654: return(0);
655: }
657: /*@
658: DMLabelGetBounds - Return the smallest and largest point in the label
660: Not collective
662: Input Parameter:
663: . label - the DMLabel
665: Output Parameters:
666: + pStart - The smallest point
667: - pEnd - The largest point + 1
669: Note: This will compute an index for the label if one does not exist.
671: Level: intermediate
673: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
674: @*/
675: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
676: {
681: if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
682: if (pStart) {
684: *pStart = label->pStart;
685: }
686: if (pEnd) {
688: *pEnd = label->pEnd;
689: }
690: return(0);
691: }
693: /*@
694: DMLabelHasValue - Determine whether a label assigns the value to any point
696: Not collective
698: Input Parameters:
699: + label - the DMLabel
700: - value - the value
702: Output Parameter:
703: . contains - Flag indicating whether the label maps this value to any point
705: Level: developer
707: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
708: @*/
709: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
710: {
711: PetscInt v;
717: DMLabelLookupStratum(label, value, &v);
718: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
719: return(0);
720: }
722: /*@
723: DMLabelHasPoint - Determine whether a label assigns a value to a point
725: Not collective
727: Input Parameters:
728: + label - the DMLabel
729: - point - the point
731: Output Parameter:
732: . contains - Flag indicating whether the label maps this point to a value
734: Note: The user must call DMLabelCreateIndex() before this function.
736: Level: developer
738: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
739: @*/
740: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
741: {
747: DMLabelMakeAllValid_Private(label);
748: #if defined(PETSC_USE_DEBUG)
749: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
750: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
751: #endif
752: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
753: return(0);
754: }
756: /*@
757: DMLabelStratumHasPoint - Return true if the stratum contains a point
759: Not collective
761: Input Parameters:
762: + label - the DMLabel
763: . value - the stratum value
764: - point - the point
766: Output Parameter:
767: . contains - true if the stratum contains the point
769: Level: intermediate
771: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
772: @*/
773: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
774: {
775: PetscInt v;
781: *contains = PETSC_FALSE;
782: DMLabelLookupStratum(label, value, &v);
783: if (v < 0) return(0);
785: if (label->validIS[v]) {
786: PetscInt i;
788: ISLocate(label->points[v], point, &i);
789: if (i >= 0) *contains = PETSC_TRUE;
790: } else {
791: PetscBool has;
793: PetscHSetIHas(label->ht[v], point, &has);
794: if (has) *contains = PETSC_TRUE;
795: }
796: return(0);
797: }
799: /*@
800: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
801: When a label is created, it is initialized to -1.
803: Not collective
805: Input parameter:
806: . label - a DMLabel object
808: Output parameter:
809: . defaultValue - the default value
811: Level: beginner
813: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
814: @*/
815: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
816: {
819: *defaultValue = label->defaultValue;
820: return(0);
821: }
823: /*@
824: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
825: When a label is created, it is initialized to -1.
827: Not collective
829: Input parameter:
830: . label - a DMLabel object
832: Output parameter:
833: . defaultValue - the default value
835: Level: beginner
837: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
838: @*/
839: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
840: {
843: label->defaultValue = defaultValue;
844: return(0);
845: }
847: /*@
848: 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())
850: Not collective
852: Input Parameters:
853: + label - the DMLabel
854: - point - the point
856: Output Parameter:
857: . value - The point value, or the default value (-1 by default)
859: Level: intermediate
861: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
862: @*/
863: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
864: {
865: PetscInt v;
871: *value = label->defaultValue;
872: for (v = 0; v < label->numStrata; ++v) {
873: if (label->validIS[v]) {
874: PetscInt i;
876: ISLocate(label->points[v], point, &i);
877: if (i >= 0) {
878: *value = label->stratumValues[v];
879: break;
880: }
881: } else {
882: PetscBool has;
884: PetscHSetIHas(label->ht[v], point, &has);
885: if (has) {
886: *value = label->stratumValues[v];
887: break;
888: }
889: }
890: }
891: return(0);
892: }
894: /*@
895: 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.
897: Not collective
899: Input Parameters:
900: + label - the DMLabel
901: . point - the point
902: - value - The point value
904: Level: intermediate
906: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
907: @*/
908: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
909: {
910: PetscInt v;
915: /* Find label value, add new entry if needed */
916: if (value == label->defaultValue) return(0);
917: DMLabelLookupAddStratum(label, value, &v);
918: /* Set key */
919: DMLabelMakeInvalid_Private(label, v);
920: PetscHSetIAdd(label->ht[v], point);
921: return(0);
922: }
924: /*@
925: DMLabelClearValue - Clear the value a label assigns to a point
927: Not collective
929: Input Parameters:
930: + label - the DMLabel
931: . point - the point
932: - value - The point value
934: Level: intermediate
936: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
937: @*/
938: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
939: {
940: PetscInt v;
945: /* Find label value */
946: DMLabelLookupStratum(label, value, &v);
947: if (v < 0) return(0);
949: if (label->bt) {
950: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
951: PetscBTClear(label->bt, point - label->pStart);
952: }
954: /* Delete key */
955: DMLabelMakeInvalid_Private(label, v);
956: PetscHSetIDel(label->ht[v], point);
957: return(0);
958: }
960: /*@
961: DMLabelInsertIS - Set all points in the IS to a value
963: Not collective
965: Input Parameters:
966: + label - the DMLabel
967: . is - the point IS
968: - value - The point value
970: Level: intermediate
972: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
973: @*/
974: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
975: {
976: PetscInt v, n, p;
977: const PetscInt *points;
978: PetscErrorCode ierr;
983: /* Find label value, add new entry if needed */
984: if (value == label->defaultValue) return(0);
985: DMLabelLookupAddStratum(label, value, &v);
986: /* Set keys */
987: DMLabelMakeInvalid_Private(label, v);
988: ISGetLocalSize(is, &n);
989: ISGetIndices(is, &points);
990: for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
991: ISRestoreIndices(is, &points);
992: return(0);
993: }
995: /*@
996: DMLabelGetNumValues - Get the number of values that the DMLabel takes
998: Not collective
1000: Input Parameter:
1001: . label - the DMLabel
1003: Output Paramater:
1004: . numValues - the number of values
1006: Level: intermediate
1008: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1009: @*/
1010: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1011: {
1015: *numValues = label->numStrata;
1016: return(0);
1017: }
1019: /*@
1020: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1022: Not collective
1024: Input Parameter:
1025: . label - the DMLabel
1027: Output Paramater:
1028: . is - the value IS
1030: Level: intermediate
1032: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1033: @*/
1034: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1035: {
1041: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1042: return(0);
1043: }
1045: /*@
1046: DMLabelHasStratum - Determine whether points exist with the given value
1048: Not collective
1050: Input Parameters:
1051: + label - the DMLabel
1052: - value - the stratum value
1054: Output Paramater:
1055: . exists - Flag saying whether points exist
1057: Level: intermediate
1059: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1060: @*/
1061: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1062: {
1063: PetscInt v;
1069: DMLabelLookupStratum(label, value, &v);
1070: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1071: return(0);
1072: }
1074: /*@
1075: DMLabelGetStratumSize - Get the size of a stratum
1077: Not collective
1079: Input Parameters:
1080: + label - the DMLabel
1081: - value - the stratum value
1083: Output Paramater:
1084: . size - The number of points in the stratum
1086: Level: intermediate
1088: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1089: @*/
1090: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1091: {
1092: PetscInt v;
1098: *size = 0;
1099: DMLabelLookupStratum(label, value, &v);
1100: if (v < 0) return(0);
1101: DMLabelMakeValid_Private(label, v);
1102: *size = label->stratumSizes[v];
1103: return(0);
1104: }
1106: /*@
1107: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1109: Not collective
1111: Input Parameters:
1112: + label - the DMLabel
1113: - value - the stratum value
1115: Output Paramaters:
1116: + start - the smallest point in the stratum
1117: - end - the largest point in the stratum
1119: Level: intermediate
1121: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1122: @*/
1123: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1124: {
1125: PetscInt v, min, max;
1132: DMLabelLookupStratum(label, value, &v);
1133: if (v < 0) return(0);
1134: DMLabelMakeValid_Private(label, v);
1135: if (label->stratumSizes[v] <= 0) return(0);
1136: ISGetMinMax(label->points[v], &min, &max);
1137: if (start) *start = min;
1138: if (end) *end = max+1;
1139: return(0);
1140: }
1142: /*@
1143: DMLabelGetStratumIS - Get an IS with the stratum points
1145: Not collective
1147: Input Parameters:
1148: + label - the DMLabel
1149: - value - the stratum value
1151: Output Paramater:
1152: . points - The stratum points
1154: Level: intermediate
1156: Notes:
1157: The output IS should be destroyed when no longer needed.
1158: Returns NULL if the stratum is empty.
1160: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1161: @*/
1162: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1163: {
1164: PetscInt v;
1170: *points = NULL;
1171: DMLabelLookupStratum(label, value, &v);
1172: if (v < 0) return(0);
1173: DMLabelMakeValid_Private(label, v);
1174: PetscObjectReference((PetscObject) label->points[v]);
1175: *points = label->points[v];
1176: return(0);
1177: }
1179: /*@
1180: DMLabelSetStratumIS - Set the stratum points using an IS
1182: Not collective
1184: Input Parameters:
1185: + label - the DMLabel
1186: . value - the stratum value
1187: - points - The stratum points
1189: Level: intermediate
1191: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1192: @*/
1193: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1194: {
1195: PetscInt v;
1201: DMLabelLookupAddStratum(label, value, &v);
1202: if (is == label->points[v]) return(0);
1203: DMLabelClearStratum(label, value);
1204: ISGetLocalSize(is, &(label->stratumSizes[v]));
1205: PetscObjectReference((PetscObject)is);
1206: ISDestroy(&(label->points[v]));
1207: label->points[v] = is;
1208: label->validIS[v] = PETSC_TRUE;
1209: PetscObjectStateIncrease((PetscObject) label);
1210: if (label->bt) {
1211: const PetscInt *points;
1212: PetscInt p;
1214: ISGetIndices(is,&points);
1215: for (p = 0; p < label->stratumSizes[v]; ++p) {
1216: const PetscInt point = points[p];
1218: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
1219: PetscBTSet(label->bt, point - label->pStart);
1220: }
1221: }
1222: return(0);
1223: }
1225: /*@
1226: DMLabelClearStratum - Remove a stratum
1228: Not collective
1230: Input Parameters:
1231: + label - the DMLabel
1232: - value - the stratum value
1234: Level: intermediate
1236: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1237: @*/
1238: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1239: {
1240: PetscInt v;
1245: DMLabelLookupStratum(label, value, &v);
1246: if (v < 0) return(0);
1247: if (label->validIS[v]) {
1248: if (label->bt) {
1249: PetscInt i;
1250: const PetscInt *points;
1252: ISGetIndices(label->points[v], &points);
1253: for (i = 0; i < label->stratumSizes[v]; ++i) {
1254: const PetscInt point = points[i];
1256: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
1257: PetscBTClear(label->bt, point - label->pStart);
1258: }
1259: ISRestoreIndices(label->points[v], &points);
1260: }
1261: label->stratumSizes[v] = 0;
1262: ISDestroy(&label->points[v]);
1263: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1264: PetscObjectSetName((PetscObject) label->points[v], "indices");
1265: PetscObjectStateIncrease((PetscObject) label);
1266: } else {
1267: PetscHSetIClear(label->ht[v]);
1268: }
1269: return(0);
1270: }
1272: /*@
1273: DMLabelFilter - Remove all points outside of [start, end)
1275: Not collective
1277: Input Parameters:
1278: + label - the DMLabel
1279: . start - the first point kept
1280: - end - one more than the last point kept
1282: Level: intermediate
1284: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1285: @*/
1286: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1287: {
1288: PetscInt v;
1293: DMLabelDestroyIndex(label);
1294: DMLabelMakeAllValid_Private(label);
1295: for (v = 0; v < label->numStrata; ++v) {
1296: ISGeneralFilter(label->points[v], start, end);
1297: }
1298: DMLabelCreateIndex(label, start, end);
1299: return(0);
1300: }
1302: /*@
1303: DMLabelPermute - Create a new label with permuted points
1305: Not collective
1307: Input Parameters:
1308: + label - the DMLabel
1309: - permutation - the point permutation
1311: Output Parameter:
1312: . labelnew - the new label containing the permuted points
1314: Level: intermediate
1316: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1317: @*/
1318: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1319: {
1320: const PetscInt *perm;
1321: PetscInt numValues, numPoints, v, q;
1322: PetscErrorCode ierr;
1327: DMLabelMakeAllValid_Private(label);
1328: DMLabelDuplicate(label, labelNew);
1329: DMLabelGetNumValues(*labelNew, &numValues);
1330: ISGetLocalSize(permutation, &numPoints);
1331: ISGetIndices(permutation, &perm);
1332: for (v = 0; v < numValues; ++v) {
1333: const PetscInt size = (*labelNew)->stratumSizes[v];
1334: const PetscInt *points;
1335: PetscInt *pointsNew;
1337: ISGetIndices((*labelNew)->points[v],&points);
1338: PetscMalloc1(size,&pointsNew);
1339: for (q = 0; q < size; ++q) {
1340: const PetscInt point = points[q];
1342: if ((point < 0) || (point >= numPoints)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [0, %D) for the remapping", point, numPoints);
1343: pointsNew[q] = perm[point];
1344: }
1345: ISRestoreIndices((*labelNew)->points[v],&points);
1346: PetscSortInt(size, pointsNew);
1347: ISDestroy(&((*labelNew)->points[v]));
1348: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1349: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1350: PetscFree(pointsNew);
1351: } else {
1352: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1353: }
1354: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1355: }
1356: ISRestoreIndices(permutation, &perm);
1357: if (label->bt) {
1358: PetscBTDestroy(&label->bt);
1359: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1360: }
1361: return(0);
1362: }
1364: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1365: {
1366: MPI_Comm comm;
1367: PetscInt s, l, nroots, nleaves, dof, offset, size;
1368: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1369: PetscSection rootSection;
1370: PetscSF labelSF;
1374: if (label) {DMLabelMakeAllValid_Private(label);}
1375: PetscObjectGetComm((PetscObject)sf, &comm);
1376: /* Build a section of stratum values per point, generate the according SF
1377: and distribute point-wise stratum values to leaves. */
1378: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1379: PetscSectionCreate(comm, &rootSection);
1380: PetscSectionSetChart(rootSection, 0, nroots);
1381: if (label) {
1382: for (s = 0; s < label->numStrata; ++s) {
1383: const PetscInt *points;
1385: ISGetIndices(label->points[s], &points);
1386: for (l = 0; l < label->stratumSizes[s]; l++) {
1387: PetscSectionGetDof(rootSection, points[l], &dof);
1388: PetscSectionSetDof(rootSection, points[l], dof+1);
1389: }
1390: ISRestoreIndices(label->points[s], &points);
1391: }
1392: }
1393: PetscSectionSetUp(rootSection);
1394: /* Create a point-wise array of stratum values */
1395: PetscSectionGetStorageSize(rootSection, &size);
1396: PetscMalloc1(size, &rootStrata);
1397: PetscCalloc1(nroots, &rootIdx);
1398: if (label) {
1399: for (s = 0; s < label->numStrata; ++s) {
1400: const PetscInt *points;
1402: ISGetIndices(label->points[s], &points);
1403: for (l = 0; l < label->stratumSizes[s]; l++) {
1404: const PetscInt p = points[l];
1405: PetscSectionGetOffset(rootSection, p, &offset);
1406: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1407: }
1408: ISRestoreIndices(label->points[s], &points);
1409: }
1410: }
1411: /* Build SF that maps label points to remote processes */
1412: PetscSectionCreate(comm, leafSection);
1413: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1414: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1415: PetscFree(remoteOffsets);
1416: /* Send the strata for each point over the derived SF */
1417: PetscSectionGetStorageSize(*leafSection, &size);
1418: PetscMalloc1(size, leafStrata);
1419: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1420: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1421: /* Clean up */
1422: PetscFree(rootStrata);
1423: PetscFree(rootIdx);
1424: PetscSectionDestroy(&rootSection);
1425: PetscSFDestroy(&labelSF);
1426: return(0);
1427: }
1429: /*@
1430: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1432: Collective on sf
1434: Input Parameters:
1435: + label - the DMLabel
1436: - sf - the map from old to new distribution
1438: Output Parameter:
1439: . labelnew - the new redistributed label
1441: Level: intermediate
1443: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1444: @*/
1445: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1446: {
1447: MPI_Comm comm;
1448: PetscSection leafSection;
1449: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1450: PetscInt *leafStrata, *strataIdx;
1451: PetscInt **points;
1452: const char *lname = NULL;
1453: char *name;
1454: PetscInt nameSize;
1455: PetscHSetI stratumHash;
1456: size_t len = 0;
1457: PetscMPIInt rank;
1462: if (label) {
1464: DMLabelMakeAllValid_Private(label);
1465: }
1466: PetscObjectGetComm((PetscObject)sf, &comm);
1467: MPI_Comm_rank(comm, &rank);
1468: /* Bcast name */
1469: if (!rank) {
1470: PetscObjectGetName((PetscObject) label, &lname);
1471: PetscStrlen(lname, &len);
1472: }
1473: nameSize = len;
1474: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1475: PetscMalloc1(nameSize+1, &name);
1476: if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1477: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1478: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1479: PetscFree(name);
1480: /* Bcast defaultValue */
1481: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1482: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1483: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1484: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1485: /* Determine received stratum values and initialise new label*/
1486: PetscHSetICreate(&stratumHash);
1487: PetscSectionGetStorageSize(leafSection, &size);
1488: for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1489: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1490: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1491: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1492: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1493: /* Turn leafStrata into indices rather than stratum values */
1494: offset = 0;
1495: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1496: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1497: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1498: PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1499: }
1500: for (p = 0; p < size; ++p) {
1501: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1502: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1503: }
1504: }
1505: /* Rebuild the point strata on the receiver */
1506: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1507: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1508: for (p=pStart; p<pEnd; p++) {
1509: PetscSectionGetDof(leafSection, p, &dof);
1510: PetscSectionGetOffset(leafSection, p, &offset);
1511: for (s=0; s<dof; s++) {
1512: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1513: }
1514: }
1515: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1516: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1517: PetscMalloc1((*labelNew)->numStrata,&points);
1518: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1519: PetscHSetICreate(&(*labelNew)->ht[s]);
1520: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1521: }
1522: /* Insert points into new strata */
1523: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1524: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1525: for (p=pStart; p<pEnd; p++) {
1526: PetscSectionGetDof(leafSection, p, &dof);
1527: PetscSectionGetOffset(leafSection, p, &offset);
1528: for (s=0; s<dof; s++) {
1529: stratum = leafStrata[offset+s];
1530: points[stratum][strataIdx[stratum]++] = p;
1531: }
1532: }
1533: for (s = 0; s < (*labelNew)->numStrata; s++) {
1534: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1535: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1536: }
1537: PetscFree(points);
1538: PetscHSetIDestroy(&stratumHash);
1539: PetscFree(leafStrata);
1540: PetscFree(strataIdx);
1541: PetscSectionDestroy(&leafSection);
1542: return(0);
1543: }
1545: /*@
1546: DMLabelGather - Gather all label values from leafs into roots
1548: Collective on sf
1550: Input Parameters:
1551: + label - the DMLabel
1552: - sf - the Star Forest point communication map
1554: Output Parameters:
1555: . labelNew - the new DMLabel with localised leaf values
1557: Level: developer
1559: Note: This is the inverse operation to DMLabelDistribute.
1561: .seealso: DMLabelDistribute()
1562: @*/
1563: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1564: {
1565: MPI_Comm comm;
1566: PetscSection rootSection;
1567: PetscSF sfLabel;
1568: PetscSFNode *rootPoints, *leafPoints;
1569: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1570: const PetscInt *rootDegree, *ilocal;
1571: PetscInt *rootStrata;
1572: const char *lname;
1573: char *name;
1574: PetscInt nameSize;
1575: size_t len = 0;
1576: PetscMPIInt rank, size;
1582: PetscObjectGetComm((PetscObject)sf, &comm);
1583: MPI_Comm_rank(comm, &rank);
1584: MPI_Comm_size(comm, &size);
1585: /* Bcast name */
1586: if (!rank) {
1587: PetscObjectGetName((PetscObject) label, &lname);
1588: PetscStrlen(lname, &len);
1589: }
1590: nameSize = len;
1591: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1592: PetscMalloc1(nameSize+1, &name);
1593: if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1594: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1595: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1596: PetscFree(name);
1597: /* Gather rank/index pairs of leaves into local roots to build
1598: an inverse, multi-rooted SF. Note that this ignores local leaf
1599: indexing due to the use of the multiSF in PetscSFGather. */
1600: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1601: PetscMalloc1(nroots, &leafPoints);
1602: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1603: for (p = 0; p < nleaves; p++) {
1604: PetscInt ilp = ilocal ? ilocal[p] : p;
1606: leafPoints[ilp].index = ilp;
1607: leafPoints[ilp].rank = rank;
1608: }
1609: PetscSFComputeDegreeBegin(sf, &rootDegree);
1610: PetscSFComputeDegreeEnd(sf, &rootDegree);
1611: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1612: PetscMalloc1(nmultiroots, &rootPoints);
1613: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1614: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1615: PetscSFCreate(comm,& sfLabel);
1616: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1617: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1618: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1619: /* Rebuild the point strata on the receiver */
1620: for (p = 0, idx = 0; p < nroots; p++) {
1621: for (d = 0; d < rootDegree[p]; d++) {
1622: PetscSectionGetDof(rootSection, idx+d, &dof);
1623: PetscSectionGetOffset(rootSection, idx+d, &offset);
1624: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1625: }
1626: idx += rootDegree[p];
1627: }
1628: PetscFree(leafPoints);
1629: PetscFree(rootStrata);
1630: PetscSectionDestroy(&rootSection);
1631: PetscSFDestroy(&sfLabel);
1632: return(0);
1633: }
1635: /*@
1636: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1638: Not collective
1640: Input Parameter:
1641: . label - the DMLabel
1643: Output Parameters:
1644: + section - the section giving offsets for each stratum
1645: - is - An IS containing all the label points
1647: Level: developer
1649: .seealso: DMLabelDistribute()
1650: @*/
1651: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1652: {
1653: IS vIS;
1654: const PetscInt *values;
1655: PetscInt *points;
1656: PetscInt nV, vS = 0, vE = 0, v, N;
1657: PetscErrorCode ierr;
1661: DMLabelGetNumValues(label, &nV);
1662: DMLabelGetValueIS(label, &vIS);
1663: ISGetIndices(vIS, &values);
1664: if (nV) {vS = values[0]; vE = values[0]+1;}
1665: for (v = 1; v < nV; ++v) {
1666: vS = PetscMin(vS, values[v]);
1667: vE = PetscMax(vE, values[v]+1);
1668: }
1669: PetscSectionCreate(PETSC_COMM_SELF, section);
1670: PetscSectionSetChart(*section, vS, vE);
1671: for (v = 0; v < nV; ++v) {
1672: PetscInt n;
1674: DMLabelGetStratumSize(label, values[v], &n);
1675: PetscSectionSetDof(*section, values[v], n);
1676: }
1677: PetscSectionSetUp(*section);
1678: PetscSectionGetStorageSize(*section, &N);
1679: PetscMalloc1(N, &points);
1680: for (v = 0; v < nV; ++v) {
1681: IS is;
1682: const PetscInt *spoints;
1683: PetscInt dof, off, p;
1685: PetscSectionGetDof(*section, values[v], &dof);
1686: PetscSectionGetOffset(*section, values[v], &off);
1687: DMLabelGetStratumIS(label, values[v], &is);
1688: ISGetIndices(is, &spoints);
1689: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1690: ISRestoreIndices(is, &spoints);
1691: ISDestroy(&is);
1692: }
1693: ISRestoreIndices(vIS, &values);
1694: ISDestroy(&vIS);
1695: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1696: return(0);
1697: }
1699: /*@
1700: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1701: the local section and an SF describing the section point overlap.
1703: Collective on sf
1705: Input Parameters:
1706: + s - The PetscSection for the local field layout
1707: . sf - The SF describing parallel layout of the section points
1708: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1709: . label - The label specifying the points
1710: - labelValue - The label stratum specifying the points
1712: Output Parameter:
1713: . gsection - The PetscSection for the global field layout
1715: Note: This gives negative sizes and offsets to points not owned by this process
1717: Level: developer
1719: .seealso: PetscSectionCreate()
1720: @*/
1721: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1722: {
1723: PetscInt *neg = NULL, *tmpOff = NULL;
1724: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1731: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1732: PetscSectionGetChart(s, &pStart, &pEnd);
1733: PetscSectionSetChart(*gsection, pStart, pEnd);
1734: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1735: if (nroots >= 0) {
1736: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1737: PetscCalloc1(nroots, &neg);
1738: if (nroots > pEnd-pStart) {
1739: PetscCalloc1(nroots, &tmpOff);
1740: } else {
1741: tmpOff = &(*gsection)->atlasDof[-pStart];
1742: }
1743: }
1744: /* Mark ghost points with negative dof */
1745: for (p = pStart; p < pEnd; ++p) {
1746: PetscInt value;
1748: DMLabelGetValue(label, p, &value);
1749: if (value != labelValue) continue;
1750: PetscSectionGetDof(s, p, &dof);
1751: PetscSectionSetDof(*gsection, p, dof);
1752: PetscSectionGetConstraintDof(s, p, &cdof);
1753: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1754: if (neg) neg[p] = -(dof+1);
1755: }
1756: PetscSectionSetUpBC(*gsection);
1757: if (nroots >= 0) {
1758: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1759: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1760: if (nroots > pEnd-pStart) {
1761: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1762: }
1763: }
1764: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1765: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1766: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1767: (*gsection)->atlasOff[p] = off;
1768: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1769: }
1770: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1771: globalOff -= off;
1772: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1773: (*gsection)->atlasOff[p] += globalOff;
1774: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1775: }
1776: /* Put in negative offsets for ghost points */
1777: if (nroots >= 0) {
1778: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1779: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1780: if (nroots > pEnd-pStart) {
1781: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1782: }
1783: }
1784: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1785: PetscFree(neg);
1786: return(0);
1787: }
1789: typedef struct _n_PetscSectionSym_Label
1790: {
1791: DMLabel label;
1792: PetscCopyMode *modes;
1793: PetscInt *sizes;
1794: const PetscInt ***perms;
1795: const PetscScalar ***rots;
1796: PetscInt (*minMaxOrients)[2];
1797: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1798: } PetscSectionSym_Label;
1800: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1801: {
1802: PetscInt i, j;
1803: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1804: PetscErrorCode ierr;
1807: for (i = 0; i <= sl->numStrata; i++) {
1808: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1809: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1810: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1811: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1812: }
1813: if (sl->perms[i]) {
1814: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1816: PetscFree(perms);
1817: }
1818: if (sl->rots[i]) {
1819: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1821: PetscFree(rots);
1822: }
1823: }
1824: }
1825: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1826: DMLabelDestroy(&sl->label);
1827: sl->numStrata = 0;
1828: return(0);
1829: }
1831: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1832: {
1836: PetscSectionSymLabelReset(sym);
1837: PetscFree(sym->data);
1838: return(0);
1839: }
1841: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1842: {
1843: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1844: PetscBool isAscii;
1845: DMLabel label = sl->label;
1846: const char *name;
1847: PetscErrorCode ierr;
1850: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1851: if (isAscii) {
1852: PetscInt i, j, k;
1853: PetscViewerFormat format;
1855: PetscViewerGetFormat(viewer,&format);
1856: if (label) {
1857: PetscViewerGetFormat(viewer,&format);
1858: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1859: PetscViewerASCIIPushTab(viewer);
1860: DMLabelView(label, viewer);
1861: PetscViewerASCIIPopTab(viewer);
1862: } else {
1863: PetscObjectGetName((PetscObject) sl->label, &name);
1864: PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);
1865: }
1866: } else {
1867: PetscViewerASCIIPrintf(viewer, "No label given\n");
1868: }
1869: PetscViewerASCIIPushTab(viewer);
1870: for (i = 0; i <= sl->numStrata; i++) {
1871: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1873: if (!(sl->perms[i] || sl->rots[i])) {
1874: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1875: } else {
1876: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1877: PetscViewerASCIIPushTab(viewer);
1878: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1879: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1880: PetscViewerASCIIPushTab(viewer);
1881: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1882: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1883: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1884: } else {
1885: PetscInt tab;
1887: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1888: PetscViewerASCIIPushTab(viewer);
1889: PetscViewerASCIIGetTab(viewer,&tab);
1890: if (sl->perms[i] && sl->perms[i][j]) {
1891: PetscViewerASCIIPrintf(viewer,"Permutation:");
1892: PetscViewerASCIISetTab(viewer,0);
1893: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1894: PetscViewerASCIIPrintf(viewer,"\n");
1895: PetscViewerASCIISetTab(viewer,tab);
1896: }
1897: if (sl->rots[i] && sl->rots[i][j]) {
1898: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1899: PetscViewerASCIISetTab(viewer,0);
1900: #if defined(PETSC_USE_COMPLEX)
1901: 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]));}
1902: #else
1903: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1904: #endif
1905: PetscViewerASCIIPrintf(viewer,"\n");
1906: PetscViewerASCIISetTab(viewer,tab);
1907: }
1908: PetscViewerASCIIPopTab(viewer);
1909: }
1910: }
1911: PetscViewerASCIIPopTab(viewer);
1912: }
1913: PetscViewerASCIIPopTab(viewer);
1914: }
1915: }
1916: PetscViewerASCIIPopTab(viewer);
1917: }
1918: return(0);
1919: }
1921: /*@
1922: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1924: Logically collective on sym
1926: Input parameters:
1927: + sym - the section symmetries
1928: - label - the DMLabel describing the types of points
1930: Level: developer:
1932: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1933: @*/
1934: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1935: {
1936: PetscSectionSym_Label *sl;
1937: PetscErrorCode ierr;
1941: sl = (PetscSectionSym_Label *) sym->data;
1942: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1943: if (label) {
1944: sl->label = label;
1945: PetscObjectReference((PetscObject) label);
1946: DMLabelGetNumValues(label,&sl->numStrata);
1947: 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);
1948: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1949: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1950: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1951: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1952: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1953: }
1954: return(0);
1955: }
1957: /*@C
1958: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1960: Logically collective on sym
1962: InputParameters:
1963: + sym - the section symmetries
1964: . stratum - the stratum value in the label that we are assigning symmetries for
1965: . size - the number of dofs for points in the stratum of the label
1966: . minOrient - the smallest orientation for a point in this stratum
1967: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1968: . mode - how sym should copy the perms and rots arrays
1969: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
1970: - 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
1972: Level: developer
1974: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1975: @*/
1976: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1977: {
1978: PetscSectionSym_Label *sl;
1979: const char *name;
1980: PetscInt i, j, k;
1981: PetscErrorCode ierr;
1985: sl = (PetscSectionSym_Label *) sym->data;
1986: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1987: for (i = 0; i <= sl->numStrata; i++) {
1988: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
1990: if (stratum == value) break;
1991: }
1992: PetscObjectGetName((PetscObject) sl->label, &name);
1993: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
1994: sl->sizes[i] = size;
1995: sl->modes[i] = mode;
1996: sl->minMaxOrients[i][0] = minOrient;
1997: sl->minMaxOrients[i][1] = maxOrient;
1998: if (mode == PETSC_COPY_VALUES) {
1999: if (perms) {
2000: PetscInt **ownPerms;
2002: PetscCalloc1(maxOrient - minOrient,&ownPerms);
2003: for (j = 0; j < maxOrient-minOrient; j++) {
2004: if (perms[j]) {
2005: PetscMalloc1(size,&ownPerms[j]);
2006: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2007: }
2008: }
2009: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2010: }
2011: if (rots) {
2012: PetscScalar **ownRots;
2014: PetscCalloc1(maxOrient - minOrient,&ownRots);
2015: for (j = 0; j < maxOrient-minOrient; j++) {
2016: if (rots[j]) {
2017: PetscMalloc1(size,&ownRots[j]);
2018: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2019: }
2020: }
2021: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2022: }
2023: } else {
2024: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2025: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2026: }
2027: return(0);
2028: }
2030: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2031: {
2032: PetscInt i, j, numStrata;
2033: PetscSectionSym_Label *sl;
2034: DMLabel label;
2035: PetscErrorCode ierr;
2038: sl = (PetscSectionSym_Label *) sym->data;
2039: numStrata = sl->numStrata;
2040: label = sl->label;
2041: for (i = 0; i < numPoints; i++) {
2042: PetscInt point = points[2*i];
2043: PetscInt ornt = points[2*i+1];
2045: for (j = 0; j < numStrata; j++) {
2046: if (label->validIS[j]) {
2047: PetscInt k;
2049: ISLocate(label->points[j],point,&k);
2050: if (k >= 0) break;
2051: } else {
2052: PetscBool has;
2054: PetscHSetIHas(label->ht[j], point, &has);
2055: if (has) break;
2056: }
2057: }
2058: if ((sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) && (ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1])) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D orientation %D not in range [%D, %D) for stratum %D",point,ornt,sl->minMaxOrients[j][0],sl->minMaxOrients[j][1],j < numStrata ? label->stratumValues[j] : label->defaultValue);
2059: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2060: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2061: }
2062: return(0);
2063: }
2065: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2066: {
2067: PetscSectionSym_Label *sl;
2068: PetscErrorCode ierr;
2071: PetscNewLog(sym,&sl);
2072: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2073: sym->ops->view = PetscSectionSymView_Label;
2074: sym->ops->destroy = PetscSectionSymDestroy_Label;
2075: sym->data = (void *) sl;
2076: return(0);
2077: }
2079: /*@
2080: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2082: Collective
2084: Input Parameters:
2085: + comm - the MPI communicator for the new symmetry
2086: - label - the label defining the strata
2088: Output Parameters:
2089: . sym - the section symmetries
2091: Level: developer
2093: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2094: @*/
2095: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2096: {
2097: PetscErrorCode ierr;
2100: DMInitializePackage();
2101: PetscSectionSymCreate(comm,sym);
2102: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2103: PetscSectionSymLabelSetLabel(*sym,label);
2104: return(0);
2105: }