Actual source code: dmlabel.c

petsc-3.10.4 2019-02-26
Report Typos and Errors
  1:  #include <petscdm.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petsc/private/isimpl.h>
  4:  #include <petscsf.h>

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

  9:   Input parameter:
 10: . name - The label name

 12:   Output parameter:
 13: . label - The DMLabel

 15:   Level: beginner

 17: .seealso: DMLabelDestroy()
 18: @*/
 19: PetscErrorCode DMLabelCreate(const char name[], DMLabel *label)
 20: {

 24:   PetscNew(label);
 25:   PetscStrallocpy(name, &(*label)->name);

 27:   (*label)->refct          = 1;
 28:   (*label)->state          = -1;
 29:   (*label)->numStrata      = 0;
 30:   (*label)->defaultValue   = -1;
 31:   (*label)->stratumValues  = NULL;
 32:   (*label)->validIS        = NULL;
 33:   (*label)->stratumSizes   = NULL;
 34:   (*label)->points         = NULL;
 35:   (*label)->ht             = NULL;
 36:   (*label)->pStart         = -1;
 37:   (*label)->pEnd           = -1;
 38:   (*label)->bt             = NULL;
 39:   return(0);
 40: }

 42: /*
 43:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 45:   Input parameter:
 46: + label - The DMLabel
 47: - v - The stratum value

 49:   Output parameter:
 50: . label - The DMLabel with stratum in sorted list format

 52:   Level: developer

 54: .seealso: DMLabelCreate()
 55: */
 56: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 57: {
 58:   PetscInt       off = 0;
 59:   PetscInt       *pointArray;

 62:   if (label->validIS[v]) return 0;
 64:   if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
 65:   PetscHSetIGetSize(label->ht[v], &label->stratumSizes[v]);

 67:   PetscMalloc1(label->stratumSizes[v], &pointArray);
 68:   PetscHSetIGetElems(label->ht[v], &off, pointArray);
 69:   if (off != label->stratumSizes[v]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of contributed points %D from value %D should be %D", off, label->stratumValues[v], label->stratumSizes[v]);
 70:   PetscHSetIClear(label->ht[v]);
 71:   PetscSortInt(label->stratumSizes[v], pointArray);
 72:   if (label->bt) {
 73:     PetscInt p;

 75:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 76:       const PetscInt point = pointArray[p];

 78:       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);
 79:       PetscBTSet(label->bt, point - label->pStart);
 80:     }
 81:   }
 82:   ISCreateGeneral(PETSC_COMM_SELF,label->stratumSizes[v],pointArray,PETSC_OWN_POINTER,&(label->points[v]));
 83:   PetscObjectSetName((PetscObject) (label->points[v]), "indices");
 84:   label->validIS[v] = PETSC_TRUE;
 85:   ++label->state;
 86:   return(0);
 87: }

 89: /*
 90:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

 92:   Input parameter:
 93: . label - The DMLabel

 95:   Output parameter:
 96: . label - The DMLabel with all strata in sorted list format

 98:   Level: developer

100: .seealso: DMLabelCreate()
101: */
102: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
103: {
104:   PetscInt       v;

108:   for (v = 0; v < label->numStrata; v++){
109:     DMLabelMakeValid_Private(label, v);
110:   }
111:   return(0);
112: }

114: /*
115:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

117:   Input parameter:
118: + label - The DMLabel
119: - v - The stratum value

121:   Output parameter:
122: . label - The DMLabel with stratum in hash format

124:   Level: developer

126: .seealso: DMLabelCreate()
127: */
128: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
129: {
130:   PetscInt                    p;
131:   const PetscInt              *points;

134:   if (!label->validIS[v]) return 0;
136:   if (v < 0 || v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeInvalid_Private\n", v);
137:   if (label->points[v]) {
138:     ISGetIndices(label->points[v],&points);
139:     for (p = 0; p < label->stratumSizes[v]; ++p) {
140:       PetscHSetIAdd(label->ht[v], points[p]);
141:     }
142:     ISRestoreIndices(label->points[v],&points);
143:     ISDestroy(&(label->points[v]));
144:   }
145:   label->validIS[v] = PETSC_FALSE;
146:   return(0);
147: }

149: PetscErrorCode DMLabelGetState(DMLabel label, PetscObjectState *state)
150: {
153:   *state = label->state;
154:   return(0);
155: }

157: PETSC_STATIC_INLINE PetscErrorCode DMLabelLookupStratum(DMLabel label, PetscInt value, PetscInt *index)
158: {
159:   PetscInt v;
161:   for (*index = -1, v = 0; v < label->numStrata; ++v)
162:     if (label->stratumValues[v] == value) { *index = v; break; }
163:   return(0);
164: }

166: static PetscErrorCode DMLabelNewStratum(DMLabel label, PetscInt value, PetscInt *index)
167: {
168:   PetscInt   v, *tmpV, *tmpS;
169:   IS         *tmpP;
170:   PetscHSetI *tmpH;
171:   PetscBool  *tmpB;

175:   PetscMalloc1((label->numStrata+1), &tmpV);
176:   PetscMalloc1((label->numStrata+1), &tmpS);
177:   PetscMalloc1((label->numStrata+1), &tmpH);
178:   PetscMalloc1((label->numStrata+1), &tmpP);
179:   PetscMalloc1((label->numStrata+1), &tmpB);
180:   for (v = 0; v < label->numStrata; ++v) {
181:     tmpV[v] = label->stratumValues[v];
182:     tmpS[v] = label->stratumSizes[v];
183:     tmpH[v] = label->ht[v];
184:     tmpP[v] = label->points[v];
185:     tmpB[v] = label->validIS[v];
186:   }
187:   tmpV[v] = value;
188:   tmpS[v] = 0;
189:   PetscHSetICreate(&tmpH[v]);
190:   ISCreateStride(PETSC_COMM_SELF,0,0,1,&tmpP[v]);
191:   tmpB[v] = PETSC_TRUE;

193:   ++label->numStrata;
194:   PetscFree(label->stratumValues);
195:   PetscFree(label->stratumSizes);
196:   PetscFree(label->ht);
197:   PetscFree(label->points);
198:   PetscFree(label->validIS);
199:   label->stratumValues = tmpV;
200:   label->stratumSizes  = tmpS;
201:   label->ht            = tmpH;
202:   label->points        = tmpP;
203:   label->validIS       = tmpB;

205:   *index = v;
206:   return(0);
207: }

209: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
210: {
211:   PetscInt       v;

215:   DMLabelLookupStratum(label, value, &v);
216:   if (v < 0) {DMLabelNewStratum(label, value, &v);}
217:   return(0);
218: }

220: /*@C
221:   DMLabelGetName - Return the name of a DMLabel object

223:   Input parameter:
224: . label - The DMLabel

226:   Output parameter:
227: . name - The label name

229:   Level: beginner

231: .seealso: DMLabelCreate()
232: @*/
233: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
234: {
237:   *name = label->name;
238:   return(0);
239: }

241: /*@C
242:   DMLabelSetName - Sets the name of a DMLabel object

244:   Input parameters:
245: + label - The DMLabel
246: - name  - The label name

248:   Level: beginner

250: .seealso: DMLabelCreate()
251: @*/
252: PetscErrorCode DMLabelSetName(DMLabel label, const char *name)
253: {

258:   PetscFree(label->name);
259:   PetscStrallocpy(name, &label->name);
260:   return(0);
261: }

263: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
264: {
265:   PetscInt       v;
266:   PetscMPIInt    rank;

270:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
271:   PetscViewerASCIIPushSynchronized(viewer);
272:   if (label) {
273:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
274:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
275:     for (v = 0; v < label->numStrata; ++v) {
276:       const PetscInt value = label->stratumValues[v];
277:       const PetscInt *points;
278:       PetscInt       p;

280:       ISGetIndices(label->points[v],&points);
281:       for (p = 0; p < label->stratumSizes[v]; ++p) {
282:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, points[p], value);
283:       }
284:       ISRestoreIndices(label->points[v],&points);
285:     }
286:   }
287:   PetscViewerFlush(viewer);
288:   PetscViewerASCIIPopSynchronized(viewer);
289:   return(0);
290: }

292: /*@C
293:   DMLabelView - View the label

295:   Input Parameters:
296: + label - The DMLabel
297: - viewer - The PetscViewer

299:   Level: intermediate

301: .seealso: DMLabelCreate(), DMLabelDestroy()
302: @*/
303: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
304: {
305:   PetscBool      iascii;

310:   if (label) {DMLabelMakeAllValid_Private(label);}
311:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
312:   if (iascii) {
313:     DMLabelView_Ascii(label, viewer);
314:   }
315:   return(0);
316: }

318: /*@
319:   DMLabelDestroy - Destroys a DMLabel

321:   Input Parameter:
322: . label - The DMLabel

324:   Level: beginner

326: .seealso: DMLabelCreate()
327: @*/
328: PetscErrorCode DMLabelDestroy(DMLabel *label)
329: {
330:   PetscInt       v;

334:   if (!(*label)) return(0);
335:   if (--(*label)->refct > 0) return(0);
336:   PetscFree((*label)->name);
337:   PetscFree((*label)->stratumValues);
338:   PetscFree((*label)->stratumSizes);
339:   for (v = 0; v < (*label)->numStrata; ++v) {
340:     PetscHSetIDestroy(&(*label)->ht[v]);
341:     ISDestroy(&(*label)->points[v]);
342:   }
343:   PetscFree((*label)->ht);
344:   PetscFree((*label)->points);
345:   PetscFree((*label)->validIS);
346:   PetscBTDestroy(&(*label)->bt);
347:   PetscFree(*label);
348:   return(0);
349: }

351: /*@
352:   DMLabelDuplicate - Duplicates a DMLabel

354:   Input Parameter:
355: . label - The DMLabel

357:   Output Parameter:
358: . labelnew - location to put new vector

360:   Level: intermediate

362: .seealso: DMLabelCreate(), DMLabelDestroy()
363: @*/
364: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
365: {
366:   PetscInt       v;

370:   DMLabelMakeAllValid_Private(label);
371:   PetscNew(labelnew);
372:   PetscStrallocpy(label->name, &(*labelnew)->name);

374:   (*labelnew)->refct        = 1;
375:   (*labelnew)->numStrata    = label->numStrata;
376:   (*labelnew)->defaultValue = label->defaultValue;
377:   if (label->numStrata) {
378:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
379:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
380:     PetscMalloc1(label->numStrata, &(*labelnew)->ht);
381:     PetscMalloc1(label->numStrata, &(*labelnew)->points);
382:     PetscMalloc1(label->numStrata, &(*labelnew)->validIS);
383:     /* Could eliminate unused space here */
384:     for (v = 0; v < label->numStrata; ++v) {
385:       PetscHSetICreate(&(*labelnew)->ht[v]);
386:       (*labelnew)->validIS[v]        = PETSC_TRUE;
387:       (*labelnew)->stratumValues[v]  = label->stratumValues[v];
388:       (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
389:       PetscObjectReference((PetscObject) (label->points[v]));
390:       (*labelnew)->points[v]         = label->points[v];
391:     }
392:   }
393:   (*labelnew)->pStart = -1;
394:   (*labelnew)->pEnd   = -1;
395:   (*labelnew)->bt     = NULL;
396:   return(0);
397: }

399: /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
400: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
401: {
402:   PetscInt       v;

406:   DMLabelDestroyIndex(label);
407:   DMLabelMakeAllValid_Private(label);
408:   label->pStart = pStart;
409:   label->pEnd   = pEnd;
410:   PetscBTCreate(pEnd - pStart, &label->bt);
411:   for (v = 0; v < label->numStrata; ++v) {
412:     const PetscInt *points;
413:     PetscInt       i;

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

419:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
420:       PetscBTSet(label->bt, point - pStart);
421:     }
422:     ISRestoreIndices(label->points[v], &points);
423:   }
424:   return(0);
425: }

427: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
428: {

432:   label->pStart = -1;
433:   label->pEnd   = -1;
434:   PetscBTDestroy(&label->bt);
435:   return(0);
436: }

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

441:   Input Parameters:
442: + label - the DMLabel
443: - value - the value

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

448:   Level: developer

450: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
451: @*/
452: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
453: {
454:   PetscInt v;

459:   DMLabelLookupStratum(label, value, &v);
460:   *contains = v < 0 ? PETSC_FALSE : PETSC_TRUE;
461:   return(0);
462: }

464: /*@
465:   DMLabelHasPoint - Determine whether a label assigns a value to a point

467:   Input Parameters:
468: + label - the DMLabel
469: - point - the point

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

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

476:   Level: developer

478: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
479: @*/
480: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
481: {

486:   DMLabelMakeAllValid_Private(label);
487: #if defined(PETSC_USE_DEBUG)
488:   if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
489:   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);
490: #endif
491:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
492:   return(0);
493: }

495: /*@
496:   DMLabelStratumHasPoint - Return true if the stratum contains a point

498:   Input Parameters:
499: + label - the DMLabel
500: . value - the stratum value
501: - point - the point

503:   Output Parameter:
504: . contains - true if the stratum contains the point

506:   Level: intermediate

508: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
509: @*/
510: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
511: {
512:   PetscInt       v;

517:   *contains = PETSC_FALSE;
518:   DMLabelLookupStratum(label, value, &v);
519:   if (v < 0) return(0);

521:   if (label->validIS[v]) {
522:     PetscInt i;

524:     ISLocate(label->points[v], point, &i);
525:     if (i >= 0) *contains = PETSC_TRUE;
526:   } else {
527:     PetscBool has;

529:     PetscHSetIHas(label->ht[v], point, &has);
530:     if (has) *contains = PETSC_TRUE;
531:   }
532:   return(0);
533: }

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

539:   Input parameter:
540: . label - a DMLabel object

542:   Output parameter:
543: . defaultValue - the default value

545:   Level: beginner

547: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
548: @*/
549: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
550: {
552:   *defaultValue = label->defaultValue;
553:   return(0);
554: }

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

560:   Input parameter:
561: . label - a DMLabel object

563:   Output parameter:
564: . defaultValue - the default value

566:   Level: beginner

568: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
569: @*/
570: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
571: {
573:   label->defaultValue = defaultValue;
574:   return(0);
575: }

577: /*@
578:   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())

580:   Input Parameters:
581: + label - the DMLabel
582: - point - the point

584:   Output Parameter:
585: . value - The point value, or -1

587:   Level: intermediate

589: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
590: @*/
591: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
592: {
593:   PetscInt       v;

598:   *value = label->defaultValue;
599:   for (v = 0; v < label->numStrata; ++v) {
600:     if (label->validIS[v]) {
601:       PetscInt i;

603:       ISLocate(label->points[v], point, &i);
604:       if (i >= 0) {
605:         *value = label->stratumValues[v];
606:         break;
607:       }
608:     } else {
609:       PetscBool has;

611:       PetscHSetIHas(label->ht[v], point, &has);
612:       if (has) {
613:         *value = label->stratumValues[v];
614:         break;
615:       }
616:     }
617:   }
618:   return(0);
619: }

621: /*@
622:   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 somethingg different), then this function will do nothing.

624:   Input Parameters:
625: + label - the DMLabel
626: . point - the point
627: - value - The point value

629:   Level: intermediate

631: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
632: @*/
633: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
634: {
635:   PetscInt                    v;
636:   PetscErrorCode              ierr;

639:   /* Find label value, add new entry if needed */
640:   if (value == label->defaultValue) return(0);
641:   DMLabelLookupStratum(label, value, &v);
642:   if (v < 0) {DMLabelNewStratum(label, value, &v);}
643:   /* Set key */
644:   DMLabelMakeInvalid_Private(label, v);
645:   PetscHSetIAdd(label->ht[v], point);
646:   return(0);
647: }

649: /*@
650:   DMLabelClearValue - Clear the value a label assigns to a point

652:   Input Parameters:
653: + label - the DMLabel
654: . point - the point
655: - value - The point value

657:   Level: intermediate

659: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
660: @*/
661: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
662: {
663:   PetscInt       v;

667:   /* Find label value */
668:   DMLabelLookupStratum(label, value, &v);
669:   if (v < 0) return(0);

671:   if (label->bt) {
672:     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);
673:     PetscBTClear(label->bt, point - label->pStart);
674:   }

676:   /* Delete key */
677:   DMLabelMakeInvalid_Private(label, v);
678:   PetscHSetIDel(label->ht[v], point);
679:   return(0);
680: }

682: /*@
683:   DMLabelInsertIS - Set all points in the IS to a value

685:   Input Parameters:
686: + label - the DMLabel
687: . is    - the point IS
688: - value - The point value

690:   Level: intermediate

692: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
693: @*/
694: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
695: {
696:   PetscInt        v, n, p;
697:   const PetscInt *points;
698:   PetscErrorCode  ierr;

702:   /* Find label value, add new entry if needed */
703:   if (value == label->defaultValue) return(0);
704:   DMLabelLookupStratum(label, value, &v);
705:   if (v < 0) {DMLabelNewStratum(label, value, &v);}
706:   /* Set keys */
707:   DMLabelMakeInvalid_Private(label, v);
708:   ISGetLocalSize(is, &n);
709:   ISGetIndices(is, &points);
710:   for (p = 0; p < n; ++p) {PetscHSetIAdd(label->ht[v], points[p]);}
711:   ISRestoreIndices(is, &points);
712:   return(0);
713: }

715: /*@
716:   DMLabelGetNumValues - Get the number of values that the DMLabel takes

718:   Input Parameter:
719: . label - the DMLabel

721:   Output Paramater:
722: . numValues - the number of values

724:   Level: intermediate

726: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
727: @*/
728: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
729: {
732:   *numValues = label->numStrata;
733:   return(0);
734: }

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

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

742:   Output Paramater:
743: . is    - the value IS

745:   Level: intermediate

747: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
748: @*/
749: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
750: {

755:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
756:   return(0);
757: }

759: /*@
760:   DMLabelHasStratum - Determine whether points exist with the given value

762:   Input Parameters:
763: + label - the DMLabel
764: - value - the stratum value

766:   Output Paramater:
767: . exists - Flag saying whether points exist

769:   Level: intermediate

771: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
772: @*/
773: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
774: {
775:   PetscInt v;

780:   DMLabelLookupStratum(label, value, &v);
781:   *exists = v < 0 ? PETSC_FALSE : PETSC_TRUE;
782:   return(0);
783: }

785: /*@
786:   DMLabelGetStratumSize - Get the size of a stratum

788:   Input Parameters:
789: + label - the DMLabel
790: - value - the stratum value

792:   Output Paramater:
793: . size - The number of points in the stratum

795:   Level: intermediate

797: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
798: @*/
799: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
800: {
801:   PetscInt       v;

806:   *size = 0;
807:   DMLabelLookupStratum(label, value, &v);
808:   if (v < 0) return(0);
809:   DMLabelMakeValid_Private(label, v);
810:   *size = label->stratumSizes[v];
811:   return(0);
812: }

814: /*@
815:   DMLabelGetStratumBounds - Get the largest and smallest point of a stratum

817:   Input Parameters:
818: + label - the DMLabel
819: - value - the stratum value

821:   Output Paramaters:
822: + start - the smallest point in the stratum
823: - end - the largest point in the stratum

825:   Level: intermediate

827: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
828: @*/
829: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
830: {
831:   PetscInt       v, min, max;

837:   DMLabelLookupStratum(label, value, &v);
838:   if (v < 0) return(0);
839:   DMLabelMakeValid_Private(label, v);
840:   if (label->stratumSizes[v] <= 0) return(0);
841:   ISGetMinMax(label->points[v], &min, &max);
842:   if (start) *start = min;
843:   if (end)   *end   = max+1;
844:   return(0);
845: }

847: /*@
848:   DMLabelGetStratumIS - Get an IS with the stratum points

850:   Input Parameters:
851: + label - the DMLabel
852: - value - the stratum value

854:   Output Paramater:
855: . points - The stratum points

857:   Level: intermediate

859: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
860: @*/
861: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
862: {
863:   PetscInt       v;

868:   *points = NULL;
869:   DMLabelLookupStratum(label, value, &v);
870:   if (v < 0) return(0);
871:   DMLabelMakeValid_Private(label, v);
872:   PetscObjectReference((PetscObject) label->points[v]);
873:   *points = label->points[v];
874:   return(0);
875: }

877: /*@
878:   DMLabelSetStratumIS - Set the stratum points using an IS

880:   Input Parameters:
881: + label - the DMLabel
882: . value - the stratum value
883: - points - The stratum points

885:   Level: intermediate

887: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
888: @*/
889: PetscErrorCode DMLabelSetStratumIS(DMLabel label, PetscInt value, IS is)
890: {
891:   PetscInt       v;

895:   DMLabelLookupStratum(label, value, &v);
896:   if (v < 0) {DMLabelNewStratum(label, value, &v);}
897:   if (is == label->points[v]) return(0);
898:   DMLabelClearStratum(label, value);
899:   ISGetLocalSize(is, &(label->stratumSizes[v]));
900:   PetscObjectReference((PetscObject)is);
901:   ISDestroy(&(label->points[v]));
902:   label->points[v] = is;
903:   label->validIS[v] = PETSC_TRUE;
904:   if (label->bt) {
905:     const PetscInt *points;
906:     PetscInt p;

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

912:       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);
913:       PetscBTSet(label->bt, point - label->pStart);
914:     }
915:   }
916:   return(0);
917: }

919: /*@
920:   DMLabelClearStratum - Remove a stratum

922:   Input Parameters:
923: + label - the DMLabel
924: - value - the stratum value

926:   Level: intermediate

928: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
929: @*/
930: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
931: {
932:   PetscInt       v;

936:   DMLabelLookupStratum(label, value, &v);
937:   if (v < 0) return(0);
938:   if (label->validIS[v]) {
939:     if (label->bt) {
940:       PetscInt       i;
941:       const PetscInt *points;

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

947:         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);
948:         PetscBTClear(label->bt, point - label->pStart);
949:       }
950:       ISRestoreIndices(label->points[v], &points);
951:     }
952:     label->stratumSizes[v] = 0;
953:     ISDestroy(&label->points[v]);
954:     ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_OWN_POINTER, &label->points[v]);
955:     PetscObjectSetName((PetscObject) label->points[v], "indices");
956:   } else {
957:     PetscHSetIClear(label->ht[v]);
958:   }
959:   return(0);
960: }

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

965:   Input Parameters:
966: + label - the DMLabel
967: . start - the first point
968: - end - the last point

970:   Level: intermediate

972: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
973: @*/
974: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
975: {
976:   PetscInt       v;

980:   DMLabelDestroyIndex(label);
981:   DMLabelMakeAllValid_Private(label);
982:   for (v = 0; v < label->numStrata; ++v) {
983:     PetscInt off, q;
984:     const PetscInt *points;
985:     PetscInt numPointsNew = 0, *pointsNew = NULL;

987:     ISGetIndices(label->points[v], &points);
988:     for (q = 0; q < label->stratumSizes[v]; ++q)
989:       if (points[q] >= start && points[q] < end)
990:         numPointsNew++;
991:     PetscMalloc1(numPointsNew, &pointsNew);
992:     for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
993:       if (points[q] >= start && points[q] < end)
994:         pointsNew[off++] = points[q];
995:     }
996:     ISRestoreIndices(label->points[v],&points);

998:     label->stratumSizes[v] = numPointsNew;
999:     ISDestroy(&label->points[v]);
1000:     ISCreateGeneral(PETSC_COMM_SELF,numPointsNew, pointsNew, PETSC_OWN_POINTER, &label->points[v]);
1001:     PetscObjectSetName((PetscObject) label->points[v], "indices");
1002:   }
1003:   DMLabelCreateIndex(label, start, end);
1004:   return(0);
1005: }

1007: /*@
1008:   DMLabelPermute - Create a new label with permuted points

1010:   Input Parameters:
1011: + label - the DMLabel
1012: - permutation - the point permutation

1014:   Output Parameter:
1015: . labelnew - the new label containing the permuted points

1017:   Level: intermediate

1019: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1020: @*/
1021: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
1022: {
1023:   const PetscInt *perm;
1024:   PetscInt        numValues, numPoints, v, q;
1025:   PetscErrorCode  ierr;

1028:   DMLabelMakeAllValid_Private(label);
1029:   DMLabelDuplicate(label, labelNew);
1030:   DMLabelGetNumValues(*labelNew, &numValues);
1031:   ISGetLocalSize(permutation, &numPoints);
1032:   ISGetIndices(permutation, &perm);
1033:   for (v = 0; v < numValues; ++v) {
1034:     const PetscInt size   = (*labelNew)->stratumSizes[v];
1035:     const PetscInt *points;
1036:     PetscInt *pointsNew;

1038:     ISGetIndices((*labelNew)->points[v],&points);
1039:     PetscMalloc1(size,&pointsNew);
1040:     for (q = 0; q < size; ++q) {
1041:       const PetscInt point = points[q];

1043:       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);
1044:       pointsNew[q] = perm[point];
1045:     }
1046:     ISRestoreIndices((*labelNew)->points[v],&points);
1047:     PetscSortInt(size, pointsNew);
1048:     ISDestroy(&((*labelNew)->points[v]));
1049:     if (size > 0 && pointsNew[size - 1] == pointsNew[0] + size - 1) {
1050:       ISCreateStride(PETSC_COMM_SELF,size,pointsNew[0],1,&((*labelNew)->points[v]));
1051:       PetscFree(pointsNew);
1052:     } else {
1053:       ISCreateGeneral(PETSC_COMM_SELF,size,pointsNew,PETSC_OWN_POINTER,&((*labelNew)->points[v]));
1054:     }
1055:     PetscObjectSetName((PetscObject) ((*labelNew)->points[v]), "indices");
1056:   }
1057:   ISRestoreIndices(permutation, &perm);
1058:   if (label->bt) {
1059:     PetscBTDestroy(&label->bt);
1060:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
1061:   }
1062:   return(0);
1063: }

1065: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
1066: {
1067:   MPI_Comm       comm;
1068:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
1069:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
1070:   PetscSection   rootSection;
1071:   PetscSF        labelSF;

1075:   if (label) {DMLabelMakeAllValid_Private(label);}
1076:   PetscObjectGetComm((PetscObject)sf, &comm);
1077:   /* Build a section of stratum values per point, generate the according SF
1078:      and distribute point-wise stratum values to leaves. */
1079:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
1080:   PetscSectionCreate(comm, &rootSection);
1081:   PetscSectionSetChart(rootSection, 0, nroots);
1082:   if (label) {
1083:     for (s = 0; s < label->numStrata; ++s) {
1084:       const PetscInt *points;

1086:       ISGetIndices(label->points[s], &points);
1087:       for (l = 0; l < label->stratumSizes[s]; l++) {
1088:         PetscSectionGetDof(rootSection, points[l], &dof);
1089:         PetscSectionSetDof(rootSection, points[l], dof+1);
1090:       }
1091:       ISRestoreIndices(label->points[s], &points);
1092:     }
1093:   }
1094:   PetscSectionSetUp(rootSection);
1095:   /* Create a point-wise array of stratum values */
1096:   PetscSectionGetStorageSize(rootSection, &size);
1097:   PetscMalloc1(size, &rootStrata);
1098:   PetscCalloc1(nroots, &rootIdx);
1099:   if (label) {
1100:     for (s = 0; s < label->numStrata; ++s) {
1101:       const PetscInt *points;

1103:       ISGetIndices(label->points[s], &points);
1104:       for (l = 0; l < label->stratumSizes[s]; l++) {
1105:         const PetscInt p = points[l];
1106:         PetscSectionGetOffset(rootSection, p, &offset);
1107:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
1108:       }
1109:       ISRestoreIndices(label->points[s], &points);
1110:     }
1111:   }
1112:   /* Build SF that maps label points to remote processes */
1113:   PetscSectionCreate(comm, leafSection);
1114:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
1115:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
1116:   PetscFree(remoteOffsets);
1117:   /* Send the strata for each point over the derived SF */
1118:   PetscSectionGetStorageSize(*leafSection, &size);
1119:   PetscMalloc1(size, leafStrata);
1120:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
1121:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
1122:   /* Clean up */
1123:   PetscFree(rootStrata);
1124:   PetscFree(rootIdx);
1125:   PetscSectionDestroy(&rootSection);
1126:   PetscSFDestroy(&labelSF);
1127:   return(0);
1128: }

1130: /*@
1131:   DMLabelDistribute - Create a new label pushed forward over the PetscSF

1133:   Input Parameters:
1134: + label - the DMLabel
1135: - sf    - the map from old to new distribution

1137:   Output Parameter:
1138: . labelnew - the new redistributed label

1140:   Level: intermediate

1142: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
1143: @*/
1144: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
1145: {
1146:   MPI_Comm       comm;
1147:   PetscSection   leafSection;
1148:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
1149:   PetscInt      *leafStrata, *strataIdx;
1150:   PetscInt     **points;
1151:   char          *name;
1152:   PetscInt       nameSize;
1153:   PetscHSetI     stratumHash;
1154:   size_t         len = 0;
1155:   PetscMPIInt    rank;

1159:   if (label) {DMLabelMakeAllValid_Private(label);}
1160:   PetscObjectGetComm((PetscObject)sf, &comm);
1161:   MPI_Comm_rank(comm, &rank);
1162:   /* Bcast name */
1163:   if (!rank) {PetscStrlen(label->name, &len);}
1164:   nameSize = len;
1165:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1166:   PetscMalloc1(nameSize+1, &name);
1167:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1168:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1169:   DMLabelCreate(name, labelNew);
1170:   PetscFree(name);
1171:   /* Bcast defaultValue */
1172:   if (!rank) (*labelNew)->defaultValue = label->defaultValue;
1173:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
1174:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
1175:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
1176:   /* Determine received stratum values and initialise new label*/
1177:   PetscHSetICreate(&stratumHash);
1178:   PetscSectionGetStorageSize(leafSection, &size);
1179:   for (p = 0; p < size; ++p) {PetscHSetIAdd(stratumHash, leafStrata[p]);}
1180:   PetscHSetIGetSize(stratumHash, &(*labelNew)->numStrata);
1181:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->validIS);
1182:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->validIS[s] = PETSC_TRUE;
1183:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
1184:   /* Turn leafStrata into indices rather than stratum values */
1185:   offset = 0;
1186:   PetscHSetIGetElems(stratumHash, &offset, (*labelNew)->stratumValues);
1187:   PetscSortInt((*labelNew)->numStrata,(*labelNew)->stratumValues);
1188:   for (p = 0; p < size; ++p) {
1189:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
1190:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
1191:     }
1192:   }
1193:   /* Rebuild the point strata on the receiver */
1194:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
1195:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1196:   for (p=pStart; p<pEnd; p++) {
1197:     PetscSectionGetDof(leafSection, p, &dof);
1198:     PetscSectionGetOffset(leafSection, p, &offset);
1199:     for (s=0; s<dof; s++) {
1200:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1201:     }
1202:   }
1203:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1204:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1205:   PetscMalloc1((*labelNew)->numStrata,&points);
1206:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1207:     PetscHSetICreate(&(*labelNew)->ht[s]);
1208:     PetscMalloc1((*labelNew)->stratumSizes[s], &(points[s]));
1209:   }
1210:   /* Insert points into new strata */
1211:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1212:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1213:   for (p=pStart; p<pEnd; p++) {
1214:     PetscSectionGetDof(leafSection, p, &dof);
1215:     PetscSectionGetOffset(leafSection, p, &offset);
1216:     for (s=0; s<dof; s++) {
1217:       stratum = leafStrata[offset+s];
1218:       points[stratum][strataIdx[stratum]++] = p;
1219:     }
1220:   }
1221:   for (s = 0; s < (*labelNew)->numStrata; s++) {
1222:     ISCreateGeneral(PETSC_COMM_SELF,(*labelNew)->stratumSizes[s],&(points[s][0]),PETSC_OWN_POINTER,&((*labelNew)->points[s]));
1223:     PetscObjectSetName((PetscObject)((*labelNew)->points[s]),"indices");
1224:   }
1225:   PetscFree(points);
1226:   PetscHSetIDestroy(&stratumHash);
1227:   PetscFree(leafStrata);
1228:   PetscFree(strataIdx);
1229:   PetscSectionDestroy(&leafSection);
1230:   return(0);
1231: }

1233: /*@
1234:   DMLabelGather - Gather all label values from leafs into roots

1236:   Input Parameters:
1237: + label - the DMLabel
1238: - sf - the Star Forest point communication map

1240:   Output Parameters:
1241: . labelNew - the new DMLabel with localised leaf values

1243:   Level: developer

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

1247: .seealso: DMLabelDistribute()
1248: @*/
1249: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1250: {
1251:   MPI_Comm       comm;
1252:   PetscSection   rootSection;
1253:   PetscSF        sfLabel;
1254:   PetscSFNode   *rootPoints, *leafPoints;
1255:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1256:   const PetscInt *rootDegree, *ilocal;
1257:   PetscInt       *rootStrata;
1258:   char          *name;
1259:   PetscInt       nameSize;
1260:   size_t         len = 0;
1261:   PetscMPIInt    rank, size;

1265:   PetscObjectGetComm((PetscObject)sf, &comm);
1266:   MPI_Comm_rank(comm, &rank);
1267:   MPI_Comm_size(comm, &size);
1268:   /* Bcast name */
1269:   if (!rank) {PetscStrlen(label->name, &len);}
1270:   nameSize = len;
1271:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1272:   PetscMalloc1(nameSize+1, &name);
1273:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1274:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1275:   DMLabelCreate(name, labelNew);
1276:   PetscFree(name);
1277:   /* Gather rank/index pairs of leaves into local roots to build
1278:      an inverse, multi-rooted SF. Note that this ignores local leaf
1279:      indexing due to the use of the multiSF in PetscSFGather. */
1280:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1281:   PetscMalloc1(nroots, &leafPoints);
1282:   for (p = 0; p < nroots; ++p) leafPoints[p].rank = leafPoints[p].index = -1;
1283:   for (p = 0; p < nleaves; p++) {
1284:     leafPoints[ilocal[p]].index = ilocal[p];
1285:     leafPoints[ilocal[p]].rank  = rank;
1286:   }
1287:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1288:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1289:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1290:   PetscMalloc1(nmultiroots, &rootPoints);
1291:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1292:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1293:   PetscSFCreate(comm,& sfLabel);
1294:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1295:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1296:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1297:   /* Rebuild the point strata on the receiver */
1298:   for (p = 0, idx = 0; p < nroots; p++) {
1299:     for (d = 0; d < rootDegree[p]; d++) {
1300:       PetscSectionGetDof(rootSection, idx+d, &dof);
1301:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1302:       for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1303:     }
1304:     idx += rootDegree[p];
1305:   }
1306:   PetscFree(leafPoints);
1307:   PetscFree(rootStrata);
1308:   PetscSectionDestroy(&rootSection);
1309:   PetscSFDestroy(&sfLabel);
1310:   return(0);
1311: }

1313: /*@
1314:   DMLabelConvertToSection - Make a PetscSection/IS pair that encodes the label

1316:   Input Parameter:
1317: . label - the DMLabel

1319:   Output Parameters:
1320: + section - the section giving offsets for each stratum
1321: - is - An IS containing all the label points

1323:   Level: developer

1325: .seealso: DMLabelDistribute()
1326: @*/
1327: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1328: {
1329:   IS              vIS;
1330:   const PetscInt *values;
1331:   PetscInt       *points;
1332:   PetscInt        nV, vS = 0, vE = 0, v, N;
1333:   PetscErrorCode  ierr;

1336:   DMLabelGetNumValues(label, &nV);
1337:   DMLabelGetValueIS(label, &vIS);
1338:   ISGetIndices(vIS, &values);
1339:   if (nV) {vS = values[0]; vE = values[0]+1;}
1340:   for (v = 1; v < nV; ++v) {
1341:     vS = PetscMin(vS, values[v]);
1342:     vE = PetscMax(vE, values[v]+1);
1343:   }
1344:   PetscSectionCreate(PETSC_COMM_SELF, section);
1345:   PetscSectionSetChart(*section, vS, vE);
1346:   for (v = 0; v < nV; ++v) {
1347:     PetscInt n;

1349:     DMLabelGetStratumSize(label, values[v], &n);
1350:     PetscSectionSetDof(*section, values[v], n);
1351:   }
1352:   PetscSectionSetUp(*section);
1353:   PetscSectionGetStorageSize(*section, &N);
1354:   PetscMalloc1(N, &points);
1355:   for (v = 0; v < nV; ++v) {
1356:     IS              is;
1357:     const PetscInt *spoints;
1358:     PetscInt        dof, off, p;

1360:     PetscSectionGetDof(*section, values[v], &dof);
1361:     PetscSectionGetOffset(*section, values[v], &off);
1362:     DMLabelGetStratumIS(label, values[v], &is);
1363:     ISGetIndices(is, &spoints);
1364:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1365:     ISRestoreIndices(is, &spoints);
1366:     ISDestroy(&is);
1367:   }
1368:   ISRestoreIndices(vIS, &values);
1369:   ISDestroy(&vIS);
1370:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1371:   return(0);
1372: }

1374: /*@
1375:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1376:   the local section and an SF describing the section point overlap.

1378:   Input Parameters:
1379:   + s - The PetscSection for the local field layout
1380:   . sf - The SF describing parallel layout of the section points
1381:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1382:   . label - The label specifying the points
1383:   - labelValue - The label stratum specifying the points

1385:   Output Parameter:
1386:   . gsection - The PetscSection for the global field layout

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

1390:   Level: developer

1392: .seealso: PetscSectionCreate()
1393: @*/
1394: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1395: {
1396:   PetscInt      *neg = NULL, *tmpOff = NULL;
1397:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1401:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1402:   PetscSectionGetChart(s, &pStart, &pEnd);
1403:   PetscSectionSetChart(*gsection, pStart, pEnd);
1404:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1405:   if (nroots >= 0) {
1406:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1407:     PetscCalloc1(nroots, &neg);
1408:     if (nroots > pEnd-pStart) {
1409:       PetscCalloc1(nroots, &tmpOff);
1410:     } else {
1411:       tmpOff = &(*gsection)->atlasDof[-pStart];
1412:     }
1413:   }
1414:   /* Mark ghost points with negative dof */
1415:   for (p = pStart; p < pEnd; ++p) {
1416:     PetscInt value;

1418:     DMLabelGetValue(label, p, &value);
1419:     if (value != labelValue) continue;
1420:     PetscSectionGetDof(s, p, &dof);
1421:     PetscSectionSetDof(*gsection, p, dof);
1422:     PetscSectionGetConstraintDof(s, p, &cdof);
1423:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1424:     if (neg) neg[p] = -(dof+1);
1425:   }
1426:   PetscSectionSetUpBC(*gsection);
1427:   if (nroots >= 0) {
1428:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1429:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1430:     if (nroots > pEnd-pStart) {
1431:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1432:     }
1433:   }
1434:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1435:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1436:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1437:     (*gsection)->atlasOff[p] = off;
1438:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1439:   }
1440:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1441:   globalOff -= off;
1442:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1443:     (*gsection)->atlasOff[p] += globalOff;
1444:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1445:   }
1446:   /* Put in negative offsets for ghost points */
1447:   if (nroots >= 0) {
1448:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1449:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1450:     if (nroots > pEnd-pStart) {
1451:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1452:     }
1453:   }
1454:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1455:   PetscFree(neg);
1456:   return(0);
1457: }

1459: typedef struct _n_PetscSectionSym_Label
1460: {
1461:   DMLabel           label;
1462:   PetscCopyMode     *modes;
1463:   PetscInt          *sizes;
1464:   const PetscInt    ***perms;
1465:   const PetscScalar ***rots;
1466:   PetscInt          (*minMaxOrients)[2];
1467:   PetscInt          numStrata; /* numStrata is only increasing, functions as a state */
1468: } PetscSectionSym_Label;

1470: static PetscErrorCode PetscSectionSymLabelReset(PetscSectionSym sym)
1471: {
1472:   PetscInt              i, j;
1473:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1474:   PetscErrorCode        ierr;

1477:   for (i = 0; i <= sl->numStrata; i++) {
1478:     if (sl->modes[i] == PETSC_OWN_POINTER || sl->modes[i] == PETSC_COPY_VALUES) {
1479:       for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1480:         if (sl->perms[i]) {PetscFree(sl->perms[i][j]);}
1481:         if (sl->rots[i]) {PetscFree(sl->rots[i][j]);}
1482:       }
1483:       if (sl->perms[i]) {
1484:         const PetscInt **perms = &sl->perms[i][sl->minMaxOrients[i][0]];

1486:         PetscFree(perms);
1487:       }
1488:       if (sl->rots[i]) {
1489:         const PetscScalar **rots = &sl->rots[i][sl->minMaxOrients[i][0]];

1491:         PetscFree(rots);
1492:       }
1493:     }
1494:   }
1495:   PetscFree5(sl->modes,sl->sizes,sl->perms,sl->rots,sl->minMaxOrients);
1496:   DMLabelDestroy(&sl->label);
1497:   sl->numStrata = 0;
1498:   return(0);
1499: }

1501: static PetscErrorCode PetscSectionSymDestroy_Label(PetscSectionSym sym)
1502: {

1506:   PetscSectionSymLabelReset(sym);
1507:   PetscFree(sym->data);
1508:   return(0);
1509: }

1511: static PetscErrorCode PetscSectionSymView_Label(PetscSectionSym sym, PetscViewer viewer)
1512: {
1513:   PetscSectionSym_Label *sl = (PetscSectionSym_Label *) sym->data;
1514:   PetscBool             isAscii;
1515:   DMLabel               label = sl->label;
1516:   PetscErrorCode        ierr;

1519:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isAscii);
1520:   if (isAscii) {
1521:     PetscInt          i, j, k;
1522:     PetscViewerFormat format;

1524:     PetscViewerGetFormat(viewer,&format);
1525:     if (label) {
1526:       PetscViewerGetFormat(viewer,&format);
1527:       if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1528:         PetscViewerASCIIPushTab(viewer);
1529:         DMLabelView(label, viewer);
1530:         PetscViewerASCIIPopTab(viewer);
1531:       } else {
1532:         PetscViewerASCIIPrintf(viewer,"  Label '%s'\n",sl->label->name);
1533:       }
1534:     } else {
1535:       PetscViewerASCIIPrintf(viewer, "No label given\n");
1536:     }
1537:     PetscViewerASCIIPushTab(viewer);
1538:     for (i = 0; i <= sl->numStrata; i++) {
1539:       PetscInt value = i < sl->numStrata ? label->stratumValues[i] : label->defaultValue;

1541:       if (!(sl->perms[i] || sl->rots[i])) {
1542:         PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point): no symmetries\n", value, sl->sizes[i]);
1543:       } else {
1544:       PetscViewerASCIIPrintf(viewer, "Symmetry for stratum value %D (%D dofs per point):\n", value, sl->sizes[i]);
1545:         PetscViewerASCIIPushTab(viewer);
1546:         PetscViewerASCIIPrintf(viewer, "Orientation range: [%D, %D)\n", sl->minMaxOrients[i][0], sl->minMaxOrients[i][1]);
1547:         if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1548:           PetscViewerASCIIPushTab(viewer);
1549:           for (j = sl->minMaxOrients[i][0]; j < sl->minMaxOrients[i][1]; j++) {
1550:             if (!((sl->perms[i] && sl->perms[i][j]) || (sl->rots[i] && sl->rots[i][j]))) {
1551:               PetscViewerASCIIPrintf(viewer, "Orientation %D: identity\n",j);
1552:             } else {
1553:               PetscInt tab;

1555:               PetscViewerASCIIPrintf(viewer, "Orientation %D:\n",j);
1556:               PetscViewerASCIIPushTab(viewer);
1557:               PetscViewerASCIIGetTab(viewer,&tab);
1558:               if (sl->perms[i] && sl->perms[i][j]) {
1559:                 PetscViewerASCIIPrintf(viewer,"Permutation:");
1560:                 PetscViewerASCIISetTab(viewer,0);
1561:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %D",sl->perms[i][j][k]);}
1562:                 PetscViewerASCIIPrintf(viewer,"\n");
1563:                 PetscViewerASCIISetTab(viewer,tab);
1564:               }
1565:               if (sl->rots[i] && sl->rots[i][j]) {
1566:                 PetscViewerASCIIPrintf(viewer,"Rotations:  ");
1567:                 PetscViewerASCIISetTab(viewer,0);
1568: #if defined(PETSC_USE_COMPLEX)
1569:                 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]));}
1570: #else
1571:                 for (k = 0; k < sl->sizes[i]; k++) {PetscViewerASCIIPrintf(viewer," %+f",sl->rots[i][j][k]);}
1572: #endif
1573:                 PetscViewerASCIIPrintf(viewer,"\n");
1574:                 PetscViewerASCIISetTab(viewer,tab);
1575:               }
1576:               PetscViewerASCIIPopTab(viewer);
1577:             }
1578:           }
1579:           PetscViewerASCIIPopTab(viewer);
1580:         }
1581:         PetscViewerASCIIPopTab(viewer);
1582:       }
1583:     }
1584:     PetscViewerASCIIPopTab(viewer);
1585:   }
1586:   return(0);
1587: }

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

1592:   Logically collective on sym

1594:   Input parameters:
1595: + sym - the section symmetries
1596: - label - the DMLabel describing the types of points

1598:   Level: developer:

1600: .seealso: PetscSectionSymLabelSetStratum(), PetscSectionSymCreateLabel(), PetscSectionGetPointSyms()
1601: @*/
1602: PetscErrorCode PetscSectionSymLabelSetLabel(PetscSectionSym sym, DMLabel label)
1603: {
1604:   PetscSectionSym_Label *sl;

1609:   sl = (PetscSectionSym_Label *) sym->data;
1610:   if (sl->label && sl->label != label) {PetscSectionSymLabelReset(sym);}
1611:   if (label) {
1612:     label->refct++;
1613:     sl->label = label;
1614:     DMLabelGetNumValues(label,&sl->numStrata);
1615:     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);
1616:     PetscMemzero((void *) sl->modes,(sl->numStrata+1)*sizeof(PetscCopyMode));
1617:     PetscMemzero((void *) sl->sizes,(sl->numStrata+1)*sizeof(PetscInt));
1618:     PetscMemzero((void *) sl->perms,(sl->numStrata+1)*sizeof(const PetscInt **));
1619:     PetscMemzero((void *) sl->rots,(sl->numStrata+1)*sizeof(const PetscScalar **));
1620:     PetscMemzero((void *) sl->minMaxOrients,(sl->numStrata+1)*sizeof(PetscInt[2]));
1621:   }
1622:   return(0);
1623: }

1625: /*@C
1626:   PetscSectionSymLabelSetStratum - set the symmetries for the orientations of a stratum

1628:   Logically collective on PetscSectionSym

1630:   InputParameters:
1631: + sys       - the section symmetries
1632: . stratum   - the stratum value in the label that we are assigning symmetries for
1633: . size      - the number of dofs for points in the stratum of the label
1634: . minOrient - the smallest orientation for a point in this stratum
1635: . maxOrient - one greater than the largest orientation for a ppoint in this stratum (i.e., orientations are in the range [minOrient, maxOrient))
1636: . mode      - how sym should copy the perms and rots arrays
1637: . perms     - NULL if there are no permutations, or (maxOrient - minOrient) permutations, one for each orientation.  A NULL permutation is the identity
1638: + 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

1640:   Level: developer

1642: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetPointSyms(), PetscSectionSymCreateLabel()
1643: @*/
1644: PetscErrorCode PetscSectionSymLabelSetStratum(PetscSectionSym sym, PetscInt stratum, PetscInt size, PetscInt minOrient, PetscInt maxOrient, PetscCopyMode mode, const PetscInt **perms, const PetscScalar **rots)
1645: {
1646:   PetscInt       i, j, k;
1647:   PetscSectionSym_Label *sl;

1652:   sl = (PetscSectionSym_Label *) sym->data;
1653:   if (!sl->label) SETERRQ(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_WRONGSTATE,"No label set yet");
1654:   for (i = 0; i <= sl->numStrata; i++) {
1655:     PetscInt value = (i < sl->numStrata) ? sl->label->stratumValues[i] : sl->label->defaultValue;

1657:     if (stratum == value) break;
1658:   }
1659:   if (i > sl->numStrata) SETERRQ2(PetscObjectComm((PetscObject)sym),PETSC_ERR_ARG_OUTOFRANGE,"Stratum %D not found in label %s\n",stratum,sl->label->name);
1660:   sl->sizes[i] = size;
1661:   sl->modes[i] = mode;
1662:   sl->minMaxOrients[i][0] = minOrient;
1663:   sl->minMaxOrients[i][1] = maxOrient;
1664:   if (mode == PETSC_COPY_VALUES) {
1665:     if (perms) {
1666:       PetscInt    **ownPerms;

1668:       PetscCalloc1(maxOrient - minOrient,&ownPerms);
1669:       for (j = 0; j < maxOrient-minOrient; j++) {
1670:         if (perms[j]) {
1671:           PetscMalloc1(size,&ownPerms[j]);
1672:           for (k = 0; k < size; k++) {ownPerms[j][k] = perms[j][k];}
1673:         }
1674:       }
1675:       sl->perms[i] = (const PetscInt **) &ownPerms[-minOrient];
1676:     }
1677:     if (rots) {
1678:       PetscScalar **ownRots;

1680:       PetscCalloc1(maxOrient - minOrient,&ownRots);
1681:       for (j = 0; j < maxOrient-minOrient; j++) {
1682:         if (rots[j]) {
1683:           PetscMalloc1(size,&ownRots[j]);
1684:           for (k = 0; k < size; k++) {ownRots[j][k] = rots[j][k];}
1685:         }
1686:       }
1687:       sl->rots[i] = (const PetscScalar **) &ownRots[-minOrient];
1688:     }
1689:   } else {
1690:     sl->perms[i] = perms ? &perms[-minOrient] : NULL;
1691:     sl->rots[i]  = rots ? &rots[-minOrient] : NULL;
1692:   }
1693:   return(0);
1694: }

1696: static PetscErrorCode PetscSectionSymGetPoints_Label(PetscSectionSym sym, PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt **perms, const PetscScalar **rots)
1697: {
1698:   PetscInt              i, j, numStrata;
1699:   PetscSectionSym_Label *sl;
1700:   DMLabel               label;
1701:   PetscErrorCode        ierr;

1704:   sl = (PetscSectionSym_Label *) sym->data;
1705:   numStrata = sl->numStrata;
1706:   label     = sl->label;
1707:   for (i = 0; i < numPoints; i++) {
1708:     PetscInt point = points[2*i];
1709:     PetscInt ornt  = points[2*i+1];

1711:     for (j = 0; j < numStrata; j++) {
1712:       if (label->validIS[j]) {
1713:         PetscInt k;

1715:         ISLocate(label->points[j],point,&k);
1716:         if (k >= 0) break;
1717:       } else {
1718:         PetscBool has;

1720:         PetscHSetIHas(label->ht[j], point, &has);
1721:         if (has) break;
1722:       }
1723:     }
1724:     if ((sl->minMaxOrients[j][1] > sl->minMaxOrients[j][0]) && (ornt < sl->minMaxOrients[j][0] || ornt >= sl->minMaxOrients[j][1])) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D orientation %D not in range [%D, %D) for stratum %D",point,ornt,sl->minMaxOrients[j][0],sl->minMaxOrients[j][1],j < numStrata ? label->stratumValues[j] : label->defaultValue);
1725:     if (perms) {perms[i] = sl->perms[j] ? sl->perms[j][ornt] : NULL;}
1726:     if (rots) {rots[i]  = sl->rots[j] ? sl->rots[j][ornt] : NULL;}
1727:   }
1728:   return(0);
1729: }

1731: PetscErrorCode PetscSectionSymCreate_Label(PetscSectionSym sym)
1732: {
1733:   PetscSectionSym_Label *sl;
1734:   PetscErrorCode        ierr;

1737:   PetscNewLog(sym,&sl);
1738:   sym->ops->getpoints = PetscSectionSymGetPoints_Label;
1739:   sym->ops->view      = PetscSectionSymView_Label;
1740:   sym->ops->destroy   = PetscSectionSymDestroy_Label;
1741:   sym->data           = (void *) sl;
1742:   return(0);
1743: }

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

1748:   Collective

1750:   Input Parameters:
1751: + comm - the MPI communicator for the new symmetry
1752: - label - the label defining the strata

1754:   Output Parameters:
1755: . sym - the section symmetries

1757:   Level: developer

1759: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
1760: @*/
1761: PetscErrorCode PetscSectionSymCreateLabel(MPI_Comm comm, DMLabel label, PetscSectionSym *sym)
1762: {
1763:   PetscErrorCode        ierr;

1766:   DMInitializePackage();
1767:   PetscSectionSymCreate(comm,sym);
1768:   PetscSectionSymSetType(*sym,PETSCSECTIONSYMLABEL);
1769:   PetscSectionSymLabelSetLabel(*sym,label);
1770:   return(0);
1771: }