Actual source code: dmlabel.c
petsc-3.9.4 2018-09-11
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 = 0;
59: PetscInt *pointArray;
62: if (label->validIS[v]) return 0;
64: 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);
65: PetscHashISize(label->ht[v], label->stratumSizes[v]);
67: PetscMalloc1(label->stratumSizes[v], &pointArray);
68: PetscHashIGetKeys(label->ht[v], &off, pointArray);
69: 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]);
70: PetscHashIClear(label->ht[v]);
71: PetscSortInt(label->stratumSizes[v], pointArray);
72: if (label->bt) {
73: PetscInt p;
75: for (p = 0; p < label->stratumSizes[v]; ++p) {
76: const PetscInt point = pointArray[p];
78: 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);
79: PetscBTSet(label->bt, point - label->pStart);
80: }
81: }
82: ISCreateGeneral(PETSC_COMM_SELF,label->stratumSizes[v],pointArray,PETSC_OWN_POINTER,&(label->points[v]));
83: PetscObjectSetName((PetscObject) (label->points[v]), "indices");
84: label->validIS[v] = PETSC_TRUE;
85: ++label->state;
86: return(0);
87: }
89: /*
90: DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format
92: Input parameter:
93: . label - The DMLabel
95: Output parameter:
96: . label - The DMLabel with all strata in sorted list format
98: Level: developer
100: .seealso: DMLabelCreate()
101: */
102: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
103: {
104: PetscInt v;
108: for (v = 0; v < label->numStrata; v++){
109: DMLabelMakeValid_Private(label, v);
110: }
111: return(0);
112: }
114: /*
115: DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format
117: Input parameter:
118: + label - The DMLabel
119: - v - The stratum value
121: Output parameter:
122: . label - The DMLabel with stratum in hash format
124: Level: developer
126: .seealso: DMLabelCreate()
127: */
128: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
129: {
130: PETSC_UNUSED PetscHashIIter ret, iter;
131: PetscInt p;
132: const PetscInt *points;
135: if (!label->validIS[v]) return 0;
137: 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);
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: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
157: {
158: PetscInt v;
160: for (*index = -1, v = 0; v < label->numStrata; ++v)
161: if (label->stratumValues[v] == value) { *index = v; break; }
162: return(0);
163: }
165: static PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
166: {
167: PetscInt v, *tmpV, *tmpS;
168: IS *tmpP;
169: PetscHashI *tmpH;
170: PetscBool *tmpB;
174: PetscMalloc1((label->numStrata+1), &tmpV);
175: PetscMalloc1((label->numStrata+1), &tmpS);
176: PetscMalloc1((label->numStrata+1), &tmpH);
177: PetscMalloc1((label->numStrata+1), &tmpP);
178: PetscMalloc1((label->numStrata+1), &tmpB);
179: for (v = 0; v < label->numStrata; ++v) {
180: tmpV[v] = label->stratumValues[v];
181: tmpS[v] = label->stratumSizes[v];
182: tmpH[v] = label->ht[v];
183: tmpP[v] = label->points[v];
184: tmpB[v] = label->validIS[v];
185: }
186: tmpV[v] = value;
187: tmpS[v] = 0;
188: PetscHashICreate(tmpH[v]);
189: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_OWN_POINTER,&tmpP[v]);
190: tmpB[v] = PETSC_TRUE;
192: ++label->numStrata;
193: PetscFree(label->stratumValues);
194: PetscFree(label->stratumSizes);
195: PetscFree(label->ht);
196: PetscFree(label->points);
197: PetscFree(label->validIS);
198: label->stratumValues = tmpV;
199: label->stratumSizes = tmpS;
200: label->ht = tmpH;
201: label->points = tmpP;
202: label->validIS = tmpB;
204: *index = v;
205: return(0);
206: }
208: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
209: {
210: PetscInt v;
214: DMLabelLookupStratum(label, value, &v);
215: if (v < 0) {DMLabelNewStratum(label, value, &v);}
216: return(0);
217: }
219: /*@C
220: DMLabelGetName - Return the name of a DMLabel object
222: Input parameter:
223: . label - The DMLabel
225: Output parameter:
226: . name - The label name
228: Level: beginner
230: .seealso: DMLabelCreate()
231: @*/
232: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
233: {
236: *name = label->name;
237: return(0);
238: }
240: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
241: {
242: PetscInt v;
243: PetscMPIInt rank;
247: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
248: PetscViewerASCIIPushSynchronized(viewer);
249: if (label) {
250: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
251: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
252: for (v = 0; v < label->numStrata; ++v) {
253: const PetscInt value = label->stratumValues[v];
254: const PetscInt *points;
255: PetscInt p;
257: ISGetIndices(label->points[v],&points);
258: for (p = 0; p < label->stratumSizes[v]; ++p) {
259: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
260: }
261: ISRestoreIndices(label->points[v],&points);
262: }
263: }
264: PetscViewerFlush(viewer);
265: PetscViewerASCIIPopSynchronized(viewer);
266: return(0);
267: }
269: /*@C
270: DMLabelView - View the label
272: Input Parameters:
273: + label - The DMLabel
274: - viewer - The PetscViewer
276: Level: intermediate
278: .seealso: DMLabelCreate(), DMLabelDestroy()
279: @*/
280: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
281: {
282: PetscBool iascii;
287: if (label) {DMLabelMakeAllValid_Private(label);}
288: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
289: if (iascii) {
290: DMLabelView_Ascii(label, viewer);
291: }
292: return(0);
293: }
295: /*@
296: DMLabelDestroy - Destroys a DMLabel
298: Input Parameter:
299: . label - The DMLabel
301: Level: beginner
303: .seealso: DMLabelCreate()
304: @*/
305: PetscErrorCode DMLabelDestroy(DMLabel *label)
306: {
307: PetscInt v;
311: if (!(*label)) return(0);
312: if (--(*label)->refct > 0) return(0);
313: PetscFree((*label)->name);
314: PetscFree((*label)->stratumValues);
315: PetscFree((*label)->stratumSizes);
316: for (v = 0; v < (*label)->numStrata; ++v) {
317: PetscHashIDestroy((*label)->ht[v]);
318: ISDestroy(&((*label)->points[v]));
319: }
320: PetscFree((*label)->ht);
321: PetscFree((*label)->points);
322: PetscFree((*label)->validIS);
323: PetscBTDestroy(&(*label)->bt);
324: PetscFree(*label);
325: return(0);
326: }
328: /*@
329: DMLabelDuplicate - Duplicates a DMLabel
331: Input Parameter:
332: . label - The DMLabel
334: Output Parameter:
335: . labelnew - location to put new vector
337: Level: intermediate
339: .seealso: DMLabelCreate(), DMLabelDestroy()
340: @*/
341: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
342: {
343: PetscInt v;
347: DMLabelMakeAllValid_Private(label);
348: PetscNew(labelnew);
349: PetscStrallocpy(label->name, &(*labelnew)->name);
351: (*labelnew)->refct = 1;
352: (*labelnew)->numStrata = label->numStrata;
353: (*labelnew)->defaultValue = label->defaultValue;
354: if (label->numStrata) {
355: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
356: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
357: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
358: PetscMalloc1(label->numStrata, &(*labelnew)->points);
359: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
360: /* Could eliminate unused space here */
361: for (v = 0; v < label->numStrata; ++v) {
362: PetscHashICreate((*labelnew)->ht[v]);
363: (*labelnew)->validIS[v] = PETSC_TRUE;
364: (*labelnew)->stratumValues[v] = label->stratumValues[v];
365: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
366: PetscObjectReference((PetscObject) (label->points[v]));
367: (*labelnew)->points[v] = label->points[v];
368: }
369: }
370: (*labelnew)->pStart = -1;
371: (*labelnew)->pEnd = -1;
372: (*labelnew)->bt = NULL;
373: return(0);
374: }
376: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
377: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
378: {
379: PetscInt v;
383: DMLabelDestroyIndex(label);
384: DMLabelMakeAllValid_Private(label);
385: label->pStart = pStart;
386: label->pEnd = pEnd;
387: PetscBTCreate(pEnd - pStart, &label->bt);
388: for (v = 0; v < label->numStrata; ++v) {
389: const PetscInt *points;
390: PetscInt i;
392: ISGetIndices(label->points[v], &points);
393: for (i = 0; i < label->stratumSizes[v]; ++i) {
394: const PetscInt point = points[i];
396: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
397: PetscBTSet(label->bt, point - pStart);
398: }
399: ISRestoreIndices(label->points[v], &points);
400: }
401: return(0);
402: }
404: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
405: {
409: label->pStart = -1;
410: label->pEnd = -1;
411: PetscBTDestroy(&label->bt);
412: return(0);
413: }
415: /*@
416: DMLabelHasValue - Determine whether a label assigns the value to any point
418: Input Parameters:
419: + label - the DMLabel
420: - value - the value
422: Output Parameter:
423: . contains - Flag indicating whether the label maps this value to any point
425: Level: developer
427: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
428: @*/
429: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
430: {
431: PetscInt v;
436: DMLabelLookupStratum(label, value, &v);
437: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
438: return(0);
439: }
441: /*@
442: DMLabelHasPoint - Determine whether a label assigns a value to a point
444: Input Parameters:
445: + label - the DMLabel
446: - point - the point
448: Output Parameter:
449: . contains - Flag indicating whether the label maps this point to a value
451: Note: The user must call DMLabelCreateIndex() before this function.
453: Level: developer
455: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
456: @*/
457: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
458: {
463: DMLabelMakeAllValid_Private(label);
464: #if defined(PETSC_USE_DEBUG)
465: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
466: 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);
467: #endif
468: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
469: return(0);
470: }
472: /*@
473: DMLabelStratumHasPoint - Return true if the stratum contains a point
475: Input Parameters:
476: + label - the DMLabel
477: . value - the stratum value
478: - point - the point
480: Output Parameter:
481: . contains - true if the stratum contains the point
483: Level: intermediate
485: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
486: @*/
487: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
488: {
489: PetscInt v;
494: *contains = PETSC_FALSE;
495: DMLabelLookupStratum(label, value, &v);
496: if (v < 0) return(0);
498: if (label->validIS[v]) {
499: PetscInt i;
501: ISLocate(label->points[v], point, &i);
502: if (i >= 0) *contains = PETSC_TRUE;
503: } else {
504: PetscBool has;
506: PetscHashIHasKey(label->ht[v], point, has);
507: if (has) *contains = PETSC_TRUE;
508: }
509: return(0);
510: }
512: /*@
513: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
514: When a label is created, it is initialized to -1.
516: Input parameter:
517: . label - a DMLabel object
519: Output parameter:
520: . defaultValue - the default value
522: Level: beginner
524: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
525: @*/
526: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
527: {
529: *defaultValue = label->defaultValue;
530: return(0);
531: }
533: /*@
534: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
535: When a label is created, it is initialized to -1.
537: Input parameter:
538: . label - a DMLabel object
540: Output parameter:
541: . defaultValue - the default value
543: Level: beginner
545: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
546: @*/
547: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
548: {
550: label->defaultValue = defaultValue;
551: return(0);
552: }
554: /*@
555: 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())
557: Input Parameters:
558: + label - the DMLabel
559: - point - the point
561: Output Parameter:
562: . value - The point value, or -1
564: Level: intermediate
566: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
567: @*/
568: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
569: {
570: PetscInt v;
575: *value = label->defaultValue;
576: for (v = 0; v < label->numStrata; ++v) {
577: if (label->validIS[v]) {
578: PetscInt i;
580: ISLocate(label->points[v], point, &i);
581: if (i >= 0) {
582: *value = label->stratumValues[v];
583: break;
584: }
585: } else {
586: PetscBool has;
588: PetscHashIHasKey(label->ht[v], point, has);
589: if (has) {
590: *value = label->stratumValues[v];
591: break;
592: }
593: }
594: }
595: return(0);
596: }
598: /*@
599: 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.
601: Input Parameters:
602: + label - the DMLabel
603: . point - the point
604: - value - The point value
606: Level: intermediate
608: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
609: @*/
610: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
611: {
612: PETSC_UNUSED PetscHashIIter iter, ret;
613: PetscInt v;
614: PetscErrorCode ierr;
617: /* Find label value, add new entry if needed */
618: if (value == label->defaultValue) return(0);
619: DMLabelLookupStratum(label, value, &v);
620: if (v < 0) {DMLabelNewStratum(label, value, &v);}
621: /* Set key */
622: DMLabelMakeInvalid_Private(label, v);
623: PetscHashIPut(label->ht[v], point, ret, iter);
624: return(0);
625: }
627: /*@
628: DMLabelClearValue - Clear the value a label assigns to a point
630: Input Parameters:
631: + label - the DMLabel
632: . point - the point
633: - value - The point value
635: Level: intermediate
637: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
638: @*/
639: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
640: {
641: PetscInt v;
645: /* Find label value */
646: DMLabelLookupStratum(label, value, &v);
647: if (v < 0) return(0);
649: if (label->bt) {
650: 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);
651: PetscBTClear(label->bt, point - label->pStart);
652: }
654: /* Delete key */
655: DMLabelMakeInvalid_Private(label, v);
656: PetscHashIDelKey(label->ht[v], point);
657: return(0);
658: }
660: /*@
661: DMLabelInsertIS - Set all points in the IS to a value
663: Input Parameters:
664: + label - the DMLabel
665: . is - the point IS
666: - value - The point value
668: Level: intermediate
670: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
671: @*/
672: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
673: {
674: PETSC_UNUSED PetscHashIIter iter, ret;
675: PetscInt v, n, p;
676: const PetscInt *points;
677: PetscErrorCode ierr;
681: /* Find label value, add new entry if needed */
682: if (value == label->defaultValue) return(0);
683: DMLabelLookupStratum(label, value, &v);
684: if (v < 0) {DMLabelNewStratum(label, value, &v);}
685: /* Set keys */
686: DMLabelMakeInvalid_Private(label, v);
687: ISGetLocalSize(is, &n);
688: ISGetIndices(is, &points);
689: for (p = 0; p < n; ++p) PetscHashIPut(label->ht[v], points[p], ret, iter);
690: ISRestoreIndices(is, &points);
691: return(0);
692: }
694: /*@
695: DMLabelGetNumValues - Get the number of values that the DMLabel takes
697: Input Parameter:
698: . label - the DMLabel
700: Output Paramater:
701: . numValues - the number of values
703: Level: intermediate
705: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
706: @*/
707: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
708: {
711: *numValues = label->numStrata;
712: return(0);
713: }
715: /*@
716: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
718: Input Parameter:
719: . label - the DMLabel
721: Output Paramater:
722: . is - the value IS
724: Level: intermediate
726: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
727: @*/
728: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
729: {
734: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
735: return(0);
736: }
738: /*@
739: DMLabelHasStratum - Determine whether points exist with the given value
741: Input Parameters:
742: + label - the DMLabel
743: - value - the stratum value
745: Output Paramater:
746: . exists - Flag saying whether points exist
748: Level: intermediate
750: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
751: @*/
752: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
753: {
754: PetscInt v;
759: DMLabelLookupStratum(label, value, &v);
760: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
761: return(0);
762: }
764: /*@
765: DMLabelGetStratumSize - Get the size of a stratum
767: Input Parameters:
768: + label - the DMLabel
769: - value - the stratum value
771: Output Paramater:
772: . size - The number of points in the stratum
774: Level: intermediate
776: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
777: @*/
778: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
779: {
780: PetscInt v;
785: *size = 0;
786: DMLabelLookupStratum(label, value, &v);
787: if (v < 0) return(0);
788: DMLabelMakeValid_Private(label, v);
789: *size = label->stratumSizes[v];
790: return(0);
791: }
793: /*@
794: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
796: Input Parameters:
797: + label - the DMLabel
798: - value - the stratum value
800: Output Paramaters:
801: + start - the smallest point in the stratum
802: - end - the largest point in the stratum
804: Level: intermediate
806: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
807: @*/
808: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
809: {
810: PetscInt v, min, max;
816: DMLabelLookupStratum(label, value, &v);
817: if (v < 0) return(0);
818: DMLabelMakeValid_Private(label, v);
819: if (label->stratumSizes[v] <= 0) return(0);
820: ISGetMinMax(label->points[v], &min, &max);
821: if (start) *start = min;
822: if (end) *end = max+1;
823: return(0);
824: }
826: /*@
827: DMLabelGetStratumIS - Get an IS with the stratum points
829: Input Parameters:
830: + label - the DMLabel
831: - value - the stratum value
833: Output Paramater:
834: . points - The stratum points
836: Level: intermediate
838: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
839: @*/
840: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
841: {
842: PetscInt v;
847: *points = NULL;
848: DMLabelLookupStratum(label, value, &v);
849: if (v < 0) return(0);
850: DMLabelMakeValid_Private(label, v);
851: PetscObjectReference((PetscObject) label->points[v]);
852: *points = label->points[v];
853: return(0);
854: }
856: /*@
857: DMLabelSetStratumIS - Set the stratum points using an IS
859: Input Parameters:
860: + label - the DMLabel
861: . value - the stratum value
862: - points - The stratum points
864: Level: intermediate
866: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
867: @*/
868: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
869: {
870: PetscInt v;
874: DMLabelLookupStratum(label, value, &v);
875: if (v < 0) {DMLabelNewStratum(label, value, &v);}
876: if (is == label->points[v]) return(0);
877: DMLabelClearStratum(label, value);
878: ISGetLocalSize(is, &(label->stratumSizes[v]));
879: PetscObjectReference((PetscObject)is);
880: ISDestroy(&(label->points[v]));
881: label->points[v] = is;
882: label->validIS[v] = PETSC_TRUE;
883: if (label->bt) {
884: const PetscInt *points;
885: PetscInt p;
887: ISGetIndices(is,&points);
888: for (p = 0; p < label->stratumSizes[v]; ++p) {
889: const PetscInt point = points[p];
891: 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);
892: PetscBTSet(label->bt, point - label->pStart);
893: }
894: }
895: return(0);
896: }
898: /*@
899: DMLabelClearStratum - Remove a stratum
901: Input Parameters:
902: + label - the DMLabel
903: - value - the stratum value
905: Level: intermediate
907: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
908: @*/
909: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
910: {
911: PetscInt v;
915: DMLabelLookupStratum(label, value, &v);
916: if (v < 0) return(0);
917: if (label->validIS[v]) {
918: if (label->bt) {
919: PetscInt i;
920: const PetscInt *points;
922: ISGetIndices(label->points[v], &points);
923: for (i = 0; i < label->stratumSizes[v]; ++i) {
924: const PetscInt point = points[i];
926: 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);
927: PetscBTClear(label->bt, point - label->pStart);
928: }
929: ISRestoreIndices(label->points[v], &points);
930: }
931: label->stratumSizes[v] = 0;
932: ISDestroy(&label->points[v]);
933: ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_OWN_POINTER, &label->points[v]);
934: PetscObjectSetName((PetscObject) label->points[v], "indices");
935: } else {
936: PetscHashIClear(label->ht[v]);
937: }
938: return(0);
939: }
941: /*@
942: DMLabelFilter - Remove all points outside of [start, end)
944: Input Parameters:
945: + label - the DMLabel
946: . start - the first point
947: - end - the last point
949: Level: intermediate
951: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
952: @*/
953: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
954: {
955: PetscInt v;
959: DMLabelDestroyIndex(label);
960: DMLabelMakeAllValid_Private(label);
961: for (v = 0; v < label->numStrata; ++v) {
962: PetscInt off, q;
963: const PetscInt *points;
964: PetscInt numPointsNew = 0, *pointsNew = NULL;
966: ISGetIndices(label->points[v], &points);
967: for (q = 0; q < label->stratumSizes[v]; ++q)
968: if (points[q] >= start && points[q] < end)
969: numPointsNew++;
970: PetscMalloc1(numPointsNew, &pointsNew);
971: for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
972: if (points[q] >= start && points[q] < end)
973: pointsNew[off++] = points[q];
974: }
975: ISRestoreIndices(label->points[v],&points);
977: label->stratumSizes[v] = numPointsNew;
978: ISDestroy(&label->points[v]);
979: ISCreateGeneral(PETSC_COMM_SELF,numPointsNew, pointsNew, PETSC_OWN_POINTER, &label->points[v]);
980: PetscObjectSetName((PetscObject) label->points[v], "indices");
981: }
982: DMLabelCreateIndex(label, start, end);
983: return(0);
984: }
986: /*@
987: DMLabelPermute - Create a new label with permuted points
989: Input Parameters:
990: + label - the DMLabel
991: - permutation - the point permutation
993: Output Parameter:
994: . labelnew - the new label containing the permuted points
996: Level: intermediate
998: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
999: @*/
1000: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1001: {
1002: const PetscInt *perm;
1003: PetscInt numValues, numPoints, v, q;
1004: PetscErrorCode ierr;
1007: DMLabelMakeAllValid_Private(label);
1008: DMLabelDuplicate(label, labelNew);
1009: DMLabelGetNumValues(*labelNew, &numValues);
1010: ISGetLocalSize(permutation, &numPoints);
1011: ISGetIndices(permutation, &perm);
1012: for (v = 0; v < numValues; ++v) {
1013: const PetscInt size = (*labelNew)->stratumSizes[v];
1014: const PetscInt *points;
1015: PetscInt *pointsNew;
1017: ISGetIndices((*labelNew)->points[v],&points);
1018: PetscMalloc1(size,&pointsNew);
1019: for (q = 0; q < size; ++q) {
1020: const PetscInt point = points[q];
1022: 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);
1023: pointsNew[q] = perm[point];
1024: }
1025: ISRestoreIndices((*labelNew)->points[v],&points);
1026: PetscSortInt(size, pointsNew);
1027: ISDestroy(&((*labelNew)->points[v]));
1028: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1029: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1030: }
1031: ISRestoreIndices(permutation, &perm);
1032: if (label->bt) {
1033: PetscBTDestroy(&label->bt);
1034: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1035: }
1036: return(0);
1037: }
1039: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1040: {
1041: MPI_Comm comm;
1042: PetscInt s, l, nroots, nleaves, dof, offset, size;
1043: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1044: PetscSection rootSection;
1045: PetscSF labelSF;
1049: if (label) {DMLabelMakeAllValid_Private(label);}
1050: PetscObjectGetComm((PetscObject)sf, &comm);
1051: /* Build a section of stratum values per point, generate the according SF
1052: and distribute point-wise stratum values to leaves. */
1053: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1054: PetscSectionCreate(comm, &rootSection);
1055: PetscSectionSetChart(rootSection, 0, nroots);
1056: if (label) {
1057: for (s = 0; s < label->numStrata; ++s) {
1058: const PetscInt *points;
1060: ISGetIndices(label->points[s], &points);
1061: for (l = 0; l < label->stratumSizes[s]; l++) {
1062: PetscSectionGetDof(rootSection, points[l], &dof);
1063: PetscSectionSetDof(rootSection, points[l], dof+1);
1064: }
1065: ISRestoreIndices(label->points[s], &points);
1066: }
1067: }
1068: PetscSectionSetUp(rootSection);
1069: /* Create a point-wise array of stratum values */
1070: PetscSectionGetStorageSize(rootSection, &size);
1071: PetscMalloc1(size, &rootStrata);
1072: PetscCalloc1(nroots, &rootIdx);
1073: if (label) {
1074: for (s = 0; s < label->numStrata; ++s) {
1075: const PetscInt *points;
1077: ISGetIndices(label->points[s], &points);
1078: for (l = 0; l < label->stratumSizes[s]; l++) {
1079: const PetscInt p = points[l];
1080: PetscSectionGetOffset(rootSection, p, &offset);
1081: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1082: }
1083: ISRestoreIndices(label->points[s], &points);
1084: }
1085: }
1086: /* Build SF that maps label points to remote processes */
1087: PetscSectionCreate(comm, leafSection);
1088: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1089: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1090: PetscFree(remoteOffsets);
1091: /* Send the strata for each point over the derived SF */
1092: PetscSectionGetStorageSize(*leafSection, &size);
1093: PetscMalloc1(size, leafStrata);
1094: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1095: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1096: /* Clean up */
1097: PetscFree(rootStrata);
1098: PetscFree(rootIdx);
1099: PetscSectionDestroy(&rootSection);
1100: PetscSFDestroy(&labelSF);
1101: return(0);
1102: }
1104: /*@
1105: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1107: Input Parameters:
1108: + label - the DMLabel
1109: - sf - the map from old to new distribution
1111: Output Parameter:
1112: . labelnew - the new redistributed label
1114: Level: intermediate
1116: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1117: @*/
1118: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1119: {
1120: MPI_Comm comm;
1121: PetscSection leafSection;
1122: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1123: PetscInt *leafStrata, *strataIdx;
1124: PetscInt **points;
1125: char *name;
1126: PetscInt nameSize;
1127: PetscHashI stratumHash;
1128: PETSC_UNUSED PetscHashIIter ret, iter;
1129: size_t len = 0;
1130: PetscMPIInt rank;
1134: if (label) {DMLabelMakeAllValid_Private(label);}
1135: PetscObjectGetComm((PetscObject)sf, &comm);
1136: MPI_Comm_rank(comm, &rank);
1137: /* Bcast name */
1138: if (!rank) {PetscStrlen(label->name, &len);}
1139: nameSize = len;
1140: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1141: PetscMalloc1(nameSize+1, &name);
1142: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1143: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1144: DMLabelCreate(name, labelNew);
1145: PetscFree(name);
1146: /* Bcast defaultValue */
1147: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1148: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1149: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1150: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1151: /* Determine received stratum values and initialise new label*/
1152: PetscHashICreate(stratumHash);
1153: PetscSectionGetStorageSize(leafSection, &size);
1154: for (p = 0; p < size; ++p) PetscHashIPut(stratumHash, leafStrata[p], ret, iter);
1155: PetscHashISize(stratumHash, (*labelNew)->numStrata);
1156: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1157: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1158: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1159: /* Turn leafStrata into indices rather than stratum values */
1160: offset = 0;
1161: PetscHashIGetKeys(stratumHash, &offset, (*labelNew)->stratumValues);
1162: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1163: for (p = 0; p < size; ++p) {
1164: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1165: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1166: }
1167: }
1168: /* Rebuild the point strata on the receiver */
1169: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1170: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1171: for (p=pStart; p<pEnd; p++) {
1172: PetscSectionGetDof(leafSection, p, &dof);
1173: PetscSectionGetOffset(leafSection, p, &offset);
1174: for (s=0; s<dof; s++) {
1175: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1176: }
1177: }
1178: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1179: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1180: PetscMalloc1((*labelNew)->numStrata,&points);
1181: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1182: PetscHashICreate((*labelNew)->ht[s]);
1183: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1184: }
1185: /* Insert points into new strata */
1186: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1187: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1188: for (p=pStart; p<pEnd; p++) {
1189: PetscSectionGetDof(leafSection, p, &dof);
1190: PetscSectionGetOffset(leafSection, p, &offset);
1191: for (s=0; s<dof; s++) {
1192: stratum = leafStrata[offset+s];
1193: points[stratum][strataIdx[stratum]++] = p;
1194: }
1195: }
1196: for (s = 0; s < (*labelNew)->numStrata; s++) {
1197: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1198: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1199: }
1200: PetscFree(points);
1201: PetscHashIDestroy(stratumHash);
1202: PetscFree(leafStrata);
1203: PetscFree(strataIdx);
1204: PetscSectionDestroy(&leafSection);
1205: return(0);
1206: }
1208: /*@
1209: DMLabelGather - Gather all label values from leafs into roots
1211: Input Parameters:
1212: + label - the DMLabel
1213: - sf - the Star Forest point communication map
1215: Output Parameters:
1216: . labelNew - the new DMLabel with localised leaf values
1218: Level: developer
1220: Note: This is the inverse operation to DMLabelDistribute.
1222: .seealso: DMLabelDistribute()
1223: @*/
1224: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1225: {
1226: MPI_Comm comm;
1227: PetscSection rootSection;
1228: PetscSF sfLabel;
1229: PetscSFNode *rootPoints, *leafPoints;
1230: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1231: const PetscInt *rootDegree, *ilocal;
1232: PetscInt *rootStrata;
1233: char *name;
1234: PetscInt nameSize;
1235: size_t len = 0;
1236: PetscMPIInt rank, size;
1240: PetscObjectGetComm((PetscObject)sf, &comm);
1241: MPI_Comm_rank(comm, &rank);
1242: MPI_Comm_size(comm, &size);
1243: /* Bcast name */
1244: if (!rank) {PetscStrlen(label->name, &len);}
1245: nameSize = len;
1246: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1247: PetscMalloc1(nameSize+1, &name);
1248: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1249: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1250: DMLabelCreate(name, labelNew);
1251: PetscFree(name);
1252: /* Gather rank/index pairs of leaves into local roots to build
1253: an inverse, multi-rooted SF. Note that this ignores local leaf
1254: indexing due to the use of the multiSF in PetscSFGather. */
1255: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1256: PetscMalloc1(nroots, &leafPoints);
1257: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1258: for (p = 0; p < nleaves; p++) {
1259: leafPoints[ilocal[p]].index = ilocal[p];
1260: leafPoints[ilocal[p]].rank = rank;
1261: }
1262: PetscSFComputeDegreeBegin(sf, &rootDegree);
1263: PetscSFComputeDegreeEnd(sf, &rootDegree);
1264: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1265: PetscMalloc1(nmultiroots, &rootPoints);
1266: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1267: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1268: PetscSFCreate(comm,& sfLabel);
1269: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1270: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1271: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1272: /* Rebuild the point strata on the receiver */
1273: for (p = 0, idx = 0; p < nroots; p++) {
1274: for (d = 0; d < rootDegree[p]; d++) {
1275: PetscSectionGetDof(rootSection, idx+d, &dof);
1276: PetscSectionGetOffset(rootSection, idx+d, &offset);
1277: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1278: }
1279: idx += rootDegree[p];
1280: }
1281: PetscFree(leafPoints);
1282: PetscFree(rootStrata);
1283: PetscSectionDestroy(&rootSection);
1284: PetscSFDestroy(&sfLabel);
1285: return(0);
1286: }
1288: /*@
1289: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1291: Input Parameter:
1292: . label - the DMLabel
1294: Output Parameters:
1295: + section - the section giving offsets for each stratum
1296: - is - An IS containing all the label points
1298: Level: developer
1300: .seealso: DMLabelDistribute()
1301: @*/
1302: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1303: {
1304: IS vIS;
1305: const PetscInt *values;
1306: PetscInt *points;
1307: PetscInt nV, vS = 0, vE = 0, v, N;
1308: PetscErrorCode ierr;
1311: DMLabelGetNumValues(label, &nV);
1312: DMLabelGetValueIS(label, &vIS);
1313: ISGetIndices(vIS, &values);
1314: if (nV) {vS = values[0]; vE = values[0]+1;}
1315: for (v = 1; v < nV; ++v) {
1316: vS = PetscMin(vS, values[v]);
1317: vE = PetscMax(vE, values[v]+1);
1318: }
1319: PetscSectionCreate(PETSC_COMM_SELF, section);
1320: PetscSectionSetChart(*section, vS, vE);
1321: for (v = 0; v < nV; ++v) {
1322: PetscInt n;
1324: DMLabelGetStratumSize(label, values[v], &n);
1325: PetscSectionSetDof(*section, values[v], n);
1326: }
1327: PetscSectionSetUp(*section);
1328: PetscSectionGetStorageSize(*section, &N);
1329: PetscMalloc1(N, &points);
1330: for (v = 0; v < nV; ++v) {
1331: IS is;
1332: const PetscInt *spoints;
1333: PetscInt dof, off, p;
1335: PetscSectionGetDof(*section, values[v], &dof);
1336: PetscSectionGetOffset(*section, values[v], &off);
1337: DMLabelGetStratumIS(label, values[v], &is);
1338: ISGetIndices(is, &spoints);
1339: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1340: ISRestoreIndices(is, &spoints);
1341: ISDestroy(&is);
1342: }
1343: ISRestoreIndices(vIS, &values);
1344: ISDestroy(&vIS);
1345: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1346: return(0);
1347: }
1349: /*@
1350: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1351: the local section and an SF describing the section point overlap.
1353: Input Parameters:
1354: + s - The PetscSection for the local field layout
1355: . sf - The SF describing parallel layout of the section points
1356: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1357: . label - The label specifying the points
1358: - labelValue - The label stratum specifying the points
1360: Output Parameter:
1361: . gsection - The PetscSection for the global field layout
1363: Note: This gives negative sizes and offsets to points not owned by this process
1365: Level: developer
1367: .seealso: PetscSectionCreate()
1368: @*/
1369: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1370: {
1371: PetscInt *neg = NULL, *tmpOff = NULL;
1372: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1376: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1377: PetscSectionGetChart(s, &pStart, &pEnd);
1378: PetscSectionSetChart(*gsection, pStart, pEnd);
1379: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1380: if (nroots >= 0) {
1381: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1382: PetscCalloc1(nroots, &neg);
1383: if (nroots > pEnd-pStart) {
1384: PetscCalloc1(nroots, &tmpOff);
1385: } else {
1386: tmpOff = &(*gsection)->atlasDof[-pStart];
1387: }
1388: }
1389: /* Mark ghost points with negative dof */
1390: for (p = pStart; p < pEnd; ++p) {
1391: PetscInt value;
1393: DMLabelGetValue(label, p, &value);
1394: if (value != labelValue) continue;
1395: PetscSectionGetDof(s, p, &dof);
1396: PetscSectionSetDof(*gsection, p, dof);
1397: PetscSectionGetConstraintDof(s, p, &cdof);
1398: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1399: if (neg) neg[p] = -(dof+1);
1400: }
1401: PetscSectionSetUpBC(*gsection);
1402: if (nroots >= 0) {
1403: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1404: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1405: if (nroots > pEnd-pStart) {
1406: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1407: }
1408: }
1409: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1410: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1411: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1412: (*gsection)->atlasOff[p] = off;
1413: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1414: }
1415: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1416: globalOff -= off;
1417: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1418: (*gsection)->atlasOff[p] += globalOff;
1419: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1420: }
1421: /* Put in negative offsets for ghost points */
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)->atlasOff[p-pStart] = tmpOff[p];}
1427: }
1428: }
1429: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1430: PetscFree(neg);
1431: return(0);
1432: }
1434: typedef struct _n_PetscSectionSym_Label
1435: {
1436: DMLabel label;
1437: PetscCopyMode *modes;
1438: PetscInt *sizes;
1439: const PetscInt ***perms;
1440: const PetscScalar ***rots;
1441: PetscInt (*minMaxOrients)[2];
1442: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1443: } PetscSectionSym_Label;
1445: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1446: {
1447: PetscInt i, j;
1448: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1449: PetscErrorCode ierr;
1452: for (i = 0; i <= sl->numStrata; i++) {
1453: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1454: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1455: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1456: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1457: }
1458: if (sl->perms[i]) {
1459: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1461: PetscFree(perms);
1462: }
1463: if (sl->rots[i]) {
1464: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1466: PetscFree(rots);
1467: }
1468: }
1469: }
1470: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1471: DMLabelDestroy(&sl->label);
1472: sl->numStrata = 0;
1473: return(0);
1474: }
1476: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1477: {
1481: PetscSectionSymLabelReset(sym);
1482: PetscFree(sym->data);
1483: return(0);
1484: }
1486: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1487: {
1488: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1489: PetscBool isAscii;
1490: DMLabel label = sl->label;
1491: PetscErrorCode ierr;
1494: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1495: if (isAscii) {
1496: PetscInt i, j, k;
1497: PetscViewerFormat format;
1499: PetscViewerGetFormat(viewer,&format);
1500: if (label) {
1501: PetscViewerGetFormat(viewer,&format);
1502: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1503: PetscViewerASCIIPushTab(viewer);
1504: DMLabelView(label, viewer);
1505: PetscViewerASCIIPopTab(viewer);
1506: } else {
1507: PetscViewerASCIIPrintf(viewer," Label '%s'\n",sl->label->name);
1508: }
1509: } else {
1510: PetscViewerASCIIPrintf(viewer, "No label given\n");
1511: }
1512: PetscViewerASCIIPushTab(viewer);
1513: for (i = 0; i <= sl->numStrata; i++) {
1514: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1516: if (!(sl->perms[i] || sl->rots[i])) {
1517: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1518: } else {
1519: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1520: PetscViewerASCIIPushTab(viewer);
1521: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1522: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1523: PetscViewerASCIIPushTab(viewer);
1524: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1525: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1526: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1527: } else {
1528: PetscInt tab;
1530: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1531: PetscViewerASCIIPushTab(viewer);
1532: PetscViewerASCIIGetTab(viewer,&tab);
1533: if (sl->perms[i] && sl->perms[i][j]) {
1534: PetscViewerASCIIPrintf(viewer,"Permutation:");
1535: PetscViewerASCIISetTab(viewer,0);
1536: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1537: PetscViewerASCIIPrintf(viewer,"\n");
1538: PetscViewerASCIISetTab(viewer,tab);
1539: }
1540: if (sl->rots[i] && sl->rots[i][j]) {
1541: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1542: PetscViewerASCIISetTab(viewer,0);
1543: #if defined(PETSC_USE_COMPLEX)
1544: 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]));}
1545: #else
1546: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1547: #endif
1548: PetscViewerASCIIPrintf(viewer,"\n");
1549: PetscViewerASCIISetTab(viewer,tab);
1550: }
1551: PetscViewerASCIIPopTab(viewer);
1552: }
1553: }
1554: PetscViewerASCIIPopTab(viewer);
1555: }
1556: PetscViewerASCIIPopTab(viewer);
1557: }
1558: }
1559: PetscViewerASCIIPopTab(viewer);
1560: }
1561: return(0);
1562: }
1564: /*@
1565: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1567: Logically collective on sym
1569: Input parameters:
1570: + sym - the section symmetries
1571: - label - the DMLabel describing the types of points
1573: Level: developer:
1575: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1576: @*/
1577: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1578: {
1579: PetscSectionSym_Label *sl;
1584: sl = (PetscSectionSym_Label *) sym->data;
1585: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1586: if (label) {
1587: label->refct++;
1588: sl->label = label;
1589: DMLabelGetNumValues(label,&sl->numStrata);
1590: 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);
1591: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1592: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1593: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1594: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1595: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1596: }
1597: return(0);
1598: }
1600: /*@C
1601: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1603: Logically collective on PetscSectionSym
1605: InputParameters:
1606: + sys - the section symmetries
1607: . stratum - the stratum value in the label that we are assigning symmetries for
1608: . size - the number of dofs for points in the stratum of the label
1609: . minOrient - the smallest orientation for a point in this stratum
1610: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1611: . mode - how sym should copy the perms and rots arrays
1612: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
1613: + 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
1615: Level: developer
1617: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1618: @*/
1619: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1620: {
1621: PetscInt i, j, k;
1622: PetscSectionSym_Label *sl;
1627: sl = (PetscSectionSym_Label *) sym->data;
1628: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1629: for (i = 0; i <= sl->numStrata; i++) {
1630: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
1632: if (stratum == value) break;
1633: }
1634: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,sl->label->name);
1635: sl->sizes[i] = size;
1636: sl->modes[i] = mode;
1637: sl->minMaxOrients[i][0] = minOrient;
1638: sl->minMaxOrients[i][1] = maxOrient;
1639: if (mode == PETSC_COPY_VALUES) {
1640: if (perms) {
1641: PetscInt **ownPerms;
1643: PetscCalloc1(maxOrient - minOrient,&ownPerms);
1644: for (j = 0; j < maxOrient-minOrient; j++) {
1645: if (perms[j]) {
1646: PetscMalloc1(size,&ownPerms[j]);
1647: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
1648: }
1649: }
1650: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
1651: }
1652: if (rots) {
1653: PetscScalar **ownRots;
1655: PetscCalloc1(maxOrient - minOrient,&ownRots);
1656: for (j = 0; j < maxOrient-minOrient; j++) {
1657: if (rots[j]) {
1658: PetscMalloc1(size,&ownRots[j]);
1659: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
1660: }
1661: }
1662: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
1663: }
1664: } else {
1665: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
1666: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
1667: }
1668: return(0);
1669: }
1671: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
1672: {
1673: PetscInt i, j, numStrata;
1674: PetscSectionSym_Label *sl;
1675: DMLabel label;
1676: PetscErrorCode ierr;
1679: sl = (PetscSectionSym_Label *) sym->data;
1680: numStrata = sl->numStrata;
1681: label = sl->label;
1682: for (i = 0; i < numPoints; i++) {
1683: PetscInt point = points[2*i];
1684: PetscInt ornt = points[2*i+1];
1686: for (j = 0; j < numStrata; j++) {
1687: if (label->validIS[j]) {
1688: PetscInt k;
1690: ISLocate(label->points[j],point,&k);
1691: if (k >= 0) break;
1692: } else {
1693: PetscBool has;
1695: PetscHashIHasKey(label->ht[j], point, has);
1696: if (has) break;
1697: }
1698: }
1699: 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);
1700: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
1701: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
1702: }
1703: return(0);
1704: }
1706: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
1707: {
1708: PetscSectionSym_Label *sl;
1709: PetscErrorCode ierr;
1712: PetscNewLog(sym,&sl);
1713: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
1714: sym->ops->view = PetscSectionSymView_Label;
1715: sym->ops->destroy = PetscSectionSymDestroy_Label;
1716: sym->data = (void *) sl;
1717: return(0);
1718: }
1720: /*@
1721: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
1723: Collective
1725: Input Parameters:
1726: + comm - the MPI communicator for the new symmetry
1727: - label - the label defining the strata
1729: Output Parameters:
1730: . sym - the section symmetries
1732: Level: developer
1734: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
1735: @*/
1736: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
1737: {
1738: PetscErrorCode ierr;
1741: DMInitializePackage();
1742: PetscSectionSymCreate(comm,sym);
1743: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
1744: PetscSectionSymLabelSetLabel(*sym,label);
1745: return(0);
1746: }