Actual source code: dmlabel.c
petsc-3.8.4 2018-03-24
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petscsf.h>
6: /*@C
7: DMLabelCreate - Create a DMLabel object, which is a multimap
9: Input parameter:
10: . name - The label name
12: Output parameter:
13: . label - The DMLabel
15: Level: beginner
17: .seealso: DMLabelDestroy()
18: @*/
19: PetscErrorCode DMLabelCreate(const char name[], DMLabel *label)
20: {
24: PetscNew(label);
25: PetscStrallocpy(name, &(*label)->name);
27: (*label)->refct = 1;
28: (*label)->state = -1;
29: (*label)->numStrata = 0;
30: (*label)->defaultValue = -1;
31: (*label)->stratumValues = NULL;
32: (*label)->validIS = NULL;
33: (*label)->stratumSizes = NULL;
34: (*label)->points = NULL;
35: (*label)->ht = NULL;
36: (*label)->pStart = -1;
37: (*label)->pEnd = -1;
38: (*label)->bt = NULL;
39: return(0);
40: }
42: /*
43: DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format
45: Input parameter:
46: + label - The DMLabel
47: - v - The stratum value
49: Output parameter:
50: . label - The DMLabel with stratum in sorted list format
52: Level: developer
54: .seealso: DMLabelCreate()
55: */
56: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
57: {
58: PetscInt off;
59: PetscInt *pointArray;
62: if (label->validIS[v]) return 0;
63: if (v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
65: PetscHashISize(label->ht[v], label->stratumSizes[v]);
67: PetscMalloc1(label->stratumSizes[v], &pointArray);
68: off = 0;
69: PetscHashIGetKeys(label->ht[v], &off, pointArray);
70: 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]);
71: PetscHashIClear(label->ht[v]);
72: PetscSortInt(label->stratumSizes[v], pointArray);
73: if (label->bt) {
74: PetscInt p;
76: for (p = 0; p < label->stratumSizes[v]; ++p) {
77: const PetscInt point = pointArray[p];
79: 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);
80: PetscBTSet(label->bt, point - label->pStart);
81: }
82: }
83: ISCreateGeneral(PETSC_COMM_SELF,label->stratumSizes[v],pointArray,PETSC_OWN_POINTER,&(label->points[v]));
84: PetscObjectSetName((PetscObject) (label->points[v]), "indices");
85: label->validIS[v] = PETSC_TRUE;
86: ++label->state;
87: return(0);
88: }
90: /*
91: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
93: Input parameter:
94: . label - The DMLabel
96: Output parameter:
97: . label - The DMLabel with all strata in sorted list format
99: Level: developer
101: .seealso: DMLabelCreate()
102: */
103: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
104: {
105: PetscInt v;
109: for (v = 0; v < label->numStrata; v++){
110: DMLabelMakeValid_Private(label, v);
111: }
112: return(0);
113: }
115: /*
116: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
118: Input parameter:
119: + label - The DMLabel
120: - v - The stratum value
122: Output parameter:
123: . label - The DMLabel with stratum in hash format
125: Level: developer
127: .seealso: DMLabelCreate()
128: */
129: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
130: {
131: PETSC_UNUSED PetscHashIIter ret, iter;
132: PetscInt p;
133: const PetscInt *points;
137: if (!label->validIS[v]) return(0);
138: if (label->points[v]) {
139: ISGetIndices(label->points[v],&points);
140: for (p = 0; p < label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], points[p], ret, iter);
141: ISRestoreIndices(label->points[v],&points);
142: ISDestroy(&(label->points[v]));
143: }
144: label->validIS[v] = PETSC_FALSE;
145: return(0);
146: }
148: PetscErrorCode DMLabelGetState(DMLabel label, PetscObjectState *state)
149: {
152: *state = label->state;
153: return(0);
154: }
156: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
157: {
158: PetscInt v, *tmpV, *tmpS;
159: IS *tmpP;
160: PetscHashI *tmpH;
161: PetscBool *tmpB;
166: for (v = 0; v < label->numStrata; v++) {
167: if (label->stratumValues[v] == value) return(0);
168: }
169: PetscMalloc1((label->numStrata+1), &tmpV);
170: PetscMalloc1((label->numStrata+1), &tmpS);
171: PetscMalloc1((label->numStrata+1), &tmpH);
172: PetscMalloc1((label->numStrata+1), &tmpP);
173: PetscMalloc1((label->numStrata+1), &tmpB);
174: for (v = 0; v < label->numStrata; ++v) {
175: tmpV[v] = label->stratumValues[v];
176: tmpS[v] = label->stratumSizes[v];
177: tmpH[v] = label->ht[v];
178: tmpP[v] = label->points[v];
179: tmpB[v] = label->validIS[v];
180: }
181: tmpV[v] = value;
182: tmpS[v] = 0;
183: PetscHashICreate(tmpH[v]);
184: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_OWN_POINTER,&tmpP[v]);
185: tmpB[v] = PETSC_TRUE;
186: ++label->numStrata;
187: PetscFree(label->stratumValues);
188: PetscFree(label->stratumSizes);
189: PetscFree(label->ht);
190: PetscFree(label->points);
191: PetscFree(label->validIS);
192: label->stratumValues = tmpV;
193: label->stratumSizes = tmpS;
194: label->ht = tmpH;
195: label->points = tmpP;
196: label->validIS = tmpB;
198: return(0);
199: }
201: /*@C
202: DMLabelGetName - Return the name of a DMLabel object
204: Input parameter:
205: . label - The DMLabel
207: Output parameter:
208: . name - The label name
210: Level: beginner
212: .seealso: DMLabelCreate()
213: @*/
214: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
215: {
218: *name = label->name;
219: return(0);
220: }
222: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
223: {
224: PetscInt v;
225: PetscMPIInt rank;
229: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
230: PetscViewerASCIIPushSynchronized(viewer);
231: if (label) {
232: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
233: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
234: for (v = 0; v < label->numStrata; ++v) {
235: const PetscInt value = label->stratumValues[v];
236: const PetscInt *points;
237: PetscInt p;
239: ISGetIndices(label->points[v],&points);
240: for (p = 0; p < label->stratumSizes[v]; ++p) {
241: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
242: }
243: ISRestoreIndices(label->points[v],&points);
244: }
245: }
246: PetscViewerFlush(viewer);
247: PetscViewerASCIIPopSynchronized(viewer);
248: return(0);
249: }
251: /*@C
252: DMLabelView - View the label
254: Input Parameters:
255: + label - The DMLabel
256: - viewer - The PetscViewer
258: Level: intermediate
260: .seealso: DMLabelCreate(), DMLabelDestroy()
261: @*/
262: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
263: {
264: PetscBool iascii;
269: if (label) {DMLabelMakeAllValid_Private(label);}
270: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
271: if (iascii) {
272: DMLabelView_Ascii(label, viewer);
273: }
274: return(0);
275: }
277: /*@
278: DMLabelDestroy - Destroys a DMLabel
280: Input Parameter:
281: . label - The DMLabel
283: Level: beginner
285: .seealso: DMLabelCreate()
286: @*/
287: PetscErrorCode DMLabelDestroy(DMLabel *label)
288: {
289: PetscInt v;
293: if (!(*label)) return(0);
294: if (--(*label)->refct > 0) return(0);
295: PetscFree((*label)->name);
296: PetscFree((*label)->stratumValues);
297: PetscFree((*label)->stratumSizes);
298: for (v = 0; v < (*label)->numStrata; ++v) {ISDestroy(&((*label)->points[v]));}
299: PetscFree((*label)->points);
300: PetscFree((*label)->validIS);
301: if ((*label)->ht) {
302: for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
303: PetscFree((*label)->ht);
304: }
305: PetscBTDestroy(&(*label)->bt);
306: PetscFree(*label);
307: return(0);
308: }
310: /*@
311: DMLabelDuplicate - Duplicates a DMLabel
313: Input Parameter:
314: . label - The DMLabel
316: Output Parameter:
317: . labelnew - location to put new vector
319: Level: intermediate
321: .seealso: DMLabelCreate(), DMLabelDestroy()
322: @*/
323: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
324: {
325: PetscInt v;
329: DMLabelMakeAllValid_Private(label);
330: PetscNew(labelnew);
331: PetscStrallocpy(label->name, &(*labelnew)->name);
333: (*labelnew)->refct = 1;
334: (*labelnew)->numStrata = label->numStrata;
335: (*labelnew)->defaultValue = label->defaultValue;
336: if (label->numStrata) {
337: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
338: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
339: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
340: PetscMalloc1(label->numStrata, &(*labelnew)->points);
341: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
342: /* Could eliminate unused space here */
343: for (v = 0; v < label->numStrata; ++v) {
344: PetscHashICreate((*labelnew)->ht[v]);
345: (*labelnew)->validIS[v] = PETSC_TRUE;
346: (*labelnew)->stratumValues[v] = label->stratumValues[v];
347: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
348: PetscObjectReference((PetscObject) (label->points[v]));
349: (*labelnew)->points[v] = label->points[v];
350: }
351: }
352: (*labelnew)->pStart = -1;
353: (*labelnew)->pEnd = -1;
354: (*labelnew)->bt = NULL;
355: return(0);
356: }
358: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
359: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
360: {
361: PetscInt v;
365: DMLabelMakeAllValid_Private(label);
366: if (label->bt) {PetscBTDestroy(&label->bt);}
367: label->pStart = pStart;
368: label->pEnd = pEnd;
369: PetscBTCreate(pEnd - pStart, &label->bt);
370: PetscBTMemzero(pEnd - pStart, label->bt);
371: for (v = 0; v < label->numStrata; ++v) {
372: const PetscInt *points;
373: PetscInt i;
375: ISGetIndices(label->points[v],&points);
376: for (i = 0; i < label->stratumSizes[v]; ++i) {
377: const PetscInt point = points[i];
379: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
380: PetscBTSet(label->bt, point - pStart);
381: }
382: ISRestoreIndices(label->points[v],&points);
383: }
384: return(0);
385: }
387: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
388: {
392: label->pStart = -1;
393: label->pEnd = -1;
394: if (label->bt) {PetscBTDestroy(&label->bt);}
395: return(0);
396: }
398: /*@
399: DMLabelHasValue - Determine whether a label assigns the value to any point
401: Input Parameters:
402: + label - the DMLabel
403: - value - the value
405: Output Parameter:
406: . contains - Flag indicating whether the label maps this value to any point
408: Level: developer
410: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
411: @*/
412: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
413: {
414: PetscInt v;
418: for (v = 0; v < label->numStrata; ++v) {
419: if (value == label->stratumValues[v]) break;
420: }
421: *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
422: return(0);
423: }
425: /*@
426: DMLabelHasPoint - Determine whether a label assigns a value to a point
428: Input Parameters:
429: + label - the DMLabel
430: - point - the point
432: Output Parameter:
433: . contains - Flag indicating whether the label maps this point to a value
435: Note: The user must call DMLabelCreateIndex() before this function.
437: Level: developer
439: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
440: @*/
441: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
442: {
447: DMLabelMakeAllValid_Private(label);
448: #if defined(PETSC_USE_DEBUG)
449: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
450: 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);
451: #endif
452: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
453: return(0);
454: }
456: /*@
457: DMLabelStratumHasPoint - Return true if the stratum contains a point
459: Input Parameters:
460: + label - the DMLabel
461: . value - the stratum value
462: - point - the point
464: Output Parameter:
465: . contains - true if the stratum contains the point
467: Level: intermediate
469: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
470: @*/
471: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
472: {
473: PetscInt v;
478: *contains = PETSC_FALSE;
479: for (v = 0; v < label->numStrata; ++v) {
480: if (label->stratumValues[v] == value) {
481: if (label->validIS[v]) {
482: PetscInt i;
484: ISLocate(label->points[v],point,&i);
485: if (i >= 0) {
486: *contains = PETSC_TRUE;
487: break;
488: }
489: } else {
490: PetscBool has;
492: PetscHashIHasKey(label->ht[v], point, has);
493: if (has) {
494: *contains = PETSC_TRUE;
495: break;
496: }
497: }
498: }
499: }
500: return(0);
501: }
503: /*@
504: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
505: When a label is created, it is initialized to -1.
507: Input parameter:
508: . label - a DMLabel object
510: Output parameter:
511: . defaultValue - the default value
513: Level: beginner
515: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
516: @*/
517: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
518: {
520: *defaultValue = label->defaultValue;
521: return(0);
522: }
524: /*@
525: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
526: When a label is created, it is initialized to -1.
528: Input parameter:
529: . label - a DMLabel object
531: Output parameter:
532: . defaultValue - the default value
534: Level: beginner
536: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
537: @*/
538: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
539: {
541: label->defaultValue = defaultValue;
542: return(0);
543: }
545: /*@
546: 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())
548: Input Parameters:
549: + label - the DMLabel
550: - point - the point
552: Output Parameter:
553: . value - The point value, or -1
555: Level: intermediate
557: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
558: @*/
559: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
560: {
561: PetscInt v;
566: *value = label->defaultValue;
567: for (v = 0; v < label->numStrata; ++v) {
568: if (label->validIS[v]) {
569: PetscInt i;
571: ISLocate(label->points[v],point,&i);
572: if (i >= 0) {
573: *value = label->stratumValues[v];
574: break;
575: }
576: } else {
577: PetscBool has;
579: PetscHashIHasKey(label->ht[v], point, has);
580: if (has) {
581: *value = label->stratumValues[v];
582: break;
583: }
584: }
585: }
586: return(0);
587: }
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: }
621: /*@
622: DMLabelClearValue - Clear the value a label assigns to a point
624: Input Parameters:
625: + label - the DMLabel
626: . point - the point
627: - value - The point value
629: Level: intermediate
631: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
632: @*/
633: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
634: {
635: PetscInt v;
639: /* Find label value */
640: for (v = 0; v < label->numStrata; ++v) {
641: if (label->stratumValues[v] == value) break;
642: }
643: if (v >= label->numStrata) return(0);
644: if (label->validIS[v]) {
645: DMLabelMakeInvalid_Private(label,v);
646: }
647: if (label->bt) {
648: 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);
649: PetscBTClear(label->bt, point - label->pStart);
650: }
651: PetscHashIDelKey(label->ht[v], point);
652: return(0);
653: }
655: /*@
656: DMLabelInsertIS - Set all points in the IS to a value
658: Input Parameters:
659: + label - the DMLabel
660: . is - the point IS
661: - value - The point value
663: Level: intermediate
665: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
666: @*/
667: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
668: {
669: const PetscInt *points;
670: PetscInt n, p;
671: PetscErrorCode ierr;
675: ISGetLocalSize(is, &n);
676: ISGetIndices(is, &points);
677: for (p = 0; p < n; ++p) {DMLabelSetValue(label, points[p], value);}
678: ISRestoreIndices(is, &points);
679: return(0);
680: }
682: /*@
683: DMLabelGetNumValues - Get the number of values that the DMLabel takes
685: Input Parameter:
686: . label - the DMLabel
688: Output Paramater:
689: . numValues - the number of values
691: Level: intermediate
693: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
694: @*/
695: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
696: {
699: *numValues = label->numStrata;
700: return(0);
701: }
703: /*@
704: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
706: Input Parameter:
707: . label - the DMLabel
709: Output Paramater:
710: . is - the value IS
712: Level: intermediate
714: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
715: @*/
716: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
717: {
722: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
723: return(0);
724: }
726: /*@
727: DMLabelHasStratum - Determine whether points exist with the given value
729: Input Parameters:
730: + label - the DMLabel
731: - value - the stratum value
733: Output Paramater:
734: . exists - Flag saying whether points exist
736: Level: intermediate
738: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
739: @*/
740: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
741: {
742: PetscInt v;
746: *exists = PETSC_FALSE;
747: for (v = 0; v < label->numStrata; ++v) {
748: if (label->stratumValues[v] == value) {
749: *exists = PETSC_TRUE;
750: break;
751: }
752: }
753: return(0);
754: }
756: /*@
757: DMLabelGetStratumSize - Get the size of a stratum
759: Input Parameters:
760: + label - the DMLabel
761: - value - the stratum value
763: Output Paramater:
764: . size - The number of points in the stratum
766: Level: intermediate
768: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
769: @*/
770: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
771: {
772: PetscInt v;
777: *size = 0;
778: for (v = 0; v < label->numStrata; ++v) {
779: if (label->stratumValues[v] == value) {
780: DMLabelMakeValid_Private(label, v);
781: *size = label->stratumSizes[v];
782: break;
783: }
784: }
785: return(0);
786: }
788: /*@
789: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
791: Input Parameters:
792: + label - the DMLabel
793: - value - the stratum value
795: Output Paramaters:
796: + start - the smallest point in the stratum
797: - end - the largest point in the stratum
799: Level: intermediate
801: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
802: @*/
803: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
804: {
805: PetscInt v;
811: for (v = 0; v < label->numStrata; ++v) {
812: PetscInt min, max;
814: if (label->stratumValues[v] != value) continue;
815: DMLabelMakeValid_Private(label, v);
816: if (label->stratumSizes[v] <= 0) break;
817: ISGetMinMax(label->points[v], &min, &max);
818: if (start) *start = min;
819: if (end) *end = max+1;
820: break;
821: }
822: return(0);
823: }
825: /*@
826: DMLabelGetStratumIS - Get an IS with the stratum points
828: Input Parameters:
829: + label - the DMLabel
830: - value - the stratum value
832: Output Paramater:
833: . points - The stratum points
835: Level: intermediate
837: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
838: @*/
839: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
840: {
841: PetscInt v;
846: *points = NULL;
847: for (v = 0; v < label->numStrata; ++v) {
848: if (label->stratumValues[v] == value) {
849: DMLabelMakeValid_Private(label, v);
850: if (label->validIS[v]) {
851: PetscObjectReference((PetscObject) label->points[v]);
852: *points = label->points[v];
853: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
854: break;
855: }
856: }
857: return(0);
858: }
860: /*@
861: DMLabelSetStratumIS - Set the stratum points using an IS
863: Input Parameters:
864: + label - the DMLabel
865: . value - the stratum value
866: - points - The stratum points
868: Level: intermediate
870: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
871: @*/
872: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
873: {
874: PetscInt v, numStrata;
878: numStrata = label->numStrata;
879: for (v = 0; v < numStrata; v++) {
880: if (label->stratumValues[v] == value) break;
881: }
882: if (v >= numStrata) {DMLabelAddStratum(label,value);}
883: if (is == label->points[v]) return(0);
884: DMLabelClearStratum(label,value);
885: ISGetLocalSize(is,&(label->stratumSizes[v]));
886: label->stratumValues[v] = value;
887: label->validIS[v] = PETSC_TRUE;
888: PetscObjectReference((PetscObject)is);
889: ISDestroy(&(label->points[v]));
890: if (label->bt) {
891: const PetscInt *points;
892: PetscInt p;
894: ISGetIndices(is,&points);
895: for (p = 0; p < label->stratumSizes[v]; ++p) {
896: const PetscInt point = points[p];
898: 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);
899: PetscBTSet(label->bt, point - label->pStart);
900: }
901: }
902: label->points[v] = is;
903: return(0);
904: }
906: /*@
907: DMLabelClearStratum - Remove a stratum
909: Input Parameters:
910: + label - the DMLabel
911: - value - the stratum value
913: Level: intermediate
915: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
916: @*/
917: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
918: {
919: PetscInt v;
923: for (v = 0; v < label->numStrata; ++v) {
924: if (label->stratumValues[v] == value) break;
925: }
926: if (v >= label->numStrata) return(0);
927: if (label->validIS[v]) {
928: if (label->bt) {
929: PetscInt i;
930: const PetscInt *points;
932: ISGetIndices(label->points[v], &points);
933: for (i = 0; i < label->stratumSizes[v]; ++i) {
934: const PetscInt point = points[i];
936: 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);
937: PetscBTClear(label->bt, point - label->pStart);
938: }
939: ISRestoreIndices(label->points[v], &points);
940: }
941: ISDestroy(&(label->points[v]));
942: label->stratumSizes[v] = 0;
943: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_OWN_POINTER,&(label->points[v]));
944: PetscObjectSetName((PetscObject) (label->points[v]), "indices");
945: } else {
946: PetscHashIClear(label->ht[v]);
947: }
948: return(0);
949: }
951: /*@
952: DMLabelFilter - Remove all points outside of [start, end)
954: Input Parameters:
955: + label - the DMLabel
956: . start - the first point
957: - end - the last point
959: Level: intermediate
961: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
962: @*/
963: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
964: {
965: PetscInt v;
969: DMLabelMakeAllValid_Private(label);
970: label->pStart = start;
971: label->pEnd = end;
972: if (label->bt) {PetscBTDestroy(&label->bt);}
973: /* Could squish offsets, but would only make sense if I reallocate the storage */
974: for (v = 0; v < label->numStrata; ++v) {
975: PetscInt off, q;
976: const PetscInt *points;
977: PetscInt *pointsNew = NULL;
979: ISGetIndices(label->points[v],&points);
980: for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
981: const PetscInt point = points[q];
983: if ((point < start) || (point >= end)) {
984: if (!pointsNew) {
985: PetscMalloc1(label->stratumSizes[v],&pointsNew);
986: PetscMemcpy(pointsNew,points,(size_t) off * sizeof(PetscInt));
987: }
988: continue;
989: }
990: if (pointsNew) {
991: pointsNew[off++] = point;
992: }
993: }
994: ISRestoreIndices(label->points[v],&points);
995: if (pointsNew) {
996: ISDestroy(&(label->points[v]));
997: ISCreateGeneral(PETSC_COMM_SELF,off,pointsNew,PETSC_OWN_POINTER,&(label->points[v]));
998: PetscObjectSetName((PetscObject) (label->points[v]), "indices");
999: }
1000: label->stratumSizes[v] = off;
1001: }
1002: DMLabelCreateIndex(label, start, end);
1003: return(0);
1004: }
1006: /*@
1007: DMLabelPermute - Create a new label with permuted points
1009: Input Parameters:
1010: + label - the DMLabel
1011: - permutation - the point permutation
1013: Output Parameter:
1014: . labelnew - the new label containing the permuted points
1016: Level: intermediate
1018: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1019: @*/
1020: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1021: {
1022: const PetscInt *perm;
1023: PetscInt numValues, numPoints, v, q;
1024: PetscErrorCode ierr;
1027: DMLabelMakeAllValid_Private(label);
1028: DMLabelDuplicate(label, labelNew);
1029: DMLabelGetNumValues(*labelNew, &numValues);
1030: ISGetLocalSize(permutation, &numPoints);
1031: ISGetIndices(permutation, &perm);
1032: for (v = 0; v < numValues; ++v) {
1033: const PetscInt size = (*labelNew)->stratumSizes[v];
1034: const PetscInt *points;
1035: PetscInt *pointsNew;
1037: ISGetIndices((*labelNew)->points[v],&points);
1038: PetscMalloc1(size,&pointsNew);
1039: for (q = 0; q < size; ++q) {
1040: const PetscInt point = points[q];
1042: 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);
1043: pointsNew[q] = perm[point];
1044: }
1045: ISRestoreIndices((*labelNew)->points[v],&points);
1046: PetscSortInt(size, pointsNew);
1047: ISDestroy(&((*labelNew)->points[v]));
1048: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1049: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1050: }
1051: ISRestoreIndices(permutation, &perm);
1052: if (label->bt) {
1053: PetscBTDestroy(&label->bt);
1054: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1055: }
1056: return(0);
1057: }
1059: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1060: {
1061: MPI_Comm comm;
1062: PetscInt s, l, nroots, nleaves, dof, offset, size;
1063: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1064: PetscSection rootSection;
1065: PetscSF labelSF;
1069: if (label) {DMLabelMakeAllValid_Private(label);}
1070: PetscObjectGetComm((PetscObject)sf, &comm);
1071: /* Build a section of stratum values per point, generate the according SF
1072: and distribute point-wise stratum values to leaves. */
1073: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1074: PetscSectionCreate(comm, &rootSection);
1075: PetscSectionSetChart(rootSection, 0, nroots);
1076: if (label) {
1077: for (s = 0; s < label->numStrata; ++s) {
1078: const PetscInt *points;
1080: ISGetIndices(label->points[s], &points);
1081: for (l = 0; l < label->stratumSizes[s]; l++) {
1082: PetscSectionGetDof(rootSection, points[l], &dof);
1083: PetscSectionSetDof(rootSection, points[l], dof+1);
1084: }
1085: ISRestoreIndices(label->points[s], &points);
1086: }
1087: }
1088: PetscSectionSetUp(rootSection);
1089: /* Create a point-wise array of stratum values */
1090: PetscSectionGetStorageSize(rootSection, &size);
1091: PetscMalloc1(size, &rootStrata);
1092: PetscCalloc1(nroots, &rootIdx);
1093: if (label) {
1094: for (s = 0; s < label->numStrata; ++s) {
1095: const PetscInt *points;
1097: ISGetIndices(label->points[s], &points);
1098: for (l = 0; l < label->stratumSizes[s]; l++) {
1099: const PetscInt p = points[l];
1100: PetscSectionGetOffset(rootSection, p, &offset);
1101: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1102: }
1103: ISRestoreIndices(label->points[s], &points);
1104: }
1105: }
1106: /* Build SF that maps label points to remote processes */
1107: PetscSectionCreate(comm, leafSection);
1108: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1109: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1110: PetscFree(remoteOffsets);
1111: /* Send the strata for each point over the derived SF */
1112: PetscSectionGetStorageSize(*leafSection, &size);
1113: PetscMalloc1(size, leafStrata);
1114: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1115: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1116: /* Clean up */
1117: PetscFree(rootStrata);
1118: PetscFree(rootIdx);
1119: PetscSectionDestroy(&rootSection);
1120: PetscSFDestroy(&labelSF);
1121: return(0);
1122: }
1124: /*@
1125: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1127: Input Parameters:
1128: + label - the DMLabel
1129: - sf - the map from old to new distribution
1131: Output Parameter:
1132: . labelnew - the new redistributed label
1134: Level: intermediate
1136: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1137: @*/
1138: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1139: {
1140: MPI_Comm comm;
1141: PetscSection leafSection;
1142: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1143: PetscInt *leafStrata, *strataIdx;
1144: PetscInt **points;
1145: char *name;
1146: PetscInt nameSize;
1147: PetscHashI stratumHash;
1148: PETSC_UNUSED PetscHashIIter ret, iter;
1149: size_t len = 0;
1150: PetscMPIInt rank;
1154: if (label) {DMLabelMakeAllValid_Private(label);}
1155: PetscObjectGetComm((PetscObject)sf, &comm);
1156: MPI_Comm_rank(comm, &rank);
1157: /* Bcast name */
1158: if (!rank) {PetscStrlen(label->name, &len);}
1159: nameSize = len;
1160: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1161: PetscMalloc1(nameSize+1, &name);
1162: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1163: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1164: DMLabelCreate(name, labelNew);
1165: PetscFree(name);
1166: /* Bcast defaultValue */
1167: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1168: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1169: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1170: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1171: /* Determine received stratum values and initialise new label*/
1172: PetscHashICreate(stratumHash);
1173: PetscSectionGetStorageSize(leafSection, &size);
1174: for (p = 0; p < size; ++p) PetscHashIPut(stratumHash, leafStrata[p], ret, iter);
1175: PetscHashISize(stratumHash, (*labelNew)->numStrata);
1176: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1177: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1178: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1179: /* Turn leafStrata into indices rather than stratum values */
1180: offset = 0;
1181: PetscHashIGetKeys(stratumHash, &offset, (*labelNew)->stratumValues);
1182: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1183: for (p = 0; p < size; ++p) {
1184: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1185: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1186: }
1187: }
1188: /* Rebuild the point strata on the receiver */
1189: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1190: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1191: for (p=pStart; p<pEnd; p++) {
1192: PetscSectionGetDof(leafSection, p, &dof);
1193: PetscSectionGetOffset(leafSection, p, &offset);
1194: for (s=0; s<dof; s++) {
1195: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1196: }
1197: }
1198: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1199: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1200: PetscMalloc1((*labelNew)->numStrata,&points);
1201: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1202: PetscHashICreate((*labelNew)->ht[s]);
1203: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1204: }
1205: /* Insert points into new strata */
1206: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1207: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1208: for (p=pStart; p<pEnd; p++) {
1209: PetscSectionGetDof(leafSection, p, &dof);
1210: PetscSectionGetOffset(leafSection, p, &offset);
1211: for (s=0; s<dof; s++) {
1212: stratum = leafStrata[offset+s];
1213: points[stratum][strataIdx[stratum]++] = p;
1214: }
1215: }
1216: for (s = 0; s < (*labelNew)->numStrata; s++) {
1217: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1218: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1219: }
1220: PetscFree(points);
1221: PetscHashIDestroy(stratumHash);
1222: PetscFree(leafStrata);
1223: PetscFree(strataIdx);
1224: PetscSectionDestroy(&leafSection);
1225: return(0);
1226: }
1228: /*@
1229: DMLabelGather - Gather all label values from leafs into roots
1231: Input Parameters:
1232: + label - the DMLabel
1233: - sf - the Star Forest point communication map
1235: Output Parameters:
1236: . labelNew - the new DMLabel with localised leaf values
1238: Level: developer
1240: Note: This is the inverse operation to DMLabelDistribute.
1242: .seealso: DMLabelDistribute()
1243: @*/
1244: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1245: {
1246: MPI_Comm comm;
1247: PetscSection rootSection;
1248: PetscSF sfLabel;
1249: PetscSFNode *rootPoints, *leafPoints;
1250: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1251: const PetscInt *rootDegree, *ilocal;
1252: PetscInt *rootStrata;
1253: char *name;
1254: PetscInt nameSize;
1255: size_t len = 0;
1256: PetscMPIInt rank, size;
1260: PetscObjectGetComm((PetscObject)sf, &comm);
1261: MPI_Comm_rank(comm, &rank);
1262: MPI_Comm_size(comm, &size);
1263: /* Bcast name */
1264: if (!rank) {PetscStrlen(label->name, &len);}
1265: nameSize = len;
1266: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1267: PetscMalloc1(nameSize+1, &name);
1268: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1269: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1270: DMLabelCreate(name, labelNew);
1271: PetscFree(name);
1272: /* Gather rank/index pairs of leaves into local roots to build
1273: an inverse, multi-rooted SF. Note that this ignores local leaf
1274: indexing due to the use of the multiSF in PetscSFGather. */
1275: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1276: PetscMalloc1(nroots, &leafPoints);
1277: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1278: for (p = 0; p < nleaves; p++) {
1279: leafPoints[ilocal[p]].index = ilocal[p];
1280: leafPoints[ilocal[p]].rank = rank;
1281: }
1282: PetscSFComputeDegreeBegin(sf, &rootDegree);
1283: PetscSFComputeDegreeEnd(sf, &rootDegree);
1284: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1285: PetscMalloc1(nmultiroots, &rootPoints);
1286: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1287: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1288: PetscSFCreate(comm,& sfLabel);
1289: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1290: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1291: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1292: /* Rebuild the point strata on the receiver */
1293: for (p = 0, idx = 0; p < nroots; p++) {
1294: for (d = 0; d < rootDegree[p]; d++) {
1295: PetscSectionGetDof(rootSection, idx+d, &dof);
1296: PetscSectionGetOffset(rootSection, idx+d, &offset);
1297: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1298: }
1299: idx += rootDegree[p];
1300: }
1301: PetscFree(leafPoints);
1302: PetscFree(rootStrata);
1303: PetscSectionDestroy(&rootSection);
1304: PetscSFDestroy(&sfLabel);
1305: return(0);
1306: }
1308: /*@
1309: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1311: Input Parameter:
1312: . label - the DMLabel
1314: Output Parameters:
1315: + section - the section giving offsets for each stratum
1316: - is - An IS containing all the label points
1318: Level: developer
1320: .seealso: DMLabelDistribute()
1321: @*/
1322: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1323: {
1324: IS vIS;
1325: const PetscInt *values;
1326: PetscInt *points;
1327: PetscInt nV, vS = 0, vE = 0, v, N;
1328: PetscErrorCode ierr;
1331: DMLabelGetNumValues(label, &nV);
1332: DMLabelGetValueIS(label, &vIS);
1333: ISGetIndices(vIS, &values);
1334: if (nV) {vS = values[0]; vE = values[0]+1;}
1335: for (v = 1; v < nV; ++v) {
1336: vS = PetscMin(vS, values[v]);
1337: vE = PetscMax(vE, values[v]+1);
1338: }
1339: PetscSectionCreate(PETSC_COMM_SELF, section);
1340: PetscSectionSetChart(*section, vS, vE);
1341: for (v = 0; v < nV; ++v) {
1342: PetscInt n;
1344: DMLabelGetStratumSize(label, values[v], &n);
1345: PetscSectionSetDof(*section, values[v], n);
1346: }
1347: PetscSectionSetUp(*section);
1348: PetscSectionGetStorageSize(*section, &N);
1349: PetscMalloc1(N, &points);
1350: for (v = 0; v < nV; ++v) {
1351: IS is;
1352: const PetscInt *spoints;
1353: PetscInt dof, off, p;
1355: PetscSectionGetDof(*section, values[v], &dof);
1356: PetscSectionGetOffset(*section, values[v], &off);
1357: DMLabelGetStratumIS(label, values[v], &is);
1358: ISGetIndices(is, &spoints);
1359: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1360: ISRestoreIndices(is, &spoints);
1361: ISDestroy(&is);
1362: }
1363: ISRestoreIndices(vIS, &values);
1364: ISDestroy(&vIS);
1365: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1366: return(0);
1367: }
1369: /*@
1370: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1371: the local section and an SF describing the section point overlap.
1373: Input Parameters:
1374: + s - The PetscSection for the local field layout
1375: . sf - The SF describing parallel layout of the section points
1376: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1377: . label - The label specifying the points
1378: - labelValue - The label stratum specifying the points
1380: Output Parameter:
1381: . gsection - The PetscSection for the global field layout
1383: Note: This gives negative sizes and offsets to points not owned by this process
1385: Level: developer
1387: .seealso: PetscSectionCreate()
1388: @*/
1389: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1390: {
1391: PetscInt *neg = NULL, *tmpOff = NULL;
1392: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1396: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1397: PetscSectionGetChart(s, &pStart, &pEnd);
1398: PetscSectionSetChart(*gsection, pStart, pEnd);
1399: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1400: if (nroots >= 0) {
1401: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1402: PetscCalloc1(nroots, &neg);
1403: if (nroots > pEnd-pStart) {
1404: PetscCalloc1(nroots, &tmpOff);
1405: } else {
1406: tmpOff = &(*gsection)->atlasDof[-pStart];
1407: }
1408: }
1409: /* Mark ghost points with negative dof */
1410: for (p = pStart; p < pEnd; ++p) {
1411: PetscInt value;
1413: DMLabelGetValue(label, p, &value);
1414: if (value != labelValue) continue;
1415: PetscSectionGetDof(s, p, &dof);
1416: PetscSectionSetDof(*gsection, p, dof);
1417: PetscSectionGetConstraintDof(s, p, &cdof);
1418: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1419: if (neg) neg[p] = -(dof+1);
1420: }
1421: PetscSectionSetUpBC(*gsection);
1422: if (nroots >= 0) {
1423: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1424: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1425: if (nroots > pEnd-pStart) {
1426: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1427: }
1428: }
1429: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1430: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1431: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1432: (*gsection)->atlasOff[p] = off;
1433: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1434: }
1435: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1436: globalOff -= off;
1437: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1438: (*gsection)->atlasOff[p] += globalOff;
1439: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1440: }
1441: /* Put in negative offsets for ghost points */
1442: if (nroots >= 0) {
1443: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1444: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1445: if (nroots > pEnd-pStart) {
1446: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1447: }
1448: }
1449: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1450: PetscFree(neg);
1451: return(0);
1452: }
1454: typedef struct _n_PetscSectionSym_Label
1455: {
1456: DMLabel label;
1457: PetscCopyMode *modes;
1458: PetscInt *sizes;
1459: const PetscInt ***perms;
1460: const PetscScalar ***rots;
1461: PetscInt (*minMaxOrients)[2];
1462: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1463: } PetscSectionSym_Label;
1465: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1466: {
1467: PetscInt i, j;
1468: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1469: PetscErrorCode ierr;
1472: for (i = 0; i <= sl->numStrata; i++) {
1473: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1474: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1475: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1476: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1477: }
1478: if (sl->perms[i]) {
1479: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1481: PetscFree(perms);
1482: }
1483: if (sl->rots[i]) {
1484: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1486: PetscFree(rots);
1487: }
1488: }
1489: }
1490: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1491: DMLabelDestroy(&sl->label);
1492: sl->numStrata = 0;
1493: return(0);
1494: }
1496: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1497: {
1501: PetscSectionSymLabelReset(sym);
1502: PetscFree(sym->data);
1503: return(0);
1504: }
1506: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1507: {
1508: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1509: PetscBool isAscii;
1510: DMLabel label = sl->label;
1511: PetscErrorCode ierr;
1514: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1515: if (isAscii) {
1516: PetscInt i, j, k;
1517: PetscViewerFormat format;
1519: PetscViewerGetFormat(viewer,&format);
1520: if (label) {
1521: PetscViewerGetFormat(viewer,&format);
1522: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1523: PetscViewerASCIIPushTab(viewer);
1524: DMLabelView(label, viewer);
1525: PetscViewerASCIIPopTab(viewer);
1526: } else {
1527: PetscViewerASCIIPrintf(viewer," Label '%s'\n",sl->label->name);
1528: }
1529: } else {
1530: PetscViewerASCIIPrintf(viewer, "No label given\n");
1531: }
1532: PetscViewerASCIIPushTab(viewer);
1533: for (i = 0; i <= sl->numStrata; i++) {
1534: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1536: if (!(sl->perms[i] || sl->rots[i])) {
1537: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1538: } else {
1539: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1540: PetscViewerASCIIPushTab(viewer);
1541: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1542: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1543: PetscViewerASCIIPushTab(viewer);
1544: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1545: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1546: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1547: } else {
1548: PetscInt tab;
1550: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1551: PetscViewerASCIIPushTab(viewer);
1552: PetscViewerASCIIGetTab(viewer,&tab);
1553: if (sl->perms[i] && sl->perms[i][j]) {
1554: PetscViewerASCIIPrintf(viewer,"Permutation:");
1555: PetscViewerASCIISetTab(viewer,0);
1556: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1557: PetscViewerASCIIPrintf(viewer,"\n");
1558: PetscViewerASCIISetTab(viewer,tab);
1559: }
1560: if (sl->rots[i] && sl->rots[i][j]) {
1561: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1562: PetscViewerASCIISetTab(viewer,0);
1563: #if defined(PETSC_USE_COMPLEX)
1564: 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]));}
1565: #else
1566: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1567: #endif
1568: PetscViewerASCIIPrintf(viewer,"\n");
1569: PetscViewerASCIISetTab(viewer,tab);
1570: }
1571: PetscViewerASCIIPopTab(viewer);
1572: }
1573: }
1574: PetscViewerASCIIPopTab(viewer);
1575: }
1576: PetscViewerASCIIPopTab(viewer);
1577: }
1578: }
1579: PetscViewerASCIIPopTab(viewer);
1580: }
1581: return(0);
1582: }
1584: /*@
1585: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1587: Logically collective on sym
1589: Input parameters:
1590: + sym - the section symmetries
1591: - label - the DMLabel describing the types of points
1593: Level: developer:
1595: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1596: @*/
1597: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1598: {
1599: PetscSectionSym_Label *sl;
1604: sl = (PetscSectionSym_Label *) sym->data;
1605: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1606: if (label) {
1607: label->refct++;
1608: sl->label = label;
1609: DMLabelGetNumValues(label,&sl->numStrata);
1610: 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);
1611: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1612: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1613: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1614: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1615: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1616: }
1617: return(0);
1618: }
1620: /*@C
1621: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1623: Logically collective on PetscSectionSym
1625: InputParameters:
1626: + sys - the section symmetries
1627: . stratum - the stratum value in the label that we are assigning symmetries for
1628: . size - the number of dofs for points in the stratum of the label
1629: . minOrient - the smallest orientation for a point in this stratum
1630: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1631: . mode - how sym should copy the perms and rots arrays
1632: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
1633: + 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
1635: Level: developer
1637: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1638: @*/
1639: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1640: {
1641: PetscInt i, j, k;
1642: PetscSectionSym_Label *sl;
1647: sl = (PetscSectionSym_Label *) sym->data;
1648: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1649: for (i = 0; i <= sl->numStrata; i++) {
1650: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
1652: if (stratum == value) break;
1653: }
1654: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,sl->label->name);
1655: sl->sizes[i] = size;
1656: sl->modes[i] = mode;
1657: sl->minMaxOrients[i][0] = minOrient;
1658: sl->minMaxOrients[i][1] = maxOrient;
1659: if (mode == PETSC_COPY_VALUES) {
1660: if (perms) {
1661: PetscInt **ownPerms;
1663: PetscCalloc1(maxOrient - minOrient,&ownPerms);
1664: for (j = 0; j < maxOrient-minOrient; j++) {
1665: if (perms[j]) {
1666: PetscMalloc1(size,&ownPerms[j]);
1667: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
1668: }
1669: }
1670: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
1671: }
1672: if (rots) {
1673: PetscScalar **ownRots;
1675: PetscCalloc1(maxOrient - minOrient,&ownRots);
1676: for (j = 0; j < maxOrient-minOrient; j++) {
1677: if (rots[j]) {
1678: PetscMalloc1(size,&ownRots[j]);
1679: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
1680: }
1681: }
1682: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
1683: }
1684: } else {
1685: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
1686: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
1687: }
1688: return(0);
1689: }
1691: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
1692: {
1693: PetscInt i, j, numStrata;
1694: PetscSectionSym_Label *sl;
1695: DMLabel label;
1696: PetscErrorCode ierr;
1699: sl = (PetscSectionSym_Label *) sym->data;
1700: numStrata = sl->numStrata;
1701: label = sl->label;
1702: for (i = 0; i < numPoints; i++) {
1703: PetscInt point = points[2*i];
1704: PetscInt ornt = points[2*i+1];
1706: for (j = 0; j < numStrata; j++) {
1707: if (label->validIS[j]) {
1708: PetscInt k;
1710: ISLocate(label->points[j],point,&k);
1711: if (k >= 0) break;
1712: } else {
1713: PetscBool has;
1715: PetscHashIHasKey(label->ht[j], point, has);
1716: if (has) break;
1717: }
1718: }
1719: 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);
1720: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
1721: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
1722: }
1723: return(0);
1724: }
1726: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
1727: {
1728: PetscSectionSym_Label *sl;
1729: PetscErrorCode ierr;
1732: PetscNewLog(sym,&sl);
1733: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
1734: sym->ops->view = PetscSectionSymView_Label;
1735: sym->ops->destroy = PetscSectionSymDestroy_Label;
1736: sym->data = (void *) sl;
1737: return(0);
1738: }
1740: /*@
1741: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
1743: Collective
1745: Input Parameters:
1746: + comm - the MPI communicator for the new symmetry
1747: - label - the label defining the strata
1749: Output Parameters:
1750: . sym - the section symmetries
1752: Level: developer
1754: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
1755: @*/
1756: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
1757: {
1758: PetscErrorCode ierr;
1761: DMInitializePackage();
1762: PetscSectionSymCreate(comm,sym);
1763: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
1764: PetscSectionSymLabelSetLabel(*sym,label);
1765: return(0);
1766: }