Actual source code: dmlabel.c
petsc-3.13.6 2020-09-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: 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 defined(PETSC_USE_DEBUG)
184: { /* Check strata hash map consistency */
185: PetscInt len, loc = -1;
186: PetscHMapIGetSize(label->hmap, &len);
187: if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
188: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
189: PetscHMapIGet(label->hmap, value, &loc);
190: } else {
191: for (v = 0; v < label->numStrata; ++v)
192: if (label->stratumValues[v] == value) {loc = v; break;}
193: }
194: if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
195: }
196: #endif
197: return(0);
198: }
200: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
201: {
202: PetscInt v;
203: PetscInt *tmpV;
204: PetscInt *tmpS;
205: PetscHSetI *tmpH, ht;
206: IS *tmpP, is;
207: PetscBool *tmpB;
208: PetscHMapI hmap = label->hmap;
212: v = label->numStrata;
213: tmpV = label->stratumValues;
214: tmpS = label->stratumSizes;
215: tmpH = label->ht;
216: tmpP = label->points;
217: tmpB = label->validIS;
218: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
219: PetscInt *oldV = tmpV;
220: PetscInt *oldS = tmpS;
221: PetscHSetI *oldH = tmpH;
222: IS *oldP = tmpP;
223: PetscBool *oldB = tmpB;
224: PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
225: PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
226: PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
227: PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
228: PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
229: PetscArraycpy(tmpV, oldV, v);
230: PetscArraycpy(tmpS, oldS, v);
231: PetscArraycpy(tmpH, oldH, v);
232: PetscArraycpy(tmpP, oldP, v);
233: PetscArraycpy(tmpB, oldB, v);
234: PetscFree(oldV);
235: PetscFree(oldS);
236: PetscFree(oldH);
237: PetscFree(oldP);
238: PetscFree(oldB);
239: }
240: label->numStrata = v+1;
241: label->stratumValues = tmpV;
242: label->stratumSizes = tmpS;
243: label->ht = tmpH;
244: label->points = tmpP;
245: label->validIS = tmpB;
246: PetscHSetICreate(&ht);
247: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
248: PetscHMapISet(hmap, value, v);
249: tmpV[v] = value;
250: tmpS[v] = 0;
251: tmpH[v] = ht;
252: tmpP[v] = is;
253: tmpB[v] = PETSC_TRUE;
254: PetscObjectStateIncrease((PetscObject) label);
255: *index = v;
256: return(0);
257: }
259: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
260: {
263: DMLabelLookupStratum(label, value, index);
264: if (*index < 0) {DMLabelNewStratum(label, value, index);}
265: return(0);
266: }
268: /*@
269: DMLabelAddStratum - Adds a new stratum value in a DMLabel
271: Input Parameter:
272: + label - The DMLabel
273: - value - The stratum value
275: Level: beginner
277: .seealso: DMLabelCreate(), DMLabelDestroy()
278: @*/
279: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
280: {
281: PetscInt v;
286: DMLabelLookupAddStratum(label, value, &v);
287: return(0);
288: }
290: /*@
291: DMLabelAddStrata - Adds new stratum values in a DMLabel
293: Not collective
295: Input Parameter:
296: + label - The DMLabel
297: . numStrata - The number of stratum values
298: - stratumValues - The stratum values
300: Level: beginner
302: .seealso: DMLabelCreate(), DMLabelDestroy()
303: @*/
304: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
305: {
306: PetscInt *values, v;
312: PetscMalloc1(numStrata, &values);
313: PetscArraycpy(values, stratumValues, numStrata);
314: PetscSortRemoveDupsInt(&numStrata, values);
315: if (!label->numStrata) { /* Fast preallocation */
316: PetscInt *tmpV;
317: PetscInt *tmpS;
318: PetscHSetI *tmpH, ht;
319: IS *tmpP, is;
320: PetscBool *tmpB;
321: PetscHMapI hmap = label->hmap;
323: PetscMalloc1(numStrata, &tmpV);
324: PetscMalloc1(numStrata, &tmpS);
325: PetscMalloc1(numStrata, &tmpH);
326: PetscMalloc1(numStrata, &tmpP);
327: PetscMalloc1(numStrata, &tmpB);
328: label->numStrata = numStrata;
329: label->stratumValues = tmpV;
330: label->stratumSizes = tmpS;
331: label->ht = tmpH;
332: label->points = tmpP;
333: label->validIS = tmpB;
334: for (v = 0; v < numStrata; ++v) {
335: PetscHSetICreate(&ht);
336: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
337: PetscHMapISet(hmap, values[v], v);
338: tmpV[v] = values[v];
339: tmpS[v] = 0;
340: tmpH[v] = ht;
341: tmpP[v] = is;
342: tmpB[v] = PETSC_TRUE;
343: }
344: PetscObjectStateIncrease((PetscObject) label);
345: } else {
346: for (v = 0; v < numStrata; ++v) {
347: DMLabelAddStratum(label, values[v]);
348: }
349: }
350: PetscFree(values);
351: return(0);
352: }
354: /*@
355: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
357: Not collective
359: Input Parameter:
360: + label - The DMLabel
361: - valueIS - Index set with stratum values
363: Level: beginner
365: .seealso: DMLabelCreate(), DMLabelDestroy()
366: @*/
367: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
368: {
369: PetscInt numStrata;
370: const PetscInt *stratumValues;
376: ISGetLocalSize(valueIS, &numStrata);
377: ISGetIndices(valueIS, &stratumValues);
378: DMLabelAddStrata(label, numStrata, stratumValues);
379: return(0);
380: }
382: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
383: {
384: PetscInt v;
385: PetscMPIInt rank;
389: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
390: PetscViewerASCIIPushSynchronized(viewer);
391: if (label) {
392: const char *name;
394: PetscObjectGetName((PetscObject) label, &name);
395: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
396: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
397: for (v = 0; v < label->numStrata; ++v) {
398: const PetscInt value = label->stratumValues[v];
399: const PetscInt *points;
400: PetscInt p;
402: ISGetIndices(label->points[v], &points);
403: for (p = 0; p < label->stratumSizes[v]; ++p) {
404: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
405: }
406: ISRestoreIndices(label->points[v],&points);
407: }
408: }
409: PetscViewerFlush(viewer);
410: PetscViewerASCIIPopSynchronized(viewer);
411: return(0);
412: }
414: /*@C
415: DMLabelView - View the label
417: Collective on viewer
419: Input Parameters:
420: + label - The DMLabel
421: - viewer - The PetscViewer
423: Level: intermediate
425: .seealso: DMLabelCreate(), DMLabelDestroy()
426: @*/
427: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
428: {
429: PetscBool iascii;
434: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
436: if (label) {DMLabelMakeAllValid_Private(label);}
437: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
438: if (iascii) {
439: DMLabelView_Ascii(label, viewer);
440: }
441: return(0);
442: }
444: /*@
445: DMLabelReset - Destroys internal data structures in a DMLabel
447: Not collective
449: Input Parameter:
450: . label - The DMLabel
452: Level: beginner
454: .seealso: DMLabelDestroy(), DMLabelCreate()
455: @*/
456: PetscErrorCode DMLabelReset(DMLabel label)
457: {
458: PetscInt v;
463: for (v = 0; v < label->numStrata; ++v) {
464: PetscHSetIDestroy(&label->ht[v]);
465: ISDestroy(&label->points[v]);
466: }
467: label->numStrata = 0;
468: PetscFree(label->stratumValues);
469: PetscFree(label->stratumSizes);
470: PetscFree(label->ht);
471: PetscFree(label->points);
472: PetscFree(label->validIS);
473: PetscHMapIReset(label->hmap);
474: label->pStart = -1;
475: label->pEnd = -1;
476: PetscBTDestroy(&label->bt);
477: return(0);
478: }
480: /*@
481: DMLabelDestroy - Destroys a DMLabel
483: Collective on label
485: Input Parameter:
486: . label - The DMLabel
488: Level: beginner
490: .seealso: DMLabelReset(), DMLabelCreate()
491: @*/
492: PetscErrorCode DMLabelDestroy(DMLabel *label)
493: {
497: if (!*label) return(0);
499: if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
500: DMLabelReset(*label);
501: PetscHMapIDestroy(&(*label)->hmap);
502: PetscHeaderDestroy(label);
503: return(0);
504: }
506: /*@
507: DMLabelDuplicate - Duplicates a DMLabel
509: Collective on label
511: Input Parameter:
512: . label - The DMLabel
514: Output Parameter:
515: . labelnew - location to put new vector
517: Level: intermediate
519: .seealso: DMLabelCreate(), DMLabelDestroy()
520: @*/
521: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
522: {
523: const char *name;
524: PetscInt v;
529: DMLabelMakeAllValid_Private(label);
530: PetscObjectGetName((PetscObject) label, &name);
531: DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);
533: (*labelnew)->numStrata = label->numStrata;
534: (*labelnew)->defaultValue = label->defaultValue;
535: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
536: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
537: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
538: PetscMalloc1(label->numStrata, &(*labelnew)->points);
539: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
540: for (v = 0; v < label->numStrata; ++v) {
541: PetscHSetICreate(&(*labelnew)->ht[v]);
542: (*labelnew)->stratumValues[v] = label->stratumValues[v];
543: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
544: PetscObjectReference((PetscObject) (label->points[v]));
545: (*labelnew)->points[v] = label->points[v];
546: (*labelnew)->validIS[v] = PETSC_TRUE;
547: }
548: PetscHMapIDestroy(&(*labelnew)->hmap);
549: PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
550: (*labelnew)->pStart = -1;
551: (*labelnew)->pEnd = -1;
552: (*labelnew)->bt = NULL;
553: return(0);
554: }
556: /*@
557: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
559: Not collective
561: Input Parameter:
562: . label - The DMLabel
564: Level: intermediate
566: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
567: @*/
568: PetscErrorCode DMLabelComputeIndex(DMLabel label)
569: {
570: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
575: DMLabelMakeAllValid_Private(label);
576: for (v = 0; v < label->numStrata; ++v) {
577: const PetscInt *points;
578: PetscInt i;
580: ISGetIndices(label->points[v], &points);
581: for (i = 0; i < label->stratumSizes[v]; ++i) {
582: const PetscInt point = points[i];
584: pStart = PetscMin(point, pStart);
585: pEnd = PetscMax(point+1, pEnd);
586: }
587: ISRestoreIndices(label->points[v], &points);
588: }
589: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
590: label->pEnd = pEnd;
591: DMLabelCreateIndex(label, label->pStart, label->pEnd);
592: return(0);
593: }
595: /*@
596: DMLabelCreateIndex - Create an index structure for membership determination
598: Not collective
600: Input Parameters:
601: + label - The DMLabel
602: . pStart - The smallest point
603: - pEnd - The largest point + 1
605: Level: intermediate
607: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
608: @*/
609: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
610: {
611: PetscInt v;
616: DMLabelDestroyIndex(label);
617: DMLabelMakeAllValid_Private(label);
618: label->pStart = pStart;
619: label->pEnd = pEnd;
620: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
621: PetscBTCreate(pEnd - pStart, &label->bt);
622: for (v = 0; v < label->numStrata; ++v) {
623: const PetscInt *points;
624: PetscInt i;
626: ISGetIndices(label->points[v], &points);
627: for (i = 0; i < label->stratumSizes[v]; ++i) {
628: const PetscInt point = points[i];
630: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
631: PetscBTSet(label->bt, point - pStart);
632: }
633: ISRestoreIndices(label->points[v], &points);
634: }
635: return(0);
636: }
638: /*@
639: DMLabelDestroyIndex - Destroy the index structure
641: Not collective
643: Input Parameter:
644: . label - the DMLabel
646: Level: intermediate
648: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
649: @*/
650: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
651: {
656: label->pStart = -1;
657: label->pEnd = -1;
658: PetscBTDestroy(&label->bt);
659: return(0);
660: }
662: /*@
663: DMLabelGetBounds - Return the smallest and largest point in the label
665: Not collective
667: Input Parameter:
668: . label - the DMLabel
670: Output Parameters:
671: + pStart - The smallest point
672: - pEnd - The largest point + 1
674: Note: This will compute an index for the label if one does not exist.
676: Level: intermediate
678: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
679: @*/
680: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
681: {
686: if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
687: if (pStart) {
689: *pStart = label->pStart;
690: }
691: if (pEnd) {
693: *pEnd = label->pEnd;
694: }
695: return(0);
696: }
698: /*@
699: DMLabelHasValue - Determine whether a label assigns the value to any point
701: Not collective
703: Input Parameters:
704: + label - the DMLabel
705: - value - the value
707: Output Parameter:
708: . contains - Flag indicating whether the label maps this value to any point
710: Level: developer
712: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
713: @*/
714: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
715: {
716: PetscInt v;
722: DMLabelLookupStratum(label, value, &v);
723: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
724: return(0);
725: }
727: /*@
728: DMLabelHasPoint - Determine whether a label assigns a value to a point
730: Not collective
732: Input Parameters:
733: + label - the DMLabel
734: - point - the point
736: Output Parameter:
737: . contains - Flag indicating whether the label maps this point to a value
739: Note: The user must call DMLabelCreateIndex() before this function.
741: Level: developer
743: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
744: @*/
745: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
746: {
752: DMLabelMakeAllValid_Private(label);
753: #if defined(PETSC_USE_DEBUG)
754: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
755: 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);
756: #endif
757: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
758: return(0);
759: }
761: /*@
762: DMLabelStratumHasPoint - Return true if the stratum contains a point
764: Not collective
766: Input Parameters:
767: + label - the DMLabel
768: . value - the stratum value
769: - point - the point
771: Output Parameter:
772: . contains - true if the stratum contains the point
774: Level: intermediate
776: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
777: @*/
778: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
779: {
780: PetscInt v;
786: *contains = PETSC_FALSE;
787: DMLabelLookupStratum(label, value, &v);
788: if (v < 0) return(0);
790: if (label->validIS[v]) {
791: PetscInt i;
793: ISLocate(label->points[v], point, &i);
794: if (i >= 0) *contains = PETSC_TRUE;
795: } else {
796: PetscBool has;
798: PetscHSetIHas(label->ht[v], point, &has);
799: if (has) *contains = PETSC_TRUE;
800: }
801: return(0);
802: }
804: /*@
805: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
806: When a label is created, it is initialized to -1.
808: Not collective
810: Input parameter:
811: . label - a DMLabel object
813: Output parameter:
814: . defaultValue - the default value
816: Level: beginner
818: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
819: @*/
820: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
821: {
824: *defaultValue = label->defaultValue;
825: return(0);
826: }
828: /*@
829: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
830: When a label is created, it is initialized to -1.
832: Not collective
834: Input parameter:
835: . label - a DMLabel object
837: Output parameter:
838: . defaultValue - the default value
840: Level: beginner
842: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
843: @*/
844: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
845: {
848: label->defaultValue = defaultValue;
849: return(0);
850: }
852: /*@
853: 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())
855: Not collective
857: Input Parameters:
858: + label - the DMLabel
859: - point - the point
861: Output Parameter:
862: . value - The point value, or the default value (-1 by default)
864: Level: intermediate
866: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
867: @*/
868: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
869: {
870: PetscInt v;
876: *value = label->defaultValue;
877: for (v = 0; v < label->numStrata; ++v) {
878: if (label->validIS[v]) {
879: PetscInt i;
881: ISLocate(label->points[v], point, &i);
882: if (i >= 0) {
883: *value = label->stratumValues[v];
884: break;
885: }
886: } else {
887: PetscBool has;
889: PetscHSetIHas(label->ht[v], point, &has);
890: if (has) {
891: *value = label->stratumValues[v];
892: break;
893: }
894: }
895: }
896: return(0);
897: }
899: /*@
900: 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.
902: Not collective
904: Input Parameters:
905: + label - the DMLabel
906: . point - the point
907: - value - The point value
909: Level: intermediate
911: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
912: @*/
913: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
914: {
915: PetscInt v;
920: /* Find label value, add new entry if needed */
921: if (value == label->defaultValue) return(0);
922: DMLabelLookupAddStratum(label, value, &v);
923: /* Set key */
924: DMLabelMakeInvalid_Private(label, v);
925: PetscHSetIAdd(label->ht[v], point);
926: return(0);
927: }
929: /*@
930: DMLabelClearValue - Clear the value a label assigns to a point
932: Not collective
934: Input Parameters:
935: + label - the DMLabel
936: . point - the point
937: - value - The point value
939: Level: intermediate
941: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
942: @*/
943: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
944: {
945: PetscInt v;
950: /* Find label value */
951: DMLabelLookupStratum(label, value, &v);
952: if (v < 0) return(0);
954: if (label->bt) {
955: 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);
956: PetscBTClear(label->bt, point - label->pStart);
957: }
959: /* Delete key */
960: DMLabelMakeInvalid_Private(label, v);
961: PetscHSetIDel(label->ht[v], point);
962: return(0);
963: }
965: /*@
966: DMLabelInsertIS - Set all points in the IS to a value
968: Not collective
970: Input Parameters:
971: + label - the DMLabel
972: . is - the point IS
973: - value - The point value
975: Level: intermediate
977: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
978: @*/
979: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
980: {
981: PetscInt v, n, p;
982: const PetscInt *points;
983: PetscErrorCode ierr;
988: /* Find label value, add new entry if needed */
989: if (value == label->defaultValue) return(0);
990: DMLabelLookupAddStratum(label, value, &v);
991: /* Set keys */
992: DMLabelMakeInvalid_Private(label, v);
993: ISGetLocalSize(is, &n);
994: ISGetIndices(is, &points);
995: for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
996: ISRestoreIndices(is, &points);
997: return(0);
998: }
1000: /*@
1001: DMLabelGetNumValues - Get the number of values that the DMLabel takes
1003: Not collective
1005: Input Parameter:
1006: . label - the DMLabel
1008: Output Paramater:
1009: . numValues - the number of values
1011: Level: intermediate
1013: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1014: @*/
1015: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1016: {
1020: *numValues = label->numStrata;
1021: return(0);
1022: }
1024: /*@
1025: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1027: Not collective
1029: Input Parameter:
1030: . label - the DMLabel
1032: Output Paramater:
1033: . is - the value IS
1035: Level: intermediate
1037: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1038: @*/
1039: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1040: {
1046: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1047: return(0);
1048: }
1050: /*@
1051: DMLabelHasStratum - Determine whether points exist with the given value
1053: Not collective
1055: Input Parameters:
1056: + label - the DMLabel
1057: - value - the stratum value
1059: Output Paramater:
1060: . exists - Flag saying whether points exist
1062: Level: intermediate
1064: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1065: @*/
1066: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1067: {
1068: PetscInt v;
1074: DMLabelLookupStratum(label, value, &v);
1075: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1076: return(0);
1077: }
1079: /*@
1080: DMLabelGetStratumSize - Get the size of a stratum
1082: Not collective
1084: Input Parameters:
1085: + label - the DMLabel
1086: - value - the stratum value
1088: Output Paramater:
1089: . size - The number of points in the stratum
1091: Level: intermediate
1093: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1094: @*/
1095: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1096: {
1097: PetscInt v;
1103: *size = 0;
1104: DMLabelLookupStratum(label, value, &v);
1105: if (v < 0) return(0);
1106: DMLabelMakeValid_Private(label, v);
1107: *size = label->stratumSizes[v];
1108: return(0);
1109: }
1111: /*@
1112: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1114: Not collective
1116: Input Parameters:
1117: + label - the DMLabel
1118: - value - the stratum value
1120: Output Paramaters:
1121: + start - the smallest point in the stratum
1122: - end - the largest point in the stratum
1124: Level: intermediate
1126: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1127: @*/
1128: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1129: {
1130: PetscInt v, min, max;
1137: DMLabelLookupStratum(label, value, &v);
1138: if (v < 0) return(0);
1139: DMLabelMakeValid_Private(label, v);
1140: if (label->stratumSizes[v] <= 0) return(0);
1141: ISGetMinMax(label->points[v], &min, &max);
1142: if (start) *start = min;
1143: if (end) *end = max+1;
1144: return(0);
1145: }
1147: /*@
1148: DMLabelGetStratumIS - Get an IS with the stratum points
1150: Not collective
1152: Input Parameters:
1153: + label - the DMLabel
1154: - value - the stratum value
1156: Output Paramater:
1157: . points - The stratum points
1159: Level: intermediate
1161: Notes:
1162: The output IS should be destroyed when no longer needed.
1163: Returns NULL if the stratum is empty.
1165: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1166: @*/
1167: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1168: {
1169: PetscInt v;
1175: *points = NULL;
1176: DMLabelLookupStratum(label, value, &v);
1177: if (v < 0) return(0);
1178: DMLabelMakeValid_Private(label, v);
1179: PetscObjectReference((PetscObject) label->points[v]);
1180: *points = label->points[v];
1181: return(0);
1182: }
1184: /*@
1185: DMLabelSetStratumIS - Set the stratum points using an IS
1187: Not collective
1189: Input Parameters:
1190: + label - the DMLabel
1191: . value - the stratum value
1192: - points - The stratum points
1194: Level: intermediate
1196: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1197: @*/
1198: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1199: {
1200: PetscInt v;
1206: DMLabelLookupAddStratum(label, value, &v);
1207: if (is == label->points[v]) return(0);
1208: DMLabelClearStratum(label, value);
1209: ISGetLocalSize(is, &(label->stratumSizes[v]));
1210: PetscObjectReference((PetscObject)is);
1211: ISDestroy(&(label->points[v]));
1212: label->points[v] = is;
1213: label->validIS[v] = PETSC_TRUE;
1214: PetscObjectStateIncrease((PetscObject) label);
1215: if (label->bt) {
1216: const PetscInt *points;
1217: PetscInt p;
1219: ISGetIndices(is,&points);
1220: for (p = 0; p < label->stratumSizes[v]; ++p) {
1221: const PetscInt point = points[p];
1223: 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);
1224: PetscBTSet(label->bt, point - label->pStart);
1225: }
1226: }
1227: return(0);
1228: }
1230: /*@
1231: DMLabelClearStratum - Remove a stratum
1233: Not collective
1235: Input Parameters:
1236: + label - the DMLabel
1237: - value - the stratum value
1239: Level: intermediate
1241: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1242: @*/
1243: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1244: {
1245: PetscInt v;
1250: DMLabelLookupStratum(label, value, &v);
1251: if (v < 0) return(0);
1252: if (label->validIS[v]) {
1253: if (label->bt) {
1254: PetscInt i;
1255: const PetscInt *points;
1257: ISGetIndices(label->points[v], &points);
1258: for (i = 0; i < label->stratumSizes[v]; ++i) {
1259: const PetscInt point = points[i];
1261: 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);
1262: PetscBTClear(label->bt, point - label->pStart);
1263: }
1264: ISRestoreIndices(label->points[v], &points);
1265: }
1266: label->stratumSizes[v] = 0;
1267: ISDestroy(&label->points[v]);
1268: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1269: PetscObjectSetName((PetscObject) label->points[v], "indices");
1270: PetscObjectStateIncrease((PetscObject) label);
1271: } else {
1272: PetscHSetIClear(label->ht[v]);
1273: }
1274: return(0);
1275: }
1277: /*@
1278: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1280: Not collective
1282: Input Parameters:
1283: + label - The DMLabel
1284: . value - The label value for all points
1285: . pStart - The first point
1286: - pEnd - A point beyond all marked points
1288: Note: The marks points are [pStart, pEnd), and only the bounds are stored.
1290: Level: intermediate
1292: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1293: @*/
1294: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1295: {
1296: IS pIS;
1300: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1301: DMLabelSetStratumIS(label, value, pIS);
1302: ISDestroy(&pIS);
1303: return(0);
1304: }
1306: /*@
1307: DMLabelFilter - Remove all points outside of [start, end)
1309: Not collective
1311: Input Parameters:
1312: + label - the DMLabel
1313: . start - the first point kept
1314: - end - one more than the last point kept
1316: Level: intermediate
1318: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1319: @*/
1320: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1321: {
1322: PetscInt v;
1327: DMLabelDestroyIndex(label);
1328: DMLabelMakeAllValid_Private(label);
1329: for (v = 0; v < label->numStrata; ++v) {
1330: ISGeneralFilter(label->points[v], start, end);
1331: }
1332: DMLabelCreateIndex(label, start, end);
1333: return(0);
1334: }
1336: /*@
1337: DMLabelPermute - Create a new label with permuted points
1339: Not collective
1341: Input Parameters:
1342: + label - the DMLabel
1343: - permutation - the point permutation
1345: Output Parameter:
1346: . labelnew - the new label containing the permuted points
1348: Level: intermediate
1350: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1351: @*/
1352: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1353: {
1354: const PetscInt *perm;
1355: PetscInt numValues, numPoints, v, q;
1356: PetscErrorCode ierr;
1361: DMLabelMakeAllValid_Private(label);
1362: DMLabelDuplicate(label, labelNew);
1363: DMLabelGetNumValues(*labelNew, &numValues);
1364: ISGetLocalSize(permutation, &numPoints);
1365: ISGetIndices(permutation, &perm);
1366: for (v = 0; v < numValues; ++v) {
1367: const PetscInt size = (*labelNew)->stratumSizes[v];
1368: const PetscInt *points;
1369: PetscInt *pointsNew;
1371: ISGetIndices((*labelNew)->points[v],&points);
1372: PetscMalloc1(size,&pointsNew);
1373: for (q = 0; q < size; ++q) {
1374: const PetscInt point = points[q];
1376: 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);
1377: pointsNew[q] = perm[point];
1378: }
1379: ISRestoreIndices((*labelNew)->points[v],&points);
1380: PetscSortInt(size, pointsNew);
1381: ISDestroy(&((*labelNew)->points[v]));
1382: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1383: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1384: PetscFree(pointsNew);
1385: } else {
1386: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1387: }
1388: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1389: }
1390: ISRestoreIndices(permutation, &perm);
1391: if (label->bt) {
1392: PetscBTDestroy(&label->bt);
1393: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1394: }
1395: return(0);
1396: }
1398: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1399: {
1400: MPI_Comm comm;
1401: PetscInt s, l, nroots, nleaves, dof, offset, size;
1402: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1403: PetscSection rootSection;
1404: PetscSF labelSF;
1408: if (label) {DMLabelMakeAllValid_Private(label);}
1409: PetscObjectGetComm((PetscObject)sf, &comm);
1410: /* Build a section of stratum values per point, generate the according SF
1411: and distribute point-wise stratum values to leaves. */
1412: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1413: PetscSectionCreate(comm, &rootSection);
1414: PetscSectionSetChart(rootSection, 0, nroots);
1415: if (label) {
1416: for (s = 0; s < label->numStrata; ++s) {
1417: const PetscInt *points;
1419: ISGetIndices(label->points[s], &points);
1420: for (l = 0; l < label->stratumSizes[s]; l++) {
1421: PetscSectionGetDof(rootSection, points[l], &dof);
1422: PetscSectionSetDof(rootSection, points[l], dof+1);
1423: }
1424: ISRestoreIndices(label->points[s], &points);
1425: }
1426: }
1427: PetscSectionSetUp(rootSection);
1428: /* Create a point-wise array of stratum values */
1429: PetscSectionGetStorageSize(rootSection, &size);
1430: PetscMalloc1(size, &rootStrata);
1431: PetscCalloc1(nroots, &rootIdx);
1432: if (label) {
1433: for (s = 0; s < label->numStrata; ++s) {
1434: const PetscInt *points;
1436: ISGetIndices(label->points[s], &points);
1437: for (l = 0; l < label->stratumSizes[s]; l++) {
1438: const PetscInt p = points[l];
1439: PetscSectionGetOffset(rootSection, p, &offset);
1440: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1441: }
1442: ISRestoreIndices(label->points[s], &points);
1443: }
1444: }
1445: /* Build SF that maps label points to remote processes */
1446: PetscSectionCreate(comm, leafSection);
1447: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1448: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1449: PetscFree(remoteOffsets);
1450: /* Send the strata for each point over the derived SF */
1451: PetscSectionGetStorageSize(*leafSection, &size);
1452: PetscMalloc1(size, leafStrata);
1453: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1454: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1455: /* Clean up */
1456: PetscFree(rootStrata);
1457: PetscFree(rootIdx);
1458: PetscSectionDestroy(&rootSection);
1459: PetscSFDestroy(&labelSF);
1460: return(0);
1461: }
1463: /*@
1464: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1466: Collective on sf
1468: Input Parameters:
1469: + label - the DMLabel
1470: - sf - the map from old to new distribution
1472: Output Parameter:
1473: . labelnew - the new redistributed label
1475: Level: intermediate
1477: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1478: @*/
1479: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1480: {
1481: MPI_Comm comm;
1482: PetscSection leafSection;
1483: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1484: PetscInt *leafStrata, *strataIdx;
1485: PetscInt **points;
1486: const char *lname = NULL;
1487: char *name;
1488: PetscInt nameSize;
1489: PetscHSetI stratumHash;
1490: size_t len = 0;
1491: PetscMPIInt rank;
1496: if (label) {
1498: DMLabelMakeAllValid_Private(label);
1499: }
1500: PetscObjectGetComm((PetscObject)sf, &comm);
1501: MPI_Comm_rank(comm, &rank);
1502: /* Bcast name */
1503: if (!rank) {
1504: PetscObjectGetName((PetscObject) label, &lname);
1505: PetscStrlen(lname, &len);
1506: }
1507: nameSize = len;
1508: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1509: PetscMalloc1(nameSize+1, &name);
1510: if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1511: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1512: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1513: PetscFree(name);
1514: /* Bcast defaultValue */
1515: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1516: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1517: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1518: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1519: /* Determine received stratum values and initialise new label*/
1520: PetscHSetICreate(&stratumHash);
1521: PetscSectionGetStorageSize(leafSection, &size);
1522: for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1523: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1524: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1525: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1526: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1527: /* Turn leafStrata into indices rather than stratum values */
1528: offset = 0;
1529: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1530: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1531: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1532: PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1533: }
1534: for (p = 0; p < size; ++p) {
1535: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1536: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1537: }
1538: }
1539: /* Rebuild the point strata on the receiver */
1540: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1541: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1542: for (p=pStart; p<pEnd; p++) {
1543: PetscSectionGetDof(leafSection, p, &dof);
1544: PetscSectionGetOffset(leafSection, p, &offset);
1545: for (s=0; s<dof; s++) {
1546: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1547: }
1548: }
1549: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1550: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1551: PetscMalloc1((*labelNew)->numStrata,&points);
1552: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1553: PetscHSetICreate(&(*labelNew)->ht[s]);
1554: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1555: }
1556: /* Insert points into new strata */
1557: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1558: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1559: for (p=pStart; p<pEnd; p++) {
1560: PetscSectionGetDof(leafSection, p, &dof);
1561: PetscSectionGetOffset(leafSection, p, &offset);
1562: for (s=0; s<dof; s++) {
1563: stratum = leafStrata[offset+s];
1564: points[stratum][strataIdx[stratum]++] = p;
1565: }
1566: }
1567: for (s = 0; s < (*labelNew)->numStrata; s++) {
1568: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1569: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1570: }
1571: PetscFree(points);
1572: PetscHSetIDestroy(&stratumHash);
1573: PetscFree(leafStrata);
1574: PetscFree(strataIdx);
1575: PetscSectionDestroy(&leafSection);
1576: return(0);
1577: }
1579: /*@
1580: DMLabelGather - Gather all label values from leafs into roots
1582: Collective on sf
1584: Input Parameters:
1585: + label - the DMLabel
1586: - sf - the Star Forest point communication map
1588: Output Parameters:
1589: . labelNew - the new DMLabel with localised leaf values
1591: Level: developer
1593: Note: This is the inverse operation to DMLabelDistribute.
1595: .seealso: DMLabelDistribute()
1596: @*/
1597: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1598: {
1599: MPI_Comm comm;
1600: PetscSection rootSection;
1601: PetscSF sfLabel;
1602: PetscSFNode *rootPoints, *leafPoints;
1603: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1604: const PetscInt *rootDegree, *ilocal;
1605: PetscInt *rootStrata;
1606: const char *lname;
1607: char *name;
1608: PetscInt nameSize;
1609: size_t len = 0;
1610: PetscMPIInt rank, size;
1616: PetscObjectGetComm((PetscObject)sf, &comm);
1617: MPI_Comm_rank(comm, &rank);
1618: MPI_Comm_size(comm, &size);
1619: /* Bcast name */
1620: if (!rank) {
1621: PetscObjectGetName((PetscObject) label, &lname);
1622: PetscStrlen(lname, &len);
1623: }
1624: nameSize = len;
1625: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1626: PetscMalloc1(nameSize+1, &name);
1627: if (!rank) {PetscArraycpy(name, lname, nameSize+1);}
1628: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1629: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1630: PetscFree(name);
1631: /* Gather rank/index pairs of leaves into local roots to build
1632: an inverse, multi-rooted SF. Note that this ignores local leaf
1633: indexing due to the use of the multiSF in PetscSFGather. */
1634: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1635: PetscMalloc1(nroots, &leafPoints);
1636: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1637: for (p = 0; p < nleaves; p++) {
1638: PetscInt ilp = ilocal ? ilocal[p] : p;
1640: leafPoints[ilp].index = ilp;
1641: leafPoints[ilp].rank = rank;
1642: }
1643: PetscSFComputeDegreeBegin(sf, &rootDegree);
1644: PetscSFComputeDegreeEnd(sf, &rootDegree);
1645: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1646: PetscMalloc1(nmultiroots, &rootPoints);
1647: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1648: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1649: PetscSFCreate(comm,& sfLabel);
1650: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1651: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1652: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1653: /* Rebuild the point strata on the receiver */
1654: for (p = 0, idx = 0; p < nroots; p++) {
1655: for (d = 0; d < rootDegree[p]; d++) {
1656: PetscSectionGetDof(rootSection, idx+d, &dof);
1657: PetscSectionGetOffset(rootSection, idx+d, &offset);
1658: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1659: }
1660: idx += rootDegree[p];
1661: }
1662: PetscFree(leafPoints);
1663: PetscFree(rootStrata);
1664: PetscSectionDestroy(&rootSection);
1665: PetscSFDestroy(&sfLabel);
1666: return(0);
1667: }
1669: /*@
1670: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1672: Not collective
1674: Input Parameter:
1675: . label - the DMLabel
1677: Output Parameters:
1678: + section - the section giving offsets for each stratum
1679: - is - An IS containing all the label points
1681: Level: developer
1683: .seealso: DMLabelDistribute()
1684: @*/
1685: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1686: {
1687: IS vIS;
1688: const PetscInt *values;
1689: PetscInt *points;
1690: PetscInt nV, vS = 0, vE = 0, v, N;
1691: PetscErrorCode ierr;
1695: DMLabelGetNumValues(label, &nV);
1696: DMLabelGetValueIS(label, &vIS);
1697: ISGetIndices(vIS, &values);
1698: if (nV) {vS = values[0]; vE = values[0]+1;}
1699: for (v = 1; v < nV; ++v) {
1700: vS = PetscMin(vS, values[v]);
1701: vE = PetscMax(vE, values[v]+1);
1702: }
1703: PetscSectionCreate(PETSC_COMM_SELF, section);
1704: PetscSectionSetChart(*section, vS, vE);
1705: for (v = 0; v < nV; ++v) {
1706: PetscInt n;
1708: DMLabelGetStratumSize(label, values[v], &n);
1709: PetscSectionSetDof(*section, values[v], n);
1710: }
1711: PetscSectionSetUp(*section);
1712: PetscSectionGetStorageSize(*section, &N);
1713: PetscMalloc1(N, &points);
1714: for (v = 0; v < nV; ++v) {
1715: IS is;
1716: const PetscInt *spoints;
1717: PetscInt dof, off, p;
1719: PetscSectionGetDof(*section, values[v], &dof);
1720: PetscSectionGetOffset(*section, values[v], &off);
1721: DMLabelGetStratumIS(label, values[v], &is);
1722: ISGetIndices(is, &spoints);
1723: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1724: ISRestoreIndices(is, &spoints);
1725: ISDestroy(&is);
1726: }
1727: ISRestoreIndices(vIS, &values);
1728: ISDestroy(&vIS);
1729: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1730: return(0);
1731: }
1733: /*@
1734: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1735: the local section and an SF describing the section point overlap.
1737: Collective on sf
1739: Input Parameters:
1740: + s - The PetscSection for the local field layout
1741: . sf - The SF describing parallel layout of the section points
1742: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1743: . label - The label specifying the points
1744: - labelValue - The label stratum specifying the points
1746: Output Parameter:
1747: . gsection - The PetscSection for the global field layout
1749: Note: This gives negative sizes and offsets to points not owned by this process
1751: Level: developer
1753: .seealso: PetscSectionCreate()
1754: @*/
1755: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1756: {
1757: PetscInt *neg = NULL, *tmpOff = NULL;
1758: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1765: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1766: PetscSectionGetChart(s, &pStart, &pEnd);
1767: PetscSectionSetChart(*gsection, pStart, pEnd);
1768: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1769: if (nroots >= 0) {
1770: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1771: PetscCalloc1(nroots, &neg);
1772: if (nroots > pEnd-pStart) {
1773: PetscCalloc1(nroots, &tmpOff);
1774: } else {
1775: tmpOff = &(*gsection)->atlasDof[-pStart];
1776: }
1777: }
1778: /* Mark ghost points with negative dof */
1779: for (p = pStart; p < pEnd; ++p) {
1780: PetscInt value;
1782: DMLabelGetValue(label, p, &value);
1783: if (value != labelValue) continue;
1784: PetscSectionGetDof(s, p, &dof);
1785: PetscSectionSetDof(*gsection, p, dof);
1786: PetscSectionGetConstraintDof(s, p, &cdof);
1787: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1788: if (neg) neg[p] = -(dof+1);
1789: }
1790: PetscSectionSetUpBC(*gsection);
1791: if (nroots >= 0) {
1792: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1793: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1794: if (nroots > pEnd-pStart) {
1795: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1796: }
1797: }
1798: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1799: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1800: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1801: (*gsection)->atlasOff[p] = off;
1802: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1803: }
1804: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1805: globalOff -= off;
1806: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1807: (*gsection)->atlasOff[p] += globalOff;
1808: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1809: }
1810: /* Put in negative offsets for ghost points */
1811: if (nroots >= 0) {
1812: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1813: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1814: if (nroots > pEnd-pStart) {
1815: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1816: }
1817: }
1818: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1819: PetscFree(neg);
1820: return(0);
1821: }
1823: typedef struct _n_PetscSectionSym_Label
1824: {
1825: DMLabel label;
1826: PetscCopyMode *modes;
1827: PetscInt *sizes;
1828: const PetscInt ***perms;
1829: const PetscScalar ***rots;
1830: PetscInt (*minMaxOrients)[2];
1831: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1832: } PetscSectionSym_Label;
1834: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1835: {
1836: PetscInt i, j;
1837: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1838: PetscErrorCode ierr;
1841: for (i = 0; i <= sl->numStrata; i++) {
1842: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1843: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1844: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1845: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1846: }
1847: if (sl->perms[i]) {
1848: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1850: PetscFree(perms);
1851: }
1852: if (sl->rots[i]) {
1853: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1855: PetscFree(rots);
1856: }
1857: }
1858: }
1859: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1860: DMLabelDestroy(&sl->label);
1861: sl->numStrata = 0;
1862: return(0);
1863: }
1865: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1866: {
1870: PetscSectionSymLabelReset(sym);
1871: PetscFree(sym->data);
1872: return(0);
1873: }
1875: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1876: {
1877: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1878: PetscBool isAscii;
1879: DMLabel label = sl->label;
1880: const char *name;
1881: PetscErrorCode ierr;
1884: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1885: if (isAscii) {
1886: PetscInt i, j, k;
1887: PetscViewerFormat format;
1889: PetscViewerGetFormat(viewer,&format);
1890: if (label) {
1891: PetscViewerGetFormat(viewer,&format);
1892: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1893: PetscViewerASCIIPushTab(viewer);
1894: DMLabelView(label, viewer);
1895: PetscViewerASCIIPopTab(viewer);
1896: } else {
1897: PetscObjectGetName((PetscObject) sl->label, &name);
1898: PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);
1899: }
1900: } else {
1901: PetscViewerASCIIPrintf(viewer, "No label given\n");
1902: }
1903: PetscViewerASCIIPushTab(viewer);
1904: for (i = 0; i <= sl->numStrata; i++) {
1905: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1907: if (!(sl->perms[i] || sl->rots[i])) {
1908: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1909: } else {
1910: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1911: PetscViewerASCIIPushTab(viewer);
1912: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1913: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1914: PetscViewerASCIIPushTab(viewer);
1915: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1916: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1917: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1918: } else {
1919: PetscInt tab;
1921: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1922: PetscViewerASCIIPushTab(viewer);
1923: PetscViewerASCIIGetTab(viewer,&tab);
1924: if (sl->perms[i] && sl->perms[i][j]) {
1925: PetscViewerASCIIPrintf(viewer,"Permutation:");
1926: PetscViewerASCIISetTab(viewer,0);
1927: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1928: PetscViewerASCIIPrintf(viewer,"\n");
1929: PetscViewerASCIISetTab(viewer,tab);
1930: }
1931: if (sl->rots[i] && sl->rots[i][j]) {
1932: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1933: PetscViewerASCIISetTab(viewer,0);
1934: #if defined(PETSC_USE_COMPLEX)
1935: 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]));}
1936: #else
1937: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1938: #endif
1939: PetscViewerASCIIPrintf(viewer,"\n");
1940: PetscViewerASCIISetTab(viewer,tab);
1941: }
1942: PetscViewerASCIIPopTab(viewer);
1943: }
1944: }
1945: PetscViewerASCIIPopTab(viewer);
1946: }
1947: PetscViewerASCIIPopTab(viewer);
1948: }
1949: }
1950: PetscViewerASCIIPopTab(viewer);
1951: }
1952: return(0);
1953: }
1955: /*@
1956: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1958: Logically collective on sym
1960: Input parameters:
1961: + sym - the section symmetries
1962: - label - the DMLabel describing the types of points
1964: Level: developer:
1966: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1967: @*/
1968: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1969: {
1970: PetscSectionSym_Label *sl;
1971: PetscErrorCode ierr;
1975: sl = (PetscSectionSym_Label *) sym->data;
1976: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1977: if (label) {
1978: sl->label = label;
1979: PetscObjectReference((PetscObject) label);
1980: DMLabelGetNumValues(label,&sl->numStrata);
1981: 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);
1982: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1983: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1984: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1985: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1986: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1987: }
1988: return(0);
1989: }
1991: /*@C
1992: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1994: Logically collective on sym
1996: InputParameters:
1997: + sym - the section symmetries
1998: . stratum - the stratum value in the label that we are assigning symmetries for
1999: . size - the number of dofs for points in the stratum of the label
2000: . minOrient - the smallest orientation for a point in this stratum
2001: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2002: . mode - how sym should copy the perms and rots arrays
2003: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2004: - 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
2006: Level: developer
2008: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2009: @*/
2010: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2011: {
2012: PetscSectionSym_Label *sl;
2013: const char *name;
2014: PetscInt i, j, k;
2015: PetscErrorCode ierr;
2019: sl = (PetscSectionSym_Label *) sym->data;
2020: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
2021: for (i = 0; i <= sl->numStrata; i++) {
2022: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2024: if (stratum == value) break;
2025: }
2026: PetscObjectGetName((PetscObject) sl->label, &name);
2027: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2028: sl->sizes[i] = size;
2029: sl->modes[i] = mode;
2030: sl->minMaxOrients[i][0] = minOrient;
2031: sl->minMaxOrients[i][1] = maxOrient;
2032: if (mode == PETSC_COPY_VALUES) {
2033: if (perms) {
2034: PetscInt **ownPerms;
2036: PetscCalloc1(maxOrient - minOrient,&ownPerms);
2037: for (j = 0; j < maxOrient-minOrient; j++) {
2038: if (perms[j]) {
2039: PetscMalloc1(size,&ownPerms[j]);
2040: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2041: }
2042: }
2043: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2044: }
2045: if (rots) {
2046: PetscScalar **ownRots;
2048: PetscCalloc1(maxOrient - minOrient,&ownRots);
2049: for (j = 0; j < maxOrient-minOrient; j++) {
2050: if (rots[j]) {
2051: PetscMalloc1(size,&ownRots[j]);
2052: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2053: }
2054: }
2055: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2056: }
2057: } else {
2058: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2059: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2060: }
2061: return(0);
2062: }
2064: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2065: {
2066: PetscInt i, j, numStrata;
2067: PetscSectionSym_Label *sl;
2068: DMLabel label;
2069: PetscErrorCode ierr;
2072: sl = (PetscSectionSym_Label *) sym->data;
2073: numStrata = sl->numStrata;
2074: label = sl->label;
2075: for (i = 0; i < numPoints; i++) {
2076: PetscInt point = points[2*i];
2077: PetscInt ornt = points[2*i+1];
2079: for (j = 0; j < numStrata; j++) {
2080: if (label->validIS[j]) {
2081: PetscInt k;
2083: ISLocate(label->points[j],point,&k);
2084: if (k >= 0) break;
2085: } else {
2086: PetscBool has;
2088: PetscHSetIHas(label->ht[j], point, &has);
2089: if (has) break;
2090: }
2091: }
2092: 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);
2093: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2094: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2095: }
2096: return(0);
2097: }
2099: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2100: {
2101: PetscSectionSym_Label *sl;
2102: PetscErrorCode ierr;
2105: PetscNewLog(sym,&sl);
2106: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2107: sym->ops->view = PetscSectionSymView_Label;
2108: sym->ops->destroy = PetscSectionSymDestroy_Label;
2109: sym->data = (void *) sl;
2110: return(0);
2111: }
2113: /*@
2114: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2116: Collective
2118: Input Parameters:
2119: + comm - the MPI communicator for the new symmetry
2120: - label - the label defining the strata
2122: Output Parameters:
2123: . sym - the section symmetries
2125: Level: developer
2127: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2128: @*/
2129: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2130: {
2131: PetscErrorCode ierr;
2134: DMInitializePackage();
2135: PetscSectionSymCreate(comm,sym);
2136: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2137: PetscSectionSymLabelSetLabel(*sym,label);
2138: return(0);
2139: }