Actual source code: dmlabel.c
petsc-3.10.5 2019-03-28
1: #include <petscdm.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petscsf.h>
6: /*@C
7: DMLabelCreate - Create a DMLabel object, which is a multimap
9: Input 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: PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
67: PetscMalloc1(label->stratumSizes[v], &pointArray);
68: PetscHSetIGetElems(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: PetscHSetIClear(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: PetscInt p;
131: const PetscInt *points;
134: if (!label->validIS[v]) return 0;
136: 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);
137: if (label->points[v]) {
138: ISGetIndices(label->points[v],&points);
139: for (p = 0; p < label->stratumSizes[v]; ++p) {
140: PetscHSetIAdd(label->ht[v], points[p]);
141: }
142: ISRestoreIndices(label->points[v],&points);
143: ISDestroy(&(label->points[v]));
144: }
145: label->validIS[v] = PETSC_FALSE;
146: return(0);
147: }
149: PetscErrorCode DMLabelGetState(DMLabel label, PetscObjectState *state)
150: {
153: *state = label->state;
154: return(0);
155: }
157: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
158: {
159: PetscInt v;
161: for (*index = -1, v = 0; v < label->numStrata; ++v)
162: if (label->stratumValues[v] == value) { *index = v; break; }
163: return(0);
164: }
166: static PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
167: {
168: PetscInt v, *tmpV, *tmpS;
169: IS *tmpP;
170: PetscHSetI *tmpH;
171: PetscBool *tmpB;
175: PetscMalloc1((label->numStrata+1), &tmpV);
176: PetscMalloc1((label->numStrata+1), &tmpS);
177: PetscMalloc1((label->numStrata+1), &tmpH);
178: PetscMalloc1((label->numStrata+1), &tmpP);
179: PetscMalloc1((label->numStrata+1), &tmpB);
180: for (v = 0; v < label->numStrata; ++v) {
181: tmpV[v] = label->stratumValues[v];
182: tmpS[v] = label->stratumSizes[v];
183: tmpH[v] = label->ht[v];
184: tmpP[v] = label->points[v];
185: tmpB[v] = label->validIS[v];
186: }
187: tmpV[v] = value;
188: tmpS[v] = 0;
189: PetscHSetICreate(&tmpH[v]);
190: ISCreateStride(PETSC_COMM_SELF,0,0,1,&tmpP[v]);
191: tmpB[v] = PETSC_TRUE;
193: ++label->numStrata;
194: PetscFree(label->stratumValues);
195: PetscFree(label->stratumSizes);
196: PetscFree(label->ht);
197: PetscFree(label->points);
198: PetscFree(label->validIS);
199: label->stratumValues = tmpV;
200: label->stratumSizes = tmpS;
201: label->ht = tmpH;
202: label->points = tmpP;
203: label->validIS = tmpB;
205: *index = v;
206: return(0);
207: }
209: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
210: {
211: PetscInt v;
215: DMLabelLookupStratum(label, value, &v);
216: if (v < 0) {DMLabelNewStratum(label, value, &v);}
217: return(0);
218: }
220: /*@C
221: DMLabelGetName - Return the name of a DMLabel object
223: Input parameter:
224: . label - The DMLabel
226: Output parameter:
227: . name - The label name
229: Level: beginner
231: .seealso: DMLabelCreate()
232: @*/
233: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
234: {
237: *name = label->name;
238: return(0);
239: }
241: /*@C
242: DMLabelSetName - Sets the name of a DMLabel object
244: Input parameters:
245: + label - The DMLabel
246: - name - The label name
248: Level: beginner
250: .seealso: DMLabelCreate()
251: @*/
252: PetscErrorCode DMLabelSetName(DMLabel label, const char *name)
253: {
258: PetscFree(label->name);
259: PetscStrallocpy(name, &label->name);
260: return(0);
261: }
263: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
264: {
265: PetscInt v;
266: PetscMPIInt rank;
270: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
271: PetscViewerASCIIPushSynchronized(viewer);
272: if (label) {
273: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
274: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
275: for (v = 0; v < label->numStrata; ++v) {
276: const PetscInt value = label->stratumValues[v];
277: const PetscInt *points;
278: PetscInt p;
280: ISGetIndices(label->points[v],&points);
281: for (p = 0; p < label->stratumSizes[v]; ++p) {
282: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
283: }
284: ISRestoreIndices(label->points[v],&points);
285: }
286: }
287: PetscViewerFlush(viewer);
288: PetscViewerASCIIPopSynchronized(viewer);
289: return(0);
290: }
292: /*@C
293: DMLabelView - View the label
295: Input Parameters:
296: + label - The DMLabel
297: - viewer - The PetscViewer
299: Level: intermediate
301: .seealso: DMLabelCreate(), DMLabelDestroy()
302: @*/
303: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
304: {
305: PetscBool iascii;
310: if (label) {DMLabelMakeAllValid_Private(label);}
311: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
312: if (iascii) {
313: DMLabelView_Ascii(label, viewer);
314: }
315: return(0);
316: }
318: /*@
319: DMLabelDestroy - Destroys a DMLabel
321: Input Parameter:
322: . label - The DMLabel
324: Level: beginner
326: .seealso: DMLabelCreate()
327: @*/
328: PetscErrorCode DMLabelDestroy(DMLabel *label)
329: {
330: PetscInt v;
334: if (!(*label)) return(0);
335: if (--(*label)->refct > 0) return(0);
336: PetscFree((*label)->name);
337: PetscFree((*label)->stratumValues);
338: PetscFree((*label)->stratumSizes);
339: for (v = 0; v < (*label)->numStrata; ++v) {
340: PetscHSetIDestroy(&(*label)->ht[v]);
341: ISDestroy(&(*label)->points[v]);
342: }
343: PetscFree((*label)->ht);
344: PetscFree((*label)->points);
345: PetscFree((*label)->validIS);
346: PetscBTDestroy(&(*label)->bt);
347: PetscFree(*label);
348: return(0);
349: }
351: /*@
352: DMLabelDuplicate - Duplicates a DMLabel
354: Input Parameter:
355: . label - The DMLabel
357: Output Parameter:
358: . labelnew - location to put new vector
360: Level: intermediate
362: .seealso: DMLabelCreate(), DMLabelDestroy()
363: @*/
364: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
365: {
366: PetscInt v;
370: DMLabelMakeAllValid_Private(label);
371: PetscNew(labelnew);
372: PetscStrallocpy(label->name, &(*labelnew)->name);
374: (*labelnew)->refct = 1;
375: (*labelnew)->numStrata = label->numStrata;
376: (*labelnew)->defaultValue = label->defaultValue;
377: if (label->numStrata) {
378: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
379: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
380: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
381: PetscMalloc1(label->numStrata, &(*labelnew)->points);
382: PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
383: /* Could eliminate unused space here */
384: for (v = 0; v < label->numStrata; ++v) {
385: PetscHSetICreate(&(*labelnew)->ht[v]);
386: (*labelnew)->validIS[v] = PETSC_TRUE;
387: (*labelnew)->stratumValues[v] = label->stratumValues[v];
388: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
389: PetscObjectReference((PetscObject) (label->points[v]));
390: (*labelnew)->points[v] = label->points[v];
391: }
392: }
393: (*labelnew)->pStart = -1;
394: (*labelnew)->pEnd = -1;
395: (*labelnew)->bt = NULL;
396: return(0);
397: }
399: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
400: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
401: {
402: PetscInt v;
406: DMLabelDestroyIndex(label);
407: DMLabelMakeAllValid_Private(label);
408: label->pStart = pStart;
409: label->pEnd = pEnd;
410: PetscBTCreate(pEnd - pStart, &label->bt);
411: for (v = 0; v < label->numStrata; ++v) {
412: const PetscInt *points;
413: PetscInt i;
415: ISGetIndices(label->points[v], &points);
416: for (i = 0; i < label->stratumSizes[v]; ++i) {
417: const PetscInt point = points[i];
419: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
420: PetscBTSet(label->bt, point - pStart);
421: }
422: ISRestoreIndices(label->points[v], &points);
423: }
424: return(0);
425: }
427: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
428: {
432: label->pStart = -1;
433: label->pEnd = -1;
434: PetscBTDestroy(&label->bt);
435: return(0);
436: }
438: /*@
439: DMLabelHasValue - Determine whether a label assigns the value to any point
441: Input Parameters:
442: + label - the DMLabel
443: - value - the value
445: Output Parameter:
446: . contains - Flag indicating whether the label maps this value to any point
448: Level: developer
450: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
451: @*/
452: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
453: {
454: PetscInt v;
459: DMLabelLookupStratum(label, value, &v);
460: *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
461: return(0);
462: }
464: /*@
465: DMLabelHasPoint - Determine whether a label assigns a value to a point
467: Input Parameters:
468: + label - the DMLabel
469: - point - the point
471: Output Parameter:
472: . contains - Flag indicating whether the label maps this point to a value
474: Note: The user must call DMLabelCreateIndex() before this function.
476: Level: developer
478: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
479: @*/
480: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
481: {
486: DMLabelMakeAllValid_Private(label);
487: #if defined(PETSC_USE_DEBUG)
488: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
489: 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);
490: #endif
491: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
492: return(0);
493: }
495: /*@
496: DMLabelStratumHasPoint - Return true if the stratum contains a point
498: Input Parameters:
499: + label - the DMLabel
500: . value - the stratum value
501: - point - the point
503: Output Parameter:
504: . contains - true if the stratum contains the point
506: Level: intermediate
508: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
509: @*/
510: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
511: {
512: PetscInt v;
517: *contains = PETSC_FALSE;
518: DMLabelLookupStratum(label, value, &v);
519: if (v < 0) return(0);
521: if (label->validIS[v]) {
522: PetscInt i;
524: ISLocate(label->points[v], point, &i);
525: if (i >= 0) *contains = PETSC_TRUE;
526: } else {
527: PetscBool has;
529: PetscHSetIHas(label->ht[v], point, &has);
530: if (has) *contains = PETSC_TRUE;
531: }
532: return(0);
533: }
535: /*@
536: DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
537: When a label is created, it is initialized to -1.
539: Input parameter:
540: . label - a DMLabel object
542: Output parameter:
543: . defaultValue - the default value
545: Level: beginner
547: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
548: @*/
549: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
550: {
552: *defaultValue = label->defaultValue;
553: return(0);
554: }
556: /*@
557: DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
558: When a label is created, it is initialized to -1.
560: Input parameter:
561: . label - a DMLabel object
563: Output parameter:
564: . defaultValue - the default value
566: Level: beginner
568: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
569: @*/
570: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
571: {
573: label->defaultValue = defaultValue;
574: return(0);
575: }
577: /*@
578: 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())
580: Input Parameters:
581: + label - the DMLabel
582: - point - the point
584: Output Parameter:
585: . value - The point value, or -1
587: Level: intermediate
589: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
590: @*/
591: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
592: {
593: PetscInt v;
598: *value = label->defaultValue;
599: for (v = 0; v < label->numStrata; ++v) {
600: if (label->validIS[v]) {
601: PetscInt i;
603: ISLocate(label->points[v], point, &i);
604: if (i >= 0) {
605: *value = label->stratumValues[v];
606: break;
607: }
608: } else {
609: PetscBool has;
611: PetscHSetIHas(label->ht[v], point, &has);
612: if (has) {
613: *value = label->stratumValues[v];
614: break;
615: }
616: }
617: }
618: return(0);
619: }
621: /*@
622: 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.
624: Input Parameters:
625: + label - the DMLabel
626: . point - the point
627: - value - The point value
629: Level: intermediate
631: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
632: @*/
633: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
634: {
635: PetscInt v;
636: PetscErrorCode ierr;
639: /* Find label value, add new entry if needed */
640: if (value == label->defaultValue) return(0);
641: DMLabelLookupStratum(label, value, &v);
642: if (v < 0) {DMLabelNewStratum(label, value, &v);}
643: /* Set key */
644: DMLabelMakeInvalid_Private(label, v);
645: PetscHSetIAdd(label->ht[v], point);
646: return(0);
647: }
649: /*@
650: DMLabelClearValue - Clear the value a label assigns to a point
652: Input Parameters:
653: + label - the DMLabel
654: . point - the point
655: - value - The point value
657: Level: intermediate
659: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
660: @*/
661: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
662: {
663: PetscInt v;
667: /* Find label value */
668: DMLabelLookupStratum(label, value, &v);
669: if (v < 0) return(0);
671: if (label->bt) {
672: 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);
673: PetscBTClear(label->bt, point - label->pStart);
674: }
676: /* Delete key */
677: DMLabelMakeInvalid_Private(label, v);
678: PetscHSetIDel(label->ht[v], point);
679: return(0);
680: }
682: /*@
683: DMLabelInsertIS - Set all points in the IS to a value
685: Input Parameters:
686: + label - the DMLabel
687: . is - the point IS
688: - value - The point value
690: Level: intermediate
692: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
693: @*/
694: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
695: {
696: PetscInt v, n, p;
697: const PetscInt *points;
698: PetscErrorCode ierr;
702: /* Find label value, add new entry if needed */
703: if (value == label->defaultValue) return(0);
704: DMLabelLookupStratum(label, value, &v);
705: if (v < 0) {DMLabelNewStratum(label, value, &v);}
706: /* Set keys */
707: DMLabelMakeInvalid_Private(label, v);
708: ISGetLocalSize(is, &n);
709: ISGetIndices(is, &points);
710: for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
711: ISRestoreIndices(is, &points);
712: return(0);
713: }
715: /*@
716: DMLabelGetNumValues - Get the number of values that the DMLabel takes
718: Input Parameter:
719: . label - the DMLabel
721: Output Paramater:
722: . numValues - the number of values
724: Level: intermediate
726: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
727: @*/
728: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
729: {
732: *numValues = label->numStrata;
733: return(0);
734: }
736: /*@
737: DMLabelGetValueIS - Get an IS of all values that the DMlabel takes
739: Input Parameter:
740: . label - the DMLabel
742: Output Paramater:
743: . is - the value IS
745: Level: intermediate
747: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
748: @*/
749: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
750: {
755: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
756: return(0);
757: }
759: /*@
760: DMLabelHasStratum - Determine whether points exist with the given value
762: Input Parameters:
763: + label - the DMLabel
764: - value - the stratum value
766: Output Paramater:
767: . exists - Flag saying whether points exist
769: Level: intermediate
771: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
772: @*/
773: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
774: {
775: PetscInt v;
780: DMLabelLookupStratum(label, value, &v);
781: *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
782: return(0);
783: }
785: /*@
786: DMLabelGetStratumSize - Get the size of a stratum
788: Input Parameters:
789: + label - the DMLabel
790: - value - the stratum value
792: Output Paramater:
793: . size - The number of points in the stratum
795: Level: intermediate
797: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
798: @*/
799: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
800: {
801: PetscInt v;
806: *size = 0;
807: DMLabelLookupStratum(label, value, &v);
808: if (v < 0) return(0);
809: DMLabelMakeValid_Private(label, v);
810: *size = label->stratumSizes[v];
811: return(0);
812: }
814: /*@
815: DMLabelGetStratumBounds - Get the largest and smallest point of a stratum
817: Input Parameters:
818: + label - the DMLabel
819: - value - the stratum value
821: Output Paramaters:
822: + start - the smallest point in the stratum
823: - end - the largest point in the stratum
825: Level: intermediate
827: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
828: @*/
829: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
830: {
831: PetscInt v, min, max;
837: DMLabelLookupStratum(label, value, &v);
838: if (v < 0) return(0);
839: DMLabelMakeValid_Private(label, v);
840: if (label->stratumSizes[v] <= 0) return(0);
841: ISGetMinMax(label->points[v], &min, &max);
842: if (start) *start = min;
843: if (end) *end = max+1;
844: return(0);
845: }
847: /*@
848: DMLabelGetStratumIS - Get an IS with the stratum points
850: Input Parameters:
851: + label - the DMLabel
852: - value - the stratum value
854: Output Paramater:
855: . points - The stratum points
857: Level: intermediate
859: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
860: @*/
861: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
862: {
863: PetscInt v;
868: *points = NULL;
869: DMLabelLookupStratum(label, value, &v);
870: if (v < 0) return(0);
871: DMLabelMakeValid_Private(label, v);
872: PetscObjectReference((PetscObject) label->points[v]);
873: *points = label->points[v];
874: return(0);
875: }
877: /*@
878: DMLabelSetStratumIS - Set the stratum points using an IS
880: Input Parameters:
881: + label - the DMLabel
882: . value - the stratum value
883: - points - The stratum points
885: Level: intermediate
887: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
888: @*/
889: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
890: {
891: PetscInt v;
895: DMLabelLookupStratum(label, value, &v);
896: if (v < 0) {DMLabelNewStratum(label, value, &v);}
897: if (is == label->points[v]) return(0);
898: DMLabelClearStratum(label, value);
899: ISGetLocalSize(is, &(label->stratumSizes[v]));
900: PetscObjectReference((PetscObject)is);
901: ISDestroy(&(label->points[v]));
902: label->points[v] = is;
903: label->validIS[v] = PETSC_TRUE;
904: if (label->bt) {
905: const PetscInt *points;
906: PetscInt p;
908: ISGetIndices(is,&points);
909: for (p = 0; p < label->stratumSizes[v]; ++p) {
910: const PetscInt point = points[p];
912: 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);
913: PetscBTSet(label->bt, point - label->pStart);
914: }
915: }
916: return(0);
917: }
919: /*@
920: DMLabelClearStratum - Remove a stratum
922: Input Parameters:
923: + label - the DMLabel
924: - value - the stratum value
926: Level: intermediate
928: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
929: @*/
930: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
931: {
932: PetscInt v;
936: DMLabelLookupStratum(label, value, &v);
937: if (v < 0) return(0);
938: if (label->validIS[v]) {
939: if (label->bt) {
940: PetscInt i;
941: const PetscInt *points;
943: ISGetIndices(label->points[v], &points);
944: for (i = 0; i < label->stratumSizes[v]; ++i) {
945: const PetscInt point = points[i];
947: 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);
948: PetscBTClear(label->bt, point - label->pStart);
949: }
950: ISRestoreIndices(label->points[v], &points);
951: }
952: label->stratumSizes[v] = 0;
953: ISDestroy(&label->points[v]);
954: ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_OWN_POINTER, &label->points[v]);
955: PetscObjectSetName((PetscObject) label->points[v], "indices");
956: } else {
957: PetscHSetIClear(label->ht[v]);
958: }
959: return(0);
960: }
962: /*@
963: DMLabelFilter - Remove all points outside of [start, end)
965: Input Parameters:
966: + label - the DMLabel
967: . start - the first point
968: - end - the last point
970: Level: intermediate
972: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
973: @*/
974: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
975: {
976: PetscInt v;
980: DMLabelDestroyIndex(label);
981: DMLabelMakeAllValid_Private(label);
982: for (v = 0; v < label->numStrata; ++v) {
983: PetscInt off, q;
984: const PetscInt *points;
985: PetscInt numPointsNew = 0, *pointsNew = NULL;
987: ISGetIndices(label->points[v], &points);
988: for (q = 0; q < label->stratumSizes[v]; ++q)
989: if (points[q] >= start && points[q] < end)
990: numPointsNew++;
991: PetscMalloc1(numPointsNew, &pointsNew);
992: for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
993: if (points[q] >= start && points[q] < end)
994: pointsNew[off++] = points[q];
995: }
996: ISRestoreIndices(label->points[v],&points);
998: label->stratumSizes[v] = numPointsNew;
999: ISDestroy(&label->points[v]);
1000: ISCreateGeneral(PETSC_COMM_SELF,numPointsNew, pointsNew, PETSC_OWN_POINTER, &label->points[v]);
1001: PetscObjectSetName((PetscObject) label->points[v], "indices");
1002: }
1003: DMLabelCreateIndex(label, start, end);
1004: return(0);
1005: }
1007: /*@
1008: DMLabelPermute - Create a new label with permuted points
1010: Input Parameters:
1011: + label - the DMLabel
1012: - permutation - the point permutation
1014: Output Parameter:
1015: . labelnew - the new label containing the permuted points
1017: Level: intermediate
1019: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1020: @*/
1021: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1022: {
1023: const PetscInt *perm;
1024: PetscInt numValues, numPoints, v, q;
1025: PetscErrorCode ierr;
1028: DMLabelMakeAllValid_Private(label);
1029: DMLabelDuplicate(label, labelNew);
1030: DMLabelGetNumValues(*labelNew, &numValues);
1031: ISGetLocalSize(permutation, &numPoints);
1032: ISGetIndices(permutation, &perm);
1033: for (v = 0; v < numValues; ++v) {
1034: const PetscInt size = (*labelNew)->stratumSizes[v];
1035: const PetscInt *points;
1036: PetscInt *pointsNew;
1038: ISGetIndices((*labelNew)->points[v],&points);
1039: PetscMalloc1(size,&pointsNew);
1040: for (q = 0; q < size; ++q) {
1041: const PetscInt point = points[q];
1043: 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);
1044: pointsNew[q] = perm[point];
1045: }
1046: ISRestoreIndices((*labelNew)->points[v],&points);
1047: PetscSortInt(size, pointsNew);
1048: ISDestroy(&((*labelNew)->points[v]));
1049: if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1050: ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1051: PetscFree(pointsNew);
1052: } else {
1053: ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1054: }
1055: PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1056: }
1057: ISRestoreIndices(permutation, &perm);
1058: if (label->bt) {
1059: PetscBTDestroy(&label->bt);
1060: DMLabelCreateIndex(label, label->pStart, label->pEnd);
1061: }
1062: return(0);
1063: }
1065: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1066: {
1067: MPI_Comm comm;
1068: PetscInt s, l, nroots, nleaves, dof, offset, size;
1069: PetscInt *remoteOffsets, *rootStrata, *rootIdx;
1070: PetscSection rootSection;
1071: PetscSF labelSF;
1075: if (label) {DMLabelMakeAllValid_Private(label);}
1076: PetscObjectGetComm((PetscObject)sf, &comm);
1077: /* Build a section of stratum values per point, generate the according SF
1078: and distribute point-wise stratum values to leaves. */
1079: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1080: PetscSectionCreate(comm, &rootSection);
1081: PetscSectionSetChart(rootSection, 0, nroots);
1082: if (label) {
1083: for (s = 0; s < label->numStrata; ++s) {
1084: const PetscInt *points;
1086: ISGetIndices(label->points[s], &points);
1087: for (l = 0; l < label->stratumSizes[s]; l++) {
1088: PetscSectionGetDof(rootSection, points[l], &dof);
1089: PetscSectionSetDof(rootSection, points[l], dof+1);
1090: }
1091: ISRestoreIndices(label->points[s], &points);
1092: }
1093: }
1094: PetscSectionSetUp(rootSection);
1095: /* Create a point-wise array of stratum values */
1096: PetscSectionGetStorageSize(rootSection, &size);
1097: PetscMalloc1(size, &rootStrata);
1098: PetscCalloc1(nroots, &rootIdx);
1099: if (label) {
1100: for (s = 0; s < label->numStrata; ++s) {
1101: const PetscInt *points;
1103: ISGetIndices(label->points[s], &points);
1104: for (l = 0; l < label->stratumSizes[s]; l++) {
1105: const PetscInt p = points[l];
1106: PetscSectionGetOffset(rootSection, p, &offset);
1107: rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1108: }
1109: ISRestoreIndices(label->points[s], &points);
1110: }
1111: }
1112: /* Build SF that maps label points to remote processes */
1113: PetscSectionCreate(comm, leafSection);
1114: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1115: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1116: PetscFree(remoteOffsets);
1117: /* Send the strata for each point over the derived SF */
1118: PetscSectionGetStorageSize(*leafSection, &size);
1119: PetscMalloc1(size, leafStrata);
1120: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1121: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1122: /* Clean up */
1123: PetscFree(rootStrata);
1124: PetscFree(rootIdx);
1125: PetscSectionDestroy(&rootSection);
1126: PetscSFDestroy(&labelSF);
1127: return(0);
1128: }
1130: /*@
1131: DMLabelDistribute - Create a new label pushed forward over the PetscSF
1133: Input Parameters:
1134: + label - the DMLabel
1135: - sf - the map from old to new distribution
1137: Output Parameter:
1138: . labelnew - the new redistributed label
1140: Level: intermediate
1142: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1143: @*/
1144: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1145: {
1146: MPI_Comm comm;
1147: PetscSection leafSection;
1148: PetscInt p, pStart, pEnd, s, size, dof, offset, stratum;
1149: PetscInt *leafStrata, *strataIdx;
1150: PetscInt **points;
1151: char *name;
1152: PetscInt nameSize;
1153: PetscHSetI stratumHash;
1154: size_t len = 0;
1155: PetscMPIInt rank;
1159: if (label) {DMLabelMakeAllValid_Private(label);}
1160: PetscObjectGetComm((PetscObject)sf, &comm);
1161: MPI_Comm_rank(comm, &rank);
1162: /* Bcast name */
1163: if (!rank) {PetscStrlen(label->name, &len);}
1164: nameSize = len;
1165: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1166: PetscMalloc1(nameSize+1, &name);
1167: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1168: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1169: DMLabelCreate(name, labelNew);
1170: PetscFree(name);
1171: /* Bcast defaultValue */
1172: if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1173: MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1174: /* Distribute stratum values over the SF and get the point mapping on the receiver */
1175: DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1176: /* Determine received stratum values and initialise new label*/
1177: PetscHSetICreate(&stratumHash);
1178: PetscSectionGetStorageSize(leafSection, &size);
1179: for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1180: PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1181: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1182: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1183: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1184: /* Turn leafStrata into indices rather than stratum values */
1185: offset = 0;
1186: PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1187: PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1188: for (p = 0; p < size; ++p) {
1189: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1190: if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1191: }
1192: }
1193: /* Rebuild the point strata on the receiver */
1194: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1195: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1196: for (p=pStart; p<pEnd; p++) {
1197: PetscSectionGetDof(leafSection, p, &dof);
1198: PetscSectionGetOffset(leafSection, p, &offset);
1199: for (s=0; s<dof; s++) {
1200: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1201: }
1202: }
1203: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1204: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1205: PetscMalloc1((*labelNew)->numStrata,&points);
1206: for (s = 0; s < (*labelNew)->numStrata; ++s) {
1207: PetscHSetICreate(&(*labelNew)->ht[s]);
1208: PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1209: }
1210: /* Insert points into new strata */
1211: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1212: PetscSectionGetChart(leafSection, &pStart, &pEnd);
1213: for (p=pStart; p<pEnd; p++) {
1214: PetscSectionGetDof(leafSection, p, &dof);
1215: PetscSectionGetOffset(leafSection, p, &offset);
1216: for (s=0; s<dof; s++) {
1217: stratum = leafStrata[offset+s];
1218: points[stratum][strataIdx[stratum]++] = p;
1219: }
1220: }
1221: for (s = 0; s < (*labelNew)->numStrata; s++) {
1222: ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1223: PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1224: }
1225: PetscFree(points);
1226: PetscHSetIDestroy(&stratumHash);
1227: PetscFree(leafStrata);
1228: PetscFree(strataIdx);
1229: PetscSectionDestroy(&leafSection);
1230: return(0);
1231: }
1233: /*@
1234: DMLabelGather - Gather all label values from leafs into roots
1236: Input Parameters:
1237: + label - the DMLabel
1238: - sf - the Star Forest point communication map
1240: Output Parameters:
1241: . labelNew - the new DMLabel with localised leaf values
1243: Level: developer
1245: Note: This is the inverse operation to DMLabelDistribute.
1247: .seealso: DMLabelDistribute()
1248: @*/
1249: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1250: {
1251: MPI_Comm comm;
1252: PetscSection rootSection;
1253: PetscSF sfLabel;
1254: PetscSFNode *rootPoints, *leafPoints;
1255: PetscInt p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1256: const PetscInt *rootDegree, *ilocal;
1257: PetscInt *rootStrata;
1258: char *name;
1259: PetscInt nameSize;
1260: size_t len = 0;
1261: PetscMPIInt rank, size;
1265: PetscObjectGetComm((PetscObject)sf, &comm);
1266: MPI_Comm_rank(comm, &rank);
1267: MPI_Comm_size(comm, &size);
1268: /* Bcast name */
1269: if (!rank) {PetscStrlen(label->name, &len);}
1270: nameSize = len;
1271: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1272: PetscMalloc1(nameSize+1, &name);
1273: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1274: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1275: DMLabelCreate(name, labelNew);
1276: PetscFree(name);
1277: /* Gather rank/index pairs of leaves into local roots to build
1278: an inverse, multi-rooted SF. Note that this ignores local leaf
1279: indexing due to the use of the multiSF in PetscSFGather. */
1280: PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1281: PetscMalloc1(nroots, &leafPoints);
1282: for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1283: for (p = 0; p < nleaves; p++) {
1284: leafPoints[ilocal[p]].index = ilocal[p];
1285: leafPoints[ilocal[p]].rank = rank;
1286: }
1287: PetscSFComputeDegreeBegin(sf, &rootDegree);
1288: PetscSFComputeDegreeEnd(sf, &rootDegree);
1289: for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1290: PetscMalloc1(nmultiroots, &rootPoints);
1291: PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1292: PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1293: PetscSFCreate(comm,& sfLabel);
1294: PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1295: /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1296: DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1297: /* Rebuild the point strata on the receiver */
1298: for (p = 0, idx = 0; p < nroots; p++) {
1299: for (d = 0; d < rootDegree[p]; d++) {
1300: PetscSectionGetDof(rootSection, idx+d, &dof);
1301: PetscSectionGetOffset(rootSection, idx+d, &offset);
1302: for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1303: }
1304: idx += rootDegree[p];
1305: }
1306: PetscFree(leafPoints);
1307: PetscFree(rootStrata);
1308: PetscSectionDestroy(&rootSection);
1309: PetscSFDestroy(&sfLabel);
1310: return(0);
1311: }
1313: /*@
1314: DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label
1316: Input Parameter:
1317: . label - the DMLabel
1319: Output Parameters:
1320: + section - the section giving offsets for each stratum
1321: - is - An IS containing all the label points
1323: Level: developer
1325: .seealso: DMLabelDistribute()
1326: @*/
1327: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1328: {
1329: IS vIS;
1330: const PetscInt *values;
1331: PetscInt *points;
1332: PetscInt nV, vS = 0, vE = 0, v, N;
1333: PetscErrorCode ierr;
1336: DMLabelGetNumValues(label, &nV);
1337: DMLabelGetValueIS(label, &vIS);
1338: ISGetIndices(vIS, &values);
1339: if (nV) {vS = values[0]; vE = values[0]+1;}
1340: for (v = 1; v < nV; ++v) {
1341: vS = PetscMin(vS, values[v]);
1342: vE = PetscMax(vE, values[v]+1);
1343: }
1344: PetscSectionCreate(PETSC_COMM_SELF, section);
1345: PetscSectionSetChart(*section, vS, vE);
1346: for (v = 0; v < nV; ++v) {
1347: PetscInt n;
1349: DMLabelGetStratumSize(label, values[v], &n);
1350: PetscSectionSetDof(*section, values[v], n);
1351: }
1352: PetscSectionSetUp(*section);
1353: PetscSectionGetStorageSize(*section, &N);
1354: PetscMalloc1(N, &points);
1355: for (v = 0; v < nV; ++v) {
1356: IS is;
1357: const PetscInt *spoints;
1358: PetscInt dof, off, p;
1360: PetscSectionGetDof(*section, values[v], &dof);
1361: PetscSectionGetOffset(*section, values[v], &off);
1362: DMLabelGetStratumIS(label, values[v], &is);
1363: ISGetIndices(is, &spoints);
1364: for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1365: ISRestoreIndices(is, &spoints);
1366: ISDestroy(&is);
1367: }
1368: ISRestoreIndices(vIS, &values);
1369: ISDestroy(&vIS);
1370: ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1371: return(0);
1372: }
1374: /*@
1375: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1376: the local section and an SF describing the section point overlap.
1378: Input Parameters:
1379: + s - The PetscSection for the local field layout
1380: . sf - The SF describing parallel layout of the section points
1381: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1382: . label - The label specifying the points
1383: - labelValue - The label stratum specifying the points
1385: Output Parameter:
1386: . gsection - The PetscSection for the global field layout
1388: Note: This gives negative sizes and offsets to points not owned by this process
1390: Level: developer
1392: .seealso: PetscSectionCreate()
1393: @*/
1394: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1395: {
1396: PetscInt *neg = NULL, *tmpOff = NULL;
1397: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
1401: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1402: PetscSectionGetChart(s, &pStart, &pEnd);
1403: PetscSectionSetChart(*gsection, pStart, pEnd);
1404: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1405: if (nroots >= 0) {
1406: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1407: PetscCalloc1(nroots, &neg);
1408: if (nroots > pEnd-pStart) {
1409: PetscCalloc1(nroots, &tmpOff);
1410: } else {
1411: tmpOff = &(*gsection)->atlasDof[-pStart];
1412: }
1413: }
1414: /* Mark ghost points with negative dof */
1415: for (p = pStart; p < pEnd; ++p) {
1416: PetscInt value;
1418: DMLabelGetValue(label, p, &value);
1419: if (value != labelValue) continue;
1420: PetscSectionGetDof(s, p, &dof);
1421: PetscSectionSetDof(*gsection, p, dof);
1422: PetscSectionGetConstraintDof(s, p, &cdof);
1423: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1424: if (neg) neg[p] = -(dof+1);
1425: }
1426: PetscSectionSetUpBC(*gsection);
1427: if (nroots >= 0) {
1428: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1429: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1430: if (nroots > pEnd-pStart) {
1431: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1432: }
1433: }
1434: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1435: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1436: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1437: (*gsection)->atlasOff[p] = off;
1438: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1439: }
1440: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1441: globalOff -= off;
1442: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1443: (*gsection)->atlasOff[p] += globalOff;
1444: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1445: }
1446: /* Put in negative offsets for ghost points */
1447: if (nroots >= 0) {
1448: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1449: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1450: if (nroots > pEnd-pStart) {
1451: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1452: }
1453: }
1454: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1455: PetscFree(neg);
1456: return(0);
1457: }
1459: typedef struct _n_PetscSectionSym_Label
1460: {
1461: DMLabel label;
1462: PetscCopyMode *modes;
1463: PetscInt *sizes;
1464: const PetscInt ***perms;
1465: const PetscScalar ***rots;
1466: PetscInt (*minMaxOrients)[2];
1467: PetscInt numStrata; /* numStrata is only increasing, functions as a state */
1468: } PetscSectionSym_Label;
1470: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1471: {
1472: PetscInt i, j;
1473: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1474: PetscErrorCode ierr;
1477: for (i = 0; i <= sl->numStrata; i++) {
1478: if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1479: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1480: if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1481: if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1482: }
1483: if (sl->perms[i]) {
1484: const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];
1486: PetscFree(perms);
1487: }
1488: if (sl->rots[i]) {
1489: const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];
1491: PetscFree(rots);
1492: }
1493: }
1494: }
1495: PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1496: DMLabelDestroy(&sl->label);
1497: sl->numStrata = 0;
1498: return(0);
1499: }
1501: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1502: {
1506: PetscSectionSymLabelReset(sym);
1507: PetscFree(sym->data);
1508: return(0);
1509: }
1511: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1512: {
1513: PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1514: PetscBool isAscii;
1515: DMLabel label = sl->label;
1516: PetscErrorCode ierr;
1519: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1520: if (isAscii) {
1521: PetscInt i, j, k;
1522: PetscViewerFormat format;
1524: PetscViewerGetFormat(viewer,&format);
1525: if (label) {
1526: PetscViewerGetFormat(viewer,&format);
1527: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1528: PetscViewerASCIIPushTab(viewer);
1529: DMLabelView(label, viewer);
1530: PetscViewerASCIIPopTab(viewer);
1531: } else {
1532: PetscViewerASCIIPrintf(viewer," Label '%s'\n",sl->label->name);
1533: }
1534: } else {
1535: PetscViewerASCIIPrintf(viewer, "No label given\n");
1536: }
1537: PetscViewerASCIIPushTab(viewer);
1538: for (i = 0; i <= sl->numStrata; i++) {
1539: PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;
1541: if (!(sl->perms[i] || sl->rots[i])) {
1542: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1543: } else {
1544: PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1545: PetscViewerASCIIPushTab(viewer);
1546: PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1547: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1548: PetscViewerASCIIPushTab(viewer);
1549: for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1550: if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1551: PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1552: } else {
1553: PetscInt tab;
1555: PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1556: PetscViewerASCIIPushTab(viewer);
1557: PetscViewerASCIIGetTab(viewer,&tab);
1558: if (sl->perms[i] && sl->perms[i][j]) {
1559: PetscViewerASCIIPrintf(viewer,"Permutation:");
1560: PetscViewerASCIISetTab(viewer,0);
1561: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1562: PetscViewerASCIIPrintf(viewer,"\n");
1563: PetscViewerASCIISetTab(viewer,tab);
1564: }
1565: if (sl->rots[i] && sl->rots[i][j]) {
1566: PetscViewerASCIIPrintf(viewer,"Rotations: ");
1567: PetscViewerASCIISetTab(viewer,0);
1568: #if defined(PETSC_USE_COMPLEX)
1569: 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]));}
1570: #else
1571: for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1572: #endif
1573: PetscViewerASCIIPrintf(viewer,"\n");
1574: PetscViewerASCIISetTab(viewer,tab);
1575: }
1576: PetscViewerASCIIPopTab(viewer);
1577: }
1578: }
1579: PetscViewerASCIIPopTab(viewer);
1580: }
1581: PetscViewerASCIIPopTab(viewer);
1582: }
1583: }
1584: PetscViewerASCIIPopTab(viewer);
1585: }
1586: return(0);
1587: }
1589: /*@
1590: PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries
1592: Logically collective on sym
1594: Input parameters:
1595: + sym - the section symmetries
1596: - label - the DMLabel describing the types of points
1598: Level: developer:
1600: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1601: @*/
1602: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1603: {
1604: PetscSectionSym_Label *sl;
1609: sl = (PetscSectionSym_Label *) sym->data;
1610: if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1611: if (label) {
1612: label->refct++;
1613: sl->label = label;
1614: DMLabelGetNumValues(label,&sl->numStrata);
1615: 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);
1616: PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1617: PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1618: PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1619: PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1620: PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1621: }
1622: return(0);
1623: }
1625: /*@C
1626: PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum
1628: Logically collective on PetscSectionSym
1630: InputParameters:
1631: + sys - the section symmetries
1632: . stratum - the stratum value in the label that we are assigning symmetries for
1633: . size - the number of dofs for points in the stratum of the label
1634: . minOrient - the smallest orientation for a point in this stratum
1635: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1636: . mode - how sym should copy the perms and rots arrays
1637: . perms - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation. A NULL permutation is the identity
1638: + 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
1640: Level: developer
1642: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1643: @*/
1644: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1645: {
1646: PetscInt i, j, k;
1647: PetscSectionSym_Label *sl;
1652: sl = (PetscSectionSym_Label *) sym->data;
1653: if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1654: for (i = 0; i <= sl->numStrata; i++) {
1655: PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;
1657: if (stratum == value) break;
1658: }
1659: if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,sl->label->name);
1660: sl->sizes[i] = size;
1661: sl->modes[i] = mode;
1662: sl->minMaxOrients[i][0] = minOrient;
1663: sl->minMaxOrients[i][1] = maxOrient;
1664: if (mode == PETSC_COPY_VALUES) {
1665: if (perms) {
1666: PetscInt **ownPerms;
1668: PetscCalloc1(maxOrient - minOrient,&ownPerms);
1669: for (j = 0; j < maxOrient-minOrient; j++) {
1670: if (perms[j]) {
1671: PetscMalloc1(size,&ownPerms[j]);
1672: for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
1673: }
1674: }
1675: sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
1676: }
1677: if (rots) {
1678: PetscScalar **ownRots;
1680: PetscCalloc1(maxOrient - minOrient,&ownRots);
1681: for (j = 0; j < maxOrient-minOrient; j++) {
1682: if (rots[j]) {
1683: PetscMalloc1(size,&ownRots[j]);
1684: for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
1685: }
1686: }
1687: sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
1688: }
1689: } else {
1690: sl->perms[i] = perms ? &perms[-minOrient] : NULL;
1691: sl->rots[i] = rots ? &rots[-minOrient] : NULL;
1692: }
1693: return(0);
1694: }
1696: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
1697: {
1698: PetscInt i, j, numStrata;
1699: PetscSectionSym_Label *sl;
1700: DMLabel label;
1701: PetscErrorCode ierr;
1704: sl = (PetscSectionSym_Label *) sym->data;
1705: numStrata = sl->numStrata;
1706: label = sl->label;
1707: for (i = 0; i < numPoints; i++) {
1708: PetscInt point = points[2*i];
1709: PetscInt ornt = points[2*i+1];
1711: for (j = 0; j < numStrata; j++) {
1712: if (label->validIS[j]) {
1713: PetscInt k;
1715: ISLocate(label->points[j],point,&k);
1716: if (k >= 0) break;
1717: } else {
1718: PetscBool has;
1720: PetscHSetIHas(label->ht[j], point, &has);
1721: if (has) break;
1722: }
1723: }
1724: 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);
1725: if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
1726: if (rots) {rots[i] = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
1727: }
1728: return(0);
1729: }
1731: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
1732: {
1733: PetscSectionSym_Label *sl;
1734: PetscErrorCode ierr;
1737: PetscNewLog(sym,&sl);
1738: sym->ops->getpoints = PetscSectionSymGetPoints_Label;
1739: sym->ops->view = PetscSectionSymView_Label;
1740: sym->ops->destroy = PetscSectionSymDestroy_Label;
1741: sym->data = (void *) sl;
1742: return(0);
1743: }
1745: /*@
1746: PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label
1748: Collective
1750: Input Parameters:
1751: + comm - the MPI communicator for the new symmetry
1752: - label - the label defining the strata
1754: Output Parameters:
1755: . sym - the section symmetries
1757: Level: developer
1759: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
1760: @*/
1761: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
1762: {
1763: PetscErrorCode ierr;
1766: DMInitializePackage();
1767: PetscSectionSymCreate(comm,sym);
1768: PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
1769: PetscSectionSymLabelSetLabel(*sym,label);
1770: return(0);
1771: }