Actual source code: plexlabel.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petscsf.h>
6: PetscErrorCode DMLabelCreate(const char name[], DMLabel *label)
7: {
11: PetscNew(label);
12: PetscStrallocpy(name, &(*label)->name);
14: (*label)->refct = 1;
15: (*label)->numStrata = 0;
16: (*label)->stratumValues = NULL;
17: (*label)->arrayValid = NULL;
18: (*label)->stratumSizes = NULL;
19: (*label)->points = NULL;
20: (*label)->ht = NULL;
21: (*label)->pStart = -1;
22: (*label)->pEnd = -1;
23: (*label)->bt = NULL;
24: return(0);
25: }
29: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
30: {
31: PetscInt off;
34: if (label->arrayValid[v]) return 0;
35: if (v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %d in DMLabelMakeValid_Private\n", v);
37: PetscHashISize(label->ht[v], label->stratumSizes[v]);
39: PetscMalloc1(label->stratumSizes[v], &label->points[v]);
40: off = 0;
41: PetscHashIGetKeys(label->ht[v], &off, &(label->points[v][0]));
42: 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]);
43: PetscHashIClear(label->ht[v]);
44: PetscSortInt(label->stratumSizes[v], label->points[v]);
45: if (label->bt) {
46: PetscInt p;
48: for (p = 0; p < label->stratumSizes[v]; ++p) {
49: const PetscInt point = label->points[v][p];
51: 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);
52: PetscBTSet(label->bt, point - label->pStart);
53: }
54: }
55: label->arrayValid[v] = PETSC_TRUE;
56: return(0);
57: }
61: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
62: {
63: PetscInt v;
67: for (v = 0; v < label->numStrata; v++){
68: DMLabelMakeValid_Private(label, v);
69: }
70: return(0);
71: }
75: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
76: {
77: PETSC_UNUSED PetscHashIIter ret, iter;
78: PetscInt p;
82: if (!label->arrayValid[v]) return(0);
83: for (p = 0; p < label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], label->points[v][p], ret, iter);
84: PetscFree(label->points[v]);
85: label->arrayValid[v] = PETSC_FALSE;
86: return(0);
87: }
91: static PetscErrorCode DMLabelAddStratum_Private(DMLabel label, PetscInt value)
92: {
93: PetscInt v, *tmpV, *tmpS, **tmpP;
94: PetscHashI *tmpH;
95: PetscBool *tmpB;
100: PetscMalloc1((label->numStrata+1), &tmpV);
101: PetscMalloc1((label->numStrata+1), &tmpS);
102: PetscMalloc1((label->numStrata+1), &tmpH);
103: PetscMalloc1((label->numStrata+1), &tmpP);
104: PetscMalloc1((label->numStrata+1), &tmpB);
105: for (v = 0; v < label->numStrata; ++v) {
106: tmpV[v] = label->stratumValues[v];
107: tmpS[v] = label->stratumSizes[v];
108: tmpH[v] = label->ht[v];
109: tmpP[v] = label->points[v];
110: tmpB[v] = label->arrayValid[v];
111: }
112: tmpV[v] = value;
113: tmpS[v] = 0;
114: PetscHashICreate(tmpH[v]);
115: tmpP[v] = NULL;
116: tmpB[v] = PETSC_TRUE;
117: ++label->numStrata;
118: PetscFree(label->stratumValues);
119: PetscFree(label->stratumSizes);
120: PetscFree(label->ht);
121: PetscFree(label->points);
122: PetscFree(label->arrayValid);
123: label->stratumValues = tmpV;
124: label->stratumSizes = tmpS;
125: label->ht = tmpH;
126: label->points = tmpP;
127: label->arrayValid = tmpB;
129: return(0);
130: }
134: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
135: {
138: *name = label->name;
139: return(0);
140: }
144: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
145: {
146: PetscInt v;
147: PetscMPIInt rank;
151: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
152: if (label) {
153: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
154: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%d, %d)\n", label->pStart, label->pEnd);}
155: for (v = 0; v < label->numStrata; ++v) {
156: const PetscInt value = label->stratumValues[v];
157: PetscInt p;
159: for (p = 0; p < label->stratumSizes[v]; ++p) {
160: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D (%D)\n", rank, label->points[v][p], value);
161: }
162: }
163: }
164: PetscViewerFlush(viewer);
165: return(0);
166: }
170: /*@C
171: DMLabelView - View the label
173: Input Parameters:
174: + label - The DMLabel
175: - viewer - The PetscViewer
177: Level: intermediate
179: .seealso: DMLabelCreate(), DMLabelDestroy()
180: @*/
181: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
182: {
183: PetscBool iascii;
188: if (label) {DMLabelMakeAllValid_Private(label);}
189: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
190: if (iascii) {
191: DMLabelView_Ascii(label, viewer);
192: } else SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not supported by this mesh object", ((PetscObject) viewer)->type_name);
193: return(0);
194: }
198: PetscErrorCode DMLabelDestroy(DMLabel *label)
199: {
200: PetscInt v;
204: if (!(*label)) return(0);
205: if (--(*label)->refct > 0) return(0);
206: PetscFree((*label)->name);
207: PetscFree((*label)->stratumValues);
208: PetscFree((*label)->stratumSizes);
209: for (v = 0; v < (*label)->numStrata; ++v) {PetscFree((*label)->points[v]);}
210: PetscFree((*label)->points);
211: PetscFree((*label)->arrayValid);
212: if ((*label)->ht) {
213: for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
214: PetscFree((*label)->ht);
215: }
216: PetscBTDestroy(&(*label)->bt);
217: PetscFree(*label);
218: return(0);
219: }
223: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
224: {
225: PetscInt v, q;
229: DMLabelMakeAllValid_Private(label);
230: PetscNew(labelnew);
231: PetscStrallocpy(label->name, &(*labelnew)->name);
233: (*labelnew)->refct = 1;
234: (*labelnew)->numStrata = label->numStrata;
235: if (label->numStrata) {
236: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
237: PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
238: PetscMalloc1(label->numStrata, &(*labelnew)->ht);
239: PetscMalloc1(label->numStrata, &(*labelnew)->points);
240: PetscMalloc1(label->numStrata, &(*labelnew)->arrayValid);
241: /* Could eliminate unused space here */
242: for (v = 0; v < label->numStrata; ++v) {
243: PetscMalloc1(label->stratumSizes[v], &(*labelnew)->points[v]);
244: PetscHashICreate((*labelnew)->ht[v]);
245: (*labelnew)->arrayValid[v] = PETSC_TRUE;
246: (*labelnew)->stratumValues[v] = label->stratumValues[v];
247: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
248: for (q = 0; q < label->stratumSizes[v]; ++q) {
249: (*labelnew)->points[v][q] = label->points[v][q];
250: }
251: }
252: }
253: (*labelnew)->pStart = -1;
254: (*labelnew)->pEnd = -1;
255: (*labelnew)->bt = NULL;
256: return(0);
257: }
261: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
262: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
263: {
264: PetscInt v;
268: DMLabelMakeAllValid_Private(label);
269: if (label->bt) {PetscBTDestroy(&label->bt);}
270: label->pStart = pStart;
271: label->pEnd = pEnd;
272: PetscBTCreate(pEnd - pStart, &label->bt);
273: PetscBTMemzero(pEnd - pStart, label->bt);
274: for (v = 0; v < label->numStrata; ++v) {
275: PetscInt i;
277: for (i = 0; i < label->stratumSizes[v]; ++i) {
278: const PetscInt point = label->points[v][i];
280: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, pStart, pEnd);
281: PetscBTSet(label->bt, point - pStart);
282: }
283: }
284: return(0);
285: }
289: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
290: {
294: label->pStart = -1;
295: label->pEnd = -1;
296: if (label->bt) {PetscBTDestroy(&label->bt);}
297: return(0);
298: }
302: /*@
303: DMLabelHasValue - Determine whether a label assigns the value to any point
305: Input Parameters:
306: + label - the DMLabel
307: - value - the value
309: Output Parameter:
310: . contains - Flag indicating whether the label maps this value to any point
312: Level: developer
314: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
315: @*/
316: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
317: {
318: PetscInt v;
322: for (v = 0; v < label->numStrata; ++v) {
323: if (value == label->stratumValues[v]) break;
324: }
325: *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
326: return(0);
327: }
331: /*@
332: DMLabelHasPoint - Determine whether a label assigns a value to a point
334: Input Parameters:
335: + label - the DMLabel
336: - point - the point
338: Output Parameter:
339: . contains - Flag indicating whether the label maps this point to a value
341: Note: The user must call DMLabelCreateIndex() before this function.
343: Level: developer
345: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
346: @*/
347: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
348: {
353: DMLabelMakeAllValid_Private(label);
354: #if defined(PETSC_USE_DEBUG)
355: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
356: 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);
357: #endif
358: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
359: return(0);
360: }
364: /*@
365: DMLabelStratumHasPoint - Return true if the stratum contains a point
367: Input Parameters:
368: + label - the DMLabel
369: . value - the stratum value
370: - point - the point
372: Output Parameter:
373: . contains - true if the stratum contains the point
375: Level: intermediate
377: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
378: @*/
379: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
380: {
381: PetscInt v;
386: *contains = PETSC_FALSE;
387: for (v = 0; v < label->numStrata; ++v) {
388: if (label->stratumValues[v] == value) {
389: if (label->arrayValid[v]) {
390: PetscInt i;
392: PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
393: if (i >= 0) {
394: *contains = PETSC_TRUE;
395: break;
396: }
397: } else {
398: PetscBool has;
400: PetscHashIHasKey(label->ht[v], point, has);
401: if (has) {
402: *contains = PETSC_TRUE;
403: break;
404: }
405: }
406: }
407: }
408: return(0);
409: }
413: /*@
414: DMLabelGetValue - Return the value a label assigns to a point, or -1
416: Input Parameters:
417: + label - the DMLabel
418: - point - the point
420: Output Parameter:
421: . value - The point value, or -1
423: Level: intermediate
425: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
426: @*/
427: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
428: {
429: PetscInt v;
434: *value = -1;
435: for (v = 0; v < label->numStrata; ++v) {
436: if (label->arrayValid[v]) {
437: PetscInt i;
439: PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
440: if (i >= 0) {
441: *value = label->stratumValues[v];
442: break;
443: }
444: } else {
445: PetscBool has;
447: PetscHashIHasKey(label->ht[v], point, has);
448: if (has) {
449: *value = label->stratumValues[v];
450: break;
451: }
452: }
453: }
454: return(0);
455: }
459: /*@
460: DMLabelSetValue - Set the value a label assigns to a point
462: Input Parameters:
463: + label - the DMLabel
464: . point - the point
465: - value - The point value
467: Level: intermediate
469: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue()
470: @*/
471: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
472: {
473: PETSC_UNUSED PetscHashIIter iter, ret;
474: PetscInt v;
475: PetscErrorCode ierr;
478: /* Find, or add, label value */
479: for (v = 0; v < label->numStrata; ++v) {
480: if (label->stratumValues[v] == value) break;
481: }
482: /* Create new table */
483: if (v >= label->numStrata) {
484: DMLabelAddStratum_Private(label, value);
485: }
486: DMLabelMakeInvalid_Private(label, v);
487: /* Set key */
488: PetscHashIPut(label->ht[v], point, ret, iter);
489: return(0);
490: }
494: /*@
495: DMLabelClearValue - Clear the value a label assigns to a point
497: Input Parameters:
498: + label - the DMLabel
499: . point - the point
500: - value - The point value
502: Level: intermediate
504: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
505: @*/
506: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
507: {
508: PetscInt v, p;
512: /* Find label value */
513: for (v = 0; v < label->numStrata; ++v) {
514: if (label->stratumValues[v] == value) break;
515: }
516: if (v >= label->numStrata) return(0);
517: if (label->arrayValid[v]) {
518: /* Check whether point exists */
519: PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &p);
520: if (p >= 0) {
521: PetscMemmove(&label->points[v][p], &label->points[v][p+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
522: --label->stratumSizes[v];
523: if (label->bt) {
524: 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);
525: PetscBTClear(label->bt, point - label->pStart);
526: }
527: }
528: } else {
529: PetscHashIDelKey(label->ht[v], point);
530: }
531: return(0);
532: }
536: /*@
537: DMLabelInsertIS - Set all points in the IS to a value
539: Input Parameters:
540: + label - the DMLabel
541: . is - the point IS
542: - value - The point value
544: Level: intermediate
546: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
547: @*/
548: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
549: {
550: const PetscInt *points;
551: PetscInt n, p;
552: PetscErrorCode ierr;
556: ISGetLocalSize(is, &n);
557: ISGetIndices(is, &points);
558: for (p = 0; p < n; ++p) {DMLabelSetValue(label, points[p], value);}
559: ISRestoreIndices(is, &points);
560: return(0);
561: }
565: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
566: {
569: *numValues = label->numStrata;
570: return(0);
571: }
575: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
576: {
581: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
582: return(0);
583: }
587: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
588: {
589: PetscInt v;
594: *size = 0;
595: for (v = 0; v < label->numStrata; ++v) {
596: if (label->stratumValues[v] == value) {
597: DMLabelMakeValid_Private(label, v);
598: *size = label->stratumSizes[v];
599: break;
600: }
601: }
602: return(0);
603: }
607: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
608: {
609: PetscInt v;
615: for (v = 0; v < label->numStrata; ++v) {
616: if (label->stratumValues[v] != value) continue;
617: DMLabelMakeValid_Private(label, v);
618: if (start) *start = label->points[v][0];
619: if (end) *end = label->points[v][label->stratumSizes[v]-1]+1;
620: break;
621: }
622: return(0);
623: }
627: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
628: {
629: PetscInt v;
634: *points = NULL;
635: for (v = 0; v < label->numStrata; ++v) {
636: if (label->stratumValues[v] == value) {
637: DMLabelMakeValid_Private(label, v);
638: if (label->arrayValid[v]) {
639: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[v][0], PETSC_COPY_VALUES, points);
640: PetscObjectSetName((PetscObject) *points, "indices");
641: } else {
642: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
643: }
644: break;
645: }
646: }
647: return(0);
648: }
652: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
653: {
654: PetscInt v;
658: for (v = 0; v < label->numStrata; ++v) {
659: if (label->stratumValues[v] == value) break;
660: }
661: if (v >= label->numStrata) return(0);
662: if (label->bt) {
663: PetscInt i;
665: for (i = 0; i < label->stratumSizes[v]; ++i) {
666: const PetscInt point = label->points[v][i];
668: 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);
669: PetscBTClear(label->bt, point - label->pStart);
670: }
671: }
672: if (label->arrayValid[v]) {
673: label->stratumSizes[v] = 0;
674: } else {
675: PetscHashIClear(label->ht[v]);
676: }
677: return(0);
678: }
682: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
683: {
684: PetscInt v;
688: DMLabelMakeAllValid_Private(label);
689: label->pStart = start;
690: label->pEnd = end;
691: if (label->bt) {PetscBTDestroy(&label->bt);}
692: /* Could squish offsets, but would only make sense if I reallocate the storage */
693: for (v = 0; v < label->numStrata; ++v) {
694: PetscInt off, q;
696: for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
697: const PetscInt point = label->points[v][q];
699: if ((point < start) || (point >= end)) continue;
700: label->points[v][off++] = point;
701: }
702: label->stratumSizes[v] = off;
703: }
704: DMLabelCreateIndex(label, start, end);
705: return(0);
706: }
710: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
711: {
712: const PetscInt *perm;
713: PetscInt numValues, numPoints, v, q;
714: PetscErrorCode ierr;
717: DMLabelMakeAllValid_Private(label);
718: DMLabelDuplicate(label, labelNew);
719: DMLabelGetNumValues(*labelNew, &numValues);
720: ISGetLocalSize(permutation, &numPoints);
721: ISGetIndices(permutation, &perm);
722: for (v = 0; v < numValues; ++v) {
723: const PetscInt size = (*labelNew)->stratumSizes[v];
725: for (q = 0; q < size; ++q) {
726: const PetscInt point = (*labelNew)->points[v][q];
728: 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);
729: (*labelNew)->points[v][q] = perm[point];
730: }
731: PetscSortInt(size, &(*labelNew)->points[v][0]);
732: }
733: ISRestoreIndices(permutation, &perm);
734: if (label->bt) {
735: PetscBTDestroy(&label->bt);
736: DMLabelCreateIndex(label, label->pStart, label->pEnd);
737: }
738: return(0);
739: }
743: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
744: {
745: MPI_Comm comm;
746: PetscSection rootSection, leafSection;
747: PetscSF labelSF;
748: PetscInt p, pStart, pEnd, l, lStart, lEnd, s, nroots, nleaves, size, dof, offset, stratum;
749: PetscInt *remoteOffsets, *rootStrata, *rootIdx, *leafStrata, *strataIdx;
750: char *name;
751: PetscInt nameSize;
752: size_t len = 0;
753: PetscMPIInt rank, numProcs;
757: if (label) {DMLabelMakeAllValid_Private(label);}
758: PetscObjectGetComm((PetscObject)sf, &comm);
759: MPI_Comm_rank(comm, &rank);
760: MPI_Comm_size(comm, &numProcs);
761: /* Bcast name */
762: if (!rank) {PetscStrlen(label->name, &len);}
763: nameSize = len;
764: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
765: PetscMalloc1(nameSize+1, &name);
766: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
767: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
768: DMLabelCreate(name, labelNew);
769: PetscFree(name);
770: /* Bcast numStrata */
771: if (!rank) (*labelNew)->numStrata = label->numStrata;
772: MPI_Bcast(&(*labelNew)->numStrata, 1, MPIU_INT, 0, comm);
773: /* Bcast stratumValues */
774: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
775: if (!rank) {PetscMemcpy((*labelNew)->stratumValues, label->stratumValues, label->numStrata * sizeof(PetscInt));}
776: MPI_Bcast((*labelNew)->stratumValues, (*labelNew)->numStrata, MPIU_INT, 0, comm);
777: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->arrayValid);
778: for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->arrayValid[s] = PETSC_TRUE;
780: /* Build a section detailing strata-per-point, distribute and build SF
781: from that and then send our points. */
782: PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
783: PetscSectionCreate(comm, &rootSection);
784: PetscSectionSetChart(rootSection, 0, nroots);
785: if (label) {
786: for (s = 0; s < label->numStrata; ++s) {
787: lStart = 0;
788: lEnd = label->stratumSizes[s];
789: for (l=lStart; l<lEnd; l++) {
790: PetscSectionGetDof(rootSection, label->points[s][l], &dof);
791: PetscSectionSetDof(rootSection, label->points[s][l], dof+1);
792: }
793: }
794: }
795: PetscSectionSetUp(rootSection);
797: /* Create a point-wise array of point strata */
798: PetscSectionGetStorageSize(rootSection, &size);
799: PetscMalloc1(size, &rootStrata);
800: PetscCalloc1(nroots, &rootIdx);
801: if (label) {
802: for (s = 0; s < label->numStrata; ++s) {
803: lStart = 0;
804: lEnd = label->stratumSizes[s];
805: for (l=lStart; l<lEnd; l++) {
806: p = label->points[s][l];
807: PetscSectionGetOffset(rootSection, p, &offset);
808: rootStrata[offset+rootIdx[p]++] = s;
809: }
810: }
811: }
813: /* Build SF that maps label points to remote processes */
814: PetscSectionCreate(comm, &leafSection);
815: PetscSFDistributeSection(sf, rootSection, &remoteOffsets, leafSection);
816: PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, leafSection, &labelSF);
818: /* Send the strata for each point over the derived SF */
819: PetscSectionGetStorageSize(leafSection, &size);
820: PetscMalloc1(size, &leafStrata);
821: PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, leafStrata);
822: PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, leafStrata);
824: /* Rebuild the point strata on the receiver */
825: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
826: PetscSectionGetChart(leafSection, &pStart, &pEnd);
827: for (p=pStart; p<pEnd; p++) {
828: PetscSectionGetDof(leafSection, p, &dof);
829: PetscSectionGetOffset(leafSection, p, &offset);
830: for (s=0; s<dof; s++) {
831: (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
832: }
833: }
834: PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
835: PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
836: for (s = 0; s < (*labelNew)->numStrata; ++s) {
837: PetscHashICreate((*labelNew)->ht[s]);
838: PetscMalloc1((*labelNew)->stratumSizes[s], &(*labelNew)->points[s]);
839: }
841: /* Insert points into new strata */
842: PetscCalloc1((*labelNew)->numStrata, &strataIdx);
843: PetscSectionGetChart(leafSection, &pStart, &pEnd);
844: for (p=pStart; p<pEnd; p++) {
845: PetscSectionGetDof(leafSection, p, &dof);
846: PetscSectionGetOffset(leafSection, p, &offset);
847: for (s=0; s<dof; s++) {
848: stratum = leafStrata[offset+s];
849: (*labelNew)->points[stratum][strataIdx[stratum]++] = p;
850: }
851: }
852: PetscFree(rootStrata);
853: PetscFree(leafStrata);
854: PetscFree(rootIdx);
855: PetscFree(strataIdx);
856: PetscSectionDestroy(&rootSection);
857: PetscSectionDestroy(&leafSection);
858: PetscSFDestroy(&labelSF);
859: return(0);
860: }
865: /*@C
866: DMPlexCreateLabel - Create a label of the given name if it does not already exist
868: Not Collective
870: Input Parameters:
871: + dm - The DMPlex object
872: - name - The label name
874: Level: intermediate
876: .keywords: mesh
877: .seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
878: @*/
879: PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
880: {
881: DM_Plex *mesh = (DM_Plex*) dm->data;
882: PlexLabel next = mesh->labels;
883: PetscBool flg = PETSC_FALSE;
889: while (next) {
890: PetscStrcmp(name, next->label->name, &flg);
891: if (flg) break;
892: next = next->next;
893: }
894: if (!flg) {
895: PlexLabel tmpLabel;
897: PetscCalloc1(1, &tmpLabel);
898: DMLabelCreate(name, &tmpLabel->label);
899: tmpLabel->output = PETSC_TRUE;
900: tmpLabel->next = mesh->labels;
901: mesh->labels = tmpLabel;
902: }
903: return(0);
904: }
908: /*@C
909: DMPlexGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
911: Not Collective
913: Input Parameters:
914: + dm - The DMPlex object
915: . name - The label name
916: - point - The mesh point
918: Output Parameter:
919: . value - The label value for this point, or -1 if the point is not in the label
921: Level: beginner
923: .keywords: mesh
924: .seealso: DMLabelGetValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
925: @*/
926: PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
927: {
928: DMLabel label;
934: DMPlexGetLabel(dm, name, &label);
935: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
936: DMLabelGetValue(label, point, value);
937: return(0);
938: }
942: /*@C
943: DMPlexSetLabelValue - Add a point to a Sieve Label with given value
945: Not Collective
947: Input Parameters:
948: + dm - The DMPlex object
949: . name - The label name
950: . point - The mesh point
951: - value - The label value for this point
953: Output Parameter:
955: Level: beginner
957: .keywords: mesh
958: .seealso: DMLabelSetValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
959: @*/
960: PetscErrorCode DMPlexSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
961: {
962: DMLabel label;
968: DMPlexGetLabel(dm, name, &label);
969: if (!label) {
970: DMPlexCreateLabel(dm, name);
971: DMPlexGetLabel(dm, name, &label);
972: }
973: DMLabelSetValue(label, point, value);
974: return(0);
975: }
979: /*@C
980: DMPlexClearLabelValue - Remove a point from a Sieve Label with given value
982: Not Collective
984: Input Parameters:
985: + dm - The DMPlex object
986: . name - The label name
987: . point - The mesh point
988: - value - The label value for this point
990: Output Parameter:
992: Level: beginner
994: .keywords: mesh
995: .seealso: DMLabelClearValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
996: @*/
997: PetscErrorCode DMPlexClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
998: {
999: DMLabel label;
1005: DMPlexGetLabel(dm, name, &label);
1006: if (!label) return(0);
1007: DMLabelClearValue(label, point, value);
1008: return(0);
1009: }
1013: /*@C
1014: DMPlexGetLabelSize - Get the number of different integer ids in a Label
1016: Not Collective
1018: Input Parameters:
1019: + dm - The DMPlex object
1020: - name - The label name
1022: Output Parameter:
1023: . size - The number of different integer ids, or 0 if the label does not exist
1025: Level: beginner
1027: .keywords: mesh
1028: .seealso: DMLabeGetNumValues(), DMPlexSetLabelValue()
1029: @*/
1030: PetscErrorCode DMPlexGetLabelSize(DM dm, const char name[], PetscInt *size)
1031: {
1032: DMLabel label;
1039: DMPlexGetLabel(dm, name, &label);
1040: *size = 0;
1041: if (!label) return(0);
1042: DMLabelGetNumValues(label, size);
1043: return(0);
1044: }
1048: /*@C
1049: DMPlexGetLabelIdIS - Get the integer ids in a label
1051: Not Collective
1053: Input Parameters:
1054: + mesh - The DMPlex object
1055: - name - The label name
1057: Output Parameter:
1058: . ids - The integer ids, or NULL if the label does not exist
1060: Level: beginner
1062: .keywords: mesh
1063: .seealso: DMLabelGetValueIS(), DMPlexGetLabelSize()
1064: @*/
1065: PetscErrorCode DMPlexGetLabelIdIS(DM dm, const char name[], IS *ids)
1066: {
1067: DMLabel label;
1074: DMPlexGetLabel(dm, name, &label);
1075: *ids = NULL;
1076: if (!label) return(0);
1077: DMLabelGetValueIS(label, ids);
1078: return(0);
1079: }
1083: /*@C
1084: DMPlexGetStratumSize - Get the number of points in a label stratum
1086: Not Collective
1088: Input Parameters:
1089: + dm - The DMPlex object
1090: . name - The label name
1091: - value - The stratum value
1093: Output Parameter:
1094: . size - The stratum size
1096: Level: beginner
1098: .keywords: mesh
1099: .seealso: DMLabelGetStratumSize(), DMPlexGetLabelSize(), DMPlexGetLabelIds()
1100: @*/
1101: PetscErrorCode DMPlexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1102: {
1103: DMLabel label;
1110: DMPlexGetLabel(dm, name, &label);
1111: *size = 0;
1112: if (!label) return(0);
1113: DMLabelGetStratumSize(label, value, size);
1114: return(0);
1115: }
1119: /*@C
1120: DMPlexGetStratumIS - Get the points in a label stratum
1122: Not Collective
1124: Input Parameters:
1125: + dm - The DMPlex object
1126: . name - The label name
1127: - value - The stratum value
1129: Output Parameter:
1130: . points - The stratum points, or NULL if the label does not exist or does not have that value
1132: Level: beginner
1134: .keywords: mesh
1135: .seealso: DMLabelGetStratumIS(), DMPlexGetStratumSize()
1136: @*/
1137: PetscErrorCode DMPlexGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
1138: {
1139: DMLabel label;
1146: DMPlexGetLabel(dm, name, &label);
1147: *points = NULL;
1148: if (!label) return(0);
1149: DMLabelGetStratumIS(label, value, points);
1150: return(0);
1151: }
1155: /*@C
1156: DMPlexClearLabelStratum - Remove all points from a stratum from a Sieve Label
1158: Not Collective
1160: Input Parameters:
1161: + dm - The DMPlex object
1162: . name - The label name
1163: - value - The label value for this point
1165: Output Parameter:
1167: Level: beginner
1169: .keywords: mesh
1170: .seealso: DMLabelClearStratum(), DMPlexSetLabelValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
1171: @*/
1172: PetscErrorCode DMPlexClearLabelStratum(DM dm, const char name[], PetscInt value)
1173: {
1174: DMLabel label;
1180: DMPlexGetLabel(dm, name, &label);
1181: if (!label) return(0);
1182: DMLabelClearStratum(label, value);
1183: return(0);
1184: }
1188: /*@
1189: DMPlexGetNumLabels - Return the number of labels defined by the mesh
1191: Not Collective
1193: Input Parameter:
1194: . dm - The DMPlex object
1196: Output Parameter:
1197: . numLabels - the number of Labels
1199: Level: intermediate
1201: .keywords: mesh
1202: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1203: @*/
1204: PetscErrorCode DMPlexGetNumLabels(DM dm, PetscInt *numLabels)
1205: {
1206: DM_Plex *mesh = (DM_Plex*) dm->data;
1207: PlexLabel next = mesh->labels;
1208: PetscInt n = 0;
1213: while (next) {++n; next = next->next;}
1214: *numLabels = n;
1215: return(0);
1216: }
1220: /*@C
1221: DMPlexGetLabelName - Return the name of nth label
1223: Not Collective
1225: Input Parameters:
1226: + dm - The DMPlex object
1227: - n - the label number
1229: Output Parameter:
1230: . name - the label name
1232: Level: intermediate
1234: .keywords: mesh
1235: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1236: @*/
1237: PetscErrorCode DMPlexGetLabelName(DM dm, PetscInt n, const char **name)
1238: {
1239: DM_Plex *mesh = (DM_Plex*) dm->data;
1240: PlexLabel next = mesh->labels;
1241: PetscInt l = 0;
1246: while (next) {
1247: if (l == n) {
1248: *name = next->label->name;
1249: return(0);
1250: }
1251: ++l;
1252: next = next->next;
1253: }
1254: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1255: }
1259: /*@C
1260: DMPlexHasLabel - Determine whether the mesh has a label of a given name
1262: Not Collective
1264: Input Parameters:
1265: + dm - The DMPlex object
1266: - name - The label name
1268: Output Parameter:
1269: . hasLabel - PETSC_TRUE if the label is present
1271: Level: intermediate
1273: .keywords: mesh
1274: .seealso: DMPlexCreateLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1275: @*/
1276: PetscErrorCode DMPlexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1277: {
1278: DM_Plex *mesh = (DM_Plex*) dm->data;
1279: PlexLabel next = mesh->labels;
1286: *hasLabel = PETSC_FALSE;
1287: while (next) {
1288: PetscStrcmp(name, next->label->name, hasLabel);
1289: if (*hasLabel) break;
1290: next = next->next;
1291: }
1292: return(0);
1293: }
1297: /*@C
1298: DMPlexGetLabel - Return the label of a given name, or NULL
1300: Not Collective
1302: Input Parameters:
1303: + dm - The DMPlex object
1304: - name - The label name
1306: Output Parameter:
1307: . label - The DMLabel, or NULL if the label is absent
1309: Level: intermediate
1311: .keywords: mesh
1312: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1313: @*/
1314: PetscErrorCode DMPlexGetLabel(DM dm, const char name[], DMLabel *label)
1315: {
1316: DM_Plex *mesh = (DM_Plex*) dm->data;
1317: PlexLabel next = mesh->labels;
1318: PetscBool hasLabel;
1325: *label = NULL;
1326: while (next) {
1327: PetscStrcmp(name, next->label->name, &hasLabel);
1328: if (hasLabel) {
1329: *label = next->label;
1330: break;
1331: }
1332: next = next->next;
1333: }
1334: return(0);
1335: }
1339: /*@C
1340: DMPlexGetLabelByNum - Return the nth label
1342: Not Collective
1344: Input Parameters:
1345: + dm - The DMPlex object
1346: - n - the label number
1348: Output Parameter:
1349: . label - the label
1351: Level: intermediate
1353: .keywords: mesh
1354: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1355: @*/
1356: PetscErrorCode DMPlexGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
1357: {
1358: DM_Plex *mesh = (DM_Plex*) dm->data;
1359: PlexLabel next = mesh->labels;
1360: PetscInt l = 0;
1365: while (next) {
1366: if (l == n) {
1367: *label = next->label;
1368: return(0);
1369: }
1370: ++l;
1371: next = next->next;
1372: }
1373: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1374: }
1378: /*@C
1379: DMPlexAddLabel - Add the label to this mesh
1381: Not Collective
1383: Input Parameters:
1384: + dm - The DMPlex object
1385: - label - The DMLabel
1387: Level: developer
1389: .keywords: mesh
1390: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1391: @*/
1392: PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
1393: {
1394: DM_Plex *mesh = (DM_Plex*) dm->data;
1395: PlexLabel tmpLabel;
1396: PetscBool hasLabel;
1401: DMPlexHasLabel(dm, label->name, &hasLabel);
1402: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
1403: PetscCalloc1(1, &tmpLabel);
1404: tmpLabel->label = label;
1405: tmpLabel->output = PETSC_TRUE;
1406: tmpLabel->next = mesh->labels;
1407: mesh->labels = tmpLabel;
1408: return(0);
1409: }
1413: /*@C
1414: DMPlexRemoveLabel - Remove the label from this mesh
1416: Not Collective
1418: Input Parameters:
1419: + dm - The DMPlex object
1420: - name - The label name
1422: Output Parameter:
1423: . label - The DMLabel, or NULL if the label is absent
1425: Level: developer
1427: .keywords: mesh
1428: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1429: @*/
1430: PetscErrorCode DMPlexRemoveLabel(DM dm, const char name[], DMLabel *label)
1431: {
1432: DM_Plex *mesh = (DM_Plex*) dm->data;
1433: PlexLabel next = mesh->labels;
1434: PlexLabel last = NULL;
1435: PetscBool hasLabel;
1440: DMPlexHasLabel(dm, name, &hasLabel);
1441: *label = NULL;
1442: if (!hasLabel) return(0);
1443: while (next) {
1444: PetscStrcmp(name, next->label->name, &hasLabel);
1445: if (hasLabel) {
1446: if (last) last->next = next->next;
1447: else mesh->labels = next->next;
1448: next->next = NULL;
1449: *label = next->label;
1450: PetscFree(next);
1451: break;
1452: }
1453: last = next;
1454: next = next->next;
1455: }
1456: return(0);
1457: }
1461: /*@C
1462: DMPlexGetLabelOutput - Get the output flag for a given label
1464: Not Collective
1466: Input Parameters:
1467: + dm - The DMPlex object
1468: - name - The label name
1470: Output Parameter:
1471: . output - The flag for output
1473: Level: developer
1475: .keywords: mesh
1476: .seealso: DMPlexSetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1477: @*/
1478: PetscErrorCode DMPlexGetLabelOutput(DM dm, const char name[], PetscBool *output)
1479: {
1480: DM_Plex *mesh = (DM_Plex*) dm->data;
1481: PlexLabel next = mesh->labels;
1488: while (next) {
1489: PetscBool flg;
1491: PetscStrcmp(name, next->label->name, &flg);
1492: if (flg) {*output = next->output; return(0);}
1493: next = next->next;
1494: }
1495: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1496: }
1500: /*@C
1501: DMPlexSetLabelOutput - Set the output flag for a given label
1503: Not Collective
1505: Input Parameters:
1506: + dm - The DMPlex object
1507: . name - The label name
1508: - output - The flag for output
1510: Level: developer
1512: .keywords: mesh
1513: .seealso: DMPlexGetLabelOutput(), DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1514: @*/
1515: PetscErrorCode DMPlexSetLabelOutput(DM dm, const char name[], PetscBool output)
1516: {
1517: DM_Plex *mesh = (DM_Plex*) dm->data;
1518: PlexLabel next = mesh->labels;
1524: while (next) {
1525: PetscBool flg;
1527: PetscStrcmp(name, next->label->name, &flg);
1528: if (flg) {next->output = output; return(0);}
1529: next = next->next;
1530: }
1531: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this mesh", name);
1532: }