Actual source code: dmlabel.c
petsc-3.14.6 2021-03-30
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: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
90: ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
91: PetscFree(pointArray);
92: } else {
93: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
94: }
95: PetscObjectSetName((PetscObject) is, "indices");
96: label->points[v] = is;
97: label->validIS[v] = PETSC_TRUE;
98: PetscObjectStateIncrease((PetscObject) label);
99: return(0);
100: }
102: /*
103: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
105: Not collective
107: Input parameter:
108: . label - The DMLabel
110: Output parameter:
111: . label - The DMLabel with all strata in sorted list format
113: Level: developer
115: .seealso: DMLabelCreate()
116: */
117: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
118: {
119: PetscInt v;
123: for (v = 0; v < label->numStrata; v++){
124: DMLabelMakeValid_Private(label, v);
125: }
126: return(0);
127: }
129: /*
130: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
132: Not collective
134: Input parameter:
135: + label - The DMLabel
136: - v - The stratum value
138: Output parameter:
139: . label - The DMLabel with stratum in hash format
141: Level: developer
143: .seealso: DMLabelCreate()
144: */
145: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
146: {
147: PetscInt p;
148: const PetscInt *points;
151: if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
153: 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);
154: if (label->points[v]) {
155: ISGetIndices(label->points[v], &points);
156: for (p = 0; p < label->stratumSizes[v]; ++p) {
157: PetscHSetIAdd(label->ht[v], points[p]);
158: }
159: ISRestoreIndices(label->points[v],&points);
160: ISDestroy(&(label->points[v]));
161: }
162: label->validIS[v] = PETSC_FALSE;
163: return(0);
164: }
166: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
167: #define DMLABEL_LOOKUP_THRESHOLD 16
168: #endif
170: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
171: {
172: PetscInt v;
176: *index = -1;
177: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
178: for (v = 0; v < label->numStrata; ++v)
179: if (label->stratumValues[v] == value) {*index = v; break;}
180: } else {
181: PetscHMapIGet(label->hmap, value, index);
182: }
183: if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
184: PetscInt len, loc = -1;
185: PetscHMapIGetSize(label->hmap, &len);
186: if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
187: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
188: PetscHMapIGet(label->hmap, value, &loc);
189: } else {
190: for (v = 0; v < label->numStrata; ++v)
191: if (label->stratumValues[v] == value) {loc = v; break;}
192: }
193: if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
194: }
195: return(0);
196: }
198: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
199: {
200: PetscInt v;
201: PetscInt *tmpV;
202: PetscInt *tmpS;
203: PetscHSetI *tmpH, ht;
204: IS *tmpP, is;
205: PetscBool *tmpB;
206: PetscHMapI hmap = label->hmap;
210: v = label->numStrata;
211: tmpV = label->stratumValues;
212: tmpS = label->stratumSizes;
213: tmpH = label->ht;
214: tmpP = label->points;
215: tmpB = label->validIS;
216: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
217: PetscInt *oldV = tmpV;
218: PetscInt *oldS = tmpS;
219: PetscHSetI *oldH = tmpH;
220: IS *oldP = tmpP;
221: PetscBool *oldB = tmpB;
222: PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
223: PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
224: PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
225: PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
226: PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
227: PetscArraycpy(tmpV, oldV, v);
228: PetscArraycpy(tmpS, oldS, v);
229: PetscArraycpy(tmpH, oldH, v);
230: PetscArraycpy(tmpP, oldP, v);
231: PetscArraycpy(tmpB, oldB, v);
232: PetscFree(oldV);
233: PetscFree(oldS);
234: PetscFree(oldH);
235: PetscFree(oldP);
236: PetscFree(oldB);
237: }
238: label->numStrata = v+1;
239: label->stratumValues = tmpV;
240: label->stratumSizes = tmpS;
241: label->ht = tmpH;
242: label->points = tmpP;
243: label->validIS = tmpB;
244: PetscHSetICreate(&ht);
245: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
246: PetscHMapISet(hmap, value, v);
247: tmpV[v] = value;
248: tmpS[v] = 0;
249: tmpH[v] = ht;
250: tmpP[v] = is;
251: tmpB[v] = PETSC_TRUE;
252: PetscObjectStateIncrease((PetscObject) label);
253: *index = v;
254: return(0);
255: }
257: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
258: {
261: DMLabelLookupStratum(label, value, index);
262: if (*index < 0) {DMLabelNewStratum(label, value, index);}
263: return(0);
264: }
266: /*@
267: DMLabelAddStratum - Adds a new stratum value in a DMLabel
269: Input Parameter:
270: + label - The DMLabel
271: - value - The stratum value
273: Level: beginner
275: .seealso: DMLabelCreate(), DMLabelDestroy()
276: @*/
277: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
278: {
279: PetscInt v;
284: DMLabelLookupAddStratum(label, value, &v);
285: return(0);
286: }
288: /*@
289: DMLabelAddStrata - Adds new stratum values in a DMLabel
291: Not collective
293: Input Parameter:
294: + label - The DMLabel
295: . numStrata - The number of stratum values
296: - stratumValues - The stratum values
298: Level: beginner
300: .seealso: DMLabelCreate(), DMLabelDestroy()
301: @*/
302: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
303: {
304: PetscInt *values, v;
310: PetscMalloc1(numStrata, &values);
311: PetscArraycpy(values, stratumValues, numStrata);
312: PetscSortRemoveDupsInt(&numStrata, values);
313: if (!label->numStrata) { /* Fast preallocation */
314: PetscInt *tmpV;
315: PetscInt *tmpS;
316: PetscHSetI *tmpH, ht;
317: IS *tmpP, is;
318: PetscBool *tmpB;
319: PetscHMapI hmap = label->hmap;
321: PetscMalloc1(numStrata, &tmpV);
322: PetscMalloc1(numStrata, &tmpS);
323: PetscMalloc1(numStrata, &tmpH);
324: PetscMalloc1(numStrata, &tmpP);
325: PetscMalloc1(numStrata, &tmpB);
326: label->numStrata = numStrata;
327: label->stratumValues = tmpV;
328: label->stratumSizes = tmpS;
329: label->ht = tmpH;
330: label->points = tmpP;
331: label->validIS = tmpB;
332: for (v = 0; v < numStrata; ++v) {
333: PetscHSetICreate(&ht);
334: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
335: PetscHMapISet(hmap, values[v], v);
336: tmpV[v] = values[v];
337: tmpS[v] = 0;
338: tmpH[v] = ht;
339: tmpP[v] = is;
340: tmpB[v] = PETSC_TRUE;
341: }
342: PetscObjectStateIncrease((PetscObject) label);
343: } else {
344: for (v = 0; v < numStrata; ++v) {
345: DMLabelAddStratum(label, values[v]);
346: }
347: }
348: PetscFree(values);
349: return(0);
350: }
352: /*@
353: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
355: Not collective
357: Input Parameter:
358: + label - The DMLabel
359: - valueIS - Index set with stratum values
361: Level: beginner
363: .seealso: DMLabelCreate(), DMLabelDestroy()
364: @*/
365: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
366: {
367: PetscInt numStrata;
368: const PetscInt *stratumValues;
374: ISGetLocalSize(valueIS, &numStrata);
375: ISGetIndices(valueIS, &stratumValues);
376: DMLabelAddStrata(label, numStrata, stratumValues);
377: return(0);
378: }
380: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
381: {
382: PetscInt v;
383: PetscMPIInt rank;
387: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
388: PetscViewerASCIIPushSynchronized(viewer);
389: if (label) {
390: const char *name;
392: PetscObjectGetName((PetscObject) label, &name);
393: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
394: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
395: for (v = 0; v < label->numStrata; ++v) {
396: const PetscInt value = label->stratumValues[v];
397: const PetscInt *points;
398: PetscInt p;
400: ISGetIndices(label->points[v], &points);
401: for (p = 0; p < label->stratumSizes[v]; ++p) {
402: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
403: }
404: ISRestoreIndices(label->points[v],&points);
405: }
406: }
407: PetscViewerFlush(viewer);
408: PetscViewerASCIIPopSynchronized(viewer);
409: return(0);
410: }
412: /*@C
413: DMLabelView - View the label
415: Collective on viewer
417: Input Parameters:
418: + label - The DMLabel
419: - viewer - The PetscViewer
421: Level: intermediate
423: .seealso: DMLabelCreate(), DMLabelDestroy()
424: @*/
425: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
426: {
427: PetscBool iascii;
432: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
434: if (label) {DMLabelMakeAllValid_Private(label);}
435: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
436: if (iascii) {
437: DMLabelView_Ascii(label, viewer);
438: }
439: return(0);
440: }
442: /*@
443: DMLabelReset - Destroys internal data structures in a DMLabel
445: Not collective
447: Input Parameter:
448: . label - The DMLabel
450: Level: beginner
452: .seealso: DMLabelDestroy(), DMLabelCreate()
453: @*/
454: PetscErrorCode DMLabelReset(DMLabel label)
455: {
456: PetscInt v;
461: for (v = 0; v < label->numStrata; ++v) {
462: PetscHSetIDestroy(&label->ht[v]);
463: ISDestroy(&label->points[v]);
464: }
465: label->numStrata = 0;
466: PetscFree(label->stratumValues);
467: PetscFree(label->stratumSizes);
468: PetscFree(label->ht);
469: PetscFree(label->points);
470: PetscFree(label->validIS);
471: PetscHMapIReset(label->hmap);
472: label->pStart = -1;
473: label->pEnd = -1;
474: PetscBTDestroy(&label->bt);
475: return(0);
476: }
478: /*@
479: DMLabelDestroy - Destroys a DMLabel
481: Collective on label
483: Input Parameter:
484: . label - The DMLabel
486: Level: beginner
488: .seealso: DMLabelReset(), DMLabelCreate()
489: @*/
490: PetscErrorCode DMLabelDestroy(DMLabel *label)
491: {
495: if (!*label) return(0);
497: if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
498: DMLabelReset(*label);
499: PetscHMapIDestroy(&(*label)->hmap);
500: PetscHeaderDestroy(label);
501: return(0);
502: }
504: /*@
505: DMLabelDuplicate - Duplicates a DMLabel
507: Collective on label
509: Input Parameter:
510: . label - The DMLabel
512: Output Parameter:
513: . labelnew - location to put new vector
515: Level: intermediate
517: .seealso: DMLabelCreate(), DMLabelDestroy()
518: @*/
519: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
520: {
521: const char *name;
522: PetscInt v;
527: DMLabelMakeAllValid_Private(label);
528: PetscObjectGetName((PetscObject) label, &name);
529: DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);
531: (*labelnew)->numStrata = label->numStrata;
532: (*labelnew)->defaultValue = label->defaultValue;
533: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
534: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
535: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
536: PetscMalloc1(label->numStrata, &(*labelnew)->points);
537: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
538: for (v = 0; v < label->numStrata; ++v) {
539: PetscHSetICreate(&(*labelnew)->ht[v]);
540: (*labelnew)->stratumValues[v] = label->stratumValues[v];
541: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
542: PetscObjectReference((PetscObject) (label->points[v]));
543: (*labelnew)->points[v] = label->points[v];
544: (*labelnew)->validIS[v] = PETSC_TRUE;
545: }
546: PetscHMapIDestroy(&(*labelnew)->hmap);
547: PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
548: (*labelnew)->pStart = -1;
549: (*labelnew)->pEnd = -1;
550: (*labelnew)->bt = NULL;
551: return(0);
552: }
554: /*@
555: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
557: Not collective
559: Input Parameter:
560: . label - The DMLabel
562: Level: intermediate
564: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
565: @*/
566: PetscErrorCode DMLabelComputeIndex(DMLabel label)
567: {
568: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
573: DMLabelMakeAllValid_Private(label);
574: for (v = 0; v < label->numStrata; ++v) {
575: const PetscInt *points;
576: PetscInt i;
578: ISGetIndices(label->points[v], &points);
579: for (i = 0; i < label->stratumSizes[v]; ++i) {
580: const PetscInt point = points[i];
582: pStart = PetscMin(point, pStart);
583: pEnd = PetscMax(point+1, pEnd);
584: }
585: ISRestoreIndices(label->points[v], &points);
586: }
587: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
588: label->pEnd = pEnd;
589: DMLabelCreateIndex(label, label->pStart, label->pEnd);
590: return(0);
591: }
593: /*@
594: DMLabelCreateIndex - Create an index structure for membership determination
596: Not collective
598: Input Parameters:
599: + label - The DMLabel
600: . pStart - The smallest point
601: - pEnd - The largest point + 1
603: Level: intermediate
605: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
606: @*/
607: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
608: {
609: PetscInt v;
614: DMLabelDestroyIndex(label);
615: DMLabelMakeAllValid_Private(label);
616: label->pStart = pStart;
617: label->pEnd = pEnd;
618: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
619: PetscBTCreate(pEnd - pStart, &label->bt);
620: for (v = 0; v < label->numStrata; ++v) {
621: const PetscInt *points;
622: PetscInt i;
624: ISGetIndices(label->points[v], &points);
625: for (i = 0; i < label->stratumSizes[v]; ++i) {
626: const PetscInt point = points[i];
628: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
629: PetscBTSet(label->bt, point - pStart);
630: }
631: ISRestoreIndices(label->points[v], &points);
632: }
633: return(0);
634: }
636: /*@
637: DMLabelDestroyIndex - Destroy the index structure
639: Not collective
641: Input Parameter:
642: . label - the DMLabel
644: Level: intermediate
646: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
647: @*/
648: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
649: {
654: label->pStart = -1;
655: label->pEnd = -1;
656: PetscBTDestroy(&label->bt);
657: return(0);
658: }
660: /*@
661: DMLabelGetBounds - Return the smallest and largest point in the label
663: Not collective
665: Input Parameter:
666: . label - the DMLabel
668: Output Parameters:
669: + pStart - The smallest point
670: - pEnd - The largest point + 1
672: Note: This will compute an index for the label if one does not exist.
674: Level: intermediate
676: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
677: @*/
678: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
679: {
684: if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
685: if (pStart) {
687: *pStart = label->pStart;
688: }
689: if (pEnd) {
691: *pEnd = label->pEnd;
692: }
693: return(0);
694: }
696: /*@
697: DMLabelHasValue - Determine whether a label assigns the value to any point
699: Not collective
701: Input Parameters:
702: + label - the DMLabel
703: - value - the value
705: Output Parameter:
706: . contains - Flag indicating whether the label maps this value to any point
708: Level: developer
710: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
711: @*/
712: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
713: {
714: PetscInt v;
720: DMLabelLookupStratum(label, value, &v);
721: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
722: return(0);
723: }
725: /*@
726: DMLabelHasPoint - Determine whether a label assigns a value to a point
728: Not collective
730: Input Parameters:
731: + label - the DMLabel
732: - point - the point
734: Output Parameter:
735: . contains - Flag indicating whether the label maps this point to a value
737: Note: The user must call DMLabelCreateIndex() before this function.
739: Level: developer
741: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
742: @*/
743: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
744: {
750: DMLabelMakeAllValid_Private(label);
751: if (PetscDefined(USE_DEBUG)) {
752: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
753: 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);
754: }
755: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
756: return(0);
757: }
759: /*@
760: DMLabelStratumHasPoint - Return true if the stratum contains a point
762: Not collective
764: Input Parameters:
765: + label - the DMLabel
766: . value - the stratum value
767: - point - the point
769: Output Parameter:
770: . contains - true if the stratum contains the point
772: Level: intermediate
774: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
775: @*/
776: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
777: {
778: PetscInt v;
784: *contains = PETSC_FALSE;
785: DMLabelLookupStratum(label, value, &v);
786: if (v < 0) return(0);
788: if (label->validIS[v]) {
789: PetscInt i;
791: ISLocate(label->points[v], point, &i);
792: if (i >= 0) *contains = PETSC_TRUE;
793: } else {
794: PetscBool has;
796: PetscHSetIHas(label->ht[v], point, &has);
797: if (has) *contains = PETSC_TRUE;
798: }
799: return(0);
800: }
802: /*@
803: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
804: When a label is created, it is initialized to -1.
806: Not collective
808: Input parameter:
809: . label - a DMLabel object
811: Output parameter:
812: . defaultValue - the default value
814: Level: beginner
816: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
817: @*/
818: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
819: {
822: *defaultValue = label->defaultValue;
823: return(0);
824: }
826: /*@
827: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
828: When a label is created, it is initialized to -1.
830: Not collective
832: Input parameter:
833: . label - a DMLabel object
835: Output parameter:
836: . defaultValue - the default value
838: Level: beginner
840: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
841: @*/
842: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
843: {
846: label->defaultValue = defaultValue;
847: return(0);
848: }
850: /*@
851: 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())
853: Not collective
855: Input Parameters:
856: + label - the DMLabel
857: - point - the point
859: Output Parameter:
860: . value - The point value, or the default value (-1 by default)
862: Level: intermediate
864: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
865: @*/
866: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
867: {
868: PetscInt v;
874: *value = label->defaultValue;
875: for (v = 0; v < label->numStrata; ++v) {
876: if (label->validIS[v]) {
877: PetscInt i;
879: ISLocate(label->points[v], point, &i);
880: if (i >= 0) {
881: *value = label->stratumValues[v];
882: break;
883: }
884: } else {
885: PetscBool has;
887: PetscHSetIHas(label->ht[v], point, &has);
888: if (has) {
889: *value = label->stratumValues[v];
890: break;
891: }
892: }
893: }
894: return(0);
895: }
897: /*@
898: 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.
900: Not collective
902: Input Parameters:
903: + label - the DMLabel
904: . point - the point
905: - value - The point value
907: Level: intermediate
909: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
910: @*/
911: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
912: {
913: PetscInt v;
918: /* Find label value, add new entry if needed */
919: if (value == label->defaultValue) return(0);
920: DMLabelLookupAddStratum(label, value, &v);
921: /* Set key */
922: DMLabelMakeInvalid_Private(label, v);
923: PetscHSetIAdd(label->ht[v], point);
924: return(0);
925: }
927: /*@
928: DMLabelClearValue - Clear the value a label assigns to a point
930: Not collective
932: Input Parameters:
933: + label - the DMLabel
934: . point - the point
935: - value - The point value
937: Level: intermediate
939: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
940: @*/
941: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
942: {
943: PetscInt v;
948: /* Find label value */
949: DMLabelLookupStratum(label, value, &v);
950: if (v < 0) return(0);
952: if (label->bt) {
953: 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);
954: PetscBTClear(label->bt, point - label->pStart);
955: }
957: /* Delete key */
958: DMLabelMakeInvalid_Private(label, v);
959: PetscHSetIDel(label->ht[v], point);
960: return(0);
961: }
963: /*@
964: DMLabelInsertIS - Set all points in the IS to a value
966: Not collective
968: Input Parameters:
969: + label - the DMLabel
970: . is - the point IS
971: - value - The point value
973: Level: intermediate
975: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
976: @*/
977: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
978: {
979: PetscInt v, n, p;
980: const PetscInt *points;
981: PetscErrorCode ierr;
986: /* Find label value, add new entry if needed */
987: if (value == label->defaultValue) return(0);
988: DMLabelLookupAddStratum(label, value, &v);
989: /* Set keys */
990: DMLabelMakeInvalid_Private(label, v);
991: ISGetLocalSize(is, &n);
992: ISGetIndices(is, &points);
993: for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
994: ISRestoreIndices(is, &points);
995: return(0);
996: }
998: /*@
999: DMLabelGetNumValues - Get the number of values that the DMLabel takes
1001: Not collective
1003: Input Parameter:
1004: . label - the DMLabel
1006: Output Paramater:
1007: . numValues - the number of values
1009: Level: intermediate
1011: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1012: @*/
1013: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1014: {
1018: *numValues = label->numStrata;
1019: return(0);
1020: }
1022: /*@
1023: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1025: Not collective
1027: Input Parameter:
1028: . label - the DMLabel
1030: Output Paramater:
1031: . is - the value IS
1033: Level: intermediate
1035: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1036: @*/
1037: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1038: {
1044: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1045: return(0);
1046: }
1048: /*@
1049: DMLabelHasStratum - Determine whether points exist with the given value
1051: Not collective
1053: Input Parameters:
1054: + label - the DMLabel
1055: - value - the stratum value
1057: Output Paramater:
1058: . exists - Flag saying whether points exist
1060: Level: intermediate
1062: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1063: @*/
1064: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1065: {
1066: PetscInt v;
1072: DMLabelLookupStratum(label, value, &v);
1073: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1074: return(0);
1075: }
1077: /*@
1078: DMLabelGetStratumSize - Get the size of a stratum
1080: Not collective
1082: Input Parameters:
1083: + label - the DMLabel
1084: - value - the stratum value
1086: Output Paramater:
1087: . size - The number of points in the stratum
1089: Level: intermediate
1091: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1092: @*/
1093: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1094: {
1095: PetscInt v;
1101: *size = 0;
1102: DMLabelLookupStratum(label, value, &v);
1103: if (v < 0) return(0);
1104: DMLabelMakeValid_Private(label, v);
1105: *size = label->stratumSizes[v];
1106: return(0);
1107: }
1109: /*@
1110: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1112: Not collective
1114: Input Parameters:
1115: + label - the DMLabel
1116: - value - the stratum value
1118: Output Paramaters:
1119: + start - the smallest point in the stratum
1120: - end - the largest point in the stratum
1122: Level: intermediate
1124: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1125: @*/
1126: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1127: {
1128: PetscInt v, min, max;
1135: DMLabelLookupStratum(label, value, &v);
1136: if (v < 0) return(0);
1137: DMLabelMakeValid_Private(label, v);
1138: if (label->stratumSizes[v] <= 0) return(0);
1139: ISGetMinMax(label->points[v], &min, &max);
1140: if (start) *start = min;
1141: if (end) *end = max+1;
1142: return(0);
1143: }
1145: /*@
1146: DMLabelGetStratumIS - Get an IS with the stratum points
1148: Not collective
1150: Input Parameters:
1151: + label - the DMLabel
1152: - value - the stratum value
1154: Output Paramater:
1155: . points - The stratum points
1157: Level: intermediate
1159: Notes:
1160: The output IS should be destroyed when no longer needed.
1161: Returns NULL if the stratum is empty.
1163: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1164: @*/
1165: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1166: {
1167: PetscInt v;
1173: *points = NULL;
1174: DMLabelLookupStratum(label, value, &v);
1175: if (v < 0) return(0);
1176: DMLabelMakeValid_Private(label, v);
1177: PetscObjectReference((PetscObject) label->points[v]);
1178: *points = label->points[v];
1179: return(0);
1180: }
1182: /*@
1183: DMLabelSetStratumIS - Set the stratum points using an IS
1185: Not collective
1187: Input Parameters:
1188: + label - the DMLabel
1189: . value - the stratum value
1190: - points - The stratum points
1192: Level: intermediate
1194: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1195: @*/
1196: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1197: {
1198: PetscInt v;
1204: DMLabelLookupAddStratum(label, value, &v);
1205: if (is == label->points[v]) return(0);
1206: DMLabelClearStratum(label, value);
1207: ISGetLocalSize(is, &(label->stratumSizes[v]));
1208: PetscObjectReference((PetscObject)is);
1209: ISDestroy(&(label->points[v]));
1210: label->points[v] = is;
1211: label->validIS[v] = PETSC_TRUE;
1212: PetscObjectStateIncrease((PetscObject) label);
1213: if (label->bt) {
1214: const PetscInt *points;
1215: PetscInt p;
1217: ISGetIndices(is,&points);
1218: for (p = 0; p < label->stratumSizes[v]; ++p) {
1219: const PetscInt point = points[p];
1221: 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);
1222: PetscBTSet(label->bt, point - label->pStart);
1223: }
1224: }
1225: return(0);
1226: }
1228: /*@
1229: DMLabelClearStratum - Remove a stratum
1231: Not collective
1233: Input Parameters:
1234: + label - the DMLabel
1235: - value - the stratum value
1237: Level: intermediate
1239: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1240: @*/
1241: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1242: {
1243: PetscInt v;
1248: DMLabelLookupStratum(label, value, &v);
1249: if (v < 0) return(0);
1250: if (label->validIS[v]) {
1251: if (label->bt) {
1252: PetscInt i;
1253: const PetscInt *points;
1255: ISGetIndices(label->points[v], &points);
1256: for (i = 0; i < label->stratumSizes[v]; ++i) {
1257: const PetscInt point = points[i];
1259: 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);
1260: PetscBTClear(label->bt, point - label->pStart);
1261: }
1262: ISRestoreIndices(label->points[v], &points);
1263: }
1264: label->stratumSizes[v] = 0;
1265: ISDestroy(&label->points[v]);
1266: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1267: PetscObjectSetName((PetscObject) label->points[v], "indices");
1268: PetscObjectStateIncrease((PetscObject) label);
1269: } else {
1270: PetscHSetIClear(label->ht[v]);
1271: }
1272: return(0);
1273: }
1275: /*@
1276: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1278: Not collective
1280: Input Parameters:
1281: + label - The DMLabel
1282: . value - The label value for all points
1283: . pStart - The first point
1284: - pEnd - A point beyond all marked points
1286: Note: The marks points are [pStart, pEnd), and only the bounds are stored.
1288: Level: intermediate
1290: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1291: @*/
1292: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1293: {
1294: IS pIS;
1298: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1299: DMLabelSetStratumIS(label, value, pIS);
1300: ISDestroy(&pIS);
1301: return(0);
1302: }
1304: /*@
1305: DMLabelFilter - Remove all points outside of [start, end)
1307: Not collective
1309: Input Parameters:
1310: + label - the DMLabel
1311: . start - the first point kept
1312: - end - one more than the last point kept
1314: Level: intermediate
1316: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1317: @*/
1318: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1319: {
1320: PetscInt v;
1325: DMLabelDestroyIndex(label);
1326: DMLabelMakeAllValid_Private(label);
1327: for (v = 0; v < label->numStrata; ++v) {
1328: ISGeneralFilter(label->points[v], start, end);
1329: }
1330: DMLabelCreateIndex(label, start, end);
1331: return(0);
1332: }
1334: /*@
1335: DMLabelPermute - Create a new label with permuted points
1337: Not collective
1339: Input Parameters:
1340: + label - the DMLabel
1341: - permutation - the point permutation
1343: Output Parameter:
1344: . labelnew - the new label containing the permuted points
1346: Level: intermediate
1348: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1349: @*/
1350: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1351: {
1352: const PetscInt *perm;
1353: PetscInt numValues, numPoints, v, q;
1354: PetscErrorCode ierr;
1359: DMLabelMakeAllValid_Private(label);
1360: DMLabelDuplicate(label, labelNew);
1361: DMLabelGetNumValues(*labelNew, &numValues);
1362: ISGetLocalSize(permutation, &numPoints);
1363: ISGetIndices(permutation, &perm);
1364: for (v = 0; v < numValues; ++v) {
1365: const PetscInt size = (*labelNew)->stratumSizes[v];
1366: const PetscInt *points;
1367: PetscInt *pointsNew;
1369: ISGetIndices((*labelNew)->points[v],&points);
1370: PetscMalloc1(size,&pointsNew);
1371: for (q = 0; q < size; ++q) {
1372: const PetscInt point = points[q];
1374: 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);
1375: pointsNew[q] = perm[point];
1376: }
1377: ISRestoreIndices((*labelNew)->points[v],&points);
1378: PetscSortInt(size, pointsNew);
1379: ISDestroy(&((*labelNew)->points[v]));
1380: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1381: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1382: PetscFree(pointsNew);
1383: } else {
1384: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1385: }
1386: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1387: }
1388: ISRestoreIndices(permutation, &perm);
1389: if (label->bt) {
1390: PetscBTDestroy(&label->bt);
1391: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1392: }
1393: return(0);
1394: }
1396: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1397: {
1398: MPI_Comm comm;
1399: PetscInt s, l, nroots, nleaves, dof, offset, size;
1400: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1401: PetscSection rootSection;
1402: PetscSF labelSF;
1406: if (label) {DMLabelMakeAllValid_Private(label);}
1407: PetscObjectGetComm((PetscObject)sf, &comm);
1408: /* Build a section of stratum values per point, generate the according SF
1409: and distribute point-wise stratum values to leaves. */
1410: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1411: PetscSectionCreate(comm, &rootSection);
1412: PetscSectionSetChart(rootSection, 0, nroots);
1413: if (label) {
1414: for (s = 0; s < label->numStrata; ++s) {
1415: const PetscInt *points;
1417: ISGetIndices(label->points[s], &points);
1418: for (l = 0; l < label->stratumSizes[s]; l++) {
1419: PetscSectionGetDof(rootSection, points[l], &dof);
1420: PetscSectionSetDof(rootSection, points[l], dof+1);
1421: }
1422: ISRestoreIndices(label->points[s], &points);
1423: }
1424: }
1425: PetscSectionSetUp(rootSection);
1426: /* Create a point-wise array of stratum values */
1427: PetscSectionGetStorageSize(rootSection, &size);
1428: PetscMalloc1(size, &rootStrata);
1429: PetscCalloc1(nroots, &rootIdx);
1430: if (label) {
1431: for (s = 0; s < label->numStrata; ++s) {
1432: const PetscInt *points;
1434: ISGetIndices(label->points[s], &points);
1435: for (l = 0; l < label->stratumSizes[s]; l++) {
1436: const PetscInt p = points[l];
1437: PetscSectionGetOffset(rootSection, p, &offset);
1438: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1439: }
1440: ISRestoreIndices(label->points[s], &points);
1441: }
1442: }
1443: /* Build SF that maps label points to remote processes */
1444: PetscSectionCreate(comm, leafSection);
1445: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1446: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1447: PetscFree(remoteOffsets);
1448: /* Send the strata for each point over the derived SF */
1449: PetscSectionGetStorageSize(*leafSection, &size);
1450: PetscMalloc1(size, leafStrata);
1451: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1452: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1453: /* Clean up */
1454: PetscFree(rootStrata);
1455: PetscFree(rootIdx);
1456: PetscSectionDestroy(&rootSection);
1457: PetscSFDestroy(&labelSF);
1458: return(0);
1459: }
1461: /*@
1462: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1464: Collective on sf
1466: Input Parameters:
1467: + label - the DMLabel
1468: - sf - the map from old to new distribution
1470: Output Parameter:
1471: . labelnew - the new redistributed label
1473: Level: intermediate
1475: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1476: @*/
1477: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1478: {
1479: MPI_Comm comm;
1480: PetscSection leafSection;
1481: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1482: PetscInt *leafStrata, *strataIdx;
1483: PetscInt **points;
1484: const char *lname = NULL;
1485: char *name;
1486: PetscInt nameSize;
1487: PetscHSetI stratumHash;
1488: size_t len = 0;
1489: PetscMPIInt rank;
1494: if (label) {
1496: DMLabelMakeAllValid_Private(label);
1497: }
1498: PetscObjectGetComm((PetscObject)sf, &comm);
1499: MPI_Comm_rank(comm, &rank);
1500: /* Bcast name */
1501: if (!rank) {
1502: PetscObjectGetName((PetscObject) label, &lname);
1503: PetscStrlen(lname, &len);
1504: }
1505: nameSize = len;
1506: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1507: PetscMalloc1(nameSize+1, &name);
1508: if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1509: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1510: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1511: PetscFree(name);
1512: /* Bcast defaultValue */
1513: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1514: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1515: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1516: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1517: /* Determine received stratum values and initialise new label*/
1518: PetscHSetICreate(&stratumHash);
1519: PetscSectionGetStorageSize(leafSection, &size);
1520: for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1521: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1522: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1523: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1524: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1525: /* Turn leafStrata into indices rather than stratum values */
1526: offset = 0;
1527: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1528: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1529: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1530: PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1531: }
1532: for (p = 0; p < size; ++p) {
1533: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1534: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1535: }
1536: }
1537: /* Rebuild the point strata on the receiver */
1538: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1539: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1540: for (p=pStart; p<pEnd; p++) {
1541: PetscSectionGetDof(leafSection, p, &dof);
1542: PetscSectionGetOffset(leafSection, p, &offset);
1543: for (s=0; s<dof; s++) {
1544: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1545: }
1546: }
1547: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1548: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1549: PetscMalloc1((*labelNew)->numStrata,&points);
1550: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1551: PetscHSetICreate(&(*labelNew)->ht[s]);
1552: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1553: }
1554: /* Insert points into new strata */
1555: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1556: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1557: for (p=pStart; p<pEnd; p++) {
1558: PetscSectionGetDof(leafSection, p, &dof);
1559: PetscSectionGetOffset(leafSection, p, &offset);
1560: for (s=0; s<dof; s++) {
1561: stratum = leafStrata[offset+s];
1562: points[stratum][strataIdx[stratum]++] = p;
1563: }
1564: }
1565: for (s = 0; s < (*labelNew)->numStrata; s++) {
1566: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1567: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1568: }
1569: PetscFree(points);
1570: PetscHSetIDestroy(&stratumHash);
1571: PetscFree(leafStrata);
1572: PetscFree(strataIdx);
1573: PetscSectionDestroy(&leafSection);
1574: return(0);
1575: }
1577: /*@
1578: DMLabelGather - Gather all label values from leafs into roots
1580: Collective on sf
1582: Input Parameters:
1583: + label - the DMLabel
1584: - sf - the Star Forest point communication map
1586: Output Parameters:
1587: . labelNew - the new DMLabel with localised leaf values
1589: Level: developer
1591: Note: This is the inverse operation to DMLabelDistribute.
1593: .seealso: DMLabelDistribute()
1594: @*/
1595: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1596: {
1597: MPI_Comm comm;
1598: PetscSection rootSection;
1599: PetscSF sfLabel;
1600: PetscSFNode *rootPoints, *leafPoints;
1601: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1602: const PetscInt *rootDegree, *ilocal;
1603: PetscInt *rootStrata;
1604: const char *lname;
1605: char *name;
1606: PetscInt nameSize;
1607: size_t len = 0;
1608: PetscMPIInt rank, size;
1614: PetscObjectGetComm((PetscObject)sf, &comm);
1615: MPI_Comm_rank(comm, &rank);
1616: MPI_Comm_size(comm, &size);
1617: /* Bcast name */
1618: if (!rank) {
1619: PetscObjectGetName((PetscObject) label, &lname);
1620: PetscStrlen(lname, &len);
1621: }
1622: nameSize = len;
1623: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1624: PetscMalloc1(nameSize+1, &name);
1625: if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1626: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1627: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1628: PetscFree(name);
1629: /* Gather rank/index pairs of leaves into local roots to build
1630: an inverse, multi-rooted SF. Note that this ignores local leaf
1631: indexing due to the use of the multiSF in PetscSFGather. */
1632: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1633: PetscMalloc1(nroots, &leafPoints);
1634: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1635: for (p = 0; p < nleaves; p++) {
1636: PetscInt ilp = ilocal ? ilocal[p] : p;
1638: leafPoints[ilp].index = ilp;
1639: leafPoints[ilp].rank = rank;
1640: }
1641: PetscSFComputeDegreeBegin(sf, &rootDegree);
1642: PetscSFComputeDegreeEnd(sf, &rootDegree);
1643: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1644: PetscMalloc1(nmultiroots, &rootPoints);
1645: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1646: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1647: PetscSFCreate(comm,& sfLabel);
1648: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1649: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1650: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1651: /* Rebuild the point strata on the receiver */
1652: for (p = 0, idx = 0; p < nroots; p++) {
1653: for (d = 0; d < rootDegree[p]; d++) {
1654: PetscSectionGetDof(rootSection, idx+d, &dof);
1655: PetscSectionGetOffset(rootSection, idx+d, &offset);
1656: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1657: }
1658: idx += rootDegree[p];
1659: }
1660: PetscFree(leafPoints);
1661: PetscFree(rootStrata);
1662: PetscSectionDestroy(&rootSection);
1663: PetscSFDestroy(&sfLabel);
1664: return(0);
1665: }
1667: /*@
1668: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1670: Not collective
1672: Input Parameter:
1673: . label - the DMLabel
1675: Output Parameters:
1676: + section - the section giving offsets for each stratum
1677: - is - An IS containing all the label points
1679: Level: developer
1681: .seealso: DMLabelDistribute()
1682: @*/
1683: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1684: {
1685: IS vIS;
1686: const PetscInt *values;
1687: PetscInt *points;
1688: PetscInt nV, vS = 0, vE = 0, v, N;
1689: PetscErrorCode ierr;
1693: DMLabelGetNumValues(label, &nV);
1694: DMLabelGetValueIS(label, &vIS);
1695: ISGetIndices(vIS, &values);
1696: if (nV) {vS = values[0]; vE = values[0]+1;}
1697: for (v = 1; v < nV; ++v) {
1698: vS = PetscMin(vS, values[v]);
1699: vE = PetscMax(vE, values[v]+1);
1700: }
1701: PetscSectionCreate(PETSC_COMM_SELF, section);
1702: PetscSectionSetChart(*section, vS, vE);
1703: for (v = 0; v < nV; ++v) {
1704: PetscInt n;
1706: DMLabelGetStratumSize(label, values[v], &n);
1707: PetscSectionSetDof(*section, values[v], n);
1708: }
1709: PetscSectionSetUp(*section);
1710: PetscSectionGetStorageSize(*section, &N);
1711: PetscMalloc1(N, &points);
1712: for (v = 0; v < nV; ++v) {
1713: IS is;
1714: const PetscInt *spoints;
1715: PetscInt dof, off, p;
1717: PetscSectionGetDof(*section, values[v], &dof);
1718: PetscSectionGetOffset(*section, values[v], &off);
1719: DMLabelGetStratumIS(label, values[v], &is);
1720: ISGetIndices(is, &spoints);
1721: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1722: ISRestoreIndices(is, &spoints);
1723: ISDestroy(&is);
1724: }
1725: ISRestoreIndices(vIS, &values);
1726: ISDestroy(&vIS);
1727: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1728: return(0);
1729: }
1731: /*@
1732: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1733: the local section and an SF describing the section point overlap.
1735: Collective on sf
1737: Input Parameters:
1738: + s - The PetscSection for the local field layout
1739: . sf - The SF describing parallel layout of the section points
1740: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1741: . label - The label specifying the points
1742: - labelValue - The label stratum specifying the points
1744: Output Parameter:
1745: . gsection - The PetscSection for the global field layout
1747: Note: This gives negative sizes and offsets to points not owned by this process
1749: Level: developer
1751: .seealso: PetscSectionCreate()
1752: @*/
1753: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1754: {
1755: PetscInt *neg = NULL, *tmpOff = NULL;
1756: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1763: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1764: PetscSectionGetChart(s, &pStart, &pEnd);
1765: PetscSectionSetChart(*gsection, pStart, pEnd);
1766: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1767: if (nroots >= 0) {
1768: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1769: PetscCalloc1(nroots, &neg);
1770: if (nroots > pEnd-pStart) {
1771: PetscCalloc1(nroots, &tmpOff);
1772: } else {
1773: tmpOff = &(*gsection)->atlasDof[-pStart];
1774: }
1775: }
1776: /* Mark ghost points with negative dof */
1777: for (p = pStart; p < pEnd; ++p) {
1778: PetscInt value;
1780: DMLabelGetValue(label, p, &value);
1781: if (value != labelValue) continue;
1782: PetscSectionGetDof(s, p, &dof);
1783: PetscSectionSetDof(*gsection, p, dof);
1784: PetscSectionGetConstraintDof(s, p, &cdof);
1785: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1786: if (neg) neg[p] = -(dof+1);
1787: }
1788: PetscSectionSetUpBC(*gsection);
1789: if (nroots >= 0) {
1790: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1791: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1792: if (nroots > pEnd-pStart) {
1793: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1794: }
1795: }
1796: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1797: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1798: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1799: (*gsection)->atlasOff[p] = off;
1800: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1801: }
1802: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1803: globalOff -= off;
1804: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1805: (*gsection)->atlasOff[p] += globalOff;
1806: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1807: }
1808: /* Put in negative offsets for ghost points */
1809: if (nroots >= 0) {
1810: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1811: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1812: if (nroots > pEnd-pStart) {
1813: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1814: }
1815: }
1816: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1817: PetscFree(neg);
1818: return(0);
1819: }
1821: typedef struct _n_PetscSectionSym_Label
1822: {
1823: DMLabel label;
1824: PetscCopyMode *modes;
1825: PetscInt *sizes;
1826: const PetscInt ***perms;
1827: const PetscScalar ***rots;
1828: PetscInt (*minMaxOrients)[2];
1829: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1830: } PetscSectionSym_Label;
1832: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1833: {
1834: PetscInt i, j;
1835: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1836: PetscErrorCode ierr;
1839: for (i = 0; i <= sl->numStrata; i++) {
1840: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1841: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1842: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1843: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1844: }
1845: if (sl->perms[i]) {
1846: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1848: PetscFree(perms);
1849: }
1850: if (sl->rots[i]) {
1851: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1853: PetscFree(rots);
1854: }
1855: }
1856: }
1857: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1858: DMLabelDestroy(&sl->label);
1859: sl->numStrata = 0;
1860: return(0);
1861: }
1863: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1864: {
1868: PetscSectionSymLabelReset(sym);
1869: PetscFree(sym->data);
1870: return(0);
1871: }
1873: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1874: {
1875: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1876: PetscBool isAscii;
1877: DMLabel label = sl->label;
1878: const char *name;
1879: PetscErrorCode ierr;
1882: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1883: if (isAscii) {
1884: PetscInt i, j, k;
1885: PetscViewerFormat format;
1887: PetscViewerGetFormat(viewer,&format);
1888: if (label) {
1889: PetscViewerGetFormat(viewer,&format);
1890: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1891: PetscViewerASCIIPushTab(viewer);
1892: DMLabelView(label, viewer);
1893: PetscViewerASCIIPopTab(viewer);
1894: } else {
1895: PetscObjectGetName((PetscObject) sl->label, &name);
1896: PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);
1897: }
1898: } else {
1899: PetscViewerASCIIPrintf(viewer, "No label given\n");
1900: }
1901: PetscViewerASCIIPushTab(viewer);
1902: for (i = 0; i <= sl->numStrata; i++) {
1903: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1905: if (!(sl->perms[i] || sl->rots[i])) {
1906: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1907: } else {
1908: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1909: PetscViewerASCIIPushTab(viewer);
1910: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1911: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1912: PetscViewerASCIIPushTab(viewer);
1913: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1914: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1915: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1916: } else {
1917: PetscInt tab;
1919: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1920: PetscViewerASCIIPushTab(viewer);
1921: PetscViewerASCIIGetTab(viewer,&tab);
1922: if (sl->perms[i] && sl->perms[i][j]) {
1923: PetscViewerASCIIPrintf(viewer,"Permutation:");
1924: PetscViewerASCIISetTab(viewer,0);
1925: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1926: PetscViewerASCIIPrintf(viewer,"\n");
1927: PetscViewerASCIISetTab(viewer,tab);
1928: }
1929: if (sl->rots[i] && sl->rots[i][j]) {
1930: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1931: PetscViewerASCIISetTab(viewer,0);
1932: #if defined(PETSC_USE_COMPLEX)
1933: 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]));}
1934: #else
1935: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1936: #endif
1937: PetscViewerASCIIPrintf(viewer,"\n");
1938: PetscViewerASCIISetTab(viewer,tab);
1939: }
1940: PetscViewerASCIIPopTab(viewer);
1941: }
1942: }
1943: PetscViewerASCIIPopTab(viewer);
1944: }
1945: PetscViewerASCIIPopTab(viewer);
1946: }
1947: }
1948: PetscViewerASCIIPopTab(viewer);
1949: }
1950: return(0);
1951: }
1953: /*@
1954: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1956: Logically collective on sym
1958: Input parameters:
1959: + sym - the section symmetries
1960: - label - the DMLabel describing the types of points
1962: Level: developer:
1964: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1965: @*/
1966: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1967: {
1968: PetscSectionSym_Label *sl;
1969: PetscErrorCode ierr;
1973: sl = (PetscSectionSym_Label *) sym->data;
1974: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1975: if (label) {
1976: sl->label = label;
1977: PetscObjectReference((PetscObject) label);
1978: DMLabelGetNumValues(label,&sl->numStrata);
1979: 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);
1980: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1981: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1982: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1983: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1984: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1985: }
1986: return(0);
1987: }
1989: /*@C
1990: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1992: Logically collective on sym
1994: InputParameters:
1995: + sym - the section symmetries
1996: . stratum - the stratum value in the label that we are assigning symmetries for
1997: . size - the number of dofs for points in the stratum of the label
1998: . minOrient - the smallest orientation for a point in this stratum
1999: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2000: . mode - how sym should copy the perms and rots arrays
2001: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2002: - 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
2004: Level: developer
2006: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2007: @*/
2008: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2009: {
2010: PetscSectionSym_Label *sl;
2011: const char *name;
2012: PetscInt i, j, k;
2013: PetscErrorCode ierr;
2017: sl = (PetscSectionSym_Label *) sym->data;
2018: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
2019: for (i = 0; i <= sl->numStrata; i++) {
2020: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2022: if (stratum == value) break;
2023: }
2024: PetscObjectGetName((PetscObject) sl->label, &name);
2025: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2026: sl->sizes[i] = size;
2027: sl->modes[i] = mode;
2028: sl->minMaxOrients[i][0] = minOrient;
2029: sl->minMaxOrients[i][1] = maxOrient;
2030: if (mode == PETSC_COPY_VALUES) {
2031: if (perms) {
2032: PetscInt **ownPerms;
2034: PetscCalloc1(maxOrient - minOrient,&ownPerms);
2035: for (j = 0; j < maxOrient-minOrient; j++) {
2036: if (perms[j]) {
2037: PetscMalloc1(size,&ownPerms[j]);
2038: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2039: }
2040: }
2041: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2042: }
2043: if (rots) {
2044: PetscScalar **ownRots;
2046: PetscCalloc1(maxOrient - minOrient,&ownRots);
2047: for (j = 0; j < maxOrient-minOrient; j++) {
2048: if (rots[j]) {
2049: PetscMalloc1(size,&ownRots[j]);
2050: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2051: }
2052: }
2053: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2054: }
2055: } else {
2056: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2057: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2058: }
2059: return(0);
2060: }
2062: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2063: {
2064: PetscInt i, j, numStrata;
2065: PetscSectionSym_Label *sl;
2066: DMLabel label;
2067: PetscErrorCode ierr;
2070: sl = (PetscSectionSym_Label *) sym->data;
2071: numStrata = sl->numStrata;
2072: label = sl->label;
2073: for (i = 0; i < numPoints; i++) {
2074: PetscInt point = points[2*i];
2075: PetscInt ornt = points[2*i+1];
2077: for (j = 0; j < numStrata; j++) {
2078: if (label->validIS[j]) {
2079: PetscInt k;
2081: ISLocate(label->points[j],point,&k);
2082: if (k >= 0) break;
2083: } else {
2084: PetscBool has;
2086: PetscHSetIHas(label->ht[j], point, &has);
2087: if (has) break;
2088: }
2089: }
2090: 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);
2091: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2092: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2093: }
2094: return(0);
2095: }
2097: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2098: {
2099: PetscSectionSym_Label *sl;
2100: PetscErrorCode ierr;
2103: PetscNewLog(sym,&sl);
2104: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2105: sym->ops->view = PetscSectionSymView_Label;
2106: sym->ops->destroy = PetscSectionSymDestroy_Label;
2107: sym->data = (void *) sl;
2108: return(0);
2109: }
2111: /*@
2112: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2114: Collective
2116: Input Parameters:
2117: + comm - the MPI communicator for the new symmetry
2118: - label - the label defining the strata
2120: Output Parameters:
2121: . sym - the section symmetries
2123: Level: developer
2125: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2126: @*/
2127: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2128: {
2129: PetscErrorCode ierr;
2132: DMInitializePackage();
2133: PetscSectionSymCreate(comm,sym);
2134: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2135: PetscSectionSymLabelSetLabel(*sym,label);
2136: return(0);
2137: }