Actual source code: plexlabel.c

petsc-dev 2014-04-18
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:   DMLabelGetValue - Return the value a label assigns to a point, or -1

314:   Input Parameters:
315: + label - the DMLabel
316: - point - the point

318:   Output Parameter:
319: . value - The point value, or -1

321:   Level: intermediate

323: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
324: @*/
325: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
326: {
327:   PetscInt       v;

332:   *value = -1;
333:   for (v = 0; v < label->numStrata; ++v) {
334:     if (label->arrayValid) {
335:       PetscInt i;

337:       PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &i);
338:       if (i >= 0) {
339:         *value = label->stratumValues[v];
340:         break;
341:       }
342:     } else {
343:       PetscBool has;

345:       PetscHashIHasKey(label->ht[v], point, has);
346:       if (has) {
347:         *value = label->stratumValues[v];
348:         break;
349:       }
350:     }
351:   }
352:   return(0);
353: }

357: /*@
358:   DMLabelSetValue - Set the value a label assigns to a point

360:   Input Parameters:
361: + label - the DMLabel
362: . point - the point
363: - value - The point value

365:   Level: intermediate

367: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue()
368: @*/
369: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
370: {
371:   PETSC_UNUSED PetscHashIIter iter, ret;
372:   PetscInt                    v;
373:   PetscErrorCode              ierr;

376:   DMLabelMakeInvalid_Private(label);
377:   /* Find, or add, label value */
378:   for (v = 0; v < label->numStrata; ++v) {
379:     if (label->stratumValues[v] == value) break;
380:   }
381:   /* Create new table */
382:   if (v >= label->numStrata) {
383:     PetscInt   *tmpV;
384:     PetscHashI *tmpH;

386:     PetscMalloc1((label->numStrata+1), &tmpV);
387:     PetscMalloc1((label->numStrata+1), &tmpH);
388:     for (v = 0; v < label->numStrata; ++v) {
389:       tmpV[v] = label->stratumValues[v];
390:       tmpH[v] = label->ht[v];
391:     }
392:     tmpV[v] = value;
393:     PetscHashICreate(tmpH[v]);
394:     ++label->numStrata;
395:     PetscFree(label->stratumValues);
396:     PetscFree(label->ht);
397:     label->stratumValues = tmpV;
398:     label->ht            = tmpH;
399:   }
400:   /* Set key */
401:   PetscHashIPut(label->ht[v], point, ret, iter);
402:   return(0);
403: }

407: /*@
408:   DMLabelClearValue - Clear the value a label assigns to a point

410:   Input Parameters:
411: + label - the DMLabel
412: . point - the point
413: - value - The point value

415:   Level: intermediate

417: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
418: @*/
419: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
420: {
421:   PetscInt       v, p;

425:   /* Find label value */
426:   for (v = 0; v < label->numStrata; ++v) {
427:     if (label->stratumValues[v] == value) break;
428:   }
429:   if (v >= label->numStrata) return(0);
430:   if (label->arrayValid) {
431:     /* Check whether point exists */
432:     PetscFindInt(point, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], &p);
433:     if (p >= 0) {
434:       PetscMemmove(&label->points[p+label->stratumOffsets[v]], &label->points[p+label->stratumOffsets[v]+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
435:       --label->stratumSizes[v];
436:       if (label->bt) {
437:         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);
438:         PetscBTClear(label->bt, point - label->pStart);
439:       }
440:     }
441:   } else {
442:     PetscHashIDelKey(label->ht[v], point);
443:   }
444:   return(0);
445: }

449: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
450: {

455:   DMLabelMakeValid_Private(label);
456:   *numValues = label->numStrata;
457:   return(0);
458: }

462: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
463: {

468:   DMLabelMakeValid_Private(label);
469:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
470:   return(0);
471: }

475: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
476: {
477:   PetscInt       v;

482:   DMLabelMakeValid_Private(label);
483:   *size = 0;
484:   for (v = 0; v < label->numStrata; ++v) {
485:     if (label->stratumValues[v] == value) {
486:       *size = label->stratumSizes[v];
487:       break;
488:     }
489:   }
490:   return(0);
491: }

495: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
496: {
497:   PetscInt       v;

503:   DMLabelMakeValid_Private(label);
504:   for (v = 0; v < label->numStrata; ++v) {
505:     if (label->stratumValues[v] != value) continue;
506:     if (start) *start = label->points[label->stratumOffsets[v]];
507:     if (end)   *end   = label->points[label->stratumOffsets[v]+label->stratumSizes[v]-1]+1;
508:     break;
509:   }
510:   return(0);
511: }

515: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
516: {
517:   PetscInt       v;

522:   DMLabelMakeValid_Private(label);
523:   *points = NULL;
524:   for (v = 0; v < label->numStrata; ++v) {
525:     if (label->stratumValues[v] == value) {
526:       if (label->arrayValid) {
527:         ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[label->stratumOffsets[v]], PETSC_COPY_VALUES, points);
528:         PetscObjectSetName((PetscObject) *points, "indices");
529:       } else {
530:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
531:       }
532:       break;
533:     }
534:   }
535:   return(0);
536: }

540: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
541: {
542:   PetscInt       v;

546:   for (v = 0; v < label->numStrata; ++v) {
547:     if (label->stratumValues[v] == value) break;
548:   }
549:   if (v >= label->numStrata) return(0);
550:   if (label->bt) {
551:     PetscInt i;

553:     for (i = 0; i < label->stratumSizes[v]; ++i) {
554:       const PetscInt point = label->points[label->stratumOffsets[v]+i];

556:       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);
557:       PetscBTClear(label->bt, point - label->pStart);
558:     }
559:   }
560:   if (label->arrayValid) {
561:     label->stratumSizes[v] = 0;
562:   } else {
563:     PetscHashIClear(label->ht[v]);
564:   }
565:   return(0);
566: }

570: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
571: {
572:   PetscInt       v;

576:   DMLabelMakeValid_Private(label);
577:   label->pStart = start;
578:   label->pEnd   = end;
579:   if (label->bt) {PetscBTDestroy(&label->bt);}
580:   /* Could squish offsets, but would only make sense if I reallocate the storage */
581:   for (v = 0; v < label->numStrata; ++v) {
582:     const PetscInt offset = label->stratumOffsets[v];
583:     const PetscInt size   = label->stratumSizes[v];
584:     PetscInt       off    = offset, q;

586:     for (q = offset; q < offset+size; ++q) {
587:       const PetscInt point = label->points[q];

589:       if ((point < start) || (point >= end)) continue;
590:       label->points[off++] = point;
591:     }
592:     label->stratumSizes[v] = off-offset;
593:   }
594:   DMLabelCreateIndex(label, start, end);
595:   return(0);
596: }

600: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
601: {
602:   const PetscInt *perm;
603:   PetscInt        numValues, numPoints, v, q;
604:   PetscErrorCode  ierr;

607:   DMLabelMakeValid_Private(label);
608:   DMLabelDuplicate(label, labelNew);
609:   DMLabelGetNumValues(*labelNew, &numValues);
610:   ISGetLocalSize(permutation, &numPoints);
611:   ISGetIndices(permutation, &perm);
612:   for (v = 0; v < numValues; ++v) {
613:     const PetscInt offset = (*labelNew)->stratumOffsets[v];
614:     const PetscInt size   = (*labelNew)->stratumSizes[v];

616:     for (q = offset; q < offset+size; ++q) {
617:       const PetscInt point = (*labelNew)->points[q];

619:       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);
620:       (*labelNew)->points[q] = perm[point];
621:     }
622:     PetscSortInt(size, &(*labelNew)->points[offset]);
623:   }
624:   ISRestoreIndices(permutation, &perm);
625:   if (label->bt) {
626:     PetscBTDestroy(&label->bt);
627:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
628:   }
629:   return(0);
630: }

634: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSection partSection, IS part, ISLocalToGlobalMapping renumbering, DMLabel *labelNew)
635: {
636:   MPI_Comm       comm = PetscObjectComm((PetscObject) partSection);
637:   PetscInt      *stratumSizes = NULL, *points = NULL, s, p;
638:   PetscMPIInt   *sendcnts = NULL, *offsets = NULL, *displs = NULL, proc;
639:   char          *name;
640:   PetscInt       nameSize;
641:   size_t         len = 0;
642:   PetscMPIInt    rank, numProcs;

646:   if (label) {DMLabelMakeValid_Private(label);}
647:   MPI_Comm_rank(comm, &rank);
648:   MPI_Comm_size(comm, &numProcs);
649:   /* Bcast name */
650:   if (!rank) {PetscStrlen(label->name, &len);}
651:   nameSize = len;
652:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
653:   PetscMalloc(nameSize+1, &name);
654:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
655:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
656:   DMLabelCreate(name, labelNew);
657:   PetscFree(name);
658:   /* Bcast numStrata */
659:   if (!rank) (*labelNew)->numStrata = label->numStrata;
660:   MPI_Bcast(&(*labelNew)->numStrata, 1, MPIU_INT, 0, comm);
661:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
662:   PetscMalloc2((*labelNew)->numStrata,&(*labelNew)->stratumSizes,(*labelNew)->numStrata+1,&(*labelNew)->stratumOffsets);
663:   /* Bcast stratumValues */
664:   if (!rank) {PetscMemcpy((*labelNew)->stratumValues, label->stratumValues, label->numStrata * sizeof(PetscInt));}
665:   MPI_Bcast((*labelNew)->stratumValues, (*labelNew)->numStrata, MPIU_INT, 0, comm);
666:   /* Find size on each process and Scatter: we use the fact that both the stratum points and partArray are sorted */
667:   if (!rank) {
668:     const PetscInt *partArray;
669:     PetscInt        proc;

671:     ISGetIndices(part, &partArray);
672:     PetscCalloc1(numProcs*label->numStrata, &stratumSizes);
673:     /* TODO We should switch to using binary search if the label is a lot smaller than partitions */
674:     for (proc = 0; proc < numProcs; ++proc) {
675:       PetscInt dof, off;

677:       PetscSectionGetDof(partSection, proc, &dof);
678:       PetscSectionGetOffset(partSection, proc, &off);
679:       for (s = 0; s < label->numStrata; ++s) {
680:         PetscInt lStart = label->stratumOffsets[s], lEnd = label->stratumOffsets[s]+label->stratumSizes[s];
681:         PetscInt pStart = off,                      pEnd = off+dof;

683:         while (pStart < pEnd && lStart < lEnd) {
684:           if (partArray[pStart] > label->points[lStart]) {
685:             ++lStart;
686:           } else if (label->points[lStart] > partArray[pStart]) {
687:             ++pStart;
688:           } else {
689:             ++stratumSizes[proc*label->numStrata+s];
690:             ++pStart; ++lStart;
691:           }
692:         }
693:       }
694:     }
695:     ISRestoreIndices(part, &partArray);
696:   }
697:   MPI_Scatter(stratumSizes, (*labelNew)->numStrata, MPIU_INT, (*labelNew)->stratumSizes, (*labelNew)->numStrata, MPIU_INT, 0, comm);
698:   /* Calculate stratumOffsets */
699:   (*labelNew)->stratumOffsets[0] = 0;
700:   for (s = 0; s < (*labelNew)->numStrata; ++s) {(*labelNew)->stratumOffsets[s+1] = (*labelNew)->stratumSizes[s] + (*labelNew)->stratumOffsets[s];}
701:   /* Pack points and Scatter */
702:   if (!rank) {
703:     const PetscInt *partArray;

705:     ISGetIndices(part, &partArray);
706:     PetscMalloc3(numProcs,&sendcnts,numProcs,&offsets,numProcs+1,&displs);
707:     displs[0] = 0;
708:     for (p = 0; p < numProcs; ++p) {
709:       sendcnts[p] = 0;
710:       for (s = 0; s < label->numStrata; ++s) sendcnts[p] += stratumSizes[p*label->numStrata+s];
711:       offsets[p]  = displs[p];
712:       displs[p+1] = displs[p] + sendcnts[p];
713:     }
714:     PetscMalloc1(displs[numProcs], &points);
715:     /* TODO We should switch to using binary search if the label is a lot smaller than partitions */
716:     for (proc = 0; proc < numProcs; ++proc) {
717:       PetscInt dof, off;

719:       PetscSectionGetDof(partSection, proc, &dof);
720:       PetscSectionGetOffset(partSection, proc, &off);
721:       for (s = 0; s < label->numStrata; ++s) {
722:         PetscInt lStart = label->stratumOffsets[s], lEnd = label->stratumOffsets[s]+label->stratumSizes[s];
723:         PetscInt pStart = off,                     pEnd = off+dof;

725:         while (pStart < pEnd && lStart < lEnd) {
726:           if (partArray[pStart] > label->points[lStart]) {
727:             ++lStart;
728:           } else if (label->points[lStart] > partArray[pStart]) {
729:             ++pStart;
730:           } else {
731:             points[offsets[proc]++] = label->points[lStart];
732:             ++pStart; ++lStart;
733:           }
734:         }
735:       }
736:     }
737:     ISRestoreIndices(part, &partArray);
738:   }
739:   PetscMalloc1((*labelNew)->stratumOffsets[(*labelNew)->numStrata], &(*labelNew)->points);
740:   MPI_Scatterv(points, sendcnts, displs, MPIU_INT, (*labelNew)->points, (*labelNew)->stratumOffsets[(*labelNew)->numStrata], MPIU_INT, 0, comm);
741:   PetscFree(points);
742:   PetscFree3(sendcnts,offsets,displs);
743:   PetscFree(stratumSizes);
744:   /* Renumber points */
745:   ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, (*labelNew)->stratumOffsets[(*labelNew)->numStrata], (*labelNew)->points, NULL, (*labelNew)->points);
746:   /* Sort points */
747:   for (s = 0; s < (*labelNew)->numStrata; ++s) {PetscSortInt((*labelNew)->stratumSizes[s], &(*labelNew)->points[(*labelNew)->stratumOffsets[s]]);}
748:   return(0);
749: }


754: /*@C
755:   DMPlexCreateLabel - Create a label of the given name if it does not already exist

757:   Not Collective

759:   Input Parameters:
760: + dm   - The DMPlex object
761: - name - The label name

763:   Level: intermediate

765: .keywords: mesh
766: .seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
767: @*/
768: PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
769: {
770:   DM_Plex        *mesh = (DM_Plex*) dm->data;
771:   DMLabel        next  = mesh->labels;
772:   PetscBool      flg   = PETSC_FALSE;

778:   while (next) {
779:     PetscStrcmp(name, next->name, &flg);
780:     if (flg) break;
781:     next = next->next;
782:   }
783:   if (!flg) {
784:     DMLabel tmpLabel = mesh->labels;

786:     DMLabelCreate(name, &mesh->labels);

788:     mesh->labels->next = tmpLabel;
789:   }
790:   return(0);
791: }

795: /*@C
796:   DMPlexGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default

798:   Not Collective

800:   Input Parameters:
801: + dm   - The DMPlex object
802: . name - The label name
803: - point - The mesh point

805:   Output Parameter:
806: . value - The label value for this point, or -1 if the point is not in the label

808:   Level: beginner

810: .keywords: mesh
811: .seealso: DMLabelGetValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
812: @*/
813: PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
814: {
815:   DMLabel        label;

821:   DMPlexGetLabel(dm, name, &label);
822:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
823:   DMLabelGetValue(label, point, value);
824:   return(0);
825: }

829: /*@C
830:   DMPlexSetLabelValue - Add a point to a Sieve Label with given value

832:   Not Collective

834:   Input Parameters:
835: + dm   - The DMPlex object
836: . name - The label name
837: . point - The mesh point
838: - value - The label value for this point

840:   Output Parameter:

842:   Level: beginner

844: .keywords: mesh
845: .seealso: DMLabelSetValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
846: @*/
847: PetscErrorCode DMPlexSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
848: {
849:   DMLabel        label;

855:   DMPlexGetLabel(dm, name, &label);
856:   if (!label) {
857:     DMPlexCreateLabel(dm, name);
858:     DMPlexGetLabel(dm, name, &label);
859:   }
860:   DMLabelSetValue(label, point, value);
861:   return(0);
862: }

866: /*@C
867:   DMPlexClearLabelValue - Remove a point from a Sieve Label with given value

869:   Not Collective

871:   Input Parameters:
872: + dm   - The DMPlex object
873: . name - The label name
874: . point - The mesh point
875: - value - The label value for this point

877:   Output Parameter:

879:   Level: beginner

881: .keywords: mesh
882: .seealso: DMLabelClearValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
883: @*/
884: PetscErrorCode DMPlexClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
885: {
886:   DMLabel        label;

892:   DMPlexGetLabel(dm, name, &label);
893:   if (!label) return(0);
894:   DMLabelClearValue(label, point, value);
895:   return(0);
896: }

900: /*@C
901:   DMPlexGetLabelSize - Get the number of different integer ids in a Label

903:   Not Collective

905:   Input Parameters:
906: + dm   - The DMPlex object
907: - name - The label name

909:   Output Parameter:
910: . size - The number of different integer ids, or 0 if the label does not exist

912:   Level: beginner

914: .keywords: mesh
915: .seealso: DMLabeGetNumValues(), DMPlexSetLabelValue()
916: @*/
917: PetscErrorCode DMPlexGetLabelSize(DM dm, const char name[], PetscInt *size)
918: {
919:   DMLabel        label;

926:   DMPlexGetLabel(dm, name, &label);
927:   *size = 0;
928:   if (!label) return(0);
929:   DMLabelGetNumValues(label, size);
930:   return(0);
931: }

935: /*@C
936:   DMPlexGetLabelIdIS - Get the integer ids in a label

938:   Not Collective

940:   Input Parameters:
941: + mesh - The DMPlex object
942: - name - The label name

944:   Output Parameter:
945: . ids - The integer ids, or NULL if the label does not exist

947:   Level: beginner

949: .keywords: mesh
950: .seealso: DMLabelGetValueIS(), DMPlexGetLabelSize()
951: @*/
952: PetscErrorCode DMPlexGetLabelIdIS(DM dm, const char name[], IS *ids)
953: {
954:   DMLabel        label;

961:   DMPlexGetLabel(dm, name, &label);
962:   *ids = NULL;
963:   if (!label) return(0);
964:   DMLabelGetValueIS(label, ids);
965:   return(0);
966: }

970: /*@C
971:   DMPlexGetStratumSize - Get the number of points in a label stratum

973:   Not Collective

975:   Input Parameters:
976: + dm - The DMPlex object
977: . name - The label name
978: - value - The stratum value

980:   Output Parameter:
981: . size - The stratum size

983:   Level: beginner

985: .keywords: mesh
986: .seealso: DMLabelGetStratumSize(), DMPlexGetLabelSize(), DMPlexGetLabelIds()
987: @*/
988: PetscErrorCode DMPlexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
989: {
990:   DMLabel        label;

997:   DMPlexGetLabel(dm, name, &label);
998:   *size = 0;
999:   if (!label) return(0);
1000:   DMLabelGetStratumSize(label, value, size);
1001:   return(0);
1002: }

1006: /*@C
1007:   DMPlexGetStratumIS - Get the points in a label stratum

1009:   Not Collective

1011:   Input Parameters:
1012: + dm - The DMPlex object
1013: . name - The label name
1014: - value - The stratum value

1016:   Output Parameter:
1017: . points - The stratum points, or NULL if the label does not exist or does not have that value

1019:   Level: beginner

1021: .keywords: mesh
1022: .seealso: DMLabelGetStratumIS(), DMPlexGetStratumSize()
1023: @*/
1024: PetscErrorCode DMPlexGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
1025: {
1026:   DMLabel        label;

1033:   DMPlexGetLabel(dm, name, &label);
1034:   *points = NULL;
1035:   if (!label) return(0);
1036:   DMLabelGetStratumIS(label, value, points);
1037:   return(0);
1038: }

1042: /*@C
1043:   DMPlexClearLabelStratum - Remove all points from a stratum from a Sieve Label

1045:   Not Collective

1047:   Input Parameters:
1048: + dm   - The DMPlex object
1049: . name - The label name
1050: - value - The label value for this point

1052:   Output Parameter:

1054:   Level: beginner

1056: .keywords: mesh
1057: .seealso: DMLabelClearStratum(), DMPlexSetLabelValue(), DMPlexGetStratumIS(), DMPlexClearLabelValue()
1058: @*/
1059: PetscErrorCode DMPlexClearLabelStratum(DM dm, const char name[], PetscInt value)
1060: {
1061:   DMLabel        label;

1067:   DMPlexGetLabel(dm, name, &label);
1068:   if (!label) return(0);
1069:   DMLabelClearStratum(label, value);
1070:   return(0);
1071: }

1075: /*@
1076:   DMPlexGetNumLabels - Return the number of labels defined by the mesh

1078:   Not Collective

1080:   Input Parameter:
1081: . dm   - The DMPlex object

1083:   Output Parameter:
1084: . numLabels - the number of Labels

1086:   Level: intermediate

1088: .keywords: mesh
1089: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1090: @*/
1091: PetscErrorCode DMPlexGetNumLabels(DM dm, PetscInt *numLabels)
1092: {
1093:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1094:   DMLabel  next  = mesh->labels;
1095:   PetscInt n     = 0;

1100:   while (next) {
1101:     ++n;
1102:     next = next->next;
1103:   }
1104:   *numLabels = n;
1105:   return(0);
1106: }

1110: /*@C
1111:   DMPlexGetLabelName - Return the name of nth label

1113:   Not Collective

1115:   Input Parameters:
1116: + dm - The DMPlex object
1117: - n  - the label number

1119:   Output Parameter:
1120: . name - the label name

1122:   Level: intermediate

1124: .keywords: mesh
1125: .seealso: DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1126: @*/
1127: PetscErrorCode DMPlexGetLabelName(DM dm, PetscInt n, const char **name)
1128: {
1129:   DM_Plex  *mesh = (DM_Plex*) dm->data;
1130:   DMLabel  next  = mesh->labels;
1131:   PetscInt l     = 0;

1136:   while (next) {
1137:     if (l == n) {
1138:       *name = next->name;
1139:       return(0);
1140:     }
1141:     ++l;
1142:     next = next->next;
1143:   }
1144:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %d does not exist in this DM", n);
1145: }

1149: /*@C
1150:   DMPlexHasLabel - Determine whether the mesh has a label of a given name

1152:   Not Collective

1154:   Input Parameters:
1155: + dm   - The DMPlex object
1156: - name - The label name

1158:   Output Parameter:
1159: . hasLabel - PETSC_TRUE if the label is present

1161:   Level: intermediate

1163: .keywords: mesh
1164: .seealso: DMPlexCreateLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1165: @*/
1166: PetscErrorCode DMPlexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1167: {
1168:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1169:   DMLabel        next  = mesh->labels;

1176:   *hasLabel = PETSC_FALSE;
1177:   while (next) {
1178:     PetscStrcmp(name, next->name, hasLabel);
1179:     if (*hasLabel) break;
1180:     next = next->next;
1181:   }
1182:   return(0);
1183: }

1187: /*@C
1188:   DMPlexGetLabel - Return the label of a given name, or NULL

1190:   Not Collective

1192:   Input Parameters:
1193: + dm   - The DMPlex object
1194: - name - The label name

1196:   Output Parameter:
1197: . label - The DMLabel, or NULL if the label is absent

1199:   Level: intermediate

1201: .keywords: mesh
1202: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1203: @*/
1204: PetscErrorCode DMPlexGetLabel(DM dm, const char name[], DMLabel *label)
1205: {
1206:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1207:   DMLabel        next  = mesh->labels;
1208:   PetscBool      hasLabel;

1215:   *label = NULL;
1216:   while (next) {
1217:     PetscStrcmp(name, next->name, &hasLabel);
1218:     if (hasLabel) {
1219:       *label = next;
1220:       break;
1221:     }
1222:     next = next->next;
1223:   }
1224:   return(0);
1225: }

1229: /*@C
1230:   DMPlexAddLabel - Add the label to this mesh

1232:   Not Collective

1234:   Input Parameters:
1235: + dm   - The DMPlex object
1236: - label - The DMLabel

1238:   Level: developer

1240: .keywords: mesh
1241: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1242: @*/
1243: PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
1244: {
1245:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1246:   PetscBool      hasLabel;

1251:   DMPlexHasLabel(dm, label->name, &hasLabel);
1252:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
1253:   label->next  = mesh->labels;
1254:   mesh->labels = label;
1255:   return(0);
1256: }

1260: /*@C
1261:   DMPlexRemoveLabel - Remove the label from this mesh

1263:   Not Collective

1265:   Input Parameters:
1266: + dm   - The DMPlex object
1267: - name - The label name

1269:   Output Parameter:
1270: . label - The DMLabel, or NULL if the label is absent

1272:   Level: developer

1274: .keywords: mesh
1275: .seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
1276: @*/
1277: PetscErrorCode DMPlexRemoveLabel(DM dm, const char name[], DMLabel *label)
1278: {
1279:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1280:   DMLabel        next  = mesh->labels;
1281:   DMLabel        last  = NULL;
1282:   PetscBool      hasLabel;

1287:   DMPlexHasLabel(dm, name, &hasLabel);
1288:   *label = NULL;
1289:   if (!hasLabel) return(0);
1290:   while (next) {
1291:     PetscStrcmp(name, next->name, &hasLabel);
1292:     if (hasLabel) {
1293:       if (last) last->next   = next->next;
1294:       else      mesh->labels = next->next;
1295:       next->next = NULL;
1296:       *label     = next;
1297:       break;
1298:     }
1299:     last = next;
1300:     next = next->next;
1301:   }
1302:   return(0);
1303: }