Actual source code: dmlabel.c
petsc-3.11.4 2019-09-28
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petscsf.h>
6: /*@C
7: DMLabelCreate - Create a DMLabel object, which is a multimap
9: Input parameters:
10: + comm - The communicator, usually PETSC_COMM_SELF
11: - name - The label name
13: Output parameter:
14: . label - The DMLabel
16: Level: beginner
18: .seealso: DMLabelDestroy()
19: @*/
20: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
21: {
26: DMInitializePackage();
28: PetscHeaderCreate(*label,DMLABEL_CLASSID,"DMLabel","DMLabel","DM",comm,DMLabelDestroy,DMLabelView);
30: (*label)->numStrata = 0;
31: (*label)->defaultValue = -1;
32: (*label)->stratumValues = NULL;
33: (*label)->validIS = NULL;
34: (*label)->stratumSizes = NULL;
35: (*label)->points = NULL;
36: (*label)->ht = NULL;
37: (*label)->pStart = -1;
38: (*label)->pEnd = -1;
39: (*label)->bt = NULL;
40: PetscHMapICreate(&(*label)->hmap);
41: PetscObjectSetName((PetscObject) *label, name);
42: return(0);
43: }
45: /*
46: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
48: Input parameter:
49: + label - The DMLabel
50: - v - The stratum value
52: Output parameter:
53: . label - The DMLabel with stratum in sorted list format
55: Level: developer
57: .seealso: DMLabelCreate()
58: */
59: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
60: {
61: PetscInt off = 0, *pointArray, p;
64: if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
66: 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);
67: PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
68: PetscMalloc1(label->stratumSizes[v], &pointArray);
69: PetscHSetIGetElems(label->ht[v], &off, pointArray);
70: PetscHSetIClear(label->ht[v]);
71: PetscSortInt(label->stratumSizes[v], pointArray);
72: if (label->bt) {
73: for (p = 0; p < label->stratumSizes[v]; ++p) {
74: const PetscInt point = pointArray[p];
75: 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);
76: PetscBTSet(label->bt, point - label->pStart);
77: }
78: }
79: ISCreateGeneral(PETSC_COMM_SELF,label->stratumSizes[v],pointArray,PETSC_OWN_POINTER,&(label->points[v]));
80: PetscObjectSetName((PetscObject) (label->points[v]), "indices");
81: label->validIS[v] = PETSC_TRUE;
82: PetscObjectStateIncrease((PetscObject) label);
83: return(0);
84: }
86: /*
87: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
89: Input parameter:
90: . label - The DMLabel
92: Output parameter:
93: . label - The DMLabel with all strata in sorted list format
95: Level: developer
97: .seealso: DMLabelCreate()
98: */
99: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
100: {
101: PetscInt v;
105: for (v = 0; v < label->numStrata; v++){
106: DMLabelMakeValid_Private(label, v);
107: }
108: return(0);
109: }
111: /*
112: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
114: Input parameter:
115: + label - The DMLabel
116: - v - The stratum value
118: Output parameter:
119: . label - The DMLabel with stratum in hash format
121: Level: developer
123: .seealso: DMLabelCreate()
124: */
125: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
126: {
127: PetscInt p;
128: const PetscInt *points;
131: if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
133: 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);
134: if (label->points[v]) {
135: ISGetIndices(label->points[v], &points);
136: for (p = 0; p < label->stratumSizes[v]; ++p) {
137: PetscHSetIAdd(label->ht[v], points[p]);
138: }
139: ISRestoreIndices(label->points[v],&points);
140: ISDestroy(&(label->points[v]));
141: }
142: label->validIS[v] = PETSC_FALSE;
143: return(0);
144: }
146: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
147: #define DMLABEL_LOOKUP_THRESHOLD 16
148: #endif
150: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
151: {
152: PetscInt v;
156: *index = -1;
157: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
158: for (v = 0; v < label->numStrata; ++v)
159: if (label->stratumValues[v] == value) {*index = v; break;}
160: } else {
161: PetscHMapIGet(label->hmap, value, index);
162: }
163: #if defined(PETSC_USE_DEBUG)
164: { /* Check strata hash map consistency */
165: PetscInt len, loc = -1;
166: PetscHMapIGetSize(label->hmap, &len);
167: if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
168: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
169: PetscHMapIGet(label->hmap, value, &loc);
170: } else {
171: for (v = 0; v < label->numStrata; ++v)
172: if (label->stratumValues[v] == value) {loc = v; break;}
173: }
174: if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
175: }
176: #endif
177: return(0);
178: }
180: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
181: {
182: PetscInt v;
183: PetscInt *tmpV;
184: PetscInt *tmpS;
185: PetscHSetI *tmpH, ht;
186: IS *tmpP, is;
187: PetscBool *tmpB;
188: PetscHMapI hmap = label->hmap;
192: v = label->numStrata;
193: tmpV = label->stratumValues;
194: tmpS = label->stratumSizes;
195: tmpH = label->ht;
196: tmpP = label->points;
197: tmpB = label->validIS;
198: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
199: PetscInt *oldV = tmpV;
200: PetscInt *oldS = tmpS;
201: PetscHSetI *oldH = tmpH;
202: IS *oldP = tmpP;
203: PetscBool *oldB = tmpB;
204: PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
205: PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
206: PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
207: PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
208: PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
209: PetscMemcpy(tmpV, oldV, v*sizeof(*tmpV));
210: PetscMemcpy(tmpS, oldS, v*sizeof(*tmpS));
211: PetscMemcpy(tmpH, oldH, v*sizeof(*tmpH));
212: PetscMemcpy(tmpP, oldP, v*sizeof(*tmpP));
213: PetscMemcpy(tmpB, oldB, v*sizeof(*tmpB));
214: PetscFree(oldV);
215: PetscFree(oldS);
216: PetscFree(oldH);
217: PetscFree(oldP);
218: PetscFree(oldB);
219: }
220: label->numStrata = v+1;
221: label->stratumValues = tmpV;
222: label->stratumSizes = tmpS;
223: label->ht = tmpH;
224: label->points = tmpP;
225: label->validIS = tmpB;
226: PetscHSetICreate(&ht);
227: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
228: PetscHMapISet(hmap, value, v);
229: tmpV[v] = value;
230: tmpS[v] = 0;
231: tmpH[v] = ht;
232: tmpP[v] = is;
233: tmpB[v] = PETSC_TRUE;
234: *index = v;
235: return(0);
236: }
238: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
239: {
242: DMLabelLookupStratum(label, value, index);
243: if (*index < 0) {DMLabelNewStratum(label, value, index);}
244: return(0);
245: }
247: /*@
248: DMLabelAddStratum - Adds a new stratum value in a DMLabel
250: Input Parameter:
251: + label - The DMLabel
252: - value - The stratum value
254: Level: beginner
256: .seealso: DMLabelCreate(), DMLabelDestroy()
257: @*/
258: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
259: {
260: PetscInt v;
265: DMLabelLookupAddStratum(label, value, &v);
266: return(0);
267: }
269: /*@
270: DMLabelAddStrata - Adds new stratum values in a DMLabel
272: Input Parameter:
273: + label - The DMLabel
274: . numStrata - The number of stratum values
275: - stratumValues - The stratum values
277: Level: beginner
279: .seealso: DMLabelCreate(), DMLabelDestroy()
280: @*/
281: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
282: {
283: PetscInt *values, v;
289: PetscMalloc1(numStrata, &values);
290: PetscMemcpy(values, stratumValues, numStrata*sizeof(PetscInt));
291: PetscSortRemoveDupsInt(&numStrata, values);
292: if (!label->numStrata) { /* Fast preallocation */
293: PetscInt *tmpV;
294: PetscInt *tmpS;
295: PetscHSetI *tmpH, ht;
296: IS *tmpP, is;
297: PetscBool *tmpB;
298: PetscHMapI hmap = label->hmap;
300: PetscMalloc1(numStrata, &tmpV);
301: PetscMalloc1(numStrata, &tmpS);
302: PetscMalloc1(numStrata, &tmpH);
303: PetscMalloc1(numStrata, &tmpP);
304: PetscMalloc1(numStrata, &tmpB);
305: label->numStrata = numStrata;
306: label->stratumValues = tmpV;
307: label->stratumSizes = tmpS;
308: label->ht = tmpH;
309: label->points = tmpP;
310: label->validIS = tmpB;
311: for (v = 0; v < numStrata; ++v) {
312: PetscHSetICreate(&ht);
313: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
314: PetscHMapISet(hmap, values[v], v);
315: tmpV[v] = values[v];
316: tmpS[v] = 0;
317: tmpH[v] = ht;
318: tmpP[v] = is;
319: tmpB[v] = PETSC_TRUE;
320: }
321: } else {
322: for (v = 0; v < numStrata; ++v) {
323: DMLabelAddStratum(label, values[v]);
324: }
325: }
326: PetscFree(values);
327: return(0);
328: }
330: /*@
331: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
333: Input Parameter:
334: + label - The DMLabel
335: - valueIS - Index set with stratum values
337: Level: beginner
339: .seealso: DMLabelCreate(), DMLabelDestroy()
340: @*/
341: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
342: {
343: PetscInt numStrata;
344: const PetscInt *stratumValues;
350: ISGetLocalSize(valueIS, &numStrata);
351: ISGetIndices(valueIS, &stratumValues);
352: DMLabelAddStrata(label, numStrata, stratumValues);
353: return(0);
354: }
356: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
357: {
358: PetscInt v;
359: PetscMPIInt rank;
363: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
364: PetscViewerASCIIPushSynchronized(viewer);
365: if (label) {
366: const char *name;
368: PetscObjectGetName((PetscObject) label, &name);
369: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
370: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
371: for (v = 0; v < label->numStrata; ++v) {
372: const PetscInt value = label->stratumValues[v];
373: const PetscInt *points;
374: PetscInt p;
376: ISGetIndices(label->points[v], &points);
377: for (p = 0; p < label->stratumSizes[v]; ++p) {
378: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
379: }
380: ISRestoreIndices(label->points[v],&points);
381: }
382: }
383: PetscViewerFlush(viewer);
384: PetscViewerASCIIPopSynchronized(viewer);
385: return(0);
386: }
388: /*@C
389: DMLabelView - View the label
391: Input Parameters:
392: + label - The DMLabel
393: - viewer - The PetscViewer
395: Level: intermediate
397: .seealso: DMLabelCreate(), DMLabelDestroy()
398: @*/
399: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
400: {
401: PetscBool iascii;
406: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
408: if (label) {DMLabelMakeAllValid_Private(label);}
409: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
410: if (iascii) {
411: DMLabelView_Ascii(label, viewer);
412: }
413: return(0);
414: }
416: /*@
417: DMLabelReset - Destroys internal data structures in a DMLabel
419: Input Parameter:
420: . label - The DMLabel
422: Level: beginner
424: .seealso: DMLabelDestroy(), DMLabelCreate()
425: @*/
426: PetscErrorCode DMLabelReset(DMLabel label)
427: {
428: PetscInt v;
433: for (v = 0; v < label->numStrata; ++v) {
434: PetscHSetIDestroy(&label->ht[v]);
435: ISDestroy(&label->points[v]);
436: }
437: label->numStrata = 0;
438: PetscFree(label->stratumValues);
439: PetscFree(label->stratumSizes);
440: PetscFree(label->ht);
441: PetscFree(label->points);
442: PetscFree(label->validIS);
443: PetscHMapIReset(label->hmap);
444: label->pStart = -1;
445: label->pEnd = -1;
446: PetscBTDestroy(&label->bt);
447: return(0);
448: }
450: /*@
451: DMLabelDestroy - Destroys a DMLabel
453: Input Parameter:
454: . label - The DMLabel
456: Level: beginner
458: .seealso: DMLabelReset(), DMLabelCreate()
459: @*/
460: PetscErrorCode DMLabelDestroy(DMLabel *label)
461: {
465: if (!*label) return(0);
467: if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
468: DMLabelReset(*label);
469: PetscHMapIDestroy(&(*label)->hmap);
470: PetscHeaderDestroy(label);
471: return(0);
472: }
474: /*@
475: DMLabelDuplicate - Duplicates a DMLabel
477: Input Parameter:
478: . label - The DMLabel
480: Output Parameter:
481: . labelnew - location to put new vector
483: Level: intermediate
485: .seealso: DMLabelCreate(), DMLabelDestroy()
486: @*/
487: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
488: {
489: const char *name;
490: PetscInt v;
495: DMLabelMakeAllValid_Private(label);
496: PetscObjectGetName((PetscObject) label, &name);
497: DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);
499: (*labelnew)->numStrata = label->numStrata;
500: (*labelnew)->defaultValue = label->defaultValue;
501: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
502: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
503: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
504: PetscMalloc1(label->numStrata, &(*labelnew)->points);
505: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
506: for (v = 0; v < label->numStrata; ++v) {
507: PetscHSetICreate(&(*labelnew)->ht[v]);
508: (*labelnew)->stratumValues[v] = label->stratumValues[v];
509: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
510: PetscObjectReference((PetscObject) (label->points[v]));
511: (*labelnew)->points[v] = label->points[v];
512: (*labelnew)->validIS[v] = PETSC_TRUE;
513: }
514: PetscHMapIDestroy(&(*labelnew)->hmap);
515: PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
516: (*labelnew)->pStart = -1;
517: (*labelnew)->pEnd = -1;
518: (*labelnew)->bt = NULL;
519: return(0);
520: }
522: /*@
523: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
525: Input Parameter:
526: . label - The DMLabel
528: Level: intermediate
530: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
531: @*/
532: PetscErrorCode DMLabelComputeIndex(DMLabel label)
533: {
534: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
539: DMLabelMakeAllValid_Private(label);
540: for (v = 0; v < label->numStrata; ++v) {
541: const PetscInt *points;
542: PetscInt i;
544: ISGetIndices(label->points[v], &points);
545: for (i = 0; i < label->stratumSizes[v]; ++i) {
546: const PetscInt point = points[i];
548: pStart = PetscMin(point, pStart);
549: pEnd = PetscMax(point+1, pEnd);
550: }
551: ISRestoreIndices(label->points[v], &points);
552: }
553: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
554: label->pEnd = pEnd;
555: DMLabelCreateIndex(label, label->pStart, label->pEnd);
556: return(0);
557: }
559: /*@
560: DMLabelCreateIndex - Create an index structure for membership determination
562: Input Parameters:
563: + label - The DMLabel
564: . pStart - The smallest point
565: - pEnd - The largest point + 1
567: Level: intermediate
569: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
570: @*/
571: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
572: {
573: PetscInt v;
578: DMLabelDestroyIndex(label);
579: DMLabelMakeAllValid_Private(label);
580: label->pStart = pStart;
581: label->pEnd = pEnd;
582: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
583: PetscBTCreate(pEnd - pStart, &label->bt);
584: for (v = 0; v < label->numStrata; ++v) {
585: const PetscInt *points;
586: PetscInt i;
588: ISGetIndices(label->points[v], &points);
589: for (i = 0; i < label->stratumSizes[v]; ++i) {
590: const PetscInt point = points[i];
592: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
593: PetscBTSet(label->bt, point - pStart);
594: }
595: ISRestoreIndices(label->points[v], &points);
596: }
597: return(0);
598: }
600: /*@
601: DMLabelDestroyIndex - Destroy the index structure
603: Input Parameter:
604: . label - the DMLabel
606: Level: intermediate
608: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
609: @*/
610: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
611: {
616: label->pStart = -1;
617: label->pEnd = -1;
618: PetscBTDestroy(&label->bt);
619: return(0);
620: }
622: /*@
623: DMLabelGetBounds - Return the smallest and largest point in the label
625: Input Parameter:
626: . label - the DMLabel
628: Output Parameters:
629: + pStart - The smallest point
630: - pEnd - The largest point + 1
632: Note: This will compute an index for the label if one does not exist.
634: Level: intermediate
636: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
637: @*/
638: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
639: {
644: if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
645: if (pStart) {
647: *pStart = label->pStart;
648: }
649: if (pEnd) {
651: *pEnd = label->pEnd;
652: }
653: return(0);
654: }
656: /*@
657: DMLabelHasValue - Determine whether a label assigns the value to any point
659: Input Parameters:
660: + label - the DMLabel
661: - value - the value
663: Output Parameter:
664: . contains - Flag indicating whether the label maps this value to any point
666: Level: developer
668: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
669: @*/
670: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
671: {
672: PetscInt v;
678: DMLabelLookupStratum(label, value, &v);
679: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
680: return(0);
681: }
683: /*@
684: DMLabelHasPoint - Determine whether a label assigns a value to a point
686: Input Parameters:
687: + label - the DMLabel
688: - point - the point
690: Output Parameter:
691: . contains - Flag indicating whether the label maps this point to a value
693: Note: The user must call DMLabelCreateIndex() before this function.
695: Level: developer
697: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
698: @*/
699: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
700: {
706: DMLabelMakeAllValid_Private(label);
707: #if defined(PETSC_USE_DEBUG)
708: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
709: 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);
710: #endif
711: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
712: return(0);
713: }
715: /*@
716: DMLabelStratumHasPoint - Return true if the stratum contains a point
718: Input Parameters:
719: + label - the DMLabel
720: . value - the stratum value
721: - point - the point
723: Output Parameter:
724: . contains - true if the stratum contains the point
726: Level: intermediate
728: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
729: @*/
730: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
731: {
732: PetscInt v;
738: *contains = PETSC_FALSE;
739: DMLabelLookupStratum(label, value, &v);
740: if (v < 0) return(0);
742: if (label->validIS[v]) {
743: PetscInt i;
745: ISLocate(label->points[v], point, &i);
746: if (i >= 0) *contains = PETSC_TRUE;
747: } else {
748: PetscBool has;
750: PetscHSetIHas(label->ht[v], point, &has);
751: if (has) *contains = PETSC_TRUE;
752: }
753: return(0);
754: }
756: /*@
757: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
758: When a label is created, it is initialized to -1.
760: Input parameter:
761: . label - a DMLabel object
763: Output parameter:
764: . defaultValue - the default value
766: Level: beginner
768: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
769: @*/
770: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
771: {
774: *defaultValue = label->defaultValue;
775: return(0);
776: }
778: /*@
779: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
780: When a label is created, it is initialized to -1.
782: Input parameter:
783: . label - a DMLabel object
785: Output parameter:
786: . defaultValue - the default value
788: Level: beginner
790: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
791: @*/
792: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
793: {
796: label->defaultValue = defaultValue;
797: return(0);
798: }
800: /*@
801: 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())
803: Input Parameters:
804: + label - the DMLabel
805: - point - the point
807: Output Parameter:
808: . value - The point value, or the default value (-1 by default)
810: Level: intermediate
812: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
813: @*/
814: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
815: {
816: PetscInt v;
822: *value = label->defaultValue;
823: for (v = 0; v < label->numStrata; ++v) {
824: if (label->validIS[v]) {
825: PetscInt i;
827: ISLocate(label->points[v], point, &i);
828: if (i >= 0) {
829: *value = label->stratumValues[v];
830: break;
831: }
832: } else {
833: PetscBool has;
835: PetscHSetIHas(label->ht[v], point, &has);
836: if (has) {
837: *value = label->stratumValues[v];
838: break;
839: }
840: }
841: }
842: return(0);
843: }
845: /*@
846: 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.
848: Input Parameters:
849: + label - the DMLabel
850: . point - the point
851: - value - The point value
853: Level: intermediate
855: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
856: @*/
857: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
858: {
859: PetscInt v;
864: /* Find label value, add new entry if needed */
865: if (value == label->defaultValue) return(0);
866: DMLabelLookupAddStratum(label, value, &v);
867: /* Set key */
868: DMLabelMakeInvalid_Private(label, v);
869: PetscHSetIAdd(label->ht[v], point);
870: return(0);
871: }
873: /*@
874: DMLabelClearValue - Clear the value a label assigns to a point
876: Input Parameters:
877: + label - the DMLabel
878: . point - the point
879: - value - The point value
881: Level: intermediate
883: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
884: @*/
885: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
886: {
887: PetscInt v;
892: /* Find label value */
893: DMLabelLookupStratum(label, value, &v);
894: if (v < 0) return(0);
896: if (label->bt) {
897: 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);
898: PetscBTClear(label->bt, point - label->pStart);
899: }
901: /* Delete key */
902: DMLabelMakeInvalid_Private(label, v);
903: PetscHSetIDel(label->ht[v], point);
904: return(0);
905: }
907: /*@
908: DMLabelInsertIS - Set all points in the IS to a value
910: Input Parameters:
911: + label - the DMLabel
912: . is - the point IS
913: - value - The point value
915: Level: intermediate
917: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
918: @*/
919: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
920: {
921: PetscInt v, n, p;
922: const PetscInt *points;
923: PetscErrorCode ierr;
928: /* Find label value, add new entry if needed */
929: if (value == label->defaultValue) return(0);
930: DMLabelLookupAddStratum(label, value, &v);
931: /* Set keys */
932: DMLabelMakeInvalid_Private(label, v);
933: ISGetLocalSize(is, &n);
934: ISGetIndices(is, &points);
935: for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
936: ISRestoreIndices(is, &points);
937: return(0);
938: }
940: /*@
941: DMLabelGetNumValues - Get the number of values that the DMLabel takes
943: Input Parameter:
944: . label - the DMLabel
946: Output Paramater:
947: . numValues - the number of values
949: Level: intermediate
951: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
952: @*/
953: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
954: {
958: *numValues = label->numStrata;
959: return(0);
960: }
962: /*@
963: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
965: Input Parameter:
966: . label - the DMLabel
968: Output Paramater:
969: . is - the value IS
971: Level: intermediate
973: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
974: @*/
975: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
976: {
982: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
983: return(0);
984: }
986: /*@
987: DMLabelHasStratum - Determine whether points exist with the given value
989: Input Parameters:
990: + label - the DMLabel
991: - value - the stratum value
993: Output Paramater:
994: . exists - Flag saying whether points exist
996: Level: intermediate
998: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
999: @*/
1000: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1001: {
1002: PetscInt v;
1008: DMLabelLookupStratum(label, value, &v);
1009: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1010: return(0);
1011: }
1013: /*@
1014: DMLabelGetStratumSize - Get the size of a stratum
1016: Input Parameters:
1017: + label - the DMLabel
1018: - value - the stratum value
1020: Output Paramater:
1021: . size - The number of points in the stratum
1023: Level: intermediate
1025: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1026: @*/
1027: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1028: {
1029: PetscInt v;
1035: *size = 0;
1036: DMLabelLookupStratum(label, value, &v);
1037: if (v < 0) return(0);
1038: DMLabelMakeValid_Private(label, v);
1039: *size = label->stratumSizes[v];
1040: return(0);
1041: }
1043: /*@
1044: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1046: Input Parameters:
1047: + label - the DMLabel
1048: - value - the stratum value
1050: Output Paramaters:
1051: + start - the smallest point in the stratum
1052: - end - the largest point in the stratum
1054: Level: intermediate
1056: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1057: @*/
1058: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1059: {
1060: PetscInt v, min, max;
1067: DMLabelLookupStratum(label, value, &v);
1068: if (v < 0) return(0);
1069: DMLabelMakeValid_Private(label, v);
1070: if (label->stratumSizes[v] <= 0) return(0);
1071: ISGetMinMax(label->points[v], &min, &max);
1072: if (start) *start = min;
1073: if (end) *end = max+1;
1074: return(0);
1075: }
1077: /*@
1078: DMLabelGetStratumIS - Get an IS with the stratum points
1080: Input Parameters:
1081: + label - the DMLabel
1082: - value - the stratum value
1084: Output Paramater:
1085: . points - The stratum points
1087: Level: intermediate
1089: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1090: @*/
1091: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1092: {
1093: PetscInt v;
1099: *points = NULL;
1100: DMLabelLookupStratum(label, value, &v);
1101: if (v < 0) return(0);
1102: DMLabelMakeValid_Private(label, v);
1103: PetscObjectReference((PetscObject) label->points[v]);
1104: *points = label->points[v];
1105: return(0);
1106: }
1108: /*@
1109: DMLabelSetStratumIS - Set the stratum points using an IS
1111: Input Parameters:
1112: + label - the DMLabel
1113: . value - the stratum value
1114: - points - The stratum points
1116: Level: intermediate
1118: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1119: @*/
1120: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1121: {
1122: PetscInt v;
1128: DMLabelLookupAddStratum(label, value, &v);
1129: if (is == label->points[v]) return(0);
1130: DMLabelClearStratum(label, value);
1131: ISGetLocalSize(is, &(label->stratumSizes[v]));
1132: PetscObjectReference((PetscObject)is);
1133: ISDestroy(&(label->points[v]));
1134: label->points[v] = is;
1135: label->validIS[v] = PETSC_TRUE;
1136: if (label->bt) {
1137: const PetscInt *points;
1138: PetscInt p;
1140: ISGetIndices(is,&points);
1141: for (p = 0; p < label->stratumSizes[v]; ++p) {
1142: const PetscInt point = points[p];
1144: 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);
1145: PetscBTSet(label->bt, point - label->pStart);
1146: }
1147: }
1148: return(0);
1149: }
1151: /*@
1152: DMLabelClearStratum - Remove a stratum
1154: Input Parameters:
1155: + label - the DMLabel
1156: - value - the stratum value
1158: Level: intermediate
1160: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1161: @*/
1162: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1163: {
1164: PetscInt v;
1169: DMLabelLookupStratum(label, value, &v);
1170: if (v < 0) return(0);
1171: if (label->validIS[v]) {
1172: if (label->bt) {
1173: PetscInt i;
1174: const PetscInt *points;
1176: ISGetIndices(label->points[v], &points);
1177: for (i = 0; i < label->stratumSizes[v]; ++i) {
1178: const PetscInt point = points[i];
1180: 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);
1181: PetscBTClear(label->bt, point - label->pStart);
1182: }
1183: ISRestoreIndices(label->points[v], &points);
1184: }
1185: label->stratumSizes[v] = 0;
1186: ISDestroy(&label->points[v]);
1187: ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_OWN_POINTER, &label->points[v]);
1188: PetscObjectSetName((PetscObject) label->points[v], "indices");
1189: } else {
1190: PetscHSetIClear(label->ht[v]);
1191: }
1192: return(0);
1193: }
1195: /*@
1196: DMLabelFilter - Remove all points outside of [start, end)
1198: Input Parameters:
1199: + label - the DMLabel
1200: . start - the first point
1201: - end - the last point
1203: Level: intermediate
1205: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1206: @*/
1207: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1208: {
1209: PetscInt v;
1214: DMLabelDestroyIndex(label);
1215: DMLabelMakeAllValid_Private(label);
1216: for (v = 0; v < label->numStrata; ++v) {
1217: PetscInt off, q;
1218: const PetscInt *points;
1219: PetscInt numPointsNew = 0, *pointsNew = NULL;
1221: ISGetIndices(label->points[v], &points);
1222: for (q = 0; q < label->stratumSizes[v]; ++q)
1223: if (points[q] >= start && points[q] < end)
1224: numPointsNew++;
1225: PetscMalloc1(numPointsNew, &pointsNew);
1226: for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
1227: if (points[q] >= start && points[q] < end)
1228: pointsNew[off++] = points[q];
1229: }
1230: ISRestoreIndices(label->points[v],&points);
1232: label->stratumSizes[v] = numPointsNew;
1233: ISDestroy(&label->points[v]);
1234: ISCreateGeneral(PETSC_COMM_SELF,numPointsNew, pointsNew, PETSC_OWN_POINTER, &label->points[v]);
1235: PetscObjectSetName((PetscObject) label->points[v], "indices");
1236: }
1237: DMLabelCreateIndex(label, start, end);
1238: return(0);
1239: }
1241: /*@
1242: DMLabelPermute - Create a new label with permuted points
1244: Input Parameters:
1245: + label - the DMLabel
1246: - permutation - the point permutation
1248: Output Parameter:
1249: . labelnew - the new label containing the permuted points
1251: Level: intermediate
1253: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1254: @*/
1255: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1256: {
1257: const PetscInt *perm;
1258: PetscInt numValues, numPoints, v, q;
1259: PetscErrorCode ierr;
1264: DMLabelMakeAllValid_Private(label);
1265: DMLabelDuplicate(label, labelNew);
1266: DMLabelGetNumValues(*labelNew, &numValues);
1267: ISGetLocalSize(permutation, &numPoints);
1268: ISGetIndices(permutation, &perm);
1269: for (v = 0; v < numValues; ++v) {
1270: const PetscInt size = (*labelNew)->stratumSizes[v];
1271: const PetscInt *points;
1272: PetscInt *pointsNew;
1274: ISGetIndices((*labelNew)->points[v],&points);
1275: PetscMalloc1(size,&pointsNew);
1276: for (q = 0; q < size; ++q) {
1277: const PetscInt point = points[q];
1279: 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);
1280: pointsNew[q] = perm[point];
1281: }
1282: ISRestoreIndices((*labelNew)->points[v],&points);
1283: PetscSortInt(size, pointsNew);
1284: ISDestroy(&((*labelNew)->points[v]));
1285: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1286: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1287: PetscFree(pointsNew);
1288: } else {
1289: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1290: }
1291: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1292: }
1293: ISRestoreIndices(permutation, &perm);
1294: if (label->bt) {
1295: PetscBTDestroy(&label->bt);
1296: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1297: }
1298: return(0);
1299: }
1301: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1302: {
1303: MPI_Comm comm;
1304: PetscInt s, l, nroots, nleaves, dof, offset, size;
1305: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1306: PetscSection rootSection;
1307: PetscSF labelSF;
1311: if (label) {DMLabelMakeAllValid_Private(label);}
1312: PetscObjectGetComm((PetscObject)sf, &comm);
1313: /* Build a section of stratum values per point, generate the according SF
1314: and distribute point-wise stratum values to leaves. */
1315: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1316: PetscSectionCreate(comm, &rootSection);
1317: PetscSectionSetChart(rootSection, 0, nroots);
1318: if (label) {
1319: for (s = 0; s < label->numStrata; ++s) {
1320: const PetscInt *points;
1322: ISGetIndices(label->points[s], &points);
1323: for (l = 0; l < label->stratumSizes[s]; l++) {
1324: PetscSectionGetDof(rootSection, points[l], &dof);
1325: PetscSectionSetDof(rootSection, points[l], dof+1);
1326: }
1327: ISRestoreIndices(label->points[s], &points);
1328: }
1329: }
1330: PetscSectionSetUp(rootSection);
1331: /* Create a point-wise array of stratum values */
1332: PetscSectionGetStorageSize(rootSection, &size);
1333: PetscMalloc1(size, &rootStrata);
1334: PetscCalloc1(nroots, &rootIdx);
1335: if (label) {
1336: for (s = 0; s < label->numStrata; ++s) {
1337: const PetscInt *points;
1339: ISGetIndices(label->points[s], &points);
1340: for (l = 0; l < label->stratumSizes[s]; l++) {
1341: const PetscInt p = points[l];
1342: PetscSectionGetOffset(rootSection, p, &offset);
1343: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1344: }
1345: ISRestoreIndices(label->points[s], &points);
1346: }
1347: }
1348: /* Build SF that maps label points to remote processes */
1349: PetscSectionCreate(comm, leafSection);
1350: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1351: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1352: PetscFree(remoteOffsets);
1353: /* Send the strata for each point over the derived SF */
1354: PetscSectionGetStorageSize(*leafSection, &size);
1355: PetscMalloc1(size, leafStrata);
1356: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1357: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1358: /* Clean up */
1359: PetscFree(rootStrata);
1360: PetscFree(rootIdx);
1361: PetscSectionDestroy(&rootSection);
1362: PetscSFDestroy(&labelSF);
1363: return(0);
1364: }
1366: /*@
1367: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1369: Input Parameters:
1370: + label - the DMLabel
1371: - sf - the map from old to new distribution
1373: Output Parameter:
1374: . labelnew - the new redistributed label
1376: Level: intermediate
1378: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1379: @*/
1380: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1381: {
1382: MPI_Comm comm;
1383: PetscSection leafSection;
1384: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1385: PetscInt *leafStrata, *strataIdx;
1386: PetscInt **points;
1387: const char *lname = NULL;
1388: char *name;
1389: PetscInt nameSize;
1390: PetscHSetI stratumHash;
1391: size_t len = 0;
1392: PetscMPIInt rank;
1397: if (label) {
1399: DMLabelMakeAllValid_Private(label);
1400: }
1401: PetscObjectGetComm((PetscObject)sf, &comm);
1402: MPI_Comm_rank(comm, &rank);
1403: /* Bcast name */
1404: if (!rank) {
1405: PetscObjectGetName((PetscObject) label, &lname);
1406: PetscStrlen(lname, &len);
1407: }
1408: nameSize = len;
1409: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1410: PetscMalloc1(nameSize+1, &name);
1411: if (!rank) {PetscMemcpy(name, lname, nameSize+1);}
1412: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1413: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1414: PetscFree(name);
1415: /* Bcast defaultValue */
1416: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1417: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1418: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1419: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1420: /* Determine received stratum values and initialise new label*/
1421: PetscHSetICreate(&stratumHash);
1422: PetscSectionGetStorageSize(leafSection, &size);
1423: for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1424: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1425: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1426: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1427: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1428: /* Turn leafStrata into indices rather than stratum values */
1429: offset = 0;
1430: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1431: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1432: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1433: PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1434: }
1435: for (p = 0; p < size; ++p) {
1436: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1437: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1438: }
1439: }
1440: /* Rebuild the point strata on the receiver */
1441: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1442: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1443: for (p=pStart; p<pEnd; p++) {
1444: PetscSectionGetDof(leafSection, p, &dof);
1445: PetscSectionGetOffset(leafSection, p, &offset);
1446: for (s=0; s<dof; s++) {
1447: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1448: }
1449: }
1450: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1451: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1452: PetscMalloc1((*labelNew)->numStrata,&points);
1453: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1454: PetscHSetICreate(&(*labelNew)->ht[s]);
1455: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1456: }
1457: /* Insert points into new strata */
1458: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1459: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1460: for (p=pStart; p<pEnd; p++) {
1461: PetscSectionGetDof(leafSection, p, &dof);
1462: PetscSectionGetOffset(leafSection, p, &offset);
1463: for (s=0; s<dof; s++) {
1464: stratum = leafStrata[offset+s];
1465: points[stratum][strataIdx[stratum]++] = p;
1466: }
1467: }
1468: for (s = 0; s < (*labelNew)->numStrata; s++) {
1469: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1470: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1471: }
1472: PetscFree(points);
1473: PetscHSetIDestroy(&stratumHash);
1474: PetscFree(leafStrata);
1475: PetscFree(strataIdx);
1476: PetscSectionDestroy(&leafSection);
1477: return(0);
1478: }
1480: /*@
1481: DMLabelGather - Gather all label values from leafs into roots
1483: Input Parameters:
1484: + label - the DMLabel
1485: - sf - the Star Forest point communication map
1487: Output Parameters:
1488: . labelNew - the new DMLabel with localised leaf values
1490: Level: developer
1492: Note: This is the inverse operation to DMLabelDistribute.
1494: .seealso: DMLabelDistribute()
1495: @*/
1496: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1497: {
1498: MPI_Comm comm;
1499: PetscSection rootSection;
1500: PetscSF sfLabel;
1501: PetscSFNode *rootPoints, *leafPoints;
1502: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1503: const PetscInt *rootDegree, *ilocal;
1504: PetscInt *rootStrata;
1505: const char *lname;
1506: char *name;
1507: PetscInt nameSize;
1508: size_t len = 0;
1509: PetscMPIInt rank, size;
1515: PetscObjectGetComm((PetscObject)sf, &comm);
1516: MPI_Comm_rank(comm, &rank);
1517: MPI_Comm_size(comm, &size);
1518: /* Bcast name */
1519: if (!rank) {
1520: PetscObjectGetName((PetscObject) label, &lname);
1521: PetscStrlen(lname, &len);
1522: }
1523: nameSize = len;
1524: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1525: PetscMalloc1(nameSize+1, &name);
1526: if (!rank) {PetscMemcpy(name, lname, nameSize+1);}
1527: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1528: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1529: PetscFree(name);
1530: /* Gather rank/index pairs of leaves into local roots to build
1531: an inverse, multi-rooted SF. Note that this ignores local leaf
1532: indexing due to the use of the multiSF in PetscSFGather. */
1533: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1534: PetscMalloc1(nroots, &leafPoints);
1535: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1536: for (p = 0; p < nleaves; p++) {
1537: PetscInt ilp = ilocal ? ilocal[p] : p;
1539: leafPoints[ilp].index = ilp;
1540: leafPoints[ilp].rank = rank;
1541: }
1542: PetscSFComputeDegreeBegin(sf, &rootDegree);
1543: PetscSFComputeDegreeEnd(sf, &rootDegree);
1544: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1545: PetscMalloc1(nmultiroots, &rootPoints);
1546: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1547: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1548: PetscSFCreate(comm,& sfLabel);
1549: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1550: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1551: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1552: /* Rebuild the point strata on the receiver */
1553: for (p = 0, idx = 0; p < nroots; p++) {
1554: for (d = 0; d < rootDegree[p]; d++) {
1555: PetscSectionGetDof(rootSection, idx+d, &dof);
1556: PetscSectionGetOffset(rootSection, idx+d, &offset);
1557: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1558: }
1559: idx += rootDegree[p];
1560: }
1561: PetscFree(leafPoints);
1562: PetscFree(rootStrata);
1563: PetscSectionDestroy(&rootSection);
1564: PetscSFDestroy(&sfLabel);
1565: return(0);
1566: }
1568: /*@
1569: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1571: Input Parameter:
1572: . label - the DMLabel
1574: Output Parameters:
1575: + section - the section giving offsets for each stratum
1576: - is - An IS containing all the label points
1578: Level: developer
1580: .seealso: DMLabelDistribute()
1581: @*/
1582: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1583: {
1584: IS vIS;
1585: const PetscInt *values;
1586: PetscInt *points;
1587: PetscInt nV, vS = 0, vE = 0, v, N;
1588: PetscErrorCode ierr;
1592: DMLabelGetNumValues(label, &nV);
1593: DMLabelGetValueIS(label, &vIS);
1594: ISGetIndices(vIS, &values);
1595: if (nV) {vS = values[0]; vE = values[0]+1;}
1596: for (v = 1; v < nV; ++v) {
1597: vS = PetscMin(vS, values[v]);
1598: vE = PetscMax(vE, values[v]+1);
1599: }
1600: PetscSectionCreate(PETSC_COMM_SELF, section);
1601: PetscSectionSetChart(*section, vS, vE);
1602: for (v = 0; v < nV; ++v) {
1603: PetscInt n;
1605: DMLabelGetStratumSize(label, values[v], &n);
1606: PetscSectionSetDof(*section, values[v], n);
1607: }
1608: PetscSectionSetUp(*section);
1609: PetscSectionGetStorageSize(*section, &N);
1610: PetscMalloc1(N, &points);
1611: for (v = 0; v < nV; ++v) {
1612: IS is;
1613: const PetscInt *spoints;
1614: PetscInt dof, off, p;
1616: PetscSectionGetDof(*section, values[v], &dof);
1617: PetscSectionGetOffset(*section, values[v], &off);
1618: DMLabelGetStratumIS(label, values[v], &is);
1619: ISGetIndices(is, &spoints);
1620: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1621: ISRestoreIndices(is, &spoints);
1622: ISDestroy(&is);
1623: }
1624: ISRestoreIndices(vIS, &values);
1625: ISDestroy(&vIS);
1626: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1627: return(0);
1628: }
1630: /*@
1631: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1632: the local section and an SF describing the section point overlap.
1634: Input Parameters:
1635: + s - The PetscSection for the local field layout
1636: . sf - The SF describing parallel layout of the section points
1637: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1638: . label - The label specifying the points
1639: - labelValue - The label stratum specifying the points
1641: Output Parameter:
1642: . gsection - The PetscSection for the global field layout
1644: Note: This gives negative sizes and offsets to points not owned by this process
1646: Level: developer
1648: .seealso: PetscSectionCreate()
1649: @*/
1650: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1651: {
1652: PetscInt *neg = NULL, *tmpOff = NULL;
1653: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1660: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1661: PetscSectionGetChart(s, &pStart, &pEnd);
1662: PetscSectionSetChart(*gsection, pStart, pEnd);
1663: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1664: if (nroots >= 0) {
1665: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1666: PetscCalloc1(nroots, &neg);
1667: if (nroots > pEnd-pStart) {
1668: PetscCalloc1(nroots, &tmpOff);
1669: } else {
1670: tmpOff = &(*gsection)->atlasDof[-pStart];
1671: }
1672: }
1673: /* Mark ghost points with negative dof */
1674: for (p = pStart; p < pEnd; ++p) {
1675: PetscInt value;
1677: DMLabelGetValue(label, p, &value);
1678: if (value != labelValue) continue;
1679: PetscSectionGetDof(s, p, &dof);
1680: PetscSectionSetDof(*gsection, p, dof);
1681: PetscSectionGetConstraintDof(s, p, &cdof);
1682: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1683: if (neg) neg[p] = -(dof+1);
1684: }
1685: PetscSectionSetUpBC(*gsection);
1686: if (nroots >= 0) {
1687: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1688: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1689: if (nroots > pEnd-pStart) {
1690: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1691: }
1692: }
1693: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1694: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1695: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1696: (*gsection)->atlasOff[p] = off;
1697: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1698: }
1699: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1700: globalOff -= off;
1701: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1702: (*gsection)->atlasOff[p] += globalOff;
1703: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1704: }
1705: /* Put in negative offsets for ghost points */
1706: if (nroots >= 0) {
1707: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1708: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1709: if (nroots > pEnd-pStart) {
1710: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1711: }
1712: }
1713: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1714: PetscFree(neg);
1715: return(0);
1716: }
1718: typedef struct _n_PetscSectionSym_Label
1719: {
1720: DMLabel label;
1721: PetscCopyMode *modes;
1722: PetscInt *sizes;
1723: const PetscInt ***perms;
1724: const PetscScalar ***rots;
1725: PetscInt (*minMaxOrients)[2];
1726: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1727: } PetscSectionSym_Label;
1729: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1730: {
1731: PetscInt i, j;
1732: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1733: PetscErrorCode ierr;
1736: for (i = 0; i <= sl->numStrata; i++) {
1737: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1738: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1739: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1740: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1741: }
1742: if (sl->perms[i]) {
1743: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1745: PetscFree(perms);
1746: }
1747: if (sl->rots[i]) {
1748: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1750: PetscFree(rots);
1751: }
1752: }
1753: }
1754: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1755: DMLabelDestroy(&sl->label);
1756: sl->numStrata = 0;
1757: return(0);
1758: }
1760: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1761: {
1765: PetscSectionSymLabelReset(sym);
1766: PetscFree(sym->data);
1767: return(0);
1768: }
1770: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1771: {
1772: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1773: PetscBool isAscii;
1774: DMLabel label = sl->label;
1775: const char *name;
1776: PetscErrorCode ierr;
1779: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1780: if (isAscii) {
1781: PetscInt i, j, k;
1782: PetscViewerFormat format;
1784: PetscViewerGetFormat(viewer,&format);
1785: if (label) {
1786: PetscViewerGetFormat(viewer,&format);
1787: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1788: PetscViewerASCIIPushTab(viewer);
1789: DMLabelView(label, viewer);
1790: PetscViewerASCIIPopTab(viewer);
1791: } else {
1792: PetscObjectGetName((PetscObject) sl->label, &name);
1793: PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);
1794: }
1795: } else {
1796: PetscViewerASCIIPrintf(viewer, "No label given\n");
1797: }
1798: PetscViewerASCIIPushTab(viewer);
1799: for (i = 0; i <= sl->numStrata; i++) {
1800: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1802: if (!(sl->perms[i] || sl->rots[i])) {
1803: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1804: } else {
1805: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1806: PetscViewerASCIIPushTab(viewer);
1807: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1808: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1809: PetscViewerASCIIPushTab(viewer);
1810: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1811: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1812: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1813: } else {
1814: PetscInt tab;
1816: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1817: PetscViewerASCIIPushTab(viewer);
1818: PetscViewerASCIIGetTab(viewer,&tab);
1819: if (sl->perms[i] && sl->perms[i][j]) {
1820: PetscViewerASCIIPrintf(viewer,"Permutation:");
1821: PetscViewerASCIISetTab(viewer,0);
1822: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1823: PetscViewerASCIIPrintf(viewer,"\n");
1824: PetscViewerASCIISetTab(viewer,tab);
1825: }
1826: if (sl->rots[i] && sl->rots[i][j]) {
1827: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1828: PetscViewerASCIISetTab(viewer,0);
1829: #if defined(PETSC_USE_COMPLEX)
1830: 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]));}
1831: #else
1832: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1833: #endif
1834: PetscViewerASCIIPrintf(viewer,"\n");
1835: PetscViewerASCIISetTab(viewer,tab);
1836: }
1837: PetscViewerASCIIPopTab(viewer);
1838: }
1839: }
1840: PetscViewerASCIIPopTab(viewer);
1841: }
1842: PetscViewerASCIIPopTab(viewer);
1843: }
1844: }
1845: PetscViewerASCIIPopTab(viewer);
1846: }
1847: return(0);
1848: }
1850: /*@
1851: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1853: Logically collective on sym
1855: Input parameters:
1856: + sym - the section symmetries
1857: - label - the DMLabel describing the types of points
1859: Level: developer:
1861: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1862: @*/
1863: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1864: {
1865: PetscSectionSym_Label *sl;
1866: PetscErrorCode ierr;
1870: sl = (PetscSectionSym_Label *) sym->data;
1871: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1872: if (label) {
1873: sl->label = label;
1874: PetscObjectReference((PetscObject) label);
1875: DMLabelGetNumValues(label,&sl->numStrata);
1876: 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);
1877: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1878: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1879: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1880: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1881: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1882: }
1883: return(0);
1884: }
1886: /*@C
1887: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1889: Logically collective on PetscSectionSym
1891: InputParameters:
1892: + sys - the section symmetries
1893: . stratum - the stratum value in the label that we are assigning symmetries for
1894: . size - the number of dofs for points in the stratum of the label
1895: . minOrient - the smallest orientation for a point in this stratum
1896: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1897: . mode - how sym should copy the perms and rots arrays
1898: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
1899: + 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
1901: Level: developer
1903: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1904: @*/
1905: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1906: {
1907: PetscSectionSym_Label *sl;
1908: const char *name;
1909: PetscInt i, j, k;
1910: PetscErrorCode ierr;
1914: sl = (PetscSectionSym_Label *) sym->data;
1915: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1916: for (i = 0; i <= sl->numStrata; i++) {
1917: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
1919: if (stratum == value) break;
1920: }
1921: PetscObjectGetName((PetscObject) sl->label, &name);
1922: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
1923: sl->sizes[i] = size;
1924: sl->modes[i] = mode;
1925: sl->minMaxOrients[i][0] = minOrient;
1926: sl->minMaxOrients[i][1] = maxOrient;
1927: if (mode == PETSC_COPY_VALUES) {
1928: if (perms) {
1929: PetscInt **ownPerms;
1931: PetscCalloc1(maxOrient - minOrient,&ownPerms);
1932: for (j = 0; j < maxOrient-minOrient; j++) {
1933: if (perms[j]) {
1934: PetscMalloc1(size,&ownPerms[j]);
1935: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
1936: }
1937: }
1938: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
1939: }
1940: if (rots) {
1941: PetscScalar **ownRots;
1943: PetscCalloc1(maxOrient - minOrient,&ownRots);
1944: for (j = 0; j < maxOrient-minOrient; j++) {
1945: if (rots[j]) {
1946: PetscMalloc1(size,&ownRots[j]);
1947: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
1948: }
1949: }
1950: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
1951: }
1952: } else {
1953: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
1954: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
1955: }
1956: return(0);
1957: }
1959: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
1960: {
1961: PetscInt i, j, numStrata;
1962: PetscSectionSym_Label *sl;
1963: DMLabel label;
1964: PetscErrorCode ierr;
1967: sl = (PetscSectionSym_Label *) sym->data;
1968: numStrata = sl->numStrata;
1969: label = sl->label;
1970: for (i = 0; i < numPoints; i++) {
1971: PetscInt point = points[2*i];
1972: PetscInt ornt = points[2*i+1];
1974: for (j = 0; j < numStrata; j++) {
1975: if (label->validIS[j]) {
1976: PetscInt k;
1978: ISLocate(label->points[j],point,&k);
1979: if (k >= 0) break;
1980: } else {
1981: PetscBool has;
1983: PetscHSetIHas(label->ht[j], point, &has);
1984: if (has) break;
1985: }
1986: }
1987: 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);
1988: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
1989: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
1990: }
1991: return(0);
1992: }
1994: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
1995: {
1996: PetscSectionSym_Label *sl;
1997: PetscErrorCode ierr;
2000: PetscNewLog(sym,&sl);
2001: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2002: sym->ops->view = PetscSectionSymView_Label;
2003: sym->ops->destroy = PetscSectionSymDestroy_Label;
2004: sym->data = (void *) sl;
2005: return(0);
2006: }
2008: /*@
2009: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2011: Collective
2013: Input Parameters:
2014: + comm - the MPI communicator for the new symmetry
2015: - label - the label defining the strata
2017: Output Parameters:
2018: . sym - the section symmetries
2020: Level: developer
2022: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2023: @*/
2024: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2025: {
2026: PetscErrorCode ierr;
2029: DMInitializePackage();
2030: PetscSectionSymCreate(comm,sym);
2031: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2032: PetscSectionSymLabelSetLabel(*sym,label);
2033: return(0);
2034: }