Actual source code: plexlabel.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  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: }