Actual source code: dmlabel.c
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/sectionimpl.h>
4: #include <petscsf.h>
5: #include <petscsection.h>
7: /*@C
8: DMLabelCreate - Create a DMLabel object, which is a multimap
10: Collective
12: Input parameters:
13: + comm - The communicator, usually PETSC_COMM_SELF
14: - name - The label name
16: Output parameter:
17: . label - The DMLabel
19: Level: beginner
21: Notes:
22: The label name is actually usual PetscObject name.
23: One can get/set it with PetscObjectGetName()/PetscObjectSetName().
25: .seealso: DMLabelDestroy()
26: @*/
27: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
28: {
33: DMInitializePackage();
35: PetscHeaderCreate(*label,DMLABEL_CLASSID,"DMLabel","DMLabel","DM",comm,DMLabelDestroy,DMLabelView);
37: (*label)->numStrata = 0;
38: (*label)->defaultValue = -1;
39: (*label)->stratumValues = NULL;
40: (*label)->validIS = NULL;
41: (*label)->stratumSizes = NULL;
42: (*label)->points = NULL;
43: (*label)->ht = NULL;
44: (*label)->pStart = -1;
45: (*label)->pEnd = -1;
46: (*label)->bt = NULL;
47: PetscHMapICreate(&(*label)->hmap);
48: PetscObjectSetName((PetscObject) *label, name);
49: return(0);
50: }
52: /*
53: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
55: Not collective
57: Input parameter:
58: + label - The DMLabel
59: - v - The stratum value
61: Output parameter:
62: . label - The DMLabel with stratum in sorted list format
64: Level: developer
66: .seealso: DMLabelCreate()
67: */
68: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
69: {
70: IS is;
71: PetscInt off = 0, *pointArray, p;
74: if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
76: if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
77: PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
78: PetscMalloc1(label->stratumSizes[v], &pointArray);
79: PetscHSetIGetElems(label->ht[v], &off, pointArray);
80: PetscHSetIClear(label->ht[v]);
81: PetscSortInt(label->stratumSizes[v], pointArray);
82: if (label->bt) {
83: for (p = 0; p < label->stratumSizes[v]; ++p) {
84: const PetscInt point = pointArray[p];
85: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
86: PetscBTSet(label->bt, point - label->pStart);
87: }
88: }
89: if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
90: ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
91: PetscFree(pointArray);
92: } else {
93: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
94: }
95: PetscObjectSetName((PetscObject) is, "indices");
96: label->points[v] = is;
97: label->validIS[v] = PETSC_TRUE;
98: PetscObjectStateIncrease((PetscObject) label);
99: return(0);
100: }
102: /*
103: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
105: Not collective
107: Input parameter:
108: . label - The DMLabel
110: Output parameter:
111: . label - The DMLabel with all strata in sorted list format
113: Level: developer
115: .seealso: DMLabelCreate()
116: */
117: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
118: {
119: PetscInt v;
123: for (v = 0; v < label->numStrata; v++) {
124: DMLabelMakeValid_Private(label, v);
125: }
126: return(0);
127: }
129: /*
130: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
132: Not collective
134: Input parameter:
135: + label - The DMLabel
136: - v - The stratum value
138: Output parameter:
139: . label - The DMLabel with stratum in hash format
141: Level: developer
143: .seealso: DMLabelCreate()
144: */
145: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
146: {
147: PetscInt p;
148: const PetscInt *points;
151: if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
153: if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeInvalid_Private\n", v);
154: if (label->points[v]) {
155: ISGetIndices(label->points[v], &points);
156: for (p = 0; p < label->stratumSizes[v]; ++p) {
157: PetscHSetIAdd(label->ht[v], points[p]);
158: }
159: ISRestoreIndices(label->points[v],&points);
160: ISDestroy(&(label->points[v]));
161: }
162: label->validIS[v] = PETSC_FALSE;
163: return(0);
164: }
166: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
167: #define DMLABEL_LOOKUP_THRESHOLD 16
168: #endif
170: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
171: {
172: PetscInt v;
176: *index = -1;
177: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
178: for (v = 0; v < label->numStrata; ++v)
179: if (label->stratumValues[v] == value) {*index = v; break;}
180: } else {
181: PetscHMapIGet(label->hmap, value, index);
182: }
183: if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
184: PetscInt len, loc = -1;
185: PetscHMapIGetSize(label->hmap, &len);
186: if (len != label->numStrata) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map size");
187: if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
188: PetscHMapIGet(label->hmap, value, &loc);
189: } else {
190: for (v = 0; v < label->numStrata; ++v)
191: if (label->stratumValues[v] == value) {loc = v; break;}
192: }
193: if (loc != *index) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent strata hash map lookup");
194: }
195: return(0);
196: }
198: PETSC_STATIC_INLINE PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
199: {
200: PetscInt v;
201: PetscInt *tmpV;
202: PetscInt *tmpS;
203: PetscHSetI *tmpH, ht;
204: IS *tmpP, is;
205: PetscBool *tmpB;
206: PetscHMapI hmap = label->hmap;
210: v = label->numStrata;
211: tmpV = label->stratumValues;
212: tmpS = label->stratumSizes;
213: tmpH = label->ht;
214: tmpP = label->points;
215: tmpB = label->validIS;
216: { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free */
217: PetscInt *oldV = tmpV;
218: PetscInt *oldS = tmpS;
219: PetscHSetI *oldH = tmpH;
220: IS *oldP = tmpP;
221: PetscBool *oldB = tmpB;
222: PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
223: PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
224: PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
225: PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
226: PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
227: PetscArraycpy(tmpV, oldV, v);
228: PetscArraycpy(tmpS, oldS, v);
229: PetscArraycpy(tmpH, oldH, v);
230: PetscArraycpy(tmpP, oldP, v);
231: PetscArraycpy(tmpB, oldB, v);
232: PetscFree(oldV);
233: PetscFree(oldS);
234: PetscFree(oldH);
235: PetscFree(oldP);
236: PetscFree(oldB);
237: }
238: label->numStrata = v+1;
239: label->stratumValues = tmpV;
240: label->stratumSizes = tmpS;
241: label->ht = tmpH;
242: label->points = tmpP;
243: label->validIS = tmpB;
244: PetscHSetICreate(&ht);
245: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
246: PetscHMapISet(hmap, value, v);
247: tmpV[v] = value;
248: tmpS[v] = 0;
249: tmpH[v] = ht;
250: tmpP[v] = is;
251: tmpB[v] = PETSC_TRUE;
252: PetscObjectStateIncrease((PetscObject) label);
253: *index = v;
254: return(0);
255: }
257: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
258: {
261: DMLabelLookupStratum(label, value, index);
262: if (*index < 0) {DMLabelNewStratum(label, value, index);}
263: return(0);
264: }
266: /*@
267: DMLabelAddStratum - Adds a new stratum value in a DMLabel
269: Input Parameters:
270: + label - The DMLabel
271: - value - The stratum value
273: Level: beginner
275: .seealso: DMLabelCreate(), DMLabelDestroy()
276: @*/
277: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
278: {
279: PetscInt v;
284: DMLabelLookupAddStratum(label, value, &v);
285: return(0);
286: }
288: /*@
289: DMLabelAddStrata - Adds new stratum values in a DMLabel
291: Not collective
293: Input Parameters:
294: + label - The DMLabel
295: . numStrata - The number of stratum values
296: - stratumValues - The stratum values
298: Level: beginner
300: .seealso: DMLabelCreate(), DMLabelDestroy()
301: @*/
302: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
303: {
304: PetscInt *values, v;
310: PetscMalloc1(numStrata, &values);
311: PetscArraycpy(values, stratumValues, numStrata);
312: PetscSortRemoveDupsInt(&numStrata, values);
313: if (!label->numStrata) { /* Fast preallocation */
314: PetscInt *tmpV;
315: PetscInt *tmpS;
316: PetscHSetI *tmpH, ht;
317: IS *tmpP, is;
318: PetscBool *tmpB;
319: PetscHMapI hmap = label->hmap;
321: PetscMalloc1(numStrata, &tmpV);
322: PetscMalloc1(numStrata, &tmpS);
323: PetscMalloc1(numStrata, &tmpH);
324: PetscMalloc1(numStrata, &tmpP);
325: PetscMalloc1(numStrata, &tmpB);
326: label->numStrata = numStrata;
327: label->stratumValues = tmpV;
328: label->stratumSizes = tmpS;
329: label->ht = tmpH;
330: label->points = tmpP;
331: label->validIS = tmpB;
332: for (v = 0; v < numStrata; ++v) {
333: PetscHSetICreate(&ht);
334: ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
335: PetscHMapISet(hmap, values[v], v);
336: tmpV[v] = values[v];
337: tmpS[v] = 0;
338: tmpH[v] = ht;
339: tmpP[v] = is;
340: tmpB[v] = PETSC_TRUE;
341: }
342: PetscObjectStateIncrease((PetscObject) label);
343: } else {
344: for (v = 0; v < numStrata; ++v) {
345: DMLabelAddStratum(label, values[v]);
346: }
347: }
348: PetscFree(values);
349: return(0);
350: }
352: /*@
353: DMLabelAddStrataIS - Adds new stratum values in a DMLabel
355: Not collective
357: Input Parameters:
358: + label - The DMLabel
359: - valueIS - Index set with stratum values
361: Level: beginner
363: .seealso: DMLabelCreate(), DMLabelDestroy()
364: @*/
365: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
366: {
367: PetscInt numStrata;
368: const PetscInt *stratumValues;
374: ISGetLocalSize(valueIS, &numStrata);
375: ISGetIndices(valueIS, &stratumValues);
376: DMLabelAddStrata(label, numStrata, stratumValues);
377: return(0);
378: }
380: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
381: {
382: PetscInt v;
383: PetscMPIInt rank;
387: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
388: PetscViewerASCIIPushSynchronized(viewer);
389: if (label) {
390: const char *name;
392: PetscObjectGetName((PetscObject) label, &name);
393: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
394: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
395: for (v = 0; v < label->numStrata; ++v) {
396: const PetscInt value = label->stratumValues[v];
397: const PetscInt *points;
398: PetscInt p;
400: ISGetIndices(label->points[v], &points);
401: for (p = 0; p < label->stratumSizes[v]; ++p) {
402: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
403: }
404: ISRestoreIndices(label->points[v],&points);
405: }
406: }
407: PetscViewerFlush(viewer);
408: PetscViewerASCIIPopSynchronized(viewer);
409: return(0);
410: }
412: /*@C
413: DMLabelView - View the label
415: Collective on viewer
417: Input Parameters:
418: + label - The DMLabel
419: - viewer - The PetscViewer
421: Level: intermediate
423: .seealso: DMLabelCreate(), DMLabelDestroy()
424: @*/
425: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
426: {
427: PetscBool iascii;
432: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);}
434: if (label) {DMLabelMakeAllValid_Private(label);}
435: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
436: if (iascii) {
437: DMLabelView_Ascii(label, viewer);
438: }
439: return(0);
440: }
442: /*@
443: DMLabelReset - Destroys internal data structures in a DMLabel
445: Not collective
447: Input Parameter:
448: . label - The DMLabel
450: Level: beginner
452: .seealso: DMLabelDestroy(), DMLabelCreate()
453: @*/
454: PetscErrorCode DMLabelReset(DMLabel label)
455: {
456: PetscInt v;
461: for (v = 0; v < label->numStrata; ++v) {
462: PetscHSetIDestroy(&label->ht[v]);
463: ISDestroy(&label->points[v]);
464: }
465: label->numStrata = 0;
466: PetscFree(label->stratumValues);
467: PetscFree(label->stratumSizes);
468: PetscFree(label->ht);
469: PetscFree(label->points);
470: PetscFree(label->validIS);
471: label->stratumValues = NULL;
472: label->stratumSizes = NULL;
473: label->ht = NULL;
474: label->points = NULL;
475: label->validIS = NULL;
476: PetscHMapIReset(label->hmap);
477: label->pStart = -1;
478: label->pEnd = -1;
479: PetscBTDestroy(&label->bt);
480: return(0);
481: }
483: /*@
484: DMLabelDestroy - Destroys a DMLabel
486: Collective on label
488: Input Parameter:
489: . label - The DMLabel
491: Level: beginner
493: .seealso: DMLabelReset(), DMLabelCreate()
494: @*/
495: PetscErrorCode DMLabelDestroy(DMLabel *label)
496: {
500: if (!*label) return(0);
502: if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return(0);}
503: DMLabelReset(*label);
504: PetscHMapIDestroy(&(*label)->hmap);
505: PetscHeaderDestroy(label);
506: return(0);
507: }
509: /*@
510: DMLabelDuplicate - Duplicates a DMLabel
512: Collective on label
514: Input Parameter:
515: . label - The DMLabel
517: Output Parameter:
518: . labelnew - location to put new vector
520: Level: intermediate
522: .seealso: DMLabelCreate(), DMLabelDestroy()
523: @*/
524: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
525: {
526: const char *name;
527: PetscInt v;
532: DMLabelMakeAllValid_Private(label);
533: PetscObjectGetName((PetscObject) label, &name);
534: DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);
536: (*labelnew)->numStrata = label->numStrata;
537: (*labelnew)->defaultValue = label->defaultValue;
538: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
539: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
540: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
541: PetscMalloc1(label->numStrata, &(*labelnew)->points);
542: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
543: for (v = 0; v < label->numStrata; ++v) {
544: PetscHSetICreate(&(*labelnew)->ht[v]);
545: (*labelnew)->stratumValues[v] = label->stratumValues[v];
546: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
547: PetscObjectReference((PetscObject) (label->points[v]));
548: (*labelnew)->points[v] = label->points[v];
549: (*labelnew)->validIS[v] = PETSC_TRUE;
550: }
551: PetscHMapIDestroy(&(*labelnew)->hmap);
552: PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
553: (*labelnew)->pStart = -1;
554: (*labelnew)->pEnd = -1;
555: (*labelnew)->bt = NULL;
556: return(0);
557: }
559: /*@
560: DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds
562: Not collective
564: Input Parameter:
565: . label - The DMLabel
567: Level: intermediate
569: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
570: @*/
571: PetscErrorCode DMLabelComputeIndex(DMLabel label)
572: {
573: PetscInt pStart = PETSC_MAX_INT, pEnd = -1, v;
578: DMLabelMakeAllValid_Private(label);
579: for (v = 0; v < label->numStrata; ++v) {
580: const PetscInt *points;
581: PetscInt i;
583: ISGetIndices(label->points[v], &points);
584: for (i = 0; i < label->stratumSizes[v]; ++i) {
585: const PetscInt point = points[i];
587: pStart = PetscMin(point, pStart);
588: pEnd = PetscMax(point+1, pEnd);
589: }
590: ISRestoreIndices(label->points[v], &points);
591: }
592: label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
593: label->pEnd = pEnd;
594: DMLabelCreateIndex(label, label->pStart, label->pEnd);
595: return(0);
596: }
598: /*@
599: DMLabelCreateIndex - Create an index structure for membership determination
601: Not collective
603: Input Parameters:
604: + label - The DMLabel
605: . pStart - The smallest point
606: - pEnd - The largest point + 1
608: Level: intermediate
610: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
611: @*/
612: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
613: {
614: PetscInt v;
619: DMLabelDestroyIndex(label);
620: DMLabelMakeAllValid_Private(label);
621: label->pStart = pStart;
622: label->pEnd = pEnd;
623: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
624: PetscBTCreate(pEnd - pStart, &label->bt);
625: for (v = 0; v < label->numStrata; ++v) {
626: const PetscInt *points;
627: PetscInt i;
629: ISGetIndices(label->points[v], &points);
630: for (i = 0; i < label->stratumSizes[v]; ++i) {
631: const PetscInt point = points[i];
633: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
634: PetscBTSet(label->bt, point - pStart);
635: }
636: ISRestoreIndices(label->points[v], &points);
637: }
638: return(0);
639: }
641: /*@
642: DMLabelDestroyIndex - Destroy the index structure
644: Not collective
646: Input Parameter:
647: . label - the DMLabel
649: Level: intermediate
651: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
652: @*/
653: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
654: {
659: label->pStart = -1;
660: label->pEnd = -1;
661: PetscBTDestroy(&label->bt);
662: return(0);
663: }
665: /*@
666: DMLabelGetBounds - Return the smallest and largest point in the label
668: Not collective
670: Input Parameter:
671: . label - the DMLabel
673: Output Parameters:
674: + pStart - The smallest point
675: - pEnd - The largest point + 1
677: Note: This will compute an index for the label if one does not exist.
679: Level: intermediate
681: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
682: @*/
683: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
684: {
689: if ((label->pStart == -1) && (label->pEnd == -1)) {DMLabelComputeIndex(label);}
690: if (pStart) {
692: *pStart = label->pStart;
693: }
694: if (pEnd) {
696: *pEnd = label->pEnd;
697: }
698: return(0);
699: }
701: /*@
702: DMLabelHasValue - Determine whether a label assigns the value to any point
704: Not collective
706: Input Parameters:
707: + label - the DMLabel
708: - value - the value
710: Output Parameter:
711: . contains - Flag indicating whether the label maps this value to any point
713: Level: developer
715: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
716: @*/
717: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
718: {
719: PetscInt v;
725: DMLabelLookupStratum(label, value, &v);
726: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
727: return(0);
728: }
730: /*@
731: DMLabelHasPoint - Determine whether a label assigns a value to a point
733: Not collective
735: Input Parameters:
736: + label - the DMLabel
737: - point - the point
739: Output Parameter:
740: . contains - Flag indicating whether the label maps this point to a value
742: Note: The user must call DMLabelCreateIndex() before this function.
744: Level: developer
746: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
747: @*/
748: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
749: {
755: DMLabelMakeAllValid_Private(label);
756: if (PetscDefined(USE_DEBUG)) {
757: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
758: 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);
759: }
760: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
761: return(0);
762: }
764: /*@
765: DMLabelStratumHasPoint - Return true if the stratum contains a point
767: Not collective
769: Input Parameters:
770: + label - the DMLabel
771: . value - the stratum value
772: - point - the point
774: Output Parameter:
775: . contains - true if the stratum contains the point
777: Level: intermediate
779: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
780: @*/
781: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
782: {
783: PetscInt v;
789: *contains = PETSC_FALSE;
790: DMLabelLookupStratum(label, value, &v);
791: if (v < 0) return(0);
793: if (label->validIS[v]) {
794: PetscInt i;
796: ISLocate(label->points[v], point, &i);
797: if (i >= 0) *contains = PETSC_TRUE;
798: } else {
799: PetscBool has;
801: PetscHSetIHas(label->ht[v], point, &has);
802: if (has) *contains = PETSC_TRUE;
803: }
804: return(0);
805: }
807: /*@
808: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
809: When a label is created, it is initialized to -1.
811: Not collective
813: Input parameter:
814: . label - a DMLabel object
816: Output parameter:
817: . defaultValue - the default value
819: Level: beginner
821: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
822: @*/
823: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
824: {
827: *defaultValue = label->defaultValue;
828: return(0);
829: }
831: /*@
832: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
833: When a label is created, it is initialized to -1.
835: Not collective
837: Input parameter:
838: . label - a DMLabel object
840: Output parameter:
841: . defaultValue - the default value
843: Level: beginner
845: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
846: @*/
847: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
848: {
851: label->defaultValue = defaultValue;
852: return(0);
853: }
855: /*@
856: 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())
858: Not collective
860: Input Parameters:
861: + label - the DMLabel
862: - point - the point
864: Output Parameter:
865: . value - The point value, or the default value (-1 by default)
867: Level: intermediate
869: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
870: @*/
871: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
872: {
873: PetscInt v;
879: *value = label->defaultValue;
880: for (v = 0; v < label->numStrata; ++v) {
881: if (label->validIS[v]) {
882: PetscInt i;
884: ISLocate(label->points[v], point, &i);
885: if (i >= 0) {
886: *value = label->stratumValues[v];
887: break;
888: }
889: } else {
890: PetscBool has;
892: PetscHSetIHas(label->ht[v], point, &has);
893: if (has) {
894: *value = label->stratumValues[v];
895: break;
896: }
897: }
898: }
899: return(0);
900: }
902: /*@
903: 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.
905: Not collective
907: Input Parameters:
908: + label - the DMLabel
909: . point - the point
910: - value - The point value
912: Level: intermediate
914: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
915: @*/
916: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
917: {
918: PetscInt v;
923: /* Find label value, add new entry if needed */
924: if (value == label->defaultValue) return(0);
925: DMLabelLookupAddStratum(label, value, &v);
926: /* Set key */
927: DMLabelMakeInvalid_Private(label, v);
928: PetscHSetIAdd(label->ht[v], point);
929: return(0);
930: }
932: /*@
933: DMLabelClearValue - Clear the value a label assigns to a point
935: Not collective
937: Input Parameters:
938: + label - the DMLabel
939: . point - the point
940: - value - The point value
942: Level: intermediate
944: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
945: @*/
946: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
947: {
948: PetscInt v;
953: /* Find label value */
954: DMLabelLookupStratum(label, value, &v);
955: if (v < 0) return(0);
957: if (label->bt) {
958: 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);
959: PetscBTClear(label->bt, point - label->pStart);
960: }
962: /* Delete key */
963: DMLabelMakeInvalid_Private(label, v);
964: PetscHSetIDel(label->ht[v], point);
965: return(0);
966: }
968: /*@
969: DMLabelInsertIS - Set all points in the IS to a value
971: Not collective
973: Input Parameters:
974: + label - the DMLabel
975: . is - the point IS
976: - value - The point value
978: Level: intermediate
980: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
981: @*/
982: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
983: {
984: PetscInt v, n, p;
985: const PetscInt *points;
986: PetscErrorCode ierr;
991: /* Find label value, add new entry if needed */
992: if (value == label->defaultValue) return(0);
993: DMLabelLookupAddStratum(label, value, &v);
994: /* Set keys */
995: DMLabelMakeInvalid_Private(label, v);
996: ISGetLocalSize(is, &n);
997: ISGetIndices(is, &points);
998: for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
999: ISRestoreIndices(is, &points);
1000: return(0);
1001: }
1003: /*@
1004: DMLabelGetNumValues - Get the number of values that the DMLabel takes
1006: Not collective
1008: Input Parameter:
1009: . label - the DMLabel
1011: Output Parameter:
1012: . numValues - the number of values
1014: Level: intermediate
1016: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1017: @*/
1018: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1019: {
1023: *numValues = label->numStrata;
1024: return(0);
1025: }
1027: /*@
1028: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
1030: Not collective
1032: Input Parameter:
1033: . label - the DMLabel
1035: Output Parameter:
1036: . is - the value IS
1038: Level: intermediate
1040: Notes:
1041: The output IS should be destroyed when no longer needed.
1043: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1044: @*/
1045: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1046: {
1052: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1053: return(0);
1054: }
1056: /*@
1057: DMLabelGetValueIndex - Get the index of a given value in the list of values for the DMlabel, or -1 if it is not present
1059: Not collective
1061: Input Parameters:
1062: + label - the DMLabel
1063: - value - the value
1065: Output Parameter:
1066: . index - the index of value in the list of values
1068: Level: intermediate
1070: .seealso: DMLabelGetValueIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1071: @*/
1072: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1073: {
1074: PetscInt v;
1079: /* Do not assume they are sorted */
1080: for (v = 0; v < label->numStrata; ++v) if (label->stratumValues[v] == value) break;
1081: if (v >= label->numStrata) *index = -1;
1082: else *index = v;
1083: return(0);
1084: }
1086: /*@
1087: DMLabelHasStratum - Determine whether points exist with the given value
1089: Not collective
1091: Input Parameters:
1092: + label - the DMLabel
1093: - value - the stratum value
1095: Output Parameter:
1096: . exists - Flag saying whether points exist
1098: Level: intermediate
1100: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1101: @*/
1102: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1103: {
1104: PetscInt v;
1110: DMLabelLookupStratum(label, value, &v);
1111: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1112: return(0);
1113: }
1115: /*@
1116: DMLabelGetStratumSize - Get the size of a stratum
1118: Not collective
1120: Input Parameters:
1121: + label - the DMLabel
1122: - value - the stratum value
1124: Output Parameter:
1125: . size - The number of points in the stratum
1127: Level: intermediate
1129: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1130: @*/
1131: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1132: {
1133: PetscInt v;
1139: *size = 0;
1140: DMLabelLookupStratum(label, value, &v);
1141: if (v < 0) return(0);
1142: DMLabelMakeValid_Private(label, v);
1143: *size = label->stratumSizes[v];
1144: return(0);
1145: }
1147: /*@
1148: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
1150: Not collective
1152: Input Parameters:
1153: + label - the DMLabel
1154: - value - the stratum value
1156: Output Parameters:
1157: + start - the smallest point in the stratum
1158: - end - the largest point in the stratum
1160: Level: intermediate
1162: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1163: @*/
1164: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1165: {
1166: PetscInt v, min, max;
1173: DMLabelLookupStratum(label, value, &v);
1174: if (v < 0) return(0);
1175: DMLabelMakeValid_Private(label, v);
1176: if (label->stratumSizes[v] <= 0) return(0);
1177: ISGetMinMax(label->points[v], &min, &max);
1178: if (start) *start = min;
1179: if (end) *end = max+1;
1180: return(0);
1181: }
1183: /*@
1184: DMLabelGetStratumIS - Get an IS with the stratum points
1186: Not collective
1188: Input Parameters:
1189: + label - the DMLabel
1190: - value - the stratum value
1192: Output Parameter:
1193: . points - The stratum points
1195: Level: intermediate
1197: Notes:
1198: The output IS should be destroyed when no longer needed.
1199: Returns NULL if the stratum is empty.
1201: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1202: @*/
1203: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1204: {
1205: PetscInt v;
1211: *points = NULL;
1212: DMLabelLookupStratum(label, value, &v);
1213: if (v < 0) return(0);
1214: DMLabelMakeValid_Private(label, v);
1215: PetscObjectReference((PetscObject) label->points[v]);
1216: *points = label->points[v];
1217: return(0);
1218: }
1220: /*@
1221: DMLabelSetStratumIS - Set the stratum points using an IS
1223: Not collective
1225: Input Parameters:
1226: + label - the DMLabel
1227: . value - the stratum value
1228: - points - The stratum points
1230: Level: intermediate
1232: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1233: @*/
1234: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1235: {
1236: PetscInt v;
1242: DMLabelLookupAddStratum(label, value, &v);
1243: if (is == label->points[v]) return(0);
1244: DMLabelClearStratum(label, value);
1245: ISGetLocalSize(is, &(label->stratumSizes[v]));
1246: PetscObjectReference((PetscObject)is);
1247: ISDestroy(&(label->points[v]));
1248: label->points[v] = is;
1249: label->validIS[v] = PETSC_TRUE;
1250: PetscObjectStateIncrease((PetscObject) label);
1251: if (label->bt) {
1252: const PetscInt *points;
1253: PetscInt p;
1255: ISGetIndices(is,&points);
1256: for (p = 0; p < label->stratumSizes[v]; ++p) {
1257: const PetscInt point = points[p];
1259: if ((point < label->pStart) || (point >= label->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, label->pStart, label->pEnd);
1260: PetscBTSet(label->bt, point - label->pStart);
1261: }
1262: }
1263: return(0);
1264: }
1266: /*@
1267: DMLabelClearStratum - Remove a stratum
1269: Not collective
1271: Input Parameters:
1272: + label - the DMLabel
1273: - value - the stratum value
1275: Level: intermediate
1277: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1278: @*/
1279: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1280: {
1281: PetscInt v;
1286: DMLabelLookupStratum(label, value, &v);
1287: if (v < 0) return(0);
1288: if (label->validIS[v]) {
1289: if (label->bt) {
1290: PetscInt i;
1291: const PetscInt *points;
1293: ISGetIndices(label->points[v], &points);
1294: for (i = 0; i < label->stratumSizes[v]; ++i) {
1295: const PetscInt point = points[i];
1297: 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);
1298: PetscBTClear(label->bt, point - label->pStart);
1299: }
1300: ISRestoreIndices(label->points[v], &points);
1301: }
1302: label->stratumSizes[v] = 0;
1303: ISDestroy(&label->points[v]);
1304: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1305: PetscObjectSetName((PetscObject) label->points[v], "indices");
1306: PetscObjectStateIncrease((PetscObject) label);
1307: } else {
1308: PetscHSetIClear(label->ht[v]);
1309: }
1310: return(0);
1311: }
1313: /*@
1314: DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value
1316: Not collective
1318: Input Parameters:
1319: + label - The DMLabel
1320: . value - The label value for all points
1321: . pStart - The first point
1322: - pEnd - A point beyond all marked points
1324: Note: The marks points are [pStart, pEnd), and only the bounds are stored.
1326: Level: intermediate
1328: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1329: @*/
1330: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1331: {
1332: IS pIS;
1336: ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1337: DMLabelSetStratumIS(label, value, pIS);
1338: ISDestroy(&pIS);
1339: return(0);
1340: }
1342: /*@
1343: DMLabelGetStratumPointIndex - Get the index of a point in a given stratum
1345: Not collective
1347: Input Parameters:
1348: + label - The DMLabel
1349: . value - The label value
1350: - p - A point with this value
1352: Output Parameter:
1353: . index - The index of this point in the stratum, or -1 if the point is not in the stratum or the stratum does not exist
1355: Level: intermediate
1357: .seealso: DMLabelGetValueIndex(), DMLabelGetStratumIS(), DMLabelCreate()
1358: @*/
1359: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1360: {
1361: const PetscInt *indices;
1362: PetscInt v;
1363: PetscErrorCode ierr;
1368: *index = -1;
1369: DMLabelLookupStratum(label, value, &v);
1370: if (v < 0) return(0);
1371: DMLabelMakeValid_Private(label, v);
1372: ISGetIndices(label->points[v], &indices);
1373: PetscFindInt(p, label->stratumSizes[v], indices, index);
1374: ISRestoreIndices(label->points[v], &indices);
1375: return(0);
1376: }
1378: /*@
1379: DMLabelFilter - Remove all points outside of [start, end)
1381: Not collective
1383: Input Parameters:
1384: + label - the DMLabel
1385: . start - the first point kept
1386: - end - one more than the last point kept
1388: Level: intermediate
1390: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1391: @*/
1392: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1393: {
1394: PetscInt v;
1399: DMLabelDestroyIndex(label);
1400: DMLabelMakeAllValid_Private(label);
1401: for (v = 0; v < label->numStrata; ++v) {
1402: ISGeneralFilter(label->points[v], start, end);
1403: }
1404: DMLabelCreateIndex(label, start, end);
1405: return(0);
1406: }
1408: /*@
1409: DMLabelPermute - Create a new label with permuted points
1411: Not collective
1413: Input Parameters:
1414: + label - the DMLabel
1415: - permutation - the point permutation
1417: Output Parameter:
1418: . labelnew - the new label containing the permuted points
1420: Level: intermediate
1422: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1423: @*/
1424: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1425: {
1426: const PetscInt *perm;
1427: PetscInt numValues, numPoints, v, q;
1428: PetscErrorCode ierr;
1433: DMLabelMakeAllValid_Private(label);
1434: DMLabelDuplicate(label, labelNew);
1435: DMLabelGetNumValues(*labelNew, &numValues);
1436: ISGetLocalSize(permutation, &numPoints);
1437: ISGetIndices(permutation, &perm);
1438: for (v = 0; v < numValues; ++v) {
1439: const PetscInt size = (*labelNew)->stratumSizes[v];
1440: const PetscInt *points;
1441: PetscInt *pointsNew;
1443: ISGetIndices((*labelNew)->points[v],&points);
1444: PetscMalloc1(size,&pointsNew);
1445: for (q = 0; q < size; ++q) {
1446: const PetscInt point = points[q];
1448: 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);
1449: pointsNew[q] = perm[point];
1450: }
1451: ISRestoreIndices((*labelNew)->points[v],&points);
1452: PetscSortInt(size, pointsNew);
1453: ISDestroy(&((*labelNew)->points[v]));
1454: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1455: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1456: PetscFree(pointsNew);
1457: } else {
1458: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1459: }
1460: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1461: }
1462: ISRestoreIndices(permutation, &perm);
1463: if (label->bt) {
1464: PetscBTDestroy(&label->bt);
1465: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1466: }
1467: return(0);
1468: }
1470: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1471: {
1472: MPI_Comm comm;
1473: PetscInt s, l, nroots, nleaves, dof, offset, size;
1474: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1475: PetscSection rootSection;
1476: PetscSF labelSF;
1480: if (label) {DMLabelMakeAllValid_Private(label);}
1481: PetscObjectGetComm((PetscObject)sf, &comm);
1482: /* Build a section of stratum values per point, generate the according SF
1483: and distribute point-wise stratum values to leaves. */
1484: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1485: PetscSectionCreate(comm, &rootSection);
1486: PetscSectionSetChart(rootSection, 0, nroots);
1487: if (label) {
1488: for (s = 0; s < label->numStrata; ++s) {
1489: const PetscInt *points;
1491: ISGetIndices(label->points[s], &points);
1492: for (l = 0; l < label->stratumSizes[s]; l++) {
1493: PetscSectionGetDof(rootSection, points[l], &dof);
1494: PetscSectionSetDof(rootSection, points[l], dof+1);
1495: }
1496: ISRestoreIndices(label->points[s], &points);
1497: }
1498: }
1499: PetscSectionSetUp(rootSection);
1500: /* Create a point-wise array of stratum values */
1501: PetscSectionGetStorageSize(rootSection, &size);
1502: PetscMalloc1(size, &rootStrata);
1503: PetscCalloc1(nroots, &rootIdx);
1504: if (label) {
1505: for (s = 0; s < label->numStrata; ++s) {
1506: const PetscInt *points;
1508: ISGetIndices(label->points[s], &points);
1509: for (l = 0; l < label->stratumSizes[s]; l++) {
1510: const PetscInt p = points[l];
1511: PetscSectionGetOffset(rootSection, p, &offset);
1512: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1513: }
1514: ISRestoreIndices(label->points[s], &points);
1515: }
1516: }
1517: /* Build SF that maps label points to remote processes */
1518: PetscSectionCreate(comm, leafSection);
1519: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1520: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1521: PetscFree(remoteOffsets);
1522: /* Send the strata for each point over the derived SF */
1523: PetscSectionGetStorageSize(*leafSection, &size);
1524: PetscMalloc1(size, leafStrata);
1525: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1526: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1527: /* Clean up */
1528: PetscFree(rootStrata);
1529: PetscFree(rootIdx);
1530: PetscSectionDestroy(&rootSection);
1531: PetscSFDestroy(&labelSF);
1532: return(0);
1533: }
1535: /*@
1536: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1538: Collective on sf
1540: Input Parameters:
1541: + label - the DMLabel
1542: - sf - the map from old to new distribution
1544: Output Parameter:
1545: . labelnew - the new redistributed label
1547: Level: intermediate
1549: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1550: @*/
1551: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1552: {
1553: MPI_Comm comm;
1554: PetscSection leafSection;
1555: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1556: PetscInt *leafStrata, *strataIdx;
1557: PetscInt **points;
1558: const char *lname = NULL;
1559: char *name;
1560: PetscInt nameSize;
1561: PetscHSetI stratumHash;
1562: size_t len = 0;
1563: PetscMPIInt rank;
1568: if (label) {
1570: DMLabelMakeAllValid_Private(label);
1571: }
1572: PetscObjectGetComm((PetscObject)sf, &comm);
1573: MPI_Comm_rank(comm, &rank);
1574: /* Bcast name */
1575: if (rank == 0) {
1576: PetscObjectGetName((PetscObject) label, &lname);
1577: PetscStrlen(lname, &len);
1578: }
1579: nameSize = len;
1580: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1581: PetscMalloc1(nameSize+1, &name);
1582: if (rank == 0) {PetscArraycpy(name, lname, nameSize+1);}
1583: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1584: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1585: PetscFree(name);
1586: /* Bcast defaultValue */
1587: if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1588: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1589: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1590: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1591: /* Determine received stratum values and initialise new label*/
1592: PetscHSetICreate(&stratumHash);
1593: PetscSectionGetStorageSize(leafSection, &size);
1594: for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1595: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1596: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1597: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1598: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1599: /* Turn leafStrata into indices rather than stratum values */
1600: offset = 0;
1601: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1602: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1603: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1604: PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1605: }
1606: for (p = 0; p < size; ++p) {
1607: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1608: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1609: }
1610: }
1611: /* Rebuild the point strata on the receiver */
1612: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1613: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1614: for (p=pStart; p<pEnd; p++) {
1615: PetscSectionGetDof(leafSection, p, &dof);
1616: PetscSectionGetOffset(leafSection, p, &offset);
1617: for (s=0; s<dof; s++) {
1618: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1619: }
1620: }
1621: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1622: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1623: PetscMalloc1((*labelNew)->numStrata,&points);
1624: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1625: PetscHSetICreate(&(*labelNew)->ht[s]);
1626: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1627: }
1628: /* Insert points into new strata */
1629: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1630: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1631: for (p=pStart; p<pEnd; p++) {
1632: PetscSectionGetDof(leafSection, p, &dof);
1633: PetscSectionGetOffset(leafSection, p, &offset);
1634: for (s=0; s<dof; s++) {
1635: stratum = leafStrata[offset+s];
1636: points[stratum][strataIdx[stratum]++] = p;
1637: }
1638: }
1639: for (s = 0; s < (*labelNew)->numStrata; s++) {
1640: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1641: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1642: }
1643: PetscFree(points);
1644: PetscHSetIDestroy(&stratumHash);
1645: PetscFree(leafStrata);
1646: PetscFree(strataIdx);
1647: PetscSectionDestroy(&leafSection);
1648: return(0);
1649: }
1651: /*@
1652: DMLabelGather - Gather all label values from leafs into roots
1654: Collective on sf
1656: Input Parameters:
1657: + label - the DMLabel
1658: - sf - the Star Forest point communication map
1660: Output Parameters:
1661: . labelNew - the new DMLabel with localised leaf values
1663: Level: developer
1665: Note: This is the inverse operation to DMLabelDistribute.
1667: .seealso: DMLabelDistribute()
1668: @*/
1669: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1670: {
1671: MPI_Comm comm;
1672: PetscSection rootSection;
1673: PetscSF sfLabel;
1674: PetscSFNode *rootPoints, *leafPoints;
1675: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1676: const PetscInt *rootDegree, *ilocal;
1677: PetscInt *rootStrata;
1678: const char *lname;
1679: char *name;
1680: PetscInt nameSize;
1681: size_t len = 0;
1682: PetscMPIInt rank, size;
1688: PetscObjectGetComm((PetscObject)sf, &comm);
1689: MPI_Comm_rank(comm, &rank);
1690: MPI_Comm_size(comm, &size);
1691: /* Bcast name */
1692: if (rank == 0) {
1693: PetscObjectGetName((PetscObject) label, &lname);
1694: PetscStrlen(lname, &len);
1695: }
1696: nameSize = len;
1697: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1698: PetscMalloc1(nameSize+1, &name);
1699: if (rank == 0) {PetscArraycpy(name, lname, nameSize+1);}
1700: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1701: DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1702: PetscFree(name);
1703: /* Gather rank/index pairs of leaves into local roots to build
1704: an inverse, multi-rooted SF. Note that this ignores local leaf
1705: indexing due to the use of the multiSF in PetscSFGather. */
1706: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1707: PetscMalloc1(nroots, &leafPoints);
1708: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1709: for (p = 0; p < nleaves; p++) {
1710: PetscInt ilp = ilocal ? ilocal[p] : p;
1712: leafPoints[ilp].index = ilp;
1713: leafPoints[ilp].rank = rank;
1714: }
1715: PetscSFComputeDegreeBegin(sf, &rootDegree);
1716: PetscSFComputeDegreeEnd(sf, &rootDegree);
1717: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1718: PetscMalloc1(nmultiroots, &rootPoints);
1719: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1720: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1721: PetscSFCreate(comm,& sfLabel);
1722: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1723: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1724: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1725: /* Rebuild the point strata on the receiver */
1726: for (p = 0, idx = 0; p < nroots; p++) {
1727: for (d = 0; d < rootDegree[p]; d++) {
1728: PetscSectionGetDof(rootSection, idx+d, &dof);
1729: PetscSectionGetOffset(rootSection, idx+d, &offset);
1730: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1731: }
1732: idx += rootDegree[p];
1733: }
1734: PetscFree(leafPoints);
1735: PetscFree(rootStrata);
1736: PetscSectionDestroy(&rootSection);
1737: PetscSFDestroy(&sfLabel);
1738: return(0);
1739: }
1741: /*@
1742: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1744: Not collective
1746: Input Parameter:
1747: . label - the DMLabel
1749: Output Parameters:
1750: + section - the section giving offsets for each stratum
1751: - is - An IS containing all the label points
1753: Level: developer
1755: .seealso: DMLabelDistribute()
1756: @*/
1757: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1758: {
1759: IS vIS;
1760: const PetscInt *values;
1761: PetscInt *points;
1762: PetscInt nV, vS = 0, vE = 0, v, N;
1763: PetscErrorCode ierr;
1767: DMLabelGetNumValues(label, &nV);
1768: DMLabelGetValueIS(label, &vIS);
1769: ISGetIndices(vIS, &values);
1770: if (nV) {vS = values[0]; vE = values[0]+1;}
1771: for (v = 1; v < nV; ++v) {
1772: vS = PetscMin(vS, values[v]);
1773: vE = PetscMax(vE, values[v]+1);
1774: }
1775: PetscSectionCreate(PETSC_COMM_SELF, section);
1776: PetscSectionSetChart(*section, vS, vE);
1777: for (v = 0; v < nV; ++v) {
1778: PetscInt n;
1780: DMLabelGetStratumSize(label, values[v], &n);
1781: PetscSectionSetDof(*section, values[v], n);
1782: }
1783: PetscSectionSetUp(*section);
1784: PetscSectionGetStorageSize(*section, &N);
1785: PetscMalloc1(N, &points);
1786: for (v = 0; v < nV; ++v) {
1787: IS is;
1788: const PetscInt *spoints;
1789: PetscInt dof, off, p;
1791: PetscSectionGetDof(*section, values[v], &dof);
1792: PetscSectionGetOffset(*section, values[v], &off);
1793: DMLabelGetStratumIS(label, values[v], &is);
1794: ISGetIndices(is, &spoints);
1795: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1796: ISRestoreIndices(is, &spoints);
1797: ISDestroy(&is);
1798: }
1799: ISRestoreIndices(vIS, &values);
1800: ISDestroy(&vIS);
1801: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1802: return(0);
1803: }
1805: /*@
1806: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1807: the local section and an SF describing the section point overlap.
1809: Collective on sf
1811: Input Parameters:
1812: + s - The PetscSection for the local field layout
1813: . sf - The SF describing parallel layout of the section points
1814: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1815: . label - The label specifying the points
1816: - labelValue - The label stratum specifying the points
1818: Output Parameter:
1819: . gsection - The PetscSection for the global field layout
1821: Note: This gives negative sizes and offsets to points not owned by this process
1823: Level: developer
1825: .seealso: PetscSectionCreate()
1826: @*/
1827: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1828: {
1829: PetscInt *neg = NULL, *tmpOff = NULL;
1830: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1837: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1838: PetscSectionGetChart(s, &pStart, &pEnd);
1839: PetscSectionSetChart(*gsection, pStart, pEnd);
1840: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1841: if (nroots >= 0) {
1842: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1843: PetscCalloc1(nroots, &neg);
1844: if (nroots > pEnd-pStart) {
1845: PetscCalloc1(nroots, &tmpOff);
1846: } else {
1847: tmpOff = &(*gsection)->atlasDof[-pStart];
1848: }
1849: }
1850: /* Mark ghost points with negative dof */
1851: for (p = pStart; p < pEnd; ++p) {
1852: PetscInt value;
1854: DMLabelGetValue(label, p, &value);
1855: if (value != labelValue) continue;
1856: PetscSectionGetDof(s, p, &dof);
1857: PetscSectionSetDof(*gsection, p, dof);
1858: PetscSectionGetConstraintDof(s, p, &cdof);
1859: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1860: if (neg) neg[p] = -(dof+1);
1861: }
1862: PetscSectionSetUpBC(*gsection);
1863: if (nroots >= 0) {
1864: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1865: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1866: if (nroots > pEnd-pStart) {
1867: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1868: }
1869: }
1870: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1871: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1872: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1873: (*gsection)->atlasOff[p] = off;
1874: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1875: }
1876: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1877: globalOff -= off;
1878: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1879: (*gsection)->atlasOff[p] += globalOff;
1880: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1881: }
1882: /* Put in negative offsets for ghost points */
1883: if (nroots >= 0) {
1884: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1885: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1886: if (nroots > pEnd-pStart) {
1887: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1888: }
1889: }
1890: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1891: PetscFree(neg);
1892: return(0);
1893: }
1895: typedef struct _n_PetscSectionSym_Label
1896: {
1897: DMLabel label;
1898: PetscCopyMode *modes;
1899: PetscInt *sizes;
1900: const PetscInt ***perms;
1901: const PetscScalar ***rots;
1902: PetscInt (*minMaxOrients)[2];
1903: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1904: } PetscSectionSym_Label;
1906: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1907: {
1908: PetscInt i, j;
1909: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1910: PetscErrorCode ierr;
1913: for (i = 0; i <= sl->numStrata; i++) {
1914: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1915: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1916: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1917: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1918: }
1919: if (sl->perms[i]) {
1920: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1922: PetscFree(perms);
1923: }
1924: if (sl->rots[i]) {
1925: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1927: PetscFree(rots);
1928: }
1929: }
1930: }
1931: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1932: DMLabelDestroy(&sl->label);
1933: sl->numStrata = 0;
1934: return(0);
1935: }
1937: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1938: {
1942: PetscSectionSymLabelReset(sym);
1943: PetscFree(sym->data);
1944: return(0);
1945: }
1947: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1948: {
1949: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1950: PetscBool isAscii;
1951: DMLabel label = sl->label;
1952: const char *name;
1953: PetscErrorCode ierr;
1956: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1957: if (isAscii) {
1958: PetscInt i, j, k;
1959: PetscViewerFormat format;
1961: PetscViewerGetFormat(viewer,&format);
1962: if (label) {
1963: PetscViewerGetFormat(viewer,&format);
1964: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1965: PetscViewerASCIIPushTab(viewer);
1966: DMLabelView(label, viewer);
1967: PetscViewerASCIIPopTab(viewer);
1968: } else {
1969: PetscObjectGetName((PetscObject) sl->label, &name);
1970: PetscViewerASCIIPrintf(viewer," Label '%s'\n",name);
1971: }
1972: } else {
1973: PetscViewerASCIIPrintf(viewer, "No label given\n");
1974: }
1975: PetscViewerASCIIPushTab(viewer);
1976: for (i = 0; i <= sl->numStrata; i++) {
1977: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1979: if (!(sl->perms[i] || sl->rots[i])) {
1980: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1981: } else {
1982: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1983: PetscViewerASCIIPushTab(viewer);
1984: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1985: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1986: PetscViewerASCIIPushTab(viewer);
1987: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1988: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1989: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1990: } else {
1991: PetscInt tab;
1993: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1994: PetscViewerASCIIPushTab(viewer);
1995: PetscViewerASCIIGetTab(viewer,&tab);
1996: if (sl->perms[i] && sl->perms[i][j]) {
1997: PetscViewerASCIIPrintf(viewer,"Permutation:");
1998: PetscViewerASCIISetTab(viewer,0);
1999: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
2000: PetscViewerASCIIPrintf(viewer,"\n");
2001: PetscViewerASCIISetTab(viewer,tab);
2002: }
2003: if (sl->rots[i] && sl->rots[i][j]) {
2004: PetscViewerASCIIPrintf(viewer,"Rotations: ");
2005: PetscViewerASCIISetTab(viewer,0);
2006: #if defined(PETSC_USE_COMPLEX)
2007: 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]));}
2008: #else
2009: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
2010: #endif
2011: PetscViewerASCIIPrintf(viewer,"\n");
2012: PetscViewerASCIISetTab(viewer,tab);
2013: }
2014: PetscViewerASCIIPopTab(viewer);
2015: }
2016: }
2017: PetscViewerASCIIPopTab(viewer);
2018: }
2019: PetscViewerASCIIPopTab(viewer);
2020: }
2021: }
2022: PetscViewerASCIIPopTab(viewer);
2023: }
2024: return(0);
2025: }
2027: /*@
2028: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
2030: Logically collective on sym
2032: Input parameters:
2033: + sym - the section symmetries
2034: - label - the DMLabel describing the types of points
2036: Level: developer:
2038: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
2039: @*/
2040: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2041: {
2042: PetscSectionSym_Label *sl;
2043: PetscErrorCode ierr;
2047: sl = (PetscSectionSym_Label *) sym->data;
2048: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
2049: if (label) {
2050: sl->label = label;
2051: PetscObjectReference((PetscObject) label);
2052: DMLabelGetNumValues(label,&sl->numStrata);
2053: 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);
2054: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
2055: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
2056: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
2057: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
2058: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
2059: }
2060: return(0);
2061: }
2063: /*@C
2064: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
2066: Logically collective on sym
2068: InputParameters:
2069: + sym - the section symmetries
2070: . stratum - the stratum value in the label that we are assigning symmetries for
2071: . size - the number of dofs for points in the stratum of the label
2072: . minOrient - the smallest orientation for a point in this stratum
2073: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2074: . mode - how sym should copy the perms and rots arrays
2075: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
2076: - 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
2078: Level: developer
2080: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2081: @*/
2082: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2083: {
2084: PetscSectionSym_Label *sl;
2085: const char *name;
2086: PetscInt i, j, k;
2087: PetscErrorCode ierr;
2091: sl = (PetscSectionSym_Label *) sym->data;
2092: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
2093: for (i = 0; i <= sl->numStrata; i++) {
2094: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
2096: if (stratum == value) break;
2097: }
2098: PetscObjectGetName((PetscObject) sl->label, &name);
2099: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,name);
2100: sl->sizes[i] = size;
2101: sl->modes[i] = mode;
2102: sl->minMaxOrients[i][0] = minOrient;
2103: sl->minMaxOrients[i][1] = maxOrient;
2104: if (mode == PETSC_COPY_VALUES) {
2105: if (perms) {
2106: PetscInt **ownPerms;
2108: PetscCalloc1(maxOrient - minOrient,&ownPerms);
2109: for (j = 0; j < maxOrient-minOrient; j++) {
2110: if (perms[j]) {
2111: PetscMalloc1(size,&ownPerms[j]);
2112: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2113: }
2114: }
2115: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2116: }
2117: if (rots) {
2118: PetscScalar **ownRots;
2120: PetscCalloc1(maxOrient - minOrient,&ownRots);
2121: for (j = 0; j < maxOrient-minOrient; j++) {
2122: if (rots[j]) {
2123: PetscMalloc1(size,&ownRots[j]);
2124: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2125: }
2126: }
2127: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2128: }
2129: } else {
2130: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2131: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
2132: }
2133: return(0);
2134: }
2136: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2137: {
2138: PetscInt i, j, numStrata;
2139: PetscSectionSym_Label *sl;
2140: DMLabel label;
2141: PetscErrorCode ierr;
2144: sl = (PetscSectionSym_Label *) sym->data;
2145: numStrata = sl->numStrata;
2146: label = sl->label;
2147: for (i = 0; i < numPoints; i++) {
2148: PetscInt point = points[2*i];
2149: PetscInt ornt = points[2*i+1];
2151: for (j = 0; j < numStrata; j++) {
2152: if (label->validIS[j]) {
2153: PetscInt k;
2155: ISLocate(label->points[j],point,&k);
2156: if (k >= 0) break;
2157: } else {
2158: PetscBool has;
2160: PetscHSetIHas(label->ht[j], point, &has);
2161: if (has) break;
2162: }
2163: }
2164: 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);
2165: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2166: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2167: }
2168: return(0);
2169: }
2171: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2172: {
2173: PetscSectionSym_Label *sl;
2174: PetscErrorCode ierr;
2177: PetscNewLog(sym,&sl);
2178: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
2179: sym->ops->view = PetscSectionSymView_Label;
2180: sym->ops->destroy = PetscSectionSymDestroy_Label;
2181: sym->data = (void *) sl;
2182: return(0);
2183: }
2185: /*@
2186: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
2188: Collective
2190: Input Parameters:
2191: + comm - the MPI communicator for the new symmetry
2192: - label - the label defining the strata
2194: Output Parameters:
2195: . sym - the section symmetries
2197: Level: developer
2199: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2200: @*/
2201: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2202: {
2203: PetscErrorCode ierr;
2206: DMInitializePackage();
2207: PetscSectionSymCreate(comm,sym);
2208: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2209: PetscSectionSymLabelSetLabel(*sym,label);
2210: return(0);
2211: }