Actual source code: plexlabel.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
5: PetscErrorCode DMLabelCreate(const char name[], DMLabel *label)
6: {
10: PetscNew(label);
11: PetscStrallocpy(name, &(*label)->name);
13: (*label)->refct = 1;
14: (*label)->numStrata = 0;
15: (*label)->stratumValues = NULL;
16: (*label)->arrayValid = PETSC_TRUE;
17: (*label)->stratumOffsets = 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: (*label)->next = NULL;
25: return(0);
26: }
30: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label)
31: {
32: PetscInt off = 0, v;
35: if (label->arrayValid) return 0;
37: PetscMalloc2(label->numStrata,&label->stratumSizes,label->numStrata+1,&label->stratumOffsets);
38: for (v = 0; v < label->numStrata; ++v) {
39: PetscInt size = 0;
40: PetscHashISize(label->ht[v], size);
41: label->stratumSizes[v] = size;
42: label->stratumOffsets[v] = off;
43: off += size;
44: }
45: label->stratumOffsets[v] = off;
46: PetscMalloc1(off, &label->points);
47: off = 0;
48: for (v = 0; v < label->numStrata; ++v) {
49: PetscHashIGetKeys(label->ht[v], &off, label->points);
50: if (off != label->stratumOffsets[v+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of contributed points %d from value %d should be %d", off-label->stratumOffsets[v], label->stratumValues[v], label->stratumOffsets[v+1]-label->stratumOffsets[v]);
51: PetscHashIDestroy(label->ht[v]);
52: PetscSortInt(label->stratumSizes[v], &label->points[label->stratumOffsets[v]]);
53: if (label->bt) {
54: PetscInt p;
56: for (p = label->stratumOffsets[v]; p < label->stratumOffsets[v]+label->stratumSizes[v]; ++p) {
57: const PetscInt point = label->points[p];
59: 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);
60: PetscBTSet(label->bt, point - label->pStart);
61: }
62: }
63: }
64: PetscFree(label->ht);
65: label->ht = NULL;
66: label->arrayValid = PETSC_TRUE;
67: return(0);
68: }
72: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label)
73: {
74: PetscInt v;
77: if (!label->arrayValid) return 0;
79: PetscMalloc1(label->numStrata, &label->ht);
80: for (v = 0; v < label->numStrata; ++v) {
81: PETSC_UNUSED PetscHashIIter ret, iter;
82: PetscInt p;
84: PetscHashICreate(label->ht[v]);
85: for (p = label->stratumOffsets[v]; p < label->stratumOffsets[v]+label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], label->points[p], ret, iter);
86: }
87: PetscFree2(label->stratumSizes,label->stratumOffsets);
88: PetscFree(label->points);
89: label->arrayValid = PETSC_FALSE;
90: return(0);
91: }
95: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
96: {
99: *name = label->name;
100: return(0);
101: }
105: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
106: {
107: PetscInt v;
108: PetscMPIInt rank;
112: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
113: if (label) {
114: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
115: if (label->bt) {PetscViewerASCIIPrintf(viewer, " Index has been calculated in [%d, %d)\n", label->pStart, label->pEnd);}
116: for (v = 0; v < label->numStrata; ++v) {
117: const PetscInt value = label->stratumValues[v];
118: PetscInt p;
120: for (p = label->stratumOffsets[v]; p < label->stratumOffsets[v]+label->stratumSizes[v]; ++p) {
121: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D (%D)\n", rank, label->points[p], value);
122: }
123: }
124: }
125: PetscViewerFlush(viewer);
126: return(0);
127: }
131: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
132: {
133: PetscBool iascii;
138: if (label) {DMLabelMakeValid_Private(label);}
139: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
140: if (iascii) {
141: DMLabelView_Ascii(label, viewer);
142: } else SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not supported by this mesh object", ((PetscObject) viewer)->type_name);
143: return(0);
144: }
148: PetscErrorCode DMLabelDestroy(DMLabel *label)
149: {
153: if (!(*label)) return(0);
154: if (--(*label)->refct > 0) return(0);
155: PetscFree((*label)->name);
156: PetscFree((*label)->stratumValues);
157: PetscFree2((*label)->stratumSizes,(*label)->stratumOffsets);
158: PetscFree((*label)->points);
159: if ((*label)->ht) {
160: PetscInt v;
161: for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
162: PetscFree((*label)->ht);
163: }
164: PetscBTDestroy(&(*label)->bt);
165: PetscFree(*label);
166: return(0);
167: }
171: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
172: {
173: PetscInt v, q;
177: DMLabelMakeValid_Private(label);
178: PetscNew(labelnew);
179: PetscStrallocpy(label->name, &(*labelnew)->name);
181: (*labelnew)->refct = 1;
182: (*labelnew)->numStrata = label->numStrata;
183: (*labelnew)->arrayValid = PETSC_TRUE;
184: if (label->numStrata) {
185: PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
186: PetscMalloc2(label->numStrata,&(*labelnew)->stratumSizes,label->numStrata+1,&(*labelnew)->stratumOffsets);
187: PetscMalloc1(label->stratumOffsets[label->numStrata], &(*labelnew)->points);
188: /* Could eliminate unused space here */
189: for (v = 0; v < label->numStrata; ++v) {
190: (*labelnew)->stratumValues[v] = label->stratumValues[v];
191: (*labelnew)->stratumOffsets[v] = label->stratumOffsets[v];
192: (*labelnew)->stratumSizes[v] = label->stratumSizes[v];
193: for (q = label->stratumOffsets[v]; q < label->stratumOffsets[v]+label->stratumSizes[v]; ++q) {
194: (*labelnew)->points[q] = label->points[q];
195: }
196: }
197: (*labelnew)->stratumOffsets[label->numStrata] = label->stratumOffsets[label->numStrata];
198: }
199: (*labelnew)->ht = NULL;
200: (*labelnew)->pStart = -1;
201: (*labelnew)->pEnd = -1;
202: (*labelnew)->bt = NULL;
203: return(0);
204: }
208: /* This can be hooked into SetValue(), ClearValue(), etc. for updating */
209: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
210: {
211: PetscInt v;
215: DMLabelMakeValid_Private(label);
216: if (label->bt) {PetscBTDestroy(&label->bt);}
217: label->pStart = pStart;
218: label->pEnd = pEnd;
219: PetscBTCreate(pEnd - pStart, &label->bt);
220: PetscBTMemzero(pEnd - pStart, label->bt);
221: for (v = 0; v < label->numStrata; ++v) {
222: PetscInt i;
224: for (i = 0; i < label->stratumSizes[v]; ++i) {
225: const PetscInt point = label->points[label->stratumOffsets[v]+i];
227: if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [%d, %d)", point, pStart, pEnd);
228: PetscBTSet(label->bt, point - pStart);
229: }
230: }
231: return(0);
232: }
236: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
237: {
241: label->pStart = -1;
242: label->pEnd = -1;
243: if (label->bt) {PetscBTDestroy(&label->bt);}
244: return(0);
245: }
249: /*@
250: DMLabelHasValue - Determine whether a label assigns the value to any point
252: Input Parameters:
253: + label - the DMLabel
254: - value - the value
256: Output Parameter:
257: . contains - Flag indicating whether the label maps this value to any point
259: Level: developer
261: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
262: @*/
263: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
264: {
265: PetscInt v;
269: for (v = 0; v < label->numStrata; ++v) {
270: if (value == label->stratumValues[v]) break;
271: }
272: *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
273: return(0);
274: }
278: /*@
279: DMLabelHasPoint - Determine whether a label assigns a value to a point
281: Input Parameters:
282: + label - the DMLabel
283: - point - the point
285: Output Parameter:
286: . contains - Flag indicating whether the label maps this point to a value
288: Note: The user must call DMLabelCreateIndex() before this function.
290: Level: developer
292: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
293: @*/
294: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
295: {
300: DMLabelMakeValid_Private(label);
301: #if defined(PETSC_USE_DEBUG)
302: if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
303: 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);
304: #endif
305: *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
306: return(0);
307: }
311: /*@
312: DMLabelStratumHasPoint - Return true if the stratum contains a point
314: Input Parameters:
315: + label - the DMLabel
316: . value - the stratum value
317: - point - the point
319: Output Parameter:
320: . contains - true if the stratum contains the point
322: Level: intermediate
324: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
325: @*/
326: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
327: {
328: PetscInt v;
333: *contains = PETSC_FALSE;
334: for (v = 0; v < label->numStrata; ++v) {
335: if (label->stratumValues[v] == value) {
336: if (label->arrayValid) {
337: PetscInt i;
339: PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &i);
340: if (i >= 0) {
341: *contains = PETSC_TRUE;
342: break;
343: }
344: } else {
345: PetscBool has;
347: PetscHashIHasKey(label->ht[v], point, has);
348: if (has) {
349: *contains = PETSC_TRUE;
350: break;
351: }
352: }
353: }
354: }
355: return(0);
356: }
360: /*@
361: DMLabelGetValue - Return the value a label assigns to a point, or -1
363: Input Parameters:
364: + label - the DMLabel
365: - point - the point
367: Output Parameter:
368: . value - The point value, or -1
370: Level: intermediate
372: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
373: @*/
374: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
375: {
376: PetscInt v;
381: *value = -1;
382: for (v = 0; v < label->numStrata; ++v) {
383: if (label->arrayValid) {
384: PetscInt i;
386: PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &i);
387: if (i >= 0) {
388: *value = label->stratumValues[v];
389: break;
390: }
391: } else {
392: PetscBool has;
394: PetscHashIHasKey(label->ht[v], point, has);
395: if (has) {
396: *value = label->stratumValues[v];
397: break;
398: }
399: }
400: }
401: return(0);
402: }
406: /*@
407: DMLabelSetValue - Set the value a label assigns to a point
409: Input Parameters:
410: + label - the DMLabel
411: . point - the point
412: - value - The point value
414: Level: intermediate
416: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue()
417: @*/
418: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
419: {
420: PETSC_UNUSED PetscHashIIter iter, ret;
421: PetscInt v;
422: PetscErrorCode ierr;
425: DMLabelMakeInvalid_Private(label);
426: /* Find, or add, label value */
427: for (v = 0; v < label->numStrata; ++v) {
428: if (label->stratumValues[v] == value) break;
429: }
430: /* Create new table */
431: if (v >= label->numStrata) {
432: PetscInt *tmpV;
433: PetscHashI *tmpH;
435: PetscMalloc1((label->numStrata+1), &tmpV);
436: PetscMalloc1((label->numStrata+1), &tmpH);
437: for (v = 0; v < label->numStrata; ++v) {
438: tmpV[v] = label->stratumValues[v];
439: tmpH[v] = label->ht[v];
440: }
441: tmpV[v] = value;
442: PetscHashICreate(tmpH[v]);
443: ++label->numStrata;
444: PetscFree(label->stratumValues);
445: PetscFree(label->ht);
446: label->stratumValues = tmpV;
447: label->ht = tmpH;
448: }
449: /* Set key */
450: PetscHashIPut(label->ht[v], point, ret, iter);
451: return(0);
452: }
456: /*@
457: DMLabelClearValue - Clear the value a label assigns to a point
459: Input Parameters:
460: + label - the DMLabel
461: . point - the point
462: - value - The point value
464: Level: intermediate
466: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
467: @*/
468: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
469: {
470: PetscInt v, p;
474: /* Find label value */
475: for (v = 0; v < label->numStrata; ++v) {
476: if (label->stratumValues[v] == value) break;
477: }
478: if (v >= label->numStrata) return(0);
479: if (label->arrayValid) {
480: /* Check whether point exists */
481: PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &p);
482: if (p >= 0) {
483: PetscMemmove(&label->points[p+label->stratumOffsets[v]], &label->points[p+label->stratumOffsets[v]+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
484: --label->stratumSizes[v];
485: if (label->bt) {
486: 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);
487: PetscBTClear(label->bt, point - label->pStart);
488: }
489: }
490: } else {
491: PetscHashIDelKey(label->ht[v], point);
492: }
493: return(0);
494: }
498: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
499: {
504: DMLabelMakeValid_Private(label);
505: *numValues = label->numStrata;
506: return(0);
507: }
511: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
512: {
517: DMLabelMakeValid_Private(label);
518: ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
519: return(0);
520: }
524: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
525: {
526: PetscInt v;
531: DMLabelMakeValid_Private(label);
532: *size = 0;
533: for (v = 0; v < label->numStrata; ++v) {
534: if (label->stratumValues[v] == value) {
535: *size = label->stratumSizes[v];
536: break;
537: }
538: }
539: return(0);
540: }
544: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
545: {
546: PetscInt v;
552: DMLabelMakeValid_Private(label);
553: for (v = 0; v < label->numStrata; ++v) {
554: if (label->stratumValues[v] != value) continue;
555: if (start) *start = label->points[label->stratumOffsets[v]];
556: if (end) *end = label->points[label->stratumOffsets[v]+label->stratumSizes[v]-1]+1;
557: break;
558: }
559: return(0);
560: }
564: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
565: {
566: PetscInt v;
571: DMLabelMakeValid_Private(label);
572: *points = NULL;
573: for (v = 0; v < label->numStrata; ++v) {
574: if (label->stratumValues[v] == value) {
575: if (label->arrayValid) {
576: ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], PETSC_COPY_VALUES, points);
577: PetscObjectSetName((PetscObject) *points, "indices");
578: } else {
579: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
580: }
581: break;
582: }
583: }
584: return(0);
585: }
589: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
590: {
591: PetscInt v;
595: for (v = 0; v < label->numStrata; ++v) {
596: if (label->stratumValues[v] == value) break;
597: }
598: if (v >= label->numStrata) return(0);
599: if (label->bt) {
600: PetscInt i;
602: for (i = 0; i < label->stratumSizes[v]; ++i) {
603: const PetscInt point = label->points[label->stratumOffsets[v]+i];
605: 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);
606: PetscBTClear(label->bt, point - label->pStart);
607: }
608: }
609: if (label->arrayValid) {
610: label->stratumSizes[v] = 0;
611: } else {
612: PetscHashIClear(label->ht[v]);
613: }
614: return(0);
615: }
619: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
620: {
621: PetscInt v;
625: DMLabelMakeValid_Private(label);
626: label->pStart = start;
627: label->pEnd = end;
628: if (label->bt) {PetscBTDestroy(&label->bt);}
629: /* Could squish offsets, but would only make sense if I reallocate the storage */
630: for (v = 0; v < label->numStrata; ++v) {
631: const PetscInt offset = label->stratumOffsets[v];
632: const PetscInt size = label->stratumSizes[v];
633: PetscInt off = offset, q;
635: for (q = offset; q < offset+size; ++q) {
636: const PetscInt point = label->points[q];
638: if ((point < start) || (point >= end)) continue;
639: label->points[off++] = point;
640: }
641: label->stratumSizes[v] = off-offset;
642: }
643: DMLabelCreateIndex(label, start, end);
644: return(0);
645: }
649: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
650: {
651: const PetscInt *perm;
652: PetscInt numValues, numPoints, v, q;
653: PetscErrorCode ierr;
656: DMLabelMakeValid_Private(label);
657: DMLabelDuplicate(label, labelNew);
658: DMLabelGetNumValues(*labelNew, &numValues);
659: ISGetLocalSize(permutation, &numPoints);
660: ISGetIndices(permutation, &perm);
661: for (v = 0; v < numValues; ++v) {
662: const PetscInt offset = (*labelNew)->stratumOffsets[v];
663: const PetscInt size = (*labelNew)->stratumSizes[v];
665: for (q = offset; q < offset+size; ++q) {
666: const PetscInt point = (*labelNew)->points[q];
668: 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);
669: (*labelNew)->points[q] = perm[point];
670: }
671: PetscSortInt(size, &(*labelNew)->points[offset]);
672: }
673: ISRestoreIndices(permutation, &perm);
674: if (label->bt) {
675: PetscBTDestroy(&label->bt);
676: DMLabelCreateIndex(label, label->pStart, label->pEnd);
677: }
678: return(0);
679: }
683: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSection partSection, IS part, ISLocalToGlobalMapping renumbering, DMLabel *labelNew)
684: {
685: MPI_Comm comm = PetscObjectComm((PetscObject) partSection);
686: PetscInt *stratumSizes = NULL, *points = NULL, s, p;
687: PetscMPIInt *sendcnts = NULL, *offsets = NULL, *displs = NULL, proc;
688: char *name;
689: PetscInt nameSize;
690: size_t len = 0;
691: PetscMPIInt rank, numProcs;
695: if (label) {DMLabelMakeValid_Private(label);}
696: MPI_Comm_rank(comm, &rank);
697: MPI_Comm_size(comm, &numProcs);
698: /* Bcast name */
699: if (!rank) {PetscStrlen(label->name, &len);}
700: nameSize = len;
701: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
702: PetscMalloc(nameSize+1, &name);
703: if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
704: MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
705: DMLabelCreate(name, labelNew);
706: PetscFree(name);
707: /* Bcast numStrata */
708: if (!rank) (*labelNew)->numStrata = label->numStrata;
709: MPI_Bcast(&(*labelNew)->numStrata, 1, MPIU_INT, 0, comm);
710: PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
711: PetscMalloc2((*labelNew)->numStrata,&(*labelNew)->stratumSizes,(*labelNew)->numStrata+1,&(*labelNew)->stratumOffsets);
712: /* Bcast stratumValues */
713: if (!rank) {PetscMemcpy((*labelNew)->stratumValues, label->stratumValues, label->numStrata * sizeof(PetscInt));}
714: MPI_Bcast((*labelNew)->stratumValues, (*labelNew)->numStrata, MPIU_INT, 0, comm);
715: /* Find size on each process and Scatter: we use the fact that both the stratum points and partArray are sorted */
716: if (!rank) {
717: const PetscInt *partArray;
718: PetscInt proc;
720: ISGetIndices(part, &partArray);
721: PetscCalloc1(numProcs*label->numStrata, &stratumSizes);
722: /* TODO We should switch to using binary search if the label is a lot smaller than partitions */
723: for (proc = 0; proc < numProcs; ++proc) {
724: PetscInt dof, off;
726: PetscSectionGetDof(partSection, proc, &dof);
727: PetscSectionGetOffset(partSection, proc, &off);
728: for (s = 0; s < label->numStrata; ++s) {
729: PetscInt lStart = label->stratumOffsets[s], lEnd = label->stratumOffsets[s]+label->stratumSizes[s];
730: PetscInt pStart = off, pEnd = off+dof;
732: while (pStart < pEnd && lStart < lEnd) {
733: if (partArray[pStart] > label->points[lStart]) {
734: ++lStart;
735: } else if (label->points[lStart] > partArray[pStart]) {
736: ++pStart;
737: } else {
738: ++stratumSizes[proc*label->numStrata+s];
739: ++pStart; ++lStart;
740: }
741: }
742: }
743: }
744: ISRestoreIndices(part, &partArray);
745: }
746: MPI_Scatter(stratumSizes, (*labelNew)->numStrata, MPIU_INT, (*labelNew)->stratumSizes, (*labelNew)->numStrata, MPIU_INT, 0, comm);
747: /* Calculate stratumOffsets */
748: (*labelNew)->stratumOffsets[0] = 0;
749: for (s = 0; s < (*labelNew)->numStrata; ++s) {(*labelNew)->stratumOffsets[s+1] = (*labelNew)->stratumSizes[s] + (*labelNew)->stratumOffsets[s];}
750: /* Pack points and Scatter */
751: if (!rank) {
752: const PetscInt *partArray;
754: ISGetIndices(part, &partArray);
755: PetscMalloc3(numProcs,&sendcnts,numProcs,&offsets,numProcs+1,&displs);
756: displs[0] = 0;
757: for (p = 0; p < numProcs; ++p) {
758: sendcnts[p] = 0;
759: for (s = 0; s < label->numStrata; ++s) sendcnts[p] += stratumSizes[p*label->numStrata+s];
760: offsets[p] = displs[p];
761: displs[p+1] = displs[p] + sendcnts[p];
762: }
763: PetscMalloc1(displs[numProcs], &points);
764: /* TODO We should switch to using binary search if the label is a lot smaller than partitions */
765: for (proc = 0; proc < numProcs; ++proc) {
766: PetscInt dof, off;
768: PetscSectionGetDof(partSection, proc, &dof);
769: PetscSectionGetOffset(partSection, proc, &off);
770: for (s = 0; s < label->numStrata; ++s) {
771: PetscInt lStart = label->stratumOffsets[s], lEnd = label->stratumOffsets[s]+label->stratumSizes[s];
772: PetscInt pStart = off, pEnd = off+dof;
774: while (pStart < pEnd && lStart < lEnd) {
775: if (partArray[pStart] > label->points[lStart]) {
776: ++lStart;
777: } else if (label->points[lStart] > partArray[pStart]) {
778: ++pStart;
779: } else {
780: points[offsets[proc]++] = label->points[lStart];
781: ++pStart; ++lStart;
782: }
783: }
784: }
785: }
786: ISRestoreIndices(part, &partArray);
787: }
788: PetscMalloc1((*labelNew)->stratumOffsets[(*labelNew)->numStrata], &(*labelNew)->points);
789: MPI_Scatterv(points, sendcnts, displs, MPIU_INT, (*labelNew)->points, (*labelNew)->stratumOffsets[(*labelNew)->numStrata], MPIU_INT, 0, comm);
790: PetscFree(points);
791: PetscFree3(sendcnts,offsets,displs);
792: PetscFree(stratumSizes);
793: /* Renumber points */
794: ISGlobalToLocalMappingApplyBlock(renumbering, IS_GTOLM_MASK, (*labelNew)->stratumOffsets[(*labelNew)->numStrata], (*labelNew)->points, NULL, (*labelNew)->points);
795: /* Sort points */
796: for (s = 0; s < (*labelNew)->numStrata; ++s) {PetscSortInt((*labelNew)->stratumSizes[s], &(*labelNew)->points[(*labelNew)->stratumOffsets[s]]);}
797: return(0);
798: }
803: /*@C
804: DMPlexCreateLabel - Create a label of the given name if it does not already exist
806: Not Collective
808: Input Parameters:
809: + dm - The DMPlex object
810: - name - The label name
812: Level: intermediate
814: .keywords: mesh
815: .seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
816: @*/
817: PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
818: {
819: DM_Plex *mesh = (DM_Plex*) dm->data;
820: DMLabel next = mesh->labels;
821: PetscBool flg = PETSC_FALSE;
827: while (next) {
828: PetscStrcmp(name, next->name, &flg);
829: if (flg) break;
830: next = next->next;
831: }
832: if (!flg) {
833: DMLabel tmpLabel = mesh->labels;
835: DMLabelCreate(name, &mesh->labels);
837: mesh->labels->next = tmpLabel;
838: }
839: return(0);
840: }
844: /*@C
845: DMPlexGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
847: Not Collective
849: Input Parameters:
850: + dm - The DMPlex object
851: . name - The label name
852: - point - The mesh point
854: Output Parameter:
855: . value - The label value for this point, or -1 if the point is not in the label
857: Level: beginner
859: .keywords: mesh
860: .seealso: DMLabelGetValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
861: @*/
862: PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
863: {
864: DMLabel label;
870: DMPlexGetLabel(dm, name, &label);
871: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
872: DMLabelGetValue(label, point, value);
873: return(0);
874: }
878: /*@C
879: DMPlexSetLabelValue - Add a point to a Sieve Label with given value
881: Not Collective
883: Input Parameters:
884: + dm - The DMPlex object
885: . name - The label name
886: . point - The mesh point
887: - value - The label value for this point
889: Output Parameter:
891: Level: beginner
893: .keywords: mesh
894: .seealso: DMLabelSetValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
895: @*/
896: PetscErrorCode DMPlexSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
897: {
898: DMLabel label;
904: DMPlexGetLabel(dm, name, &label);
905: if (!label) {
906: DMPlexCreateLabel(dm, name);
907: DMPlexGetLabel(dm, name, &label);
908: }
909: DMLabelSetValue(label, point, value);
910: return(0);
911: }
915: /*@C
916: DMPlexClearLabelValue - Remove a point from a Sieve Label with given value
918: Not Collective
920: Input Parameters:
921: + dm - The DMPlex object
922: . name - The label name
923: . point - The mesh point
924: - value - The label value for this point
926: Output Parameter:
928: Level: beginner
930: .keywords: mesh
931: .seealso: DMLabelClearValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
932: @*/
933: PetscErrorCode DMPlexClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
934: {
935: DMLabel label;
941: DMPlexGetLabel(dm, name, &label);
942: if (!label) return(0);
943: DMLabelClearValue(label, point, value);
944: return(0);
945: }
949: /*@C
950: DMPlexGetLabelSize - Get the number of different integer ids in a Label
952: Not Collective
954: Input Parameters:
955: + dm - The DMPlex object
956: - name - The label name
958: Output Parameter:
959: . size - The number of different integer ids, or 0 if the label does not exist
961: Level: beginner
963: .keywords: mesh
964: .seealso: DMLabeGetNumValues(), DMPlexSetLabelValue()
965: @*/
966: PetscErrorCode DMPlexGetLabelSize(DM dm, const char name[], PetscInt *size)
967: {
968: DMLabel label;
975: DMPlexGetLabel(dm, name, &label);
976: *size = 0;
977: if (!label) return(0);
978: DMLabelGetNumValues(label, size);
979: return(0);
980: }
984: /*@C
985: DMPlexGetLabelIdIS - Get the integer ids in a label
987: Not Collective
989: Input Parameters:
990: + mesh - The DMPlex object
991: - name - The label name
993: Output Parameter:
994: . ids - The integer ids, or NULL if the label does not exist
996: Level: beginner
998: .keywords: mesh
999: .seealso: DMLabelGetValueIS(), DMPlexGetLabelSize()
1000: @*/
1001: PetscErrorCode DMPlexGetLabelIdIS(DM dm, const char name[], IS *ids)
1002: {
1003: DMLabel label;
1010: DMPlexGetLabel(dm, name, &label);
1011: *ids = NULL;
1012: if (!label) return(0);
1013: DMLabelGetValueIS(label, ids);
1014: return(0);
1015: }
1019: /*@C
1020: DMPlexGetStratumSize - Get the number of points in a label stratum
1022: Not Collective
1024: Input Parameters:
1025: + dm - The DMPlex object
1026: . name - The label name
1027: - value - The stratum value
1029: Output Parameter:
1030: . size - The stratum size
1032: Level: beginner
1034: .keywords: mesh
1035: .seealso: DMLabelGetStratumSize(), DMPlexGetLabelSize(), DMPlexGetLabelIds()
1036: @*/
1037: PetscErrorCode DMPlexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1038: {
1039: DMLabel label;
1046: DMPlexGetLabel(dm, name, &label);
1047: *size = 0;
1048: if (!label) return(0);
1049: DMLabelGetStratumSize(label, value, size);
1050: return(0);
1051: }
1055: /*@C
1056: DMPlexGetStratumIS - Get the points in a label stratum
1058: Not Collective
1060: Input Parameters:
1061: + dm - The DMPlex object
1062: . name - The label name
1063: - value - The stratum value
1065: Output Parameter:
1066: . points - The stratum points, or NULL if the label does not exist or does not have that value
1068: Level: beginner
1070: .keywords: mesh
1071: .seealso: DMLabelGetStratumIS(), DMPlexGetStratumSize()
1072: @*/
1073: PetscErrorCode DMPlexGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
1074: {
1075: DMLabel label;
1082: DMPlexGetLabel(dm, name, &label);
1083: *points = NULL;
1084: if (!label) return(0);
1085: DMLabelGetStratumIS(label, value, points);
1086: return(0);
1087: }
1091: /*@C
1092: DMPlexClearLabelStratum - Remove all points from a stratum from a Sieve Label
1094: Not Collective
1096: Input Parameters:
1097: + dm - The DMPlex object
1098: . name - The label name
1099: - value - The label value for this point
1101: Output Parameter:
1103: Level: beginner
1105: .keywords: mesh
1106: .seealso: DMLabelClearStratum(), DMPlexSetLabelValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
1107: @*/
1108: PetscErrorCode DMPlexClearLabelStratum(DM dm, const char name[], PetscInt value)
1109: {
1110: DMLabel label;
1116: DMPlexGetLabel(dm, name, &label);
1117: if (!label) return(0);
1118: DMLabelClearStratum(label, value);
1119: return(0);
1120: }
1124: /*@
1125: DMPlexGetNumLabels - Return the number of labels defined by the mesh
1127: Not Collective
1129: Input Parameter:
1130: . dm - The DMPlex object
1132: Output Parameter:
1133: . numLabels - the number of Labels
1135: Level: intermediate
1137: .keywords: mesh
1138: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1139: @*/
1140: PetscErrorCode DMPlexGetNumLabels(DM dm, PetscInt *numLabels)
1141: {
1142: DM_Plex *mesh = (DM_Plex*) dm->data;
1143: DMLabel next = mesh->labels;
1144: PetscInt n = 0;
1149: while (next) {
1150: ++n;
1151: next = next->next;
1152: }
1153: *numLabels = n;
1154: return(0);
1155: }
1159: /*@C
1160: DMPlexGetLabelName - Return the name of nth label
1162: Not Collective
1164: Input Parameters:
1165: + dm - The DMPlex object
1166: - n - the label number
1168: Output Parameter:
1169: . name - the label name
1171: Level: intermediate
1173: .keywords: mesh
1174: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1175: @*/
1176: PetscErrorCode DMPlexGetLabelName(DM dm, PetscInt n, const char **name)
1177: {
1178: DM_Plex *mesh = (DM_Plex*) dm->data;
1179: DMLabel next = mesh->labels;
1180: PetscInt l = 0;
1185: while (next) {
1186: if (l == n) {
1187: *name = next->name;
1188: return(0);
1189: }
1190: ++l;
1191: next = next->next;
1192: }
1193: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1194: }
1198: /*@C
1199: DMPlexHasLabel - Determine whether the mesh has a label of a given name
1201: Not Collective
1203: Input Parameters:
1204: + dm - The DMPlex object
1205: - name - The label name
1207: Output Parameter:
1208: . hasLabel - PETSC_TRUE if the label is present
1210: Level: intermediate
1212: .keywords: mesh
1213: .seealso: DMPlexCreateLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1214: @*/
1215: PetscErrorCode DMPlexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1216: {
1217: DM_Plex *mesh = (DM_Plex*) dm->data;
1218: DMLabel next = mesh->labels;
1225: *hasLabel = PETSC_FALSE;
1226: while (next) {
1227: PetscStrcmp(name, next->name, hasLabel);
1228: if (*hasLabel) break;
1229: next = next->next;
1230: }
1231: return(0);
1232: }
1236: /*@C
1237: DMPlexGetLabel - Return the label of a given name, or NULL
1239: Not Collective
1241: Input Parameters:
1242: + dm - The DMPlex object
1243: - name - The label name
1245: Output Parameter:
1246: . label - The DMLabel, or NULL if the label is absent
1248: Level: intermediate
1250: .keywords: mesh
1251: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1252: @*/
1253: PetscErrorCode DMPlexGetLabel(DM dm, const char name[], DMLabel *label)
1254: {
1255: DM_Plex *mesh = (DM_Plex*) dm->data;
1256: DMLabel next = mesh->labels;
1257: PetscBool hasLabel;
1264: *label = NULL;
1265: while (next) {
1266: PetscStrcmp(name, next->name, &hasLabel);
1267: if (hasLabel) {
1268: *label = next;
1269: break;
1270: }
1271: next = next->next;
1272: }
1273: return(0);
1274: }
1278: /*@C
1279: DMPlexAddLabel - Add the label to this mesh
1281: Not Collective
1283: Input Parameters:
1284: + dm - The DMPlex object
1285: - label - The DMLabel
1287: Level: developer
1289: .keywords: mesh
1290: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1291: @*/
1292: PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
1293: {
1294: DM_Plex *mesh = (DM_Plex*) dm->data;
1295: PetscBool hasLabel;
1300: DMPlexHasLabel(dm, label->name, &hasLabel);
1301: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
1302: label->next = mesh->labels;
1303: mesh->labels = label;
1304: return(0);
1305: }
1309: /*@C
1310: DMPlexRemoveLabel - Remove the label from this mesh
1312: Not Collective
1314: Input Parameters:
1315: + dm - The DMPlex object
1316: - name - The label name
1318: Output Parameter:
1319: . label - The DMLabel, or NULL if the label is absent
1321: Level: developer
1323: .keywords: mesh
1324: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1325: @*/
1326: PetscErrorCode DMPlexRemoveLabel(DM dm, const char name[], DMLabel *label)
1327: {
1328: DM_Plex *mesh = (DM_Plex*) dm->data;
1329: DMLabel next = mesh->labels;
1330: DMLabel last = NULL;
1331: PetscBool hasLabel;
1336: DMPlexHasLabel(dm, name, &hasLabel);
1337: *label = NULL;
1338: if (!hasLabel) return(0);
1339: while (next) {
1340: PetscStrcmp(name, next->name, &hasLabel);
1341: if (hasLabel) {
1342: if (last) last->next = next->next;
1343: else mesh->labels = next->next;
1344: next->next = NULL;
1345: *label = next;
1346: break;
1347: }
1348: last = next;
1349: next = next->next;
1350: }
1351: return(0);
1352: }