Actual source code: dmlabel.c

  1: #include <petscdm.h>
  2: #include <petsc/private/dmlabelimpl.h>
  3: #include <petsc/private/sectionimpl.h>
  4: #include <petscsf.h>
  5: #include <petscsection.h>

  7: /*@C
  8:   DMLabelCreate - Create a DMLabel object, which is a multimap

 10:   Collective

 12:   Input parameters:
 13: + comm - The communicator, usually PETSC_COMM_SELF
 14: - name - The label name

 16:   Output parameter:
 17: . label - The DMLabel

 19:   Level: beginner

 21:   Notes:
 22:   The label name is actually usual PetscObject name.
 23:   One can get/set it with PetscObjectGetName()/PetscObjectSetName().

 25: .seealso: DMLabelDestroy()
 26: @*/
 27: PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
 28: {
 30:   DMInitializePackage();

 32:   PetscHeaderCreate(*label,DMLABEL_CLASSID,"DMLabel","DMLabel","DM",comm,DMLabelDestroy,DMLabelView);

 34:   (*label)->numStrata      = 0;
 35:   (*label)->defaultValue   = -1;
 36:   (*label)->stratumValues  = NULL;
 37:   (*label)->validIS        = NULL;
 38:   (*label)->stratumSizes   = NULL;
 39:   (*label)->points         = NULL;
 40:   (*label)->ht             = NULL;
 41:   (*label)->pStart         = -1;
 42:   (*label)->pEnd           = -1;
 43:   (*label)->bt             = NULL;
 44:   PetscHMapICreate(&(*label)->hmap);
 45:   PetscObjectSetName((PetscObject) *label, name);
 46:   return 0;
 47: }

 49: /*
 50:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 52:   Not collective

 54:   Input parameter:
 55: + label - The DMLabel
 56: - v - The stratum value

 58:   Output parameter:
 59: . label - The DMLabel with stratum in sorted list format

 61:   Level: developer

 63: .seealso: DMLabelCreate()
 64: */
 65: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 66: {
 67:   IS             is;
 68:   PetscInt       off = 0, *pointArray, p;

 70:   if (PetscLikely(v >= 0 && v < label->numStrata) && label->validIS[v]) return 0;
 72:   PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);
 73:   PetscMalloc1(label->stratumSizes[v], &pointArray);
 74:   PetscHSetIGetElems(label->ht[v], &off, pointArray);
 75:   PetscHSetIClear(label->ht[v]);
 76:   PetscSortInt(label->stratumSizes[v], pointArray);
 77:   if (label->bt) {
 78:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 79:       const PetscInt point = pointArray[p];
 81:       PetscBTSet(label->bt, point - label->pStart);
 82:     }
 83:   }
 84:   if (label->stratumSizes[v] > 0 && pointArray[label->stratumSizes[v]-1] == pointArray[0] + label->stratumSizes[v]-1) {
 85:     ISCreateStride(PETSC_COMM_SELF, label->stratumSizes[v], pointArray[0], 1, &is);
 86:     PetscFree(pointArray);
 87:   } else {
 88:     ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], pointArray, PETSC_OWN_POINTER, &is);
 89:   }
 90:   PetscObjectSetName((PetscObject) is, "indices");
 91:   label->points[v]  = is;
 92:   label->validIS[v] = PETSC_TRUE;
 93:   PetscObjectStateIncrease((PetscObject) label);
 94:   return 0;
 95: }

 97: /*
 98:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

100:   Not collective

102:   Input parameter:
103: . label - The DMLabel

105:   Output parameter:
106: . label - The DMLabel with all strata in sorted list format

108:   Level: developer

110: .seealso: DMLabelCreate()
111: */
112: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
113: {
114:   PetscInt       v;

116:   for (v = 0; v < label->numStrata; v++) {
117:     DMLabelMakeValid_Private(label, v);
118:   }
119:   return 0;
120: }

122: /*
123:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

125:   Not collective

127:   Input parameter:
128: + label - The DMLabel
129: - v - The stratum value

131:   Output parameter:
132: . label - The DMLabel with stratum in hash format

134:   Level: developer

136: .seealso: DMLabelCreate()
137: */
138: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
139: {
140:   PetscInt       p;
141:   const PetscInt *points;

143:   if (PetscLikely(v >= 0 && v < label->numStrata) && !label->validIS[v]) return 0;
145:   if (label->points[v]) {
146:     ISGetIndices(label->points[v], &points);
147:     for (p = 0; p < label->stratumSizes[v]; ++p) {
148:       PetscHSetIAdd(label->ht[v], points[p]);
149:     }
150:     ISRestoreIndices(label->points[v],&points);
151:     ISDestroy(&(label->points[v]));
152:   }
153:   label->validIS[v] = PETSC_FALSE;
154:   return 0;
155: }

157: #if !defined(DMLABEL_LOOKUP_THRESHOLD)
158: #define DMLABEL_LOOKUP_THRESHOLD 16
159: #endif

161: static inline PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
162: {
163:   PetscInt       v;

165:   *index = -1;
166:   if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
167:     for (v = 0; v < label->numStrata; ++v)
168:       if (label->stratumValues[v] == value) {*index = v; break;}
169:   } else {
170:     PetscHMapIGet(label->hmap, value, index);
171:   }
172:   if (PetscDefined(USE_DEBUG)) { /* Check strata hash map consistency */
173:     PetscInt len, loc = -1;
174:     PetscHMapIGetSize(label->hmap, &len);
176:     if (label->numStrata <= DMLABEL_LOOKUP_THRESHOLD) {
177:       PetscHMapIGet(label->hmap, value, &loc);
178:     } else {
179:       for (v = 0; v < label->numStrata; ++v)
180:         if (label->stratumValues[v] == value) {loc = v; break;}
181:     }
183:   }
184:   return 0;
185: }

187: static inline PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
188: {
189:   PetscInt       v;
190:   PetscInt      *tmpV;
191:   PetscInt      *tmpS;
192:   PetscHSetI    *tmpH, ht;
193:   IS            *tmpP, is;
194:   PetscBool     *tmpB;
195:   PetscHMapI     hmap = label->hmap;

197:   v    = label->numStrata;
198:   tmpV = label->stratumValues;
199:   tmpS = label->stratumSizes;
200:   tmpH = label->ht;
201:   tmpP = label->points;
202:   tmpB = label->validIS;
203:   { /* TODO: PetscRealloc() is broken, use malloc+memcpy+free  */
204:     PetscInt   *oldV = tmpV;
205:     PetscInt   *oldS = tmpS;
206:     PetscHSetI *oldH = tmpH;
207:     IS         *oldP = tmpP;
208:     PetscBool  *oldB = tmpB;
209:     PetscMalloc((v+1)*sizeof(*tmpV), &tmpV);
210:     PetscMalloc((v+1)*sizeof(*tmpS), &tmpS);
211:     PetscMalloc((v+1)*sizeof(*tmpH), &tmpH);
212:     PetscMalloc((v+1)*sizeof(*tmpP), &tmpP);
213:     PetscMalloc((v+1)*sizeof(*tmpB), &tmpB);
214:     PetscArraycpy(tmpV, oldV, v);
215:     PetscArraycpy(tmpS, oldS, v);
216:     PetscArraycpy(tmpH, oldH, v);
217:     PetscArraycpy(tmpP, oldP, v);
218:     PetscArraycpy(tmpB, oldB, v);
219:     PetscFree(oldV);
220:     PetscFree(oldS);
221:     PetscFree(oldH);
222:     PetscFree(oldP);
223:     PetscFree(oldB);
224:   }
225:   label->numStrata     = v+1;
226:   label->stratumValues = tmpV;
227:   label->stratumSizes  = tmpS;
228:   label->ht            = tmpH;
229:   label->points        = tmpP;
230:   label->validIS       = tmpB;
231:   PetscHSetICreate(&ht);
232:   ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
233:   PetscHMapISet(hmap, value, v);
234:   tmpV[v] = value;
235:   tmpS[v] = 0;
236:   tmpH[v] = ht;
237:   tmpP[v] = is;
238:   tmpB[v] = PETSC_TRUE;
239:   PetscObjectStateIncrease((PetscObject) label);
240:   *index = v;
241:   return 0;
242: }

244: static inline PetscErrorCode DMLabelLookupAddStratum(DMLabel label, PetscInt value, PetscInt *index)
245: {
246:   DMLabelLookupStratum(label, value, index);
247:   if (*index < 0) DMLabelNewStratum(label, value, index);
248:   return 0;
249: }

251: static inline PetscErrorCode DMLabelGetStratumSize_Private(DMLabel label, PetscInt v, PetscInt *size)
252: {
253:   *size = 0;
254:   if (v < 0) return 0;
255:   if (label->validIS[v]) {
256:     *size = label->stratumSizes[v];
257:   } else {
258:     PetscHSetIGetSize(label->ht[v], size);
259:   }
260:   return 0;
261: }

263: /*@
264:   DMLabelAddStratum - Adds a new stratum value in a DMLabel

266:   Input Parameters:
267: + label - The DMLabel
268: - value - The stratum value

270:   Level: beginner

272: .seealso:  DMLabelCreate(), DMLabelDestroy()
273: @*/
274: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
275: {
276:   PetscInt       v;

279:   DMLabelLookupAddStratum(label, value, &v);
280:   return 0;
281: }

283: /*@
284:   DMLabelAddStrata - Adds new stratum values in a DMLabel

286:   Not collective

288:   Input Parameters:
289: + label - The DMLabel
290: . numStrata - The number of stratum values
291: - stratumValues - The stratum values

293:   Level: beginner

295: .seealso:  DMLabelCreate(), DMLabelDestroy()
296: @*/
297: PetscErrorCode DMLabelAddStrata(DMLabel label, PetscInt numStrata, const PetscInt stratumValues[])
298: {
299:   PetscInt       *values, v;

303:   PetscMalloc1(numStrata, &values);
304:   PetscArraycpy(values, stratumValues, numStrata);
305:   PetscSortRemoveDupsInt(&numStrata, values);
306:   if (!label->numStrata) { /* Fast preallocation */
307:     PetscInt   *tmpV;
308:     PetscInt   *tmpS;
309:     PetscHSetI *tmpH, ht;
310:     IS         *tmpP, is;
311:     PetscBool  *tmpB;
312:     PetscHMapI  hmap = label->hmap;

314:     PetscMalloc1(numStrata, &tmpV);
315:     PetscMalloc1(numStrata, &tmpS);
316:     PetscMalloc1(numStrata, &tmpH);
317:     PetscMalloc1(numStrata, &tmpP);
318:     PetscMalloc1(numStrata, &tmpB);
319:     label->numStrata     = numStrata;
320:     label->stratumValues = tmpV;
321:     label->stratumSizes  = tmpS;
322:     label->ht            = tmpH;
323:     label->points        = tmpP;
324:     label->validIS       = tmpB;
325:     for (v = 0; v < numStrata; ++v) {
326:       PetscHSetICreate(&ht);
327:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&is);
328:       PetscHMapISet(hmap, values[v], v);
329:       tmpV[v] = values[v];
330:       tmpS[v] = 0;
331:       tmpH[v] = ht;
332:       tmpP[v] = is;
333:       tmpB[v] = PETSC_TRUE;
334:     }
335:     PetscObjectStateIncrease((PetscObject) label);
336:   } else {
337:     for (v = 0; v < numStrata; ++v) {
338:       DMLabelAddStratum(label, values[v]);
339:     }
340:   }
341:   PetscFree(values);
342:   return 0;
343: }

345: /*@
346:   DMLabelAddStrataIS - Adds new stratum values in a DMLabel

348:   Not collective

350:   Input Parameters:
351: + label - The DMLabel
352: - valueIS - Index set with stratum values

354:   Level: beginner

356: .seealso:  DMLabelCreate(), DMLabelDestroy()
357: @*/
358: PetscErrorCode DMLabelAddStrataIS(DMLabel label, IS valueIS)
359: {
360:   PetscInt       numStrata;
361:   const PetscInt *stratumValues;

365:   ISGetLocalSize(valueIS, &numStrata);
366:   ISGetIndices(valueIS, &stratumValues);
367:   DMLabelAddStrata(label, numStrata, stratumValues);
368:   return 0;
369: }

371: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
372: {
373:   PetscInt       v;
374:   PetscMPIInt    rank;

376:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
377:   PetscViewerASCIIPushSynchronized(viewer);
378:   if (label) {
379:     const char *name;

381:     PetscObjectGetName((PetscObject) label, &name);
382:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
383:     if (label->bt) PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);
384:     for (v = 0; v < label->numStrata; ++v) {
385:       const PetscInt value = label->stratumValues[v];
386:       const PetscInt *points;
387:       PetscInt       p;

389:       ISGetIndices(label->points[v], &points);
390:       for (p = 0; p < label->stratumSizes[v]; ++p) {
391:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
392:       }
393:       ISRestoreIndices(label->points[v],&points);
394:     }
395:   }
396:   PetscViewerFlush(viewer);
397:   PetscViewerASCIIPopSynchronized(viewer);
398:   return 0;
399: }

401: /*@C
402:   DMLabelView - View the label

404:   Collective on viewer

406:   Input Parameters:
407: + label - The DMLabel
408: - viewer - The PetscViewer

410:   Level: intermediate

412: .seealso: DMLabelCreate(), DMLabelDestroy()
413: @*/
414: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
415: {
416:   PetscBool      iascii;

419:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)label), &viewer);
421:   if (label) DMLabelMakeAllValid_Private(label);
422:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
423:   if (iascii) {
424:     DMLabelView_Ascii(label, viewer);
425:   }
426:   return 0;
427: }

429: /*@
430:   DMLabelReset - Destroys internal data structures in a DMLabel

432:   Not collective

434:   Input Parameter:
435: . label - The DMLabel

437:   Level: beginner

439: .seealso: DMLabelDestroy(), DMLabelCreate()
440: @*/
441: PetscErrorCode DMLabelReset(DMLabel label)
442: {
443:   PetscInt       v;

446:   for (v = 0; v < label->numStrata; ++v) {
447:     PetscHSetIDestroy(&label->ht[v]);
448:     ISDestroy(&label->points[v]);
449:   }
450:   label->numStrata = 0;
451:   PetscFree(label->stratumValues);
452:   PetscFree(label->stratumSizes);
453:   PetscFree(label->ht);
454:   PetscFree(label->points);
455:   PetscFree(label->validIS);
456:   label->stratumValues = NULL;
457:   label->stratumSizes  = NULL;
458:   label->ht            = NULL;
459:   label->points        = NULL;
460:   label->validIS       = NULL;
461:   PetscHMapIReset(label->hmap);
462:   label->pStart = -1;
463:   label->pEnd   = -1;
464:   PetscBTDestroy(&label->bt);
465:   return 0;
466: }

468: /*@
469:   DMLabelDestroy - Destroys a DMLabel

471:   Collective on label

473:   Input Parameter:
474: . label - The DMLabel

476:   Level: beginner

478: .seealso: DMLabelReset(), DMLabelCreate()
479: @*/
480: PetscErrorCode DMLabelDestroy(DMLabel *label)
481: {
482:   if (!*label) return 0;
484:   if (--((PetscObject)(*label))->refct > 0) {*label = NULL; return 0;}
485:   DMLabelReset(*label);
486:   PetscHMapIDestroy(&(*label)->hmap);
487:   PetscHeaderDestroy(label);
488:   return 0;
489: }

491: /*@
492:   DMLabelDuplicate - Duplicates a DMLabel

494:   Collective on label

496:   Input Parameter:
497: . label - The DMLabel

499:   Output Parameter:
500: . labelnew - location to put new vector

502:   Level: intermediate

504: .seealso: DMLabelCreate(), DMLabelDestroy()
505: @*/
506: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
507: {
508:   const char    *name;
509:   PetscInt       v;

512:   DMLabelMakeAllValid_Private(label);
513:   PetscObjectGetName((PetscObject) label, &name);
514:   DMLabelCreate(PetscObjectComm((PetscObject) label), name, labelnew);

516:   (*labelnew)->numStrata    = label->numStrata;
517:   (*labelnew)->defaultValue = label->defaultValue;
518:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
519:   PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
520:   PetscMalloc1(label->numStrata, &(*labelnew)->ht);
521:   PetscMalloc1(label->numStrata, &(*labelnew)->points);
522:   PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
523:   for (v = 0; v < label->numStrata; ++v) {
524:     PetscHSetICreate(&(*labelnew)->ht[v]);
525:     (*labelnew)->stratumValues[v]  = label->stratumValues[v];
526:     (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
527:     PetscObjectReference((PetscObject) (label->points[v]));
528:     (*labelnew)->points[v]         = label->points[v];
529:     (*labelnew)->validIS[v]        = PETSC_TRUE;
530:   }
531:   PetscHMapIDestroy(&(*labelnew)->hmap);
532:   PetscHMapIDuplicate(label->hmap,&(*labelnew)->hmap);
533:   (*labelnew)->pStart = -1;
534:   (*labelnew)->pEnd   = -1;
535:   (*labelnew)->bt     = NULL;
536:   return 0;
537: }

539: /*@C
540:   DMLabelCompare - Compare two DMLabel objects

542:   Collective on comm

544:   Input Parameters:
545: + comm - Comm over which to compare labels
546: . l0 - First DMLabel
547: - l1 - Second DMLabel

549:   Output Parameters
550: + equal   - (Optional) Flag whether the two labels are equal
551: - message - (Optional) Message describing the difference

553:   Level: intermediate

555:   Notes:
556:   The output flag equal is the same on all processes.
557:   If it is passed as NULL and difference is found, an error is thrown on all processes.
558:   Make sure to pass NULL on all processes.

560:   The output message is set independently on each rank.
561:   It is set to NULL if no difference was found on the current rank. It must be freed by user.
562:   If message is passed as NULL and difference is found, the difference description is printed to stderr in synchronized manner.
563:   Make sure to pass NULL on all processes.

565:   For the comparison, we ignore the order of stratum values, and strata with no points.

567:   The communicator needs to be specified because currently DMLabel can live on PETSC_COMM_SELF even if the underlying DM is parallel.

569:   Fortran Notes:
570:   This function is currently not available from Fortran.

572: .seealso: DMCompareLabels(), DMLabelGetNumValues(), DMLabelGetDefaultValue(), DMLabelGetNonEmptyStratumValuesIS(), DMLabelGetStratumIS()
573: @*/
574: PetscErrorCode DMLabelCompare(MPI_Comm comm, DMLabel l0, DMLabel l1, PetscBool *equal, char **message)
575: {
576:   const char     *name0, *name1;
577:   char            msg[PETSC_MAX_PATH_LEN] = "";
578:   PetscBool       eq;
579:   PetscMPIInt     rank;

585:   MPI_Comm_rank(comm, &rank);
586:   PetscObjectGetName((PetscObject)l0, &name0);
587:   PetscObjectGetName((PetscObject)l1, &name1);
588:   {
589:     PetscInt v0, v1;

591:     DMLabelGetDefaultValue(l0, &v0);
592:     DMLabelGetDefaultValue(l1, &v1);
593:     eq = (PetscBool) (v0 == v1);
594:     if (!eq) {
595:       PetscSNPrintf(msg, sizeof(msg), "Default value of DMLabel l0 \"%s\" = %D != %D = Default value of DMLabel l1 \"%s\"", name0, v0, v1, name1);
596:     }
597:     MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
598:     if (!eq) goto finish;
599:   }
600:   {
601:     IS              is0, is1;

603:     DMLabelGetNonEmptyStratumValuesIS(l0, &is0);
604:     DMLabelGetNonEmptyStratumValuesIS(l1, &is1);
605:     ISEqual(is0, is1, &eq);
606:     ISDestroy(&is0);
607:     ISDestroy(&is1);
608:     if (!eq) {
609:       PetscSNPrintf(msg, sizeof(msg), "Stratum values in DMLabel l0 \"%s\" are different than in DMLabel l1 \"%s\"", name0, name1);
610:     }
611:     MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
612:     if (!eq) goto finish;
613:   }
614:   {
615:     PetscInt i, nValues;

617:     DMLabelGetNumValues(l0, &nValues);
618:     for (i=0; i<nValues; i++) {
619:       const PetscInt  v = l0->stratumValues[i];
620:       PetscInt        n;
621:       IS              is0, is1;

623:       DMLabelGetStratumSize_Private(l0, i, &n);
624:       if (!n) continue;
625:       DMLabelGetStratumIS(l0, v, &is0);
626:       DMLabelGetStratumIS(l1, v, &is1);
627:       ISEqualUnsorted(is0, is1, &eq);
628:       ISDestroy(&is0);
629:       ISDestroy(&is1);
630:       if (!eq) {
631:         PetscSNPrintf(msg, sizeof(msg), "Stratum #%D with value %D contains different points in DMLabel l0 \"%s\" and DMLabel l1 \"%s\"", i, v, name0, name1);
632:         break;
633:       }
634:     }
635:     MPI_Allreduce(MPI_IN_PLACE, &eq, 1, MPIU_BOOL, MPI_LAND, comm);
636:   }
637: finish:
638:   /* If message output arg not set, print to stderr */
639:   if (message) {
640:     *message = NULL;
641:     if (msg[0]) {
642:       PetscStrallocpy(msg, message);
643:     }
644:   } else {
645:     if (msg[0]) {
646:       PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] %s\n", rank, msg);
647:     }
648:     PetscSynchronizedFlush(comm, PETSC_STDERR);
649:   }
650:   /* If same output arg not ser and labels are not equal, throw error */
651:   if (equal) *equal = eq;
653:   return 0;
654: }

656: /*@
657:   DMLabelComputeIndex - Create an index structure for membership determination, automatically determining the bounds

659:   Not collective

661:   Input Parameter:
662: . label  - The DMLabel

664:   Level: intermediate

666: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
667: @*/
668: PetscErrorCode DMLabelComputeIndex(DMLabel label)
669: {
670:   PetscInt       pStart = PETSC_MAX_INT, pEnd = -1, v;

673:   DMLabelMakeAllValid_Private(label);
674:   for (v = 0; v < label->numStrata; ++v) {
675:     const PetscInt *points;
676:     PetscInt       i;

678:     ISGetIndices(label->points[v], &points);
679:     for (i = 0; i < label->stratumSizes[v]; ++i) {
680:       const PetscInt point = points[i];

682:       pStart = PetscMin(point,   pStart);
683:       pEnd   = PetscMax(point+1, pEnd);
684:     }
685:     ISRestoreIndices(label->points[v], &points);
686:   }
687:   label->pStart = pStart == PETSC_MAX_INT ? -1 : pStart;
688:   label->pEnd   = pEnd;
689:   DMLabelCreateIndex(label, label->pStart, label->pEnd);
690:   return 0;
691: }

693: /*@
694:   DMLabelCreateIndex - Create an index structure for membership determination

696:   Not collective

698:   Input Parameters:
699: + label  - The DMLabel
700: . pStart - The smallest point
701: - pEnd   - The largest point + 1

703:   Level: intermediate

705: .seealso: DMLabelHasPoint(), DMLabelComputeIndex(), DMLabelDestroyIndex(), DMLabelGetValue(), DMLabelSetValue()
706: @*/
707: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
708: {
709:   PetscInt       v;

712:   DMLabelDestroyIndex(label);
713:   DMLabelMakeAllValid_Private(label);
714:   label->pStart = pStart;
715:   label->pEnd   = pEnd;
716:   /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
717:   PetscBTCreate(pEnd - pStart, &label->bt);
718:   for (v = 0; v < label->numStrata; ++v) {
719:     const PetscInt *points;
720:     PetscInt       i;

722:     ISGetIndices(label->points[v], &points);
723:     for (i = 0; i < label->stratumSizes[v]; ++i) {
724:       const PetscInt point = points[i];

727:       PetscBTSet(label->bt, point - pStart);
728:     }
729:     ISRestoreIndices(label->points[v], &points);
730:   }
731:   return 0;
732: }

734: /*@
735:   DMLabelDestroyIndex - Destroy the index structure

737:   Not collective

739:   Input Parameter:
740: . label - the DMLabel

742:   Level: intermediate

744: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
745: @*/
746: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
747: {
749:   label->pStart = -1;
750:   label->pEnd   = -1;
751:   PetscBTDestroy(&label->bt);
752:   return 0;
753: }

755: /*@
756:   DMLabelGetBounds - Return the smallest and largest point in the label

758:   Not collective

760:   Input Parameter:
761: . label - the DMLabel

763:   Output Parameters:
764: + pStart - The smallest point
765: - pEnd   - The largest point + 1

767:   Note: This will compute an index for the label if one does not exist.

769:   Level: intermediate

771: .seealso: DMLabelHasPoint(), DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
772: @*/
773: PetscErrorCode DMLabelGetBounds(DMLabel label, PetscInt *pStart, PetscInt *pEnd)
774: {
776:   if ((label->pStart == -1) && (label->pEnd == -1)) DMLabelComputeIndex(label);
777:   if (pStart) {
779:     *pStart = label->pStart;
780:   }
781:   if (pEnd) {
783:     *pEnd = label->pEnd;
784:   }
785:   return 0;
786: }

788: /*@
789:   DMLabelHasValue - Determine whether a label assigns the value to any point

791:   Not collective

793:   Input Parameters:
794: + label - the DMLabel
795: - value - the value

797:   Output Parameter:
798: . contains - Flag indicating whether the label maps this value to any point

800:   Level: developer

802: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
803: @*/
804: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
805: {
806:   PetscInt v;

810:   DMLabelLookupStratum(label, value, &v);
811:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
812:   return 0;
813: }

815: /*@
816:   DMLabelHasPoint - Determine whether a label assigns a value to a point

818:   Not collective

820:   Input Parameters:
821: + label - the DMLabel
822: - point - the point

824:   Output Parameter:
825: . contains - Flag indicating whether the label maps this point to a value

827:   Note: The user must call DMLabelCreateIndex() before this function.

829:   Level: developer

831: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
832: @*/
833: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
834: {
838:   DMLabelMakeAllValid_Private(label);
839:   if (PetscDefined(USE_DEBUG)) {
842:   }
843:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
844:   return 0;
845: }

847: /*@
848:   DMLabelStratumHasPoint - Return true if the stratum contains a point

850:   Not collective

852:   Input Parameters:
853: + label - the DMLabel
854: . value - the stratum value
855: - point - the point

857:   Output Parameter:
858: . contains - true if the stratum contains the point

860:   Level: intermediate

862: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
863: @*/
864: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
865: {
866:   PetscInt       v;

871:   *contains = PETSC_FALSE;
872:   DMLabelLookupStratum(label, value, &v);
873:   if (v < 0) return 0;

875:   if (label->validIS[v]) {
876:     PetscInt i;

878:     ISLocate(label->points[v], point, &i);
879:     if (i >= 0) *contains = PETSC_TRUE;
880:   } else {
881:     PetscBool has;

883:     PetscHSetIHas(label->ht[v], point, &has);
884:     if (has) *contains = PETSC_TRUE;
885:   }
886:   return 0;
887: }

889: /*@
890:   DMLabelGetDefaultValue - Get the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
891:   When a label is created, it is initialized to -1.

893:   Not collective

895:   Input parameter:
896: . label - a DMLabel object

898:   Output parameter:
899: . defaultValue - the default value

901:   Level: beginner

903: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
904: @*/
905: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
906: {
908:   *defaultValue = label->defaultValue;
909:   return 0;
910: }

912: /*@
913:   DMLabelSetDefaultValue - Set the default value returned by DMLabelGetValue() if a point has not been explicitly given a value.
914:   When a label is created, it is initialized to -1.

916:   Not collective

918:   Input parameter:
919: . label - a DMLabel object

921:   Output parameter:
922: . defaultValue - the default value

924:   Level: beginner

926: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
927: @*/
928: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
929: {
931:   label->defaultValue = defaultValue;
932:   return 0;
933: }

935: /*@
936:   DMLabelGetValue - Return the value a label assigns to a point, or the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue())

938:   Not collective

940:   Input Parameters:
941: + label - the DMLabel
942: - point - the point

944:   Output Parameter:
945: . value - The point value, or the default value (-1 by default)

947:   Level: intermediate

949: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
950: @*/
951: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
952: {
953:   PetscInt       v;

958:   *value = label->defaultValue;
959:   for (v = 0; v < label->numStrata; ++v) {
960:     if (label->validIS[v]) {
961:       PetscInt i;

963:       ISLocate(label->points[v], point, &i);
964:       if (i >= 0) {
965:         *value = label->stratumValues[v];
966:         break;
967:       }
968:     } else {
969:       PetscBool has;

971:       PetscHSetIHas(label->ht[v], point, &has);
972:       if (has) {
973:         *value = label->stratumValues[v];
974:         break;
975:       }
976:     }
977:   }
978:   return 0;
979: }

981: /*@
982:   DMLabelSetValue - Set the value a label assigns to a point.  If the value is the same as the label's default value (which is initially -1, and can be changed with DMLabelSetDefaultValue() to something different), then this function will do nothing.

984:   Not collective

986:   Input Parameters:
987: + label - the DMLabel
988: . point - the point
989: - value - The point value

991:   Level: intermediate

993: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
994: @*/
995: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
996: {
997:   PetscInt       v;

1000:   /* Find label value, add new entry if needed */
1001:   if (value == label->defaultValue) return 0;
1002:   DMLabelLookupAddStratum(label, value, &v);
1003:   /* Set key */
1004:   DMLabelMakeInvalid_Private(label, v);
1005:   PetscHSetIAdd(label->ht[v], point);
1006:   return 0;
1007: }

1009: /*@
1010:   DMLabelClearValue - Clear the value a label assigns to a point

1012:   Not collective

1014:   Input Parameters:
1015: + label - the DMLabel
1016: . point - the point
1017: - value - The point value

1019:   Level: intermediate

1021: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
1022: @*/
1023: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
1024: {
1025:   PetscInt       v;

1028:   /* Find label value */
1029:   DMLabelLookupStratum(label, value, &v);
1030:   if (v < 0) return 0;

1032:   if (label->bt) {
1034:     PetscBTClear(label->bt, point - label->pStart);
1035:   }

1037:   /* Delete key */
1038:   DMLabelMakeInvalid_Private(label, v);
1039:   PetscHSetIDel(label->ht[v], point);
1040:   return 0;
1041: }

1043: /*@
1044:   DMLabelInsertIS - Set all points in the IS to a value

1046:   Not collective

1048:   Input Parameters:
1049: + label - the DMLabel
1050: . is    - the point IS
1051: - value - The point value

1053:   Level: intermediate

1055: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1056: @*/
1057: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
1058: {
1059:   PetscInt        v, n, p;
1060:   const PetscInt *points;

1064:   /* Find label value, add new entry if needed */
1065:   if (value == label->defaultValue) return 0;
1066:   DMLabelLookupAddStratum(label, value, &v);
1067:   /* Set keys */
1068:   DMLabelMakeInvalid_Private(label, v);
1069:   ISGetLocalSize(is, &n);
1070:   ISGetIndices(is, &points);
1071:   for (p = 0; p < n; ++p) PetscHSetIAdd(label->ht[v], points[p]);
1072:   ISRestoreIndices(is, &points);
1073:   return 0;
1074: }

1076: /*@
1077:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

1079:   Not collective

1081:   Input Parameter:
1082: . label - the DMLabel

1084:   Output Parameter:
1085: . numValues - the number of values

1087:   Level: intermediate

1089: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1090: @*/
1091: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
1092: {
1095:   *numValues = label->numStrata;
1096:   return 0;
1097: }

1099: /*@
1100:   DMLabelGetValueIS - Get an IS of all values that the DMlabel takes

1102:   Not collective

1104:   Input Parameter:
1105: . label - the DMLabel

1107:   Output Parameter:
1108: . is    - the value IS

1110:   Level: intermediate

1112:   Notes:
1113:   The output IS should be destroyed when no longer needed.
1114:   Strata which are allocated but empty [DMLabelGetStratumSize() yields 0] are counted.
1115:   If you need to count only nonempty strata, use DMLabelGetNonEmptyStratumValuesIS().

1117: .seealso: DMLabelGetNonEmptyStratumValuesIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1118: @*/
1119: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
1120: {
1123:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1124:   return 0;
1125: }

1127: /*@
1128:   DMLabelGetNonEmptyStratumValuesIS - Get an IS of all values that the DMlabel takes

1130:   Not collective

1132:   Input Parameter:
1133: . label - the DMLabel

1135:   Output Paramater:
1136: . is    - the value IS

1138:   Level: intermediate

1140:   Notes:
1141:   The output IS should be destroyed when no longer needed.
1142:   This is similar to DMLabelGetValueIS() but counts only nonempty strata.

1144: .seealso: DMLabelGetValueIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1145: @*/
1146: PetscErrorCode DMLabelGetNonEmptyStratumValuesIS(DMLabel label, IS *values)
1147: {
1148:   PetscInt        i, j;
1149:   PetscInt       *valuesArr;

1153:   PetscMalloc1(label->numStrata, &valuesArr);
1154:   for (i = 0, j = 0; i < label->numStrata; i++) {
1155:     PetscInt        n;

1157:     DMLabelGetStratumSize_Private(label, i, &n);
1158:     if (n) valuesArr[j++] = label->stratumValues[i];
1159:   }
1160:   if (j == label->numStrata) {
1161:     ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
1162:   } else {
1163:     ISCreateGeneral(PETSC_COMM_SELF, j, valuesArr, PETSC_COPY_VALUES, values);
1164:   }
1165:   PetscFree(valuesArr);
1166:   return 0;
1167: }

1169: /*@
1170:   DMLabelGetValueIndex - Get the index of a given value in the list of values for the DMlabel, or -1 if it is not present

1172:   Not collective

1174:   Input Parameters:
1175: + label - the DMLabel
1176: - value - the value

1178:   Output Parameter:
1179: . index - the index of value in the list of values

1181:   Level: intermediate

1183: .seealso: DMLabelGetValueIS(), DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1184: @*/
1185: PetscErrorCode DMLabelGetValueIndex(DMLabel label, PetscInt value, PetscInt *index)
1186: {
1187:   PetscInt v;

1191:   /* Do not assume they are sorted */
1192:   for (v = 0; v < label->numStrata; ++v) if (label->stratumValues[v] == value) break;
1193:   if (v >= label->numStrata) *index = -1;
1194:   else                       *index = v;
1195:   return 0;
1196: }

1198: /*@
1199:   DMLabelHasStratum - Determine whether points exist with the given value

1201:   Not collective

1203:   Input Parameters:
1204: + label - the DMLabel
1205: - value - the stratum value

1207:   Output Parameter:
1208: . exists - Flag saying whether points exist

1210:   Level: intermediate

1212: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1213: @*/
1214: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
1215: {
1216:   PetscInt       v;

1220:   DMLabelLookupStratum(label, value, &v);
1221:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
1222:   return 0;
1223: }

1225: /*@
1226:   DMLabelGetStratumSize - Get the size of a stratum

1228:   Not collective

1230:   Input Parameters:
1231: + label - the DMLabel
1232: - value - the stratum value

1234:   Output Parameter:
1235: . size - The number of points in the stratum

1237:   Level: intermediate

1239: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1240: @*/
1241: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
1242: {
1243:   PetscInt       v;

1247:   DMLabelLookupStratum(label, value, &v);
1248:   DMLabelGetStratumSize_Private(label, v, size);
1249:   return 0;
1250: }

1252: /*@
1253:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

1255:   Not collective

1257:   Input Parameters:
1258: + label - the DMLabel
1259: - value - the stratum value

1261:   Output Parameters:
1262: + start - the smallest point in the stratum
1263: - end - the largest point in the stratum

1265:   Level: intermediate

1267: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1268: @*/
1269: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
1270: {
1271:   PetscInt       v, min, max;

1276:   DMLabelLookupStratum(label, value, &v);
1277:   if (v < 0) return 0;
1278:   DMLabelMakeValid_Private(label, v);
1279:   if (label->stratumSizes[v] <= 0) return 0;
1280:   ISGetMinMax(label->points[v], &min, &max);
1281:   if (start) *start = min;
1282:   if (end)   *end   = max+1;
1283:   return 0;
1284: }

1286: /*@
1287:   DMLabelGetStratumIS - Get an IS with the stratum points

1289:   Not collective

1291:   Input Parameters:
1292: + label - the DMLabel
1293: - value - the stratum value

1295:   Output Parameter:
1296: . points - The stratum points

1298:   Level: intermediate

1300:   Notes:
1301:   The output IS should be destroyed when no longer needed.
1302:   Returns NULL if the stratum is empty.

1304: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1305: @*/
1306: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
1307: {
1308:   PetscInt       v;

1312:   *points = NULL;
1313:   DMLabelLookupStratum(label, value, &v);
1314:   if (v < 0) return 0;
1315:   DMLabelMakeValid_Private(label, v);
1316:   PetscObjectReference((PetscObject) label->points[v]);
1317:   *points = label->points[v];
1318:   return 0;
1319: }

1321: /*@
1322:   DMLabelSetStratumIS - Set the stratum points using an IS

1324:   Not collective

1326:   Input Parameters:
1327: + label - the DMLabel
1328: . value - the stratum value
1329: - points - The stratum points

1331:   Level: intermediate

1333: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1334: @*/
1335: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
1336: {
1337:   PetscInt       v;

1341:   DMLabelLookupAddStratum(label, value, &v);
1342:   if (is == label->points[v]) return 0;
1343:   DMLabelClearStratum(label, value);
1344:   ISGetLocalSize(is, &(label->stratumSizes[v]));
1345:   PetscObjectReference((PetscObject)is);
1346:   ISDestroy(&(label->points[v]));
1347:   label->points[v]  = is;
1348:   label->validIS[v] = PETSC_TRUE;
1349:   PetscObjectStateIncrease((PetscObject) label);
1350:   if (label->bt) {
1351:     const PetscInt *points;
1352:     PetscInt p;

1354:     ISGetIndices(is,&points);
1355:     for (p = 0; p < label->stratumSizes[v]; ++p) {
1356:       const PetscInt point = points[p];

1359:       PetscBTSet(label->bt, point - label->pStart);
1360:     }
1361:   }
1362:   return 0;
1363: }

1365: /*@
1366:   DMLabelClearStratum - Remove a stratum

1368:   Not collective

1370:   Input Parameters:
1371: + label - the DMLabel
1372: - value - the stratum value

1374:   Level: intermediate

1376: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1377: @*/
1378: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
1379: {
1380:   PetscInt       v;

1383:   DMLabelLookupStratum(label, value, &v);
1384:   if (v < 0) return 0;
1385:   if (label->validIS[v]) {
1386:     if (label->bt) {
1387:       PetscInt       i;
1388:       const PetscInt *points;

1390:       ISGetIndices(label->points[v], &points);
1391:       for (i = 0; i < label->stratumSizes[v]; ++i) {
1392:         const PetscInt point = points[i];

1395:         PetscBTClear(label->bt, point - label->pStart);
1396:       }
1397:       ISRestoreIndices(label->points[v], &points);
1398:     }
1399:     label->stratumSizes[v] = 0;
1400:     ISDestroy(&label->points[v]);
1401:     ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &label->points[v]);
1402:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1403:     PetscObjectStateIncrease((PetscObject) label);
1404:   } else {
1405:     PetscHSetIClear(label->ht[v]);
1406:   }
1407:   return 0;
1408: }

1410: /*@
1411:   DMLabelSetStratumBounds - Efficiently give a contiguous set of points a given label value

1413:   Not collective

1415:   Input Parameters:
1416: + label  - The DMLabel
1417: . value  - The label value for all points
1418: . pStart - The first point
1419: - pEnd   - A point beyond all marked points

1421:   Note: The marks points are [pStart, pEnd), and only the bounds are stored.

1423:   Level: intermediate

1425: .seealso: DMLabelCreate(), DMLabelSetStratumIS(), DMLabelGetStratumIS()
1426: @*/
1427: PetscErrorCode DMLabelSetStratumBounds(DMLabel label, PetscInt value, PetscInt pStart, PetscInt pEnd)
1428: {
1429:   IS             pIS;

1431:   ISCreateStride(PETSC_COMM_SELF, pEnd - pStart, pStart, 1, &pIS);
1432:   DMLabelSetStratumIS(label, value, pIS);
1433:   ISDestroy(&pIS);
1434:   return 0;
1435: }

1437: /*@
1438:   DMLabelGetStratumPointIndex - Get the index of a point in a given stratum

1440:   Not collective

1442:   Input Parameters:
1443: + label  - The DMLabel
1444: . value  - The label value
1445: - p      - A point with this value

1447:   Output Parameter:
1448: . index  - The index of this point in the stratum, or -1 if the point is not in the stratum or the stratum does not exist

1450:   Level: intermediate

1452: .seealso: DMLabelGetValueIndex(), DMLabelGetStratumIS(), DMLabelCreate()
1453: @*/
1454: PetscErrorCode DMLabelGetStratumPointIndex(DMLabel label, PetscInt value, PetscInt p, PetscInt *index)
1455: {
1456:   const PetscInt *indices;
1457:   PetscInt        v;

1461:   *index = -1;
1462:   DMLabelLookupStratum(label, value, &v);
1463:   if (v < 0) return 0;
1464:   DMLabelMakeValid_Private(label, v);
1465:   ISGetIndices(label->points[v], &indices);
1466:   PetscFindInt(p, label->stratumSizes[v], indices, index);
1467:   ISRestoreIndices(label->points[v], &indices);
1468:   return 0;
1469: }

1471: /*@
1472:   DMLabelFilter - Remove all points outside of [start, end)

1474:   Not collective

1476:   Input Parameters:
1477: + label - the DMLabel
1478: . start - the first point kept
1479: - end - one more than the last point kept

1481:   Level: intermediate

1483: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1484: @*/
1485: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
1486: {
1487:   PetscInt       v;

1490:   DMLabelDestroyIndex(label);
1491:   DMLabelMakeAllValid_Private(label);
1492:   for (v = 0; v < label->numStrata; ++v) {
1493:     ISGeneralFilter(label->points[v], start, end);
1494:   }
1495:   DMLabelCreateIndex(label, start, end);
1496:   return 0;
1497: }

1499: /*@
1500:   DMLabelPermute - Create a new label with permuted points

1502:   Not collective

1504:   Input Parameters:
1505: + label - the DMLabel
1506: - permutation - the point permutation

1508:   Output Parameter:
1509: . labelnew - the new label containing the permuted points

1511:   Level: intermediate

1513: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1514: @*/
1515: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1516: {
1517:   const PetscInt *perm;
1518:   PetscInt        numValues, numPoints, v, q;

1522:   DMLabelMakeAllValid_Private(label);
1523:   DMLabelDuplicate(label, labelNew);
1524:   DMLabelGetNumValues(*labelNew, &numValues);
1525:   ISGetLocalSize(permutation, &numPoints);
1526:   ISGetIndices(permutation, &perm);
1527:   for (v = 0; v < numValues; ++v) {
1528:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1529:     const PetscInt *points;
1530:     PetscInt *pointsNew;

1532:     ISGetIndices((*labelNew)->points[v],&points);
1533:     PetscMalloc1(size,&pointsNew);
1534:     for (q = 0; q < size; ++q) {
1535:       const PetscInt point = points[q];

1538:       pointsNew[q] = perm[point];
1539:     }
1540:     ISRestoreIndices((*labelNew)->points[v],&points);
1541:     PetscSortInt(size, pointsNew);
1542:     ISDestroy(&((*labelNew)->points[v]));
1543:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1544:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1545:       PetscFree(pointsNew);
1546:     } else {
1547:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1548:     }
1549:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1550:   }
1551:   ISRestoreIndices(permutation, &perm);
1552:   if (label->bt) {
1553:     PetscBTDestroy(&label->bt);
1554:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1555:   }
1556:   return 0;
1557: }

1559: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1560: {
1561:   MPI_Comm       comm;
1562:   PetscInt       s, l, nroots, nleaves, offset, size;
1563:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1564:   PetscSection   rootSection;
1565:   PetscSF        labelSF;

1567:   if (label) DMLabelMakeAllValid_Private(label);
1568:   PetscObjectGetComm((PetscObject)sf, &comm);
1569:   /* Build a section of stratum values per point, generate the according SF
1570:      and distribute point-wise stratum values to leaves. */
1571:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1572:   PetscSectionCreate(comm, &rootSection);
1573:   PetscSectionSetChart(rootSection, 0, nroots);
1574:   if (label) {
1575:     for (s = 0; s < label->numStrata; ++s) {
1576:       const PetscInt *points;

1578:       ISGetIndices(label->points[s], &points);
1579:       for (l = 0; l < label->stratumSizes[s]; l++) {
1580:         PetscSectionAddDof(rootSection, points[l], 1);
1581:       }
1582:       ISRestoreIndices(label->points[s], &points);
1583:     }
1584:   }
1585:   PetscSectionSetUp(rootSection);
1586:   /* Create a point-wise array of stratum values */
1587:   PetscSectionGetStorageSize(rootSection, &size);
1588:   PetscMalloc1(size, &rootStrata);
1589:   PetscCalloc1(nroots, &rootIdx);
1590:   if (label) {
1591:     for (s = 0; s < label->numStrata; ++s) {
1592:       const PetscInt *points;

1594:       ISGetIndices(label->points[s], &points);
1595:       for (l = 0; l < label->stratumSizes[s]; l++) {
1596:         const PetscInt p = points[l];
1597:         PetscSectionGetOffset(rootSection, p, &offset);
1598:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1599:       }
1600:       ISRestoreIndices(label->points[s], &points);
1601:     }
1602:   }
1603:   /* Build SF that maps label points to remote processes */
1604:   PetscSectionCreate(comm, leafSection);
1605:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1606:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1607:   PetscFree(remoteOffsets);
1608:   /* Send the strata for each point over the derived SF */
1609:   PetscSectionGetStorageSize(*leafSection, &size);
1610:   PetscMalloc1(size, leafStrata);
1611:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1612:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata,MPI_REPLACE);
1613:   /* Clean up */
1614:   PetscFree(rootStrata);
1615:   PetscFree(rootIdx);
1616:   PetscSectionDestroy(&rootSection);
1617:   PetscSFDestroy(&labelSF);
1618:   return 0;
1619: }

1621: /*@
1622:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1624:   Collective on sf

1626:   Input Parameters:
1627: + label - the DMLabel
1628: - sf    - the map from old to new distribution

1630:   Output Parameter:
1631: . labelnew - the new redistributed label

1633:   Level: intermediate

1635: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1636: @*/
1637: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1638: {
1639:   MPI_Comm       comm;
1640:   PetscSection   leafSection;
1641:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1642:   PetscInt      *leafStrata, *strataIdx;
1643:   PetscInt     **points;
1644:   const char    *lname = NULL;
1645:   char          *name;
1646:   PetscInt       nameSize;
1647:   PetscHSetI     stratumHash;
1648:   size_t         len = 0;
1649:   PetscMPIInt    rank;

1652:   if (label) {
1654:     DMLabelMakeAllValid_Private(label);
1655:   }
1656:   PetscObjectGetComm((PetscObject)sf, &comm);
1657:   MPI_Comm_rank(comm, &rank);
1658:   /* Bcast name */
1659:   if (rank == 0) {
1660:     PetscObjectGetName((PetscObject) label, &lname);
1661:     PetscStrlen(lname, &len);
1662:   }
1663:   nameSize = len;
1664:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1665:   PetscMalloc1(nameSize+1, &name);
1666:   if (rank == 0) PetscArraycpy(name, lname, nameSize+1);
1667:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1668:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1669:   PetscFree(name);
1670:   /* Bcast defaultValue */
1671:   if (rank == 0) (*labelNew)->defaultValue = label->defaultValue;
1672:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1673:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1674:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1675:   /* Determine received stratum values and initialise new label*/
1676:   PetscHSetICreate(&stratumHash);
1677:   PetscSectionGetStorageSize(leafSection, &size);
1678:   for (p = 0; p < size; ++p) PetscHSetIAdd(stratumHash, leafStrata[p]);
1679:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1680:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1681:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1682:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1683:   /* Turn leafStrata into indices rather than stratum values */
1684:   offset = 0;
1685:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1686:   PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1687:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1688:     PetscHMapISet((*labelNew)->hmap, (*labelNew)->stratumValues[s], s);
1689:   }
1690:   for (p = 0; p < size; ++p) {
1691:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1692:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1693:     }
1694:   }
1695:   /* Rebuild the point strata on the receiver */
1696:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1697:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1698:   for (p=pStart; p<pEnd; p++) {
1699:     PetscSectionGetDof(leafSection, p, &dof);
1700:     PetscSectionGetOffset(leafSection, p, &offset);
1701:     for (s=0; s<dof; s++) {
1702:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1703:     }
1704:   }
1705:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1706:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1707:   PetscMalloc1((*labelNew)->numStrata,&points);
1708:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1709:     PetscHSetICreate(&(*labelNew)->ht[s]);
1710:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1711:   }
1712:   /* Insert points into new strata */
1713:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1714:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1715:   for (p=pStart; p<pEnd; p++) {
1716:     PetscSectionGetDof(leafSection, p, &dof);
1717:     PetscSectionGetOffset(leafSection, p, &offset);
1718:     for (s=0; s<dof; s++) {
1719:       stratum = leafStrata[offset+s];
1720:       points[stratum][strataIdx[stratum]++] = p;
1721:     }
1722:   }
1723:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1724:     ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1725:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1726:   }
1727:   PetscFree(points);
1728:   PetscHSetIDestroy(&stratumHash);
1729:   PetscFree(leafStrata);
1730:   PetscFree(strataIdx);
1731:   PetscSectionDestroy(&leafSection);
1732:   return 0;
1733: }

1735: /*@
1736:   DMLabelGather - Gather all label values from leafs into roots

1738:   Collective on sf

1740:   Input Parameters:
1741: + label - the DMLabel
1742: - sf - the Star Forest point communication map

1744:   Output Parameters:
1745: . labelNew - the new DMLabel with localised leaf values

1747:   Level: developer

1749:   Note: This is the inverse operation to DMLabelDistribute.

1751: .seealso: DMLabelDistribute()
1752: @*/
1753: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1754: {
1755:   MPI_Comm       comm;
1756:   PetscSection   rootSection;
1757:   PetscSF        sfLabel;
1758:   PetscSFNode   *rootPoints, *leafPoints;
1759:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1760:   const PetscInt *rootDegree, *ilocal;
1761:   PetscInt       *rootStrata;
1762:   const char    *lname;
1763:   char          *name;
1764:   PetscInt       nameSize;
1765:   size_t         len = 0;
1766:   PetscMPIInt    rank, size;

1770:   PetscObjectGetComm((PetscObject)sf, &comm);
1771:   MPI_Comm_rank(comm, &rank);
1772:   MPI_Comm_size(comm, &size);
1773:   /* Bcast name */
1774:   if (rank == 0) {
1775:     PetscObjectGetName((PetscObject) label, &lname);
1776:     PetscStrlen(lname, &len);
1777:   }
1778:   nameSize = len;
1779:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1780:   PetscMalloc1(nameSize+1, &name);
1781:   if (rank == 0) PetscArraycpy(name, lname, nameSize+1);
1782:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1783:   DMLabelCreate(PETSC_COMM_SELF, name, labelNew);
1784:   PetscFree(name);
1785:   /* Gather rank/index pairs of leaves into local roots to build
1786:      an inverse, multi-rooted SF. Note that this ignores local leaf
1787:      indexing due to the use of the multiSF in PetscSFGather. */
1788:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1789:   PetscMalloc1(nroots, &leafPoints);
1790:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1791:   for (p = 0; p < nleaves; p++) {
1792:     PetscInt ilp = ilocal ? ilocal[p] : p;

1794:     leafPoints[ilp].index = ilp;
1795:     leafPoints[ilp].rank  = rank;
1796:   }
1797:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1798:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1799:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1800:   PetscMalloc1(nmultiroots, &rootPoints);
1801:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1802:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1803:   PetscSFCreate(comm,& sfLabel);
1804:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1805:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1806:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1807:   /* Rebuild the point strata on the receiver */
1808:   for (p = 0, idx = 0; p < nroots; p++) {
1809:     for (d = 0; d < rootDegree[p]; d++) {
1810:       PetscSectionGetDof(rootSection, idx+d, &dof);
1811:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1812:       for (s = 0; s < dof; s++) DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);
1813:     }
1814:     idx += rootDegree[p];
1815:   }
1816:   PetscFree(leafPoints);
1817:   PetscFree(rootStrata);
1818:   PetscSectionDestroy(&rootSection);
1819:   PetscSFDestroy(&sfLabel);
1820:   return 0;
1821: }

1823: /*@
1824:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1826:   Not collective

1828:   Input Parameter:
1829: . label - the DMLabel

1831:   Output Parameters:
1832: + section - the section giving offsets for each stratum
1833: - is - An IS containing all the label points

1835:   Level: developer

1837: .seealso: DMLabelDistribute()
1838: @*/
1839: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1840: {
1841:   IS              vIS;
1842:   const PetscInt *values;
1843:   PetscInt       *points;
1844:   PetscInt        nV, vS = 0, vE = 0, v, N;

1847:   DMLabelGetNumValues(label, &nV);
1848:   DMLabelGetValueIS(label, &vIS);
1849:   ISGetIndices(vIS, &values);
1850:   if (nV) {vS = values[0]; vE = values[0]+1;}
1851:   for (v = 1; v < nV; ++v) {
1852:     vS = PetscMin(vS, values[v]);
1853:     vE = PetscMax(vE, values[v]+1);
1854:   }
1855:   PetscSectionCreate(PETSC_COMM_SELF, section);
1856:   PetscSectionSetChart(*section, vS, vE);
1857:   for (v = 0; v < nV; ++v) {
1858:     PetscInt n;

1860:     DMLabelGetStratumSize(label, values[v], &n);
1861:     PetscSectionSetDof(*section, values[v], n);
1862:   }
1863:   PetscSectionSetUp(*section);
1864:   PetscSectionGetStorageSize(*section, &N);
1865:   PetscMalloc1(N, &points);
1866:   for (v = 0; v < nV; ++v) {
1867:     IS              is;
1868:     const PetscInt *spoints;
1869:     PetscInt        dof, off, p;

1871:     PetscSectionGetDof(*section, values[v], &dof);
1872:     PetscSectionGetOffset(*section, values[v], &off);
1873:     DMLabelGetStratumIS(label, values[v], &is);
1874:     ISGetIndices(is, &spoints);
1875:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1876:     ISRestoreIndices(is, &spoints);
1877:     ISDestroy(&is);
1878:   }
1879:   ISRestoreIndices(vIS, &values);
1880:   ISDestroy(&vIS);
1881:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1882:   return 0;
1883: }

1885: /*@
1886:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1887:   the local section and an SF describing the section point overlap.

1889:   Collective on sf

1891:   Input Parameters:
1892: + s - The PetscSection for the local field layout
1893: . sf - The SF describing parallel layout of the section points
1894: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1895: . label - The label specifying the points
1896: - labelValue - The label stratum specifying the points

1898:   Output Parameter:
1899: . gsection - The PetscSection for the global field layout

1901:   Note: This gives negative sizes and offsets to points not owned by this process

1903:   Level: developer

1905: .seealso: PetscSectionCreate()
1906: @*/
1907: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1908: {
1909:   PetscInt      *neg = NULL, *tmpOff = NULL;
1910:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1915:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1916:   PetscSectionGetChart(s, &pStart, &pEnd);
1917:   PetscSectionSetChart(*gsection, pStart, pEnd);
1918:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1919:   if (nroots >= 0) {
1921:     PetscCalloc1(nroots, &neg);
1922:     if (nroots > pEnd-pStart) {
1923:       PetscCalloc1(nroots, &tmpOff);
1924:     } else {
1925:       tmpOff = &(*gsection)->atlasDof[-pStart];
1926:     }
1927:   }
1928:   /* Mark ghost points with negative dof */
1929:   for (p = pStart; p < pEnd; ++p) {
1930:     PetscInt value;

1932:     DMLabelGetValue(label, p, &value);
1933:     if (value != labelValue) continue;
1934:     PetscSectionGetDof(s, p, &dof);
1935:     PetscSectionSetDof(*gsection, p, dof);
1936:     PetscSectionGetConstraintDof(s, p, &cdof);
1937:     if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
1938:     if (neg) neg[p] = -(dof+1);
1939:   }
1940:   PetscSectionSetUpBC(*gsection);
1941:   if (nroots >= 0) {
1942:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1943:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1944:     if (nroots > pEnd-pStart) {
1945:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1946:     }
1947:   }
1948:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1949:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1950:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1951:     (*gsection)->atlasOff[p] = off;
1952:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1953:   }
1954:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1955:   globalOff -= off;
1956:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1957:     (*gsection)->atlasOff[p] += globalOff;
1958:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1959:   }
1960:   /* Put in negative offsets for ghost points */
1961:   if (nroots >= 0) {
1962:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1963:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1964:     if (nroots > pEnd-pStart) {
1965:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1966:     }
1967:   }
1968:   if (nroots >= 0 && nroots > pEnd-pStart) PetscFree(tmpOff);
1969:   PetscFree(neg);
1970:   return 0;
1971: }

1973: typedef struct _n_PetscSectionSym_Label
1974: {
1975:   DMLabel           label;
1976:   PetscCopyMode     *modes;
1977:   PetscInt          *sizes;
1978:   const PetscInt    ***perms;
1979:   const PetscScalar ***rots;
1980:   PetscInt          (*minMaxOrients)[2];
1981:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1982: } PetscSectionSym_Label;

1984: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1985: {
1986:   PetscInt              i, j;
1987:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;

1989:   for (i = 0; i <= sl->numStrata; i++) {
1990:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1991:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1992:         if (sl->perms[i]) PetscFree(sl->perms[i][j]);
1993:         if (sl->rots[i]) PetscFree(sl->rots[i][j]);
1994:       }
1995:       if (sl->perms[i]) {
1996:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1998:         PetscFree(perms);
1999:       }
2000:       if (sl->rots[i]) {
2001:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

2003:         PetscFree(rots);
2004:       }
2005:     }
2006:   }
2007:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
2008:   DMLabelDestroy(&sl->label);
2009:   sl->numStrata = 0;
2010:   return 0;
2011: }

2013: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
2014: {
2015:   PetscSectionSymLabelReset(sym);
2016:   PetscFree(sym->data);
2017:   return 0;
2018: }

2020: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
2021: {
2022:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
2023:   PetscBool             isAscii;
2024:   DMLabel               label = sl->label;
2025:   const char           *name;

2027:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
2028:   if (isAscii) {
2029:     PetscInt          i, j, k;
2030:     PetscViewerFormat format;

2032:     PetscViewerGetFormat(viewer,&format);
2033:     if (label) {
2034:       PetscViewerGetFormat(viewer,&format);
2035:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2036:         PetscViewerASCIIPushTab(viewer);
2037:         DMLabelView(label, viewer);
2038:         PetscViewerASCIIPopTab(viewer);
2039:       } else {
2040:         PetscObjectGetName((PetscObject) sl->label, &name);
2041:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",name);
2042:       }
2043:     } else {
2044:       PetscViewerASCIIPrintf(viewer, "No label given\n");
2045:     }
2046:     PetscViewerASCIIPushTab(viewer);
2047:     for (i = 0; i <= sl->numStrata; i++) {
2048:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

2050:       if (!(sl->perms[i] || sl->rots[i])) {
2051:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
2052:       } else {
2053:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
2054:         PetscViewerASCIIPushTab(viewer);
2055:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
2056:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2057:           PetscViewerASCIIPushTab(viewer);
2058:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
2059:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
2060:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
2061:             } else {
2062:               PetscInt tab;

2064:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
2065:               PetscViewerASCIIPushTab(viewer);
2066:               PetscViewerASCIIGetTab(viewer,&tab);
2067:               if (sl->perms[i] && sl->perms[i][j]) {
2068:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
2069:                 PetscViewerASCIISetTab(viewer,0);
2070:                 for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);
2071:                 PetscViewerASCIIPrintf(viewer,"\n");
2072:                 PetscViewerASCIISetTab(viewer,tab);
2073:               }
2074:               if (sl->rots[i] && sl->rots[i][j]) {
2075:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
2076:                 PetscViewerASCIISetTab(viewer,0);
2077: #if defined(PETSC_USE_COMPLEX)
2078:                 for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer," %+f+i*%+f",PetscRealPart(sl->rots[i][j][k]),PetscImaginaryPart(sl->rots[i][j][k]));
2079: #else
2080:                 for (k = 0; k < sl->sizes[i]; k++) PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);
2081: #endif
2082:                 PetscViewerASCIIPrintf(viewer,"\n");
2083:                 PetscViewerASCIISetTab(viewer,tab);
2084:               }
2085:               PetscViewerASCIIPopTab(viewer);
2086:             }
2087:           }
2088:           PetscViewerASCIIPopTab(viewer);
2089:         }
2090:         PetscViewerASCIIPopTab(viewer);
2091:       }
2092:     }
2093:     PetscViewerASCIIPopTab(viewer);
2094:   }
2095:   return 0;
2096: }

2098: /*@
2099:   PetscSectionSymLabelSetLabel - set the label whose strata will define the points that receive symmetries

2101:   Logically collective on sym

2103:   Input parameters:
2104: + sym - the section symmetries
2105: - label - the DMLabel describing the types of points

2107:   Level: developer:

2109: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
2110: @*/
2111: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
2112: {
2113:   PetscSectionSym_Label *sl;

2116:   sl = (PetscSectionSym_Label *) sym->data;
2117:   if (sl->label && sl->label != label) PetscSectionSymLabelReset(sym);
2118:   if (label) {
2119:     sl->label = label;
2120:     PetscObjectReference((PetscObject) label);
2121:     DMLabelGetNumValues(label,&sl->numStrata);
2122:     PetscMalloc5(sl->numStrata+1,&sl->modes,sl->numStrata+1,&sl->sizes,sl->numStrata+1,&sl->perms,sl->numStrata+1,&sl->rots,sl->numStrata+1,&sl->minMaxOrients);
2123:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
2124:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
2125:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
2126:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
2127:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
2128:   }
2129:   return 0;
2130: }

2132: /*@C
2133:   PetscSectionSymLabelGetStratum - get the symmetries for the orientations of a stratum

2135:   Logically collective on sym

2137:   Input Parameters:
2138: + sym       - the section symmetries
2139: - stratum   - the stratum value in the label that we are assigning symmetries for

2141:   Output Parameters:
2142: + size      - the number of dofs for points in the stratum of the label
2143: . minOrient - the smallest orientation for a point in this stratum
2144: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2145: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
2146: - rots      - NULL if there are no rotations, or (maxOrient - minOrient) sets of rotations, one for each orientation.  A NULL set of orientations is the identity

2148:   Level: developer

2150: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2151: @*/
2152: PetscErrorCode PetscSectionSymLabelGetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt *size, PetscInt *minOrient, PetscInt *maxOrient, const PetscInt ***perms, const PetscScalar ***rots)
2153: {
2154:   PetscSectionSym_Label *sl;
2155:   const char            *name;
2156:   PetscInt               i;

2159:   sl = (PetscSectionSym_Label *) sym->data;
2161:   for (i = 0; i <= sl->numStrata; i++) {
2162:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

2164:     if (stratum == value) break;
2165:   }
2166:   PetscObjectGetName((PetscObject) sl->label, &name);
2173:   return 0;
2174: }

2176: /*@C
2177:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

2179:   Logically collective on sym

2181:   InputParameters:
2182: + sym       - the section symmetries
2183: . stratum   - the stratum value in the label that we are assigning symmetries for
2184: . size      - the number of dofs for points in the stratum of the label
2185: . minOrient - the smallest orientation for a point in this stratum
2186: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
2187: . mode      - how sym should copy the perms and rots arrays
2188: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
2189: - rots      - NULL if there are no rotations, or (maxOrient - minOrient) sets of rotations, one for each orientation.  A NULL set of orientations is the identity

2191:   Level: developer

2193: .seealso: PetscSectionSymLabelGetStratum(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
2194: @*/
2195: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
2196: {
2197:   PetscSectionSym_Label *sl;
2198:   const char            *name;
2199:   PetscInt               i, j, k;

2202:   sl = (PetscSectionSym_Label *) sym->data;
2204:   for (i = 0; i <= sl->numStrata; i++) {
2205:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

2207:     if (stratum == value) break;
2208:   }
2209:   PetscObjectGetName((PetscObject) sl->label, &name);
2211:   sl->sizes[i] = size;
2212:   sl->modes[i] = mode;
2213:   sl->minMaxOrients[i][0] = minOrient;
2214:   sl->minMaxOrients[i][1] = maxOrient;
2215:   if (mode == PETSC_COPY_VALUES) {
2216:     if (perms) {
2217:       PetscInt    **ownPerms;

2219:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
2220:       for (j = 0; j < maxOrient-minOrient; j++) {
2221:         if (perms[j]) {
2222:           PetscMalloc1(size,&ownPerms[j]);
2223:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
2224:         }
2225:       }
2226:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
2227:     }
2228:     if (rots) {
2229:       PetscScalar **ownRots;

2231:       PetscCalloc1(maxOrient - minOrient,&ownRots);
2232:       for (j = 0; j < maxOrient-minOrient; j++) {
2233:         if (rots[j]) {
2234:           PetscMalloc1(size,&ownRots[j]);
2235:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
2236:         }
2237:       }
2238:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
2239:     }
2240:   } else {
2241:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
2242:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
2243:   }
2244:   return 0;
2245: }

2247: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
2248: {
2249:   PetscInt              i, j, numStrata;
2250:   PetscSectionSym_Label *sl;
2251:   DMLabel               label;

2253:   sl = (PetscSectionSym_Label *) sym->data;
2254:   numStrata = sl->numStrata;
2255:   label     = sl->label;
2256:   for (i = 0; i < numPoints; i++) {
2257:     PetscInt point = points[2*i];
2258:     PetscInt ornt  = points[2*i+1];

2260:     for (j = 0; j < numStrata; j++) {
2261:       if (label->validIS[j]) {
2262:         PetscInt k;

2264:         ISLocate(label->points[j],point,&k);
2265:         if (k >= 0) break;
2266:       } else {
2267:         PetscBool has;

2269:         PetscHSetIHas(label->ht[j], point, &has);
2270:         if (has) break;
2271:       }
2272:     }
2274:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
2275:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
2276:   }
2277:   return 0;
2278: }

2280: static PetscErrorCode PetscSectionSymCopy_Label(PetscSectionSym sym, PetscSectionSym nsym)
2281: {
2282:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) nsym->data;
2283:   IS                     valIS;
2284:   const PetscInt        *values;
2285:   PetscInt               Nv, v;

2287:   DMLabelGetNumValues(sl->label, &Nv);
2288:   DMLabelGetValueIS(sl->label, &valIS);
2289:   ISGetIndices(valIS, &values);
2290:   for (v = 0; v < Nv; ++v) {
2291:     const PetscInt      val = values[v];
2292:     PetscInt            size, minOrient, maxOrient;
2293:     const PetscInt    **perms;
2294:     const PetscScalar **rots;

2296:     PetscSectionSymLabelGetStratum(sym,  val, &size, &minOrient, &maxOrient, &perms, &rots);
2297:     PetscSectionSymLabelSetStratum(nsym, val,  size,  minOrient,  maxOrient, PETSC_COPY_VALUES, perms, rots);
2298:   }
2299:   ISDestroy(&valIS);
2300:   return 0;
2301: }

2303: static PetscErrorCode PetscSectionSymDistribute_Label(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
2304: {
2305:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
2306:   DMLabel                dlabel;

2308:   DMLabelDistribute(sl->label, migrationSF, &dlabel);
2309:   PetscSectionSymCreateLabel(PetscObjectComm((PetscObject) sym), dlabel, dsym);
2310:   DMLabelDestroy(&dlabel);
2311:   PetscSectionSymCopy(sym, *dsym);
2312:   return 0;
2313: }

2315: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
2316: {
2317:   PetscSectionSym_Label *sl;

2319:   PetscNewLog(sym,&sl);
2320:   sym->ops->getpoints  = PetscSectionSymGetPoints_Label;
2321:   sym->ops->distribute = PetscSectionSymDistribute_Label;
2322:   sym->ops->copy       = PetscSectionSymCopy_Label;
2323:   sym->ops->view       = PetscSectionSymView_Label;
2324:   sym->ops->destroy    = PetscSectionSymDestroy_Label;
2325:   sym->data            = (void *) sl;
2326:   return 0;
2327: }

2329: /*@
2330:   PetscSectionSymCreateLabel - Create a section symmetry that assigns one symmetry to each stratum of a label

2332:   Collective

2334:   Input Parameters:
2335: + comm - the MPI communicator for the new symmetry
2336: - label - the label defining the strata

2338:   Output Parameters:
2339: . sym - the section symmetries

2341:   Level: developer

2343: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
2344: @*/
2345: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
2346: {
2347:   DMInitializePackage();
2348:   PetscSectionSymCreate(comm,sym);
2349:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
2350:   PetscSectionSymLabelSetLabel(*sym,label);
2351:   return 0;
2352: }