Actual source code: dmlabel.c
petsc-3.7.7 2017-09-25
1: #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/
2: #include <petsc/private/isimpl.h> /*I "petscis.h" I*/
3: #include <petscsf.h>
7: /*@C
8: DMLabelCreate - Create a DMLabel object, which is a multimap
10: Input parameter:
11: . name - The label name
13: Output parameter:
14: . label - The DMLabel
16: Level: beginner
18: .seealso: DMLabelDestroy()
19: @*/
20: PetscErrorCode DMLabelCreate(const char name[], DMLabel *label)
21: {
25: PetscNew(label);
26: PetscStrallocpy(name, &(*label)->name);
28: (*label)->refct = 1;
29: (*label)->state = -1;
30: (*label)->numStrata = 0;
31: (*label)->defaultValue = -1;
32: (*label)->stratumValues = NULL;
33: (*label)->arrayValid = 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: return(0);
41: }
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;
64: if (label->arrayValid[v]) return 0;
65: if (v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
67: PetscHashISize(label->ht[v], label->stratumSizes[v]);
69: PetscMalloc1(label->stratumSizes[v], &label->points[v]);
70: off = 0;
71: PetscHashIGetKeys(label->ht[v], &off, &(label->points[v][0]));
72: if (off != label->stratumSizes[v]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of contributed points %D from value %D should be %D", off, label->stratumValues[v], label->stratumSizes[v]);
73: PetscHashIClear(label->ht[v]);
74: PetscSortInt(label->stratumSizes[v], label->points[v]);
75: if (label->bt) {
76: PetscInt p;
78: for (p = 0; p < label->stratumSizes[v]; ++p) {
79: const PetscInt point = label->points[v][p];
81: 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);
82: PetscBTSet(label->bt, point - label->pStart);
83: }
84: }
85: label->arrayValid[v] = PETSC_TRUE;
86: ++label->state;
87: return(0);
88: }
92: /*
93: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
95: Input parameter:
96: . label - The DMLabel
98: Output parameter:
99: . label - The DMLabel with all strata in sorted list format
101: Level: developer
103: .seealso: DMLabelCreate()
104: */
105: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
106: {
107: PetscInt v;
111: for (v = 0; v < label->numStrata; v++){
112: DMLabelMakeValid_Private(label, v);
113: }
114: return(0);
115: }
119: /*
120: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
122: Input parameter:
123: + label - The DMLabel
124: - v - The stratum value
126: Output parameter:
127: . label - The DMLabel with stratum in hash format
129: Level: developer
131: .seealso: DMLabelCreate()
132: */
133: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
134: {
135: PETSC_UNUSED PetscHashIIter ret, iter;
136: PetscInt p;
140: if (!label->arrayValid[v]) return(0);
141: for (p = 0; p < label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], label->points[v][p], ret, iter);
142: PetscFree(label->points[v]);
143: label->arrayValid[v] = PETSC_FALSE;
144: return(0);
145: }
149: PetscErrorCode DMLabelGetState(DMLabel label, PetscObjectState *state)
150: {
153: *state = label->state;
154: return(0);
155: }
159: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
160: {
161: PetscInt v, *tmpV, *tmpS, **tmpP;
162: PetscHashI *tmpH;
163: PetscBool *tmpB;
168: PetscMalloc1((label->numStrata+1), &tmpV);
169: PetscMalloc1((label->numStrata+1), &tmpS);
170: PetscMalloc1((label->numStrata+1), &tmpH);
171: PetscMalloc1((label->numStrata+1), &tmpP);
172: PetscMalloc1((label->numStrata+1), &tmpB);
173: for (v = 0; v < label->numStrata; ++v) {
174: tmpV[v] = label->stratumValues[v];
175: tmpS[v] = label->stratumSizes[v];
176: tmpH[v] = label->ht[v];
177: tmpP[v] = label->points[v];
178: tmpB[v] = label->arrayValid[v];
179: }
180: tmpV[v] = value;
181: tmpS[v] = 0;
182: PetscHashICreate(tmpH[v]);
183: tmpP[v] = NULL;
184: tmpB[v] = PETSC_TRUE;
185: ++label->numStrata;
186: PetscFree(label->stratumValues);
187: PetscFree(label->stratumSizes);
188: PetscFree(label->ht);
189: PetscFree(label->points);
190: PetscFree(label->arrayValid);
191: label->stratumValues = tmpV;
192: label->stratumSizes = tmpS;
193: label->ht = tmpH;
194: label->points = tmpP;
195: label->arrayValid = tmpB;
197: return(0);
198: }
202: /*@C
203: DMLabelGetName - Return the name of a DMLabel object
205: Input parameter:
206: . label - The DMLabel
208: Output parameter:
209: . name - The label name
211: Level: beginner
213: .seealso: DMLabelCreate()
214: @*/
215: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
216: {
219: *name = label->name;
220: return(0);
221: }
225: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
226: {
227: PetscInt v;
228: PetscMPIInt rank;
232: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
233: PetscViewerASCIIPushSynchronized(viewer);
234: if (label) {
235: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
236: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
237: for (v = 0; v < label->numStrata; ++v) {
238: const PetscInt value = label->stratumValues[v];
239: PetscInt p;
241: for (p = 0; p < label->stratumSizes[v]; ++p) {
242: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, label->points[v][p], value);
243: }
244: }
245: }
246: PetscViewerFlush(viewer);
247: PetscViewerASCIIPopSynchronized(viewer);
248: return(0);
249: }
253: /*@C
254: DMLabelView - View the label
256: Input Parameters:
257: + label - The DMLabel
258: - viewer - The PetscViewer
260: Level: intermediate
262: .seealso: DMLabelCreate(), DMLabelDestroy()
263: @*/
264: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
265: {
266: PetscBool iascii;
271: if (label) {DMLabelMakeAllValid_Private(label);}
272: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
273: if (iascii) {
274: DMLabelView_Ascii(label, viewer);
275: }
276: return(0);
277: }
281: PetscErrorCode DMLabelDestroy(DMLabel *label)
282: {
283: PetscInt v;
287: if (!(*label)) return(0);
288: if (--(*label)->refct > 0) return(0);
289: PetscFree((*label)->name);
290: PetscFree((*label)->stratumValues);
291: PetscFree((*label)->stratumSizes);
292: for (v = 0; v < (*label)->numStrata; ++v) {PetscFree((*label)->points[v]);}
293: PetscFree((*label)->points);
294: PetscFree((*label)->arrayValid);
295: if ((*label)->ht) {
296: for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
297: PetscFree((*label)->ht);
298: }
299: PetscBTDestroy(&(*label)->bt);
300: PetscFree(*label);
301: return(0);
302: }
306: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
307: {
308: PetscInt v, q;
312: DMLabelMakeAllValid_Private(label);
313: PetscNew(labelnew);
314: PetscStrallocpy(label->name, &(*labelnew)->name);
316: (*labelnew)->refct = 1;
317: (*labelnew)->numStrata = label->numStrata;
318: (*labelnew)->defaultValue = label->defaultValue;
319: if (label->numStrata) {
320: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
321: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
322: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
323: PetscMalloc1(label->numStrata, &(*labelnew)->points);
324: PetscMalloc1(label->numStrata, &(*labelnew)->arrayValid);
325: /* Could eliminate unused space here */
326: for (v = 0; v < label->numStrata; ++v) {
327: PetscMalloc1(label->stratumSizes[v], &(*labelnew)->points[v]);
328: PetscHashICreate((*labelnew)->ht[v]);
329: (*labelnew)->arrayValid[v] = PETSC_TRUE;
330: (*labelnew)->stratumValues[v] = label->stratumValues[v];
331: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
332: for (q = 0; q < label->stratumSizes[v]; ++q) {
333: (*labelnew)->points[v][q] = label->points[v][q];
334: }
335: }
336: }
337: (*labelnew)->pStart = -1;
338: (*labelnew)->pEnd = -1;
339: (*labelnew)->bt = NULL;
340: return(0);
341: }
345: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
346: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
347: {
348: PetscInt v;
352: DMLabelMakeAllValid_Private(label);
353: if (label->bt) {PetscBTDestroy(&label->bt);}
354: label->pStart = pStart;
355: label->pEnd = pEnd;
356: PetscBTCreate(pEnd - pStart, &label->bt);
357: PetscBTMemzero(pEnd - pStart, label->bt);
358: for (v = 0; v < label->numStrata; ++v) {
359: PetscInt i;
361: for (i = 0; i < label->stratumSizes[v]; ++i) {
362: const PetscInt point = label->points[v][i];
364: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
365: PetscBTSet(label->bt, point - pStart);
366: }
367: }
368: return(0);
369: }
373: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
374: {
378: label->pStart = -1;
379: label->pEnd = -1;
380: if (label->bt) {PetscBTDestroy(&label->bt);}
381: return(0);
382: }
386: /*@
387: DMLabelHasValue - Determine whether a label assigns the value to any point
389: Input Parameters:
390: + label - the DMLabel
391: - value - the value
393: Output Parameter:
394: . contains - Flag indicating whether the label maps this value to any point
396: Level: developer
398: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
399: @*/
400: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
401: {
402: PetscInt v;
406: for (v = 0; v < label->numStrata; ++v) {
407: if (value == label->stratumValues[v]) break;
408: }
409: *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
410: return(0);
411: }
415: /*@
416: DMLabelHasPoint - Determine whether a label assigns a value to a point
418: Input Parameters:
419: + label - the DMLabel
420: - point - the point
422: Output Parameter:
423: . contains - Flag indicating whether the label maps this point to a value
425: Note: The user must call DMLabelCreateIndex() before this function.
427: Level: developer
429: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
430: @*/
431: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
432: {
437: DMLabelMakeAllValid_Private(label);
438: #if defined(PETSC_USE_DEBUG)
439: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
440: 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);
441: #endif
442: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
443: return(0);
444: }
448: /*@
449: DMLabelStratumHasPoint - Return true if the stratum contains a point
451: Input Parameters:
452: + label - the DMLabel
453: . value - the stratum value
454: - point - the point
456: Output Parameter:
457: . contains - true if the stratum contains the point
459: Level: intermediate
461: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
462: @*/
463: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
464: {
465: PetscInt v;
470: *contains = PETSC_FALSE;
471: for (v = 0; v < label->numStrata; ++v) {
472: if (label->stratumValues[v] == value) {
473: if (label->arrayValid[v]) {
474: PetscInt i;
476: PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
477: if (i >= 0) {
478: *contains = PETSC_TRUE;
479: break;
480: }
481: } else {
482: PetscBool has;
484: PetscHashIHasKey(label->ht[v], point, has);
485: if (has) {
486: *contains = PETSC_TRUE;
487: break;
488: }
489: }
490: }
491: }
492: return(0);
493: }
497: /*
498: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
499: When a label is created, it is initialized to -1.
501: Input parameter:
502: . label - a DMLabel object
504: Output parameter:
505: . defaultValue - the default value
507: Level: beginner
509: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
510: */
511: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
512: {
514: *defaultValue = label->defaultValue;
515: return(0);
516: }
520: /*
521: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
522: When a label is created, it is initialized to -1.
524: Input parameter:
525: . label - a DMLabel object
527: Output parameter:
528: . defaultValue - the default value
530: Level: beginner
532: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
533: */
534: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
535: {
537: label->defaultValue = defaultValue;
538: return(0);
539: }
543: /*@
544: 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())
546: Input Parameters:
547: + label - the DMLabel
548: - point - the point
550: Output Parameter:
551: . value - The point value, or -1
553: Level: intermediate
555: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
556: @*/
557: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
558: {
559: PetscInt v;
564: *value = label->defaultValue;
565: for (v = 0; v < label->numStrata; ++v) {
566: if (label->arrayValid[v]) {
567: PetscInt i;
569: PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
570: if (i >= 0) {
571: *value = label->stratumValues[v];
572: break;
573: }
574: } else {
575: PetscBool has;
577: PetscHashIHasKey(label->ht[v], point, has);
578: if (has) {
579: *value = label->stratumValues[v];
580: break;
581: }
582: }
583: }
584: return(0);
585: }
589: /*@
590: 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 somethingg different), then this function will do nothing.
592: Input Parameters:
593: + label - the DMLabel
594: . point - the point
595: - value - The point value
597: Level: intermediate
599: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
600: @*/
601: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
602: {
603: PETSC_UNUSED PetscHashIIter iter, ret;
604: PetscInt v;
605: PetscErrorCode ierr;
608: /* Find, or add, label value */
609: if (value == label->defaultValue) return(0);
610: for (v = 0; v < label->numStrata; ++v) {
611: if (label->stratumValues[v] == value) break;
612: }
613: /* Create new table */
614: if (v >= label->numStrata) {DMLabelAddStratum(label, value);}
615: DMLabelMakeInvalid_Private(label, v);
616: /* Set key */
617: PetscHashIPut(label->ht[v], point, ret, iter);
618: return(0);
619: }
623: /*@
624: DMLabelClearValue - Clear the value a label assigns to a point
626: Input Parameters:
627: + label - the DMLabel
628: . point - the point
629: - value - The point value
631: Level: intermediate
633: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
634: @*/
635: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
636: {
637: PetscInt v, p;
641: /* Find label value */
642: for (v = 0; v < label->numStrata; ++v) {
643: if (label->stratumValues[v] == value) break;
644: }
645: if (v >= label->numStrata) return(0);
646: if (label->arrayValid[v]) {
647: /* Check whether point exists */
648: PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &p);
649: if (p >= 0) {
650: PetscMemmove(&label->points[v][p], &label->points[v][p+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
651: --label->stratumSizes[v];
652: if (label->bt) {
653: 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);
654: PetscBTClear(label->bt, point - label->pStart);
655: }
656: }
657: } else {
658: PetscHashIDelKey(label->ht[v], point);
659: }
660: return(0);
661: }
665: /*@
666: DMLabelInsertIS - Set all points in the IS to a value
668: Input Parameters:
669: + label - the DMLabel
670: . is - the point IS
671: - value - The point value
673: Level: intermediate
675: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
676: @*/
677: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
678: {
679: const PetscInt *points;
680: PetscInt n, p;
681: PetscErrorCode ierr;
685: ISGetLocalSize(is, &n);
686: ISGetIndices(is, &points);
687: for (p = 0; p < n; ++p) {DMLabelSetValue(label, points[p], value);}
688: ISRestoreIndices(is, &points);
689: return(0);
690: }
694: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
695: {
698: *numValues = label->numStrata;
699: return(0);
700: }
704: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
705: {
710: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
711: return(0);
712: }
716: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
717: {
718: PetscInt v;
722: *exists = PETSC_FALSE;
723: for (v = 0; v < label->numStrata; ++v) {
724: if (label->stratumValues[v] == value) {
725: *exists = PETSC_TRUE;
726: break;
727: }
728: }
729: return(0);
730: }
734: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
735: {
736: PetscInt v;
741: *size = 0;
742: for (v = 0; v < label->numStrata; ++v) {
743: if (label->stratumValues[v] == value) {
744: DMLabelMakeValid_Private(label, v);
745: *size = label->stratumSizes[v];
746: break;
747: }
748: }
749: return(0);
750: }
754: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
755: {
756: PetscInt v;
762: for (v = 0; v < label->numStrata; ++v) {
763: if (label->stratumValues[v] != value) continue;
764: DMLabelMakeValid_Private(label, v);
765: if (label->stratumSizes[v] <= 0) break;
766: if (start) *start = label->points[v][0];
767: if (end) *end = label->points[v][label->stratumSizes[v]-1]+1;
768: break;
769: }
770: return(0);
771: }
775: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
776: {
777: PetscInt v;
782: *points = NULL;
783: for (v = 0; v < label->numStrata; ++v) {
784: if (label->stratumValues[v] == value) {
785: DMLabelMakeValid_Private(label, v);
786: if (label->arrayValid[v]) {
787: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[v][0], PETSC_COPY_VALUES, points);
788: PetscObjectSetName((PetscObject) *points, "indices");
789: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
790: break;
791: }
792: }
793: return(0);
794: }
798: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
799: {
800: PetscInt v;
804: for (v = 0; v < label->numStrata; ++v) {
805: if (label->stratumValues[v] == value) break;
806: }
807: if (v >= label->numStrata) return(0);
808: if (label->bt) {
809: PetscInt i;
811: for (i = 0; i < label->stratumSizes[v]; ++i) {
812: const PetscInt point = label->points[v][i];
814: 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);
815: PetscBTClear(label->bt, point - label->pStart);
816: }
817: }
818: if (label->arrayValid[v]) {
819: label->stratumSizes[v] = 0;
820: } else {
821: PetscHashIClear(label->ht[v]);
822: }
823: return(0);
824: }
828: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
829: {
830: PetscInt v;
834: DMLabelMakeAllValid_Private(label);
835: label->pStart = start;
836: label->pEnd = end;
837: if (label->bt) {PetscBTDestroy(&label->bt);}
838: /* Could squish offsets, but would only make sense if I reallocate the storage */
839: for (v = 0; v < label->numStrata; ++v) {
840: PetscInt off, q;
842: for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
843: const PetscInt point = label->points[v][q];
845: if ((point < start) || (point >= end)) continue;
846: label->points[v][off++] = point;
847: }
848: label->stratumSizes[v] = off;
849: }
850: DMLabelCreateIndex(label, start, end);
851: return(0);
852: }
856: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
857: {
858: const PetscInt *perm;
859: PetscInt numValues, numPoints, v, q;
860: PetscErrorCode ierr;
863: DMLabelMakeAllValid_Private(label);
864: DMLabelDuplicate(label, labelNew);
865: DMLabelGetNumValues(*labelNew, &numValues);
866: ISGetLocalSize(permutation, &numPoints);
867: ISGetIndices(permutation, &perm);
868: for (v = 0; v < numValues; ++v) {
869: const PetscInt size = (*labelNew)->stratumSizes[v];
871: for (q = 0; q < size; ++q) {
872: const PetscInt point = (*labelNew)->points[v][q];
874: 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);
875: (*labelNew)->points[v][q] = perm[point];
876: }
877: PetscSortInt(size, &(*labelNew)->points[v][0]);
878: }
879: ISRestoreIndices(permutation, &perm);
880: if (label->bt) {
881: PetscBTDestroy(&label->bt);
882: DMLabelCreateIndex(label, label->pStart, label->pEnd);
883: }
884: return(0);
885: }
889: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
890: {
891: MPI_Comm comm;
892: PetscInt s, l, nroots, nleaves, dof, offset, size;
893: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
894: PetscSection rootSection;
895: PetscSF labelSF;
899: if (label) {DMLabelMakeAllValid_Private(label);}
900: PetscObjectGetComm((PetscObject)sf, &comm);
901: /* Build a section of stratum values per point, generate the according SF
902: and distribute point-wise stratum values to leaves. */
903: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
904: PetscSectionCreate(comm, &rootSection);
905: PetscSectionSetChart(rootSection, 0, nroots);
906: if (label) {
907: for (s = 0; s < label->numStrata; ++s) {
908: for (l = 0; l < label->stratumSizes[s]; l++) {
909: PetscSectionGetDof(rootSection, label->points[s][l], &dof);
910: PetscSectionSetDof(rootSection, label->points[s][l], dof+1);
911: }
912: }
913: }
914: PetscSectionSetUp(rootSection);
915: /* Create a point-wise array of stratum values */
916: PetscSectionGetStorageSize(rootSection, &size);
917: PetscMalloc1(size, &rootStrata);
918: PetscCalloc1(nroots, &rootIdx);
919: if (label) {
920: for (s = 0; s < label->numStrata; ++s) {
921: for (l = 0; l < label->stratumSizes[s]; l++) {
922: const PetscInt p = label->points[s][l];
923: PetscSectionGetOffset(rootSection, p, &offset);
924: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
925: }
926: }
927: }
928: /* Build SF that maps label points to remote processes */
929: PetscSectionCreate(comm, leafSection);
930: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
931: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
932: PetscFree(remoteOffsets);
933: /* Send the strata for each point over the derived SF */
934: PetscSectionGetStorageSize(*leafSection, &size);
935: PetscMalloc1(size, leafStrata);
936: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
937: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
938: /* Clean up */
939: PetscFree(rootStrata);
940: PetscFree(rootIdx);
941: PetscSectionDestroy(&rootSection);
942: PetscSFDestroy(&labelSF);
943: return(0);
944: }
948: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
949: {
950: MPI_Comm comm;
951: PetscSection leafSection;
952: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
953: PetscInt *leafStrata, *strataIdx;
954: char *name;
955: PetscInt nameSize;
956: PetscHashI stratumHash;
957: PETSC_UNUSED PetscHashIIter ret, iter;
958: size_t len = 0;
959: PetscMPIInt rank;
963: if (label) {DMLabelMakeAllValid_Private(label);}
964: PetscObjectGetComm((PetscObject)sf, &comm);
965: MPI_Comm_rank(comm, &rank);
966: /* Bcast name */
967: if (!rank) {PetscStrlen(label->name, &len);}
968: nameSize = len;
969: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
970: PetscMalloc1(nameSize+1, &name);
971: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
972: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
973: DMLabelCreate(name, labelNew);
974: PetscFree(name);
975: /* Bcast defaultValue */
976: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
977: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
978: /* Distribute stratum values over the SF and get the point mapping on the receiver */
979: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
980: /* Determine received stratum values and initialise new label*/
981: PetscHashICreate(stratumHash);
982: PetscSectionGetStorageSize(leafSection, &size);
983: for (p = 0; p < size; ++p) PetscHashIPut(stratumHash, leafStrata[p], ret, iter);
984: PetscHashISize(stratumHash, (*labelNew)->numStrata);
985: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->arrayValid);
986: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->arrayValid[s] = PETSC_TRUE;
987: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
988: /* Turn leafStrata into indices rather than stratum values */
989: offset = 0;
990: PetscHashIGetKeys(stratumHash, &offset, (*labelNew)->stratumValues);
991: for (p = 0; p < size; ++p) {
992: for (s = 0; s < (*labelNew)->numStrata; ++s) {
993: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
994: }
995: }
996: /* Rebuild the point strata on the receiver */
997: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
998: PetscSectionGetChart(leafSection, &pStart, &pEnd);
999: for (p=pStart; p<pEnd; p++) {
1000: PetscSectionGetDof(leafSection, p, &dof);
1001: PetscSectionGetOffset(leafSection, p, &offset);
1002: for (s=0; s<dof; s++) {
1003: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1004: }
1005: }
1006: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1007: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1008: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1009: PetscHashICreate((*labelNew)->ht[s]);
1010: PetscMalloc1((*labelNew)->stratumSizes[s], &(*labelNew)->points[s]);
1011: }
1012: /* Insert points into new strata */
1013: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1014: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1015: for (p=pStart; p<pEnd; p++) {
1016: PetscSectionGetDof(leafSection, p, &dof);
1017: PetscSectionGetOffset(leafSection, p, &offset);
1018: for (s=0; s<dof; s++) {
1019: stratum = leafStrata[offset+s];
1020: (*labelNew)->points[stratum][strataIdx[stratum]++] = p;
1021: }
1022: }
1023: PetscHashIDestroy(stratumHash);
1024: PetscFree(leafStrata);
1025: PetscFree(strataIdx);
1026: PetscSectionDestroy(&leafSection);
1027: return(0);
1028: }
1032: /*@
1033: DMLabelGather - Gather all label values from leafs into roots
1035: Input Parameters:
1036: + label - the DMLabel
1037: . point - the Star Forest point communication map
1039: Input Parameters:
1040: + label - the new DMLabel with localised leaf values
1042: Level: developer
1044: Note: This is the inverse operation to DMLabelDistribute.
1046: .seealso: DMLabelDistribute()
1047: @*/
1048: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1049: {
1050: MPI_Comm comm;
1051: PetscSection rootSection;
1052: PetscSF sfLabel;
1053: PetscSFNode *rootPoints, *leafPoints;
1054: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1055: const PetscInt *rootDegree, *ilocal;
1056: PetscInt *rootStrata;
1057: char *name;
1058: PetscInt nameSize;
1059: size_t len = 0;
1060: PetscMPIInt rank, numProcs;
1064: PetscObjectGetComm((PetscObject)sf, &comm);
1065: MPI_Comm_rank(comm, &rank);
1066: MPI_Comm_size(comm, &numProcs);
1067: /* Bcast name */
1068: if (!rank) {PetscStrlen(label->name, &len);}
1069: nameSize = len;
1070: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1071: PetscMalloc1(nameSize+1, &name);
1072: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1073: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1074: DMLabelCreate(name, labelNew);
1075: PetscFree(name);
1076: /* Gather rank/index pairs of leaves into local roots to build
1077: an inverse, multi-rooted SF. Note that this ignores local leaf
1078: indexing due to the use of the multiSF in PetscSFGather. */
1079: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1080: PetscMalloc1(nleaves, &leafPoints);
1081: for (p = 0; p < nleaves; p++) {
1082: leafPoints[p].index = ilocal[p];
1083: leafPoints[p].rank = rank;
1084: }
1085: PetscSFComputeDegreeBegin(sf, &rootDegree);
1086: PetscSFComputeDegreeEnd(sf, &rootDegree);
1087: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1088: PetscMalloc1(nmultiroots, &rootPoints);
1089: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1090: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1091: PetscSFCreate(comm,& sfLabel);
1092: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1093: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1094: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1095: /* Rebuild the point strata on the receiver */
1096: for (p = 0, idx = 0; p < nroots; p++) {
1097: for (d = 0; d < rootDegree[p]; d++) {
1098: PetscSectionGetDof(rootSection, idx+d, &dof);
1099: PetscSectionGetOffset(rootSection, idx+d, &offset);
1100: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1101: }
1102: idx += rootDegree[p];
1103: }
1104: PetscFree(leafPoints);
1105: PetscFree(rootStrata);
1106: PetscSectionDestroy(&rootSection);
1107: PetscSFDestroy(&sfLabel);
1108: return(0);
1109: }
1113: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1114: {
1115: IS vIS;
1116: const PetscInt *values;
1117: PetscInt *points;
1118: PetscInt nV, vS = 0, vE = 0, v, N;
1119: PetscErrorCode ierr;
1122: DMLabelGetNumValues(label, &nV);
1123: DMLabelGetValueIS(label, &vIS);
1124: ISGetIndices(vIS, &values);
1125: if (nV) {vS = values[0]; vE = values[0]+1;}
1126: for (v = 1; v < nV; ++v) {
1127: vS = PetscMin(vS, values[v]);
1128: vE = PetscMax(vE, values[v]+1);
1129: }
1130: PetscSectionCreate(PETSC_COMM_SELF, section);
1131: PetscSectionSetChart(*section, vS, vE);
1132: for (v = 0; v < nV; ++v) {
1133: PetscInt n;
1135: DMLabelGetStratumSize(label, values[v], &n);
1136: PetscSectionSetDof(*section, values[v], n);
1137: }
1138: PetscSectionSetUp(*section);
1139: PetscSectionGetStorageSize(*section, &N);
1140: PetscMalloc1(N, &points);
1141: for (v = 0; v < nV; ++v) {
1142: IS is;
1143: const PetscInt *spoints;
1144: PetscInt dof, off, p;
1146: PetscSectionGetDof(*section, values[v], &dof);
1147: PetscSectionGetOffset(*section, values[v], &off);
1148: DMLabelGetStratumIS(label, values[v], &is);
1149: ISGetIndices(is, &spoints);
1150: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1151: ISRestoreIndices(is, &spoints);
1152: ISDestroy(&is);
1153: }
1154: ISRestoreIndices(vIS, &values);
1155: ISDestroy(&vIS);
1156: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1157: return(0);
1158: }
1162: /*@C
1163: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1164: the local section and an SF describing the section point overlap.
1166: Input Parameters:
1167: + s - The PetscSection for the local field layout
1168: . sf - The SF describing parallel layout of the section points
1169: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1170: . label - The label specifying the points
1171: - labelValue - The label stratum specifying the points
1173: Output Parameter:
1174: . gsection - The PetscSection for the global field layout
1176: Note: This gives negative sizes and offsets to points not owned by this process
1178: Level: developer
1180: .seealso: PetscSectionCreate()
1181: @*/
1182: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1183: {
1184: PetscInt *neg = NULL, *tmpOff = NULL;
1185: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1189: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1190: PetscSectionGetChart(s, &pStart, &pEnd);
1191: PetscSectionSetChart(*gsection, pStart, pEnd);
1192: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1193: if (nroots >= 0) {
1194: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1195: PetscCalloc1(nroots, &neg);
1196: if (nroots > pEnd-pStart) {
1197: PetscCalloc1(nroots, &tmpOff);
1198: } else {
1199: tmpOff = &(*gsection)->atlasDof[-pStart];
1200: }
1201: }
1202: /* Mark ghost points with negative dof */
1203: for (p = pStart; p < pEnd; ++p) {
1204: PetscInt value;
1206: DMLabelGetValue(label, p, &value);
1207: if (value != labelValue) continue;
1208: PetscSectionGetDof(s, p, &dof);
1209: PetscSectionSetDof(*gsection, p, dof);
1210: PetscSectionGetConstraintDof(s, p, &cdof);
1211: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1212: if (neg) neg[p] = -(dof+1);
1213: }
1214: PetscSectionSetUpBC(*gsection);
1215: if (nroots >= 0) {
1216: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1217: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1218: if (nroots > pEnd-pStart) {
1219: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1220: }
1221: }
1222: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1223: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1224: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1225: (*gsection)->atlasOff[p] = off;
1226: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1227: }
1228: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1229: globalOff -= off;
1230: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1231: (*gsection)->atlasOff[p] += globalOff;
1232: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1233: }
1234: /* Put in negative offsets for ghost points */
1235: if (nroots >= 0) {
1236: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1237: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1238: if (nroots > pEnd-pStart) {
1239: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1240: }
1241: }
1242: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1243: PetscFree(neg);
1244: return(0);
1245: }