Actual source code: dmlabel.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmlabelimpl.h>   /*I      "petscdmlabel.h"   I*/
  2: #include <petsc/private/isimpl.h>        /*I      "petscis.h"        I*/
  3: #include <petscsf.h>

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

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

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

 16:   Level: beginner

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

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

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

 45: /*
 46:   DMLabelMakeValid_Private - Transfer stratum data from the hash format to the sorted list format

 48:   Input parameter:
 49: + label - The DMLabel
 50: - v - The stratum value

 52:   Output parameter:
 53: . label - The DMLabel with stratum in sorted list format

 55:   Level: developer

 57: .seealso: DMLabelCreate()
 58: */
 59: static PetscErrorCode DMLabelMakeValid_Private(DMLabel label, PetscInt v)
 60: {
 61:   PetscInt       off;

 64:   if (label->arrayValid[v]) return 0;
 65:   if (v >= label->numStrata) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Trying to access invalid stratum %D in DMLabelMakeValid_Private\n", v);
 67:   PetscHashISize(label->ht[v], label->stratumSizes[v]);

 69:   PetscMalloc1(label->stratumSizes[v], &label->points[v]);
 70:   off = 0;
 71:   PetscHashIGetKeys(label->ht[v], &off, &(label->points[v][0]));
 72:   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]);
 73:   PetscHashIClear(label->ht[v]);
 74:   PetscSortInt(label->stratumSizes[v], label->points[v]);
 75:   if (label->bt) {
 76:     PetscInt p;

 78:     for (p = 0; p < label->stratumSizes[v]; ++p) {
 79:       const PetscInt point = label->points[v][p];

 81:       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);
 82:       PetscBTSet(label->bt, point - label->pStart);
 83:     }
 84:   }
 85:   label->arrayValid[v] = PETSC_TRUE;
 86:   ++label->state;
 87:   return(0);
 88: }

 92: /*
 93:   DMLabelMakeAllValid_Private - Transfer all strata from the hash format to the sorted list format

 95:   Input parameter:
 96: . label - The DMLabel

 98:   Output parameter:
 99: . label - The DMLabel with all strata in sorted list format

101:   Level: developer

103: .seealso: DMLabelCreate()
104: */
105: static PetscErrorCode DMLabelMakeAllValid_Private(DMLabel label)
106: {
107:   PetscInt       v;

111:   for (v = 0; v < label->numStrata; v++){
112:     DMLabelMakeValid_Private(label, v);
113:   }
114:   return(0);
115: }

119: /*
120:   DMLabelMakeInvalid_Private - Transfer stratum data from the sorted list format to the hash format

122:   Input parameter:
123: + label - The DMLabel
124: - v - The stratum value

126:   Output parameter:
127: . label - The DMLabel with stratum in hash format

129:   Level: developer

131: .seealso: DMLabelCreate()
132: */
133: static PetscErrorCode DMLabelMakeInvalid_Private(DMLabel label, PetscInt v)
134: {
135:   PETSC_UNUSED PetscHashIIter ret, iter;
136:   PetscInt                    p;

140:   if (!label->arrayValid[v]) return(0);
141:   for (p = 0; p < label->stratumSizes[v]; ++p) PetscHashIPut(label->ht[v], label->points[v][p], ret, iter);
142:   PetscFree(label->points[v]);
143:   label->arrayValid[v] = PETSC_FALSE;
144:   return(0);
145: }

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

159: PetscErrorCode DMLabelAddStratum(DMLabel label, PetscInt value)
160: {
161:   PetscInt    v, *tmpV, *tmpS, **tmpP;
162:   PetscHashI *tmpH;
163:   PetscBool  *tmpB;


168:   PetscMalloc1((label->numStrata+1), &tmpV);
169:   PetscMalloc1((label->numStrata+1), &tmpS);
170:   PetscMalloc1((label->numStrata+1), &tmpH);
171:   PetscMalloc1((label->numStrata+1), &tmpP);
172:   PetscMalloc1((label->numStrata+1), &tmpB);
173:   for (v = 0; v < label->numStrata; ++v) {
174:     tmpV[v] = label->stratumValues[v];
175:     tmpS[v] = label->stratumSizes[v];
176:     tmpH[v] = label->ht[v];
177:     tmpP[v] = label->points[v];
178:     tmpB[v] = label->arrayValid[v];
179:   }
180:   tmpV[v] = value;
181:   tmpS[v] = 0;
182:   PetscHashICreate(tmpH[v]);
183:   tmpP[v] = NULL;
184:   tmpB[v] = PETSC_TRUE;
185:   ++label->numStrata;
186:   PetscFree(label->stratumValues);
187:   PetscFree(label->stratumSizes);
188:   PetscFree(label->ht);
189:   PetscFree(label->points);
190:   PetscFree(label->arrayValid);
191:   label->stratumValues = tmpV;
192:   label->stratumSizes  = tmpS;
193:   label->ht            = tmpH;
194:   label->points        = tmpP;
195:   label->arrayValid    = tmpB;

197:   return(0);
198: }

202: /*@C
203:   DMLabelGetName - Return the name of a DMLabel object

205:   Input parameter:
206: . label - The DMLabel

208:   Output parameter:
209: . name - The label name

211:   Level: beginner

213: .seealso: DMLabelCreate()
214: @*/
215: PetscErrorCode DMLabelGetName(DMLabel label, const char **name)
216: {
219:   *name = label->name;
220:   return(0);
221: }

225: static PetscErrorCode DMLabelView_Ascii(DMLabel label, PetscViewer viewer)
226: {
227:   PetscInt       v;
228:   PetscMPIInt    rank;

232:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
233:   PetscViewerASCIIPushSynchronized(viewer);
234:   if (label) {
235:     PetscViewerASCIIPrintf(viewer, "Label '%s':\n", label->name);
236:     if (label->bt) {PetscViewerASCIIPrintf(viewer, "  Index has been calculated in [%D, %D)\n", label->pStart, label->pEnd);}
237:     for (v = 0; v < label->numStrata; ++v) {
238:       const PetscInt value = label->stratumValues[v];
239:       PetscInt       p;

241:       for (p = 0; p < label->stratumSizes[v]; ++p) {
242:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D (%D)\n", rank, label->points[v][p], value);
243:       }
244:     }
245:   }
246:   PetscViewerFlush(viewer);
247:   PetscViewerASCIIPopSynchronized(viewer);
248:   return(0);
249: }

253: /*@C
254:   DMLabelView - View the label

256:   Input Parameters:
257: + label - The DMLabel
258: - viewer - The PetscViewer

260:   Level: intermediate

262: .seealso: DMLabelCreate(), DMLabelDestroy()
263: @*/
264: PetscErrorCode DMLabelView(DMLabel label, PetscViewer viewer)
265: {
266:   PetscBool      iascii;

271:   if (label) {DMLabelMakeAllValid_Private(label);}
272:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
273:   if (iascii) {
274:     DMLabelView_Ascii(label, viewer);
275:   }
276:   return(0);
277: }

281: PetscErrorCode DMLabelDestroy(DMLabel *label)
282: {
283:   PetscInt       v;

287:   if (!(*label)) return(0);
288:   if (--(*label)->refct > 0) return(0);
289:   PetscFree((*label)->name);
290:   PetscFree((*label)->stratumValues);
291:   PetscFree((*label)->stratumSizes);
292:   for (v = 0; v < (*label)->numStrata; ++v) {PetscFree((*label)->points[v]);}
293:   PetscFree((*label)->points);
294:   PetscFree((*label)->arrayValid);
295:   if ((*label)->ht) {
296:     for (v = 0; v < (*label)->numStrata; ++v) {PetscHashIDestroy((*label)->ht[v]);}
297:     PetscFree((*label)->ht);
298:   }
299:   PetscBTDestroy(&(*label)->bt);
300:   PetscFree(*label);
301:   return(0);
302: }

306: PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
307: {
308:   PetscInt       v, q;

312:   DMLabelMakeAllValid_Private(label);
313:   PetscNew(labelnew);
314:   PetscStrallocpy(label->name, &(*labelnew)->name);

316:   (*labelnew)->refct        = 1;
317:   (*labelnew)->numStrata    = label->numStrata;
318:   (*labelnew)->defaultValue = label->defaultValue;
319:   if (label->numStrata) {
320:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);
321:     PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);
322:     PetscMalloc1(label->numStrata, &(*labelnew)->ht);
323:     PetscMalloc1(label->numStrata, &(*labelnew)->points);
324:     PetscMalloc1(label->numStrata, &(*labelnew)->arrayValid);
325:     /* Could eliminate unused space here */
326:     for (v = 0; v < label->numStrata; ++v) {
327:       PetscMalloc1(label->stratumSizes[v], &(*labelnew)->points[v]);
328:       PetscHashICreate((*labelnew)->ht[v]);
329:       (*labelnew)->arrayValid[v]     = PETSC_TRUE;
330:       (*labelnew)->stratumValues[v]  = label->stratumValues[v];
331:       (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
332:       for (q = 0; q < label->stratumSizes[v]; ++q) {
333:         (*labelnew)->points[v][q] = label->points[v][q];
334:       }
335:     }
336:   }
337:   (*labelnew)->pStart = -1;
338:   (*labelnew)->pEnd   = -1;
339:   (*labelnew)->bt     = NULL;
340:   return(0);
341: }

345: /* This can be hooked into SetValue(),  ClearValue(), etc. for updating */
346: PetscErrorCode DMLabelCreateIndex(DMLabel label, PetscInt pStart, PetscInt pEnd)
347: {
348:   PetscInt       v;

352:   DMLabelMakeAllValid_Private(label);
353:   if (label->bt) {PetscBTDestroy(&label->bt);}
354:   label->pStart = pStart;
355:   label->pEnd   = pEnd;
356:   PetscBTCreate(pEnd - pStart, &label->bt);
357:   PetscBTMemzero(pEnd - pStart, label->bt);
358:   for (v = 0; v < label->numStrata; ++v) {
359:     PetscInt i;

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

364:       if ((point < pStart) || (point >= pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %D is not in [%D, %D)", point, pStart, pEnd);
365:       PetscBTSet(label->bt, point - pStart);
366:     }
367:   }
368:   return(0);
369: }

373: PetscErrorCode DMLabelDestroyIndex(DMLabel label)
374: {

378:   label->pStart = -1;
379:   label->pEnd   = -1;
380:   if (label->bt) {PetscBTDestroy(&label->bt);}
381:   return(0);
382: }

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

389:   Input Parameters:
390: + label - the DMLabel
391: - value - the value

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

396:   Level: developer

398: .seealso: DMLabelHasPoint(), DMLabelGetValue(), DMLabelSetValue()
399: @*/
400: PetscErrorCode DMLabelHasValue(DMLabel label, PetscInt value, PetscBool *contains)
401: {
402:   PetscInt v;

406:   for (v = 0; v < label->numStrata; ++v) {
407:     if (value == label->stratumValues[v]) break;
408:   }
409:   *contains = (v < label->numStrata ? PETSC_TRUE : PETSC_FALSE);
410:   return(0);
411: }

415: /*@
416:   DMLabelHasPoint - Determine whether a label assigns a value to a point

418:   Input Parameters:
419: + label - the DMLabel
420: - point - the point

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

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

427:   Level: developer

429: .seealso: DMLabelCreateIndex(), DMLabelGetValue(), DMLabelSetValue()
430: @*/
431: PetscErrorCode DMLabelHasPoint(DMLabel label, PetscInt point, PetscBool *contains)
432: {

437:   DMLabelMakeAllValid_Private(label);
438: #if defined(PETSC_USE_DEBUG)
439:   if (!label->bt) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call DMLabelCreateIndex() before DMLabelHasPoint()");
440:   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);
441: #endif
442:   *contains = PetscBTLookup(label->bt, point - label->pStart) ? PETSC_TRUE : PETSC_FALSE;
443:   return(0);
444: }

448: /*@
449:   DMLabelStratumHasPoint - Return true if the stratum contains a point

451:   Input Parameters:
452: + label - the DMLabel
453: . value - the stratum value
454: - point - the point

456:   Output Parameter:
457: . contains - true if the stratum contains the point

459:   Level: intermediate

461: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue()
462: @*/
463: PetscErrorCode DMLabelStratumHasPoint(DMLabel label, PetscInt value, PetscInt point, PetscBool *contains)
464: {
465:   PetscInt       v;

470:   *contains = PETSC_FALSE;
471:   for (v = 0; v < label->numStrata; ++v) {
472:     if (label->stratumValues[v] == value) {
473:       if (label->arrayValid[v]) {
474:         PetscInt i;

476:         PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
477:         if (i >= 0) {
478:           *contains = PETSC_TRUE;
479:           break;
480:         }
481:       } else {
482:         PetscBool has;

484:         PetscHashIHasKey(label->ht[v], point, has);
485:         if (has) {
486:           *contains = PETSC_TRUE;
487:           break;
488:         }
489:       }
490:     }
491:   }
492:   return(0);
493: }

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

501:   Input parameter:
502: . label - a DMLabel object

504:   Output parameter:
505: . defaultValue - the default value

507:   Level: beginner

509: .seealso: DMLabelSetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
510: */
511: PetscErrorCode DMLabelGetDefaultValue(DMLabel label, PetscInt *defaultValue)
512: {
514:   *defaultValue = label->defaultValue;
515:   return(0);
516: }

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

524:   Input parameter:
525: . label - a DMLabel object

527:   Output parameter:
528: . defaultValue - the default value

530:   Level: beginner

532: .seealso: DMLabelGetDefaultValue(), DMLabelGetValue(), DMLabelSetValue()
533: */
534: PetscErrorCode DMLabelSetDefaultValue(DMLabel label, PetscInt defaultValue)
535: {
537:   label->defaultValue = defaultValue;
538:   return(0);
539: }

543: /*@
544:   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())

546:   Input Parameters:
547: + label - the DMLabel
548: - point - the point

550:   Output Parameter:
551: . value - The point value, or -1

553:   Level: intermediate

555: .seealso: DMLabelCreate(), DMLabelSetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
556: @*/
557: PetscErrorCode DMLabelGetValue(DMLabel label, PetscInt point, PetscInt *value)
558: {
559:   PetscInt       v;

564:   *value = label->defaultValue;
565:   for (v = 0; v < label->numStrata; ++v) {
566:     if (label->arrayValid[v]) {
567:       PetscInt i;

569:       PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &i);
570:       if (i >= 0) {
571:         *value = label->stratumValues[v];
572:         break;
573:       }
574:     } else {
575:       PetscBool has;

577:       PetscHashIHasKey(label->ht[v], point, has);
578:       if (has) {
579:         *value = label->stratumValues[v];
580:         break;
581:       }
582:     }
583:   }
584:   return(0);
585: }

589: /*@
590:   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.

592:   Input Parameters:
593: + label - the DMLabel
594: . point - the point
595: - value - The point value

597:   Level: intermediate

599: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelClearValue(), DMLabelGetDefaultValue(), DMLabelSetDefaultValue()
600: @*/
601: PetscErrorCode DMLabelSetValue(DMLabel label, PetscInt point, PetscInt value)
602: {
603:   PETSC_UNUSED PetscHashIIter iter, ret;
604:   PetscInt                    v;
605:   PetscErrorCode              ierr;

608:   /* Find, or add, label value */
609:   if (value == label->defaultValue) return(0);
610:   for (v = 0; v < label->numStrata; ++v) {
611:     if (label->stratumValues[v] == value) break;
612:   }
613:   /* Create new table */
614:   if (v >= label->numStrata) {DMLabelAddStratum(label, value);}
615:   DMLabelMakeInvalid_Private(label, v);
616:   /* Set key */
617:   PetscHashIPut(label->ht[v], point, ret, iter);
618:   return(0);
619: }

623: /*@
624:   DMLabelClearValue - Clear the value a label assigns to a point

626:   Input Parameters:
627: + label - the DMLabel
628: . point - the point
629: - value - The point value

631:   Level: intermediate

633: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue()
634: @*/
635: PetscErrorCode DMLabelClearValue(DMLabel label, PetscInt point, PetscInt value)
636: {
637:   PetscInt       v, p;

641:   /* Find label value */
642:   for (v = 0; v < label->numStrata; ++v) {
643:     if (label->stratumValues[v] == value) break;
644:   }
645:   if (v >= label->numStrata) return(0);
646:   if (label->arrayValid[v]) {
647:     /* Check whether point exists */
648:     PetscFindInt(point, label->stratumSizes[v], &label->points[v][0], &p);
649:     if (p >= 0) {
650:       PetscMemmove(&label->points[v][p], &label->points[v][p+1], (label->stratumSizes[v]-p-1) * sizeof(PetscInt));
651:       --label->stratumSizes[v];
652:       if (label->bt) {
653:         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);
654:         PetscBTClear(label->bt, point - label->pStart);
655:       }
656:     }
657:   } else {
658:     PetscHashIDelKey(label->ht[v], point);
659:   }
660:   return(0);
661: }

665: /*@
666:   DMLabelInsertIS - Set all points in the IS to a value

668:   Input Parameters:
669: + label - the DMLabel
670: . is    - the point IS
671: - value - The point value

673:   Level: intermediate

675: .seealso: DMLabelCreate(), DMLabelGetValue(), DMLabelSetValue(), DMLabelClearValue()
676: @*/
677: PetscErrorCode DMLabelInsertIS(DMLabel label, IS is, PetscInt value)
678: {
679:   const PetscInt *points;
680:   PetscInt        n, p;
681:   PetscErrorCode  ierr;

685:   ISGetLocalSize(is, &n);
686:   ISGetIndices(is, &points);
687:   for (p = 0; p < n; ++p) {DMLabelSetValue(label, points[p], value);}
688:   ISRestoreIndices(is, &points);
689:   return(0);
690: }

694: PetscErrorCode DMLabelGetNumValues(DMLabel label, PetscInt *numValues)
695: {
698:   *numValues = label->numStrata;
699:   return(0);
700: }

704: PetscErrorCode DMLabelGetValueIS(DMLabel label, IS *values)
705: {

710:   ISCreateGeneral(PETSC_COMM_SELF, label->numStrata, label->stratumValues, PETSC_USE_POINTER, values);
711:   return(0);
712: }

716: PetscErrorCode DMLabelHasStratum(DMLabel label, PetscInt value, PetscBool *exists)
717: {
718:   PetscInt v;

722:   *exists = PETSC_FALSE;
723:   for (v = 0; v < label->numStrata; ++v) {
724:     if (label->stratumValues[v] == value) {
725:       *exists = PETSC_TRUE;
726:       break;
727:     }
728:   }
729:   return(0);
730: }

734: PetscErrorCode DMLabelGetStratumSize(DMLabel label, PetscInt value, PetscInt *size)
735: {
736:   PetscInt       v;

741:   *size = 0;
742:   for (v = 0; v < label->numStrata; ++v) {
743:     if (label->stratumValues[v] == value) {
744:       DMLabelMakeValid_Private(label, v);
745:       *size = label->stratumSizes[v];
746:       break;
747:     }
748:   }
749:   return(0);
750: }

754: PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *start, PetscInt *end)
755: {
756:   PetscInt       v;

762:   for (v = 0; v < label->numStrata; ++v) {
763:     if (label->stratumValues[v] != value) continue;
764:     DMLabelMakeValid_Private(label, v);
765:     if (label->stratumSizes[v]  <= 0)     break;
766:     if (start) *start = label->points[v][0];
767:     if (end)   *end   = label->points[v][label->stratumSizes[v]-1]+1;
768:     break;
769:   }
770:   return(0);
771: }

775: PetscErrorCode DMLabelGetStratumIS(DMLabel label, PetscInt value, IS *points)
776: {
777:   PetscInt       v;

782:   *points = NULL;
783:   for (v = 0; v < label->numStrata; ++v) {
784:     if (label->stratumValues[v] == value) {
785:       DMLabelMakeValid_Private(label, v);
786:       if (label->arrayValid[v]) {
787:         ISCreateGeneral(PETSC_COMM_SELF, label->stratumSizes[v], &label->points[v][0], PETSC_COPY_VALUES, points);
788:         PetscObjectSetName((PetscObject) *points, "indices");
789:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Need to implement this to speedup Stratify");
790:       break;
791:     }
792:   }
793:   return(0);
794: }

798: PetscErrorCode DMLabelClearStratum(DMLabel label, PetscInt value)
799: {
800:   PetscInt       v;

804:   for (v = 0; v < label->numStrata; ++v) {
805:     if (label->stratumValues[v] == value) break;
806:   }
807:   if (v >= label->numStrata) return(0);
808:   if (label->bt) {
809:     PetscInt i;

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

814:       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);
815:       PetscBTClear(label->bt, point - label->pStart);
816:     }
817:   }
818:   if (label->arrayValid[v]) {
819:     label->stratumSizes[v] = 0;
820:   } else {
821:     PetscHashIClear(label->ht[v]);
822:   }
823:   return(0);
824: }

828: PetscErrorCode DMLabelFilter(DMLabel label, PetscInt start, PetscInt end)
829: {
830:   PetscInt       v;

834:   DMLabelMakeAllValid_Private(label);
835:   label->pStart = start;
836:   label->pEnd   = end;
837:   if (label->bt) {PetscBTDestroy(&label->bt);}
838:   /* Could squish offsets, but would only make sense if I reallocate the storage */
839:   for (v = 0; v < label->numStrata; ++v) {
840:     PetscInt off, q;

842:     for (off = 0, q = 0; q < label->stratumSizes[v]; ++q) {
843:       const PetscInt point = label->points[v][q];

845:       if ((point < start) || (point >= end)) continue;
846:       label->points[v][off++] = point;
847:     }
848:     label->stratumSizes[v] = off;
849:   }
850:   DMLabelCreateIndex(label, start, end);
851:   return(0);
852: }

856: PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
857: {
858:   const PetscInt *perm;
859:   PetscInt        numValues, numPoints, v, q;
860:   PetscErrorCode  ierr;

863:   DMLabelMakeAllValid_Private(label);
864:   DMLabelDuplicate(label, labelNew);
865:   DMLabelGetNumValues(*labelNew, &numValues);
866:   ISGetLocalSize(permutation, &numPoints);
867:   ISGetIndices(permutation, &perm);
868:   for (v = 0; v < numValues; ++v) {
869:     const PetscInt size   = (*labelNew)->stratumSizes[v];

871:     for (q = 0; q < size; ++q) {
872:       const PetscInt point = (*labelNew)->points[v][q];

874:       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);
875:       (*labelNew)->points[v][q] = perm[point];
876:     }
877:     PetscSortInt(size, &(*labelNew)->points[v][0]);
878:   }
879:   ISRestoreIndices(permutation, &perm);
880:   if (label->bt) {
881:     PetscBTDestroy(&label->bt);
882:     DMLabelCreateIndex(label, label->pStart, label->pEnd);
883:   }
884:   return(0);
885: }

889: PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
890: {
891:   MPI_Comm       comm;
892:   PetscInt       s, l, nroots, nleaves, dof, offset, size;
893:   PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
894:   PetscSection   rootSection;
895:   PetscSF        labelSF;

899:   if (label) {DMLabelMakeAllValid_Private(label);}
900:   PetscObjectGetComm((PetscObject)sf, &comm);
901:   /* Build a section of stratum values per point, generate the according SF
902:      and distribute point-wise stratum values to leaves. */
903:   PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);
904:   PetscSectionCreate(comm, &rootSection);
905:   PetscSectionSetChart(rootSection, 0, nroots);
906:   if (label) {
907:     for (s = 0; s < label->numStrata; ++s) {
908:       for (l = 0; l < label->stratumSizes[s]; l++) {
909:         PetscSectionGetDof(rootSection, label->points[s][l], &dof);
910:         PetscSectionSetDof(rootSection, label->points[s][l], dof+1);
911:       }
912:     }
913:   }
914:   PetscSectionSetUp(rootSection);
915:   /* Create a point-wise array of stratum values */
916:   PetscSectionGetStorageSize(rootSection, &size);
917:   PetscMalloc1(size, &rootStrata);
918:   PetscCalloc1(nroots, &rootIdx);
919:   if (label) {
920:     for (s = 0; s < label->numStrata; ++s) {
921:       for (l = 0; l < label->stratumSizes[s]; l++) {
922:         const PetscInt p = label->points[s][l];
923:         PetscSectionGetOffset(rootSection, p, &offset);
924:         rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
925:       }
926:     }
927:   }
928:   /* Build SF that maps label points to remote processes */
929:   PetscSectionCreate(comm, leafSection);
930:   PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);
931:   PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);
932:   PetscFree(remoteOffsets);
933:   /* Send the strata for each point over the derived SF */
934:   PetscSectionGetStorageSize(*leafSection, &size);
935:   PetscMalloc1(size, leafStrata);
936:   PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);
937:   PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);
938:   /* Clean up */
939:   PetscFree(rootStrata);
940:   PetscFree(rootIdx);
941:   PetscSectionDestroy(&rootSection);
942:   PetscSFDestroy(&labelSF);
943:   return(0);
944: }

948: PetscErrorCode DMLabelDistribute(DMLabel label, PetscSF sf, DMLabel *labelNew)
949: {
950:   MPI_Comm       comm;
951:   PetscSection   leafSection;
952:   PetscInt       p, pStart, pEnd, s, size, dof, offset, stratum;
953:   PetscInt      *leafStrata, *strataIdx;
954:   char          *name;
955:   PetscInt       nameSize;
956:   PetscHashI     stratumHash;
957:   PETSC_UNUSED   PetscHashIIter ret, iter;
958:   size_t         len = 0;
959:   PetscMPIInt    rank;

963:   if (label) {DMLabelMakeAllValid_Private(label);}
964:   PetscObjectGetComm((PetscObject)sf, &comm);
965:   MPI_Comm_rank(comm, &rank);
966:   /* Bcast name */
967:   if (!rank) {PetscStrlen(label->name, &len);}
968:   nameSize = len;
969:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
970:   PetscMalloc1(nameSize+1, &name);
971:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
972:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
973:   DMLabelCreate(name, labelNew);
974:   PetscFree(name);
975:   /* Bcast defaultValue */
976:   if (!rank) (*labelNew)->defaultValue = label->defaultValue;
977:   MPI_Bcast(&(*labelNew)->defaultValue, 1, MPIU_INT, 0, comm);
978:   /* Distribute stratum values over the SF and get the point mapping on the receiver */
979:   DMLabelDistribute_Internal(label, sf, &leafSection, &leafStrata);
980:   /* Determine received stratum values and initialise new label*/
981:   PetscHashICreate(stratumHash);
982:   PetscSectionGetStorageSize(leafSection, &size);
983:   for (p = 0; p < size; ++p) PetscHashIPut(stratumHash, leafStrata[p], ret, iter);
984:   PetscHashISize(stratumHash, (*labelNew)->numStrata);
985:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->arrayValid);
986:   for (s = 0; s < (*labelNew)->numStrata; ++s) (*labelNew)->arrayValid[s] = PETSC_TRUE;
987:   PetscMalloc1((*labelNew)->numStrata, &(*labelNew)->stratumValues);
988:   /* Turn leafStrata into indices rather than stratum values */
989:   offset = 0;
990:   PetscHashIGetKeys(stratumHash, &offset, (*labelNew)->stratumValues);
991:   for (p = 0; p < size; ++p) {
992:     for (s = 0; s < (*labelNew)->numStrata; ++s) {
993:       if (leafStrata[p] == (*labelNew)->stratumValues[s]) {leafStrata[p] = s; break;}
994:     }
995:   }
996:   /* Rebuild the point strata on the receiver */
997:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->stratumSizes);
998:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
999:   for (p=pStart; p<pEnd; p++) {
1000:     PetscSectionGetDof(leafSection, p, &dof);
1001:     PetscSectionGetOffset(leafSection, p, &offset);
1002:     for (s=0; s<dof; s++) {
1003:       (*labelNew)->stratumSizes[leafStrata[offset+s]]++;
1004:     }
1005:   }
1006:   PetscCalloc1((*labelNew)->numStrata,&(*labelNew)->ht);
1007:   PetscMalloc1((*labelNew)->numStrata,&(*labelNew)->points);
1008:   for (s = 0; s < (*labelNew)->numStrata; ++s) {
1009:     PetscHashICreate((*labelNew)->ht[s]);
1010:     PetscMalloc1((*labelNew)->stratumSizes[s], &(*labelNew)->points[s]);
1011:   }
1012:   /* Insert points into new strata */
1013:   PetscCalloc1((*labelNew)->numStrata, &strataIdx);
1014:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1015:   for (p=pStart; p<pEnd; p++) {
1016:     PetscSectionGetDof(leafSection, p, &dof);
1017:     PetscSectionGetOffset(leafSection, p, &offset);
1018:     for (s=0; s<dof; s++) {
1019:       stratum = leafStrata[offset+s];
1020:       (*labelNew)->points[stratum][strataIdx[stratum]++] = p;
1021:     }
1022:   }
1023:   PetscHashIDestroy(stratumHash);
1024:   PetscFree(leafStrata);
1025:   PetscFree(strataIdx);
1026:   PetscSectionDestroy(&leafSection);
1027:   return(0);
1028: }

1032: /*@
1033:   DMLabelGather - Gather all label values from leafs into roots

1035:   Input Parameters:
1036: + label - the DMLabel
1037: . point - the Star Forest point communication map

1039:   Input Parameters:
1040: + label - the new DMLabel with localised leaf values

1042:   Level: developer

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

1046: .seealso: DMLabelDistribute()
1047: @*/
1048: PetscErrorCode DMLabelGather(DMLabel label, PetscSF sf, DMLabel *labelNew)
1049: {
1050:   MPI_Comm       comm;
1051:   PetscSection   rootSection;
1052:   PetscSF        sfLabel;
1053:   PetscSFNode   *rootPoints, *leafPoints;
1054:   PetscInt       p, s, d, nroots, nleaves, nmultiroots, idx, dof, offset;
1055:   const PetscInt *rootDegree, *ilocal;
1056:   PetscInt       *rootStrata;
1057:   char          *name;
1058:   PetscInt       nameSize;
1059:   size_t         len = 0;
1060:   PetscMPIInt    rank, numProcs;

1064:   PetscObjectGetComm((PetscObject)sf, &comm);
1065:   MPI_Comm_rank(comm, &rank);
1066:   MPI_Comm_size(comm, &numProcs);
1067:   /* Bcast name */
1068:   if (!rank) {PetscStrlen(label->name, &len);}
1069:   nameSize = len;
1070:   MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
1071:   PetscMalloc1(nameSize+1, &name);
1072:   if (!rank) {PetscMemcpy(name, label->name, nameSize+1);}
1073:   MPI_Bcast(name, nameSize+1, MPI_CHAR, 0, comm);
1074:   DMLabelCreate(name, labelNew);
1075:   PetscFree(name);
1076:   /* Gather rank/index pairs of leaves into local roots to build
1077:      an inverse, multi-rooted SF. Note that this ignores local leaf
1078:      indexing due to the use of the multiSF in PetscSFGather. */
1079:   PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, NULL);
1080:   PetscMalloc1(nleaves, &leafPoints);
1081:   for (p = 0; p < nleaves; p++) {
1082:     leafPoints[p].index = ilocal[p];
1083:     leafPoints[p].rank = rank;
1084:   }
1085:   PetscSFComputeDegreeBegin(sf, &rootDegree);
1086:   PetscSFComputeDegreeEnd(sf, &rootDegree);
1087:   for (p = 0, nmultiroots = 0; p < nroots; ++p) nmultiroots += rootDegree[p];
1088:   PetscMalloc1(nmultiroots, &rootPoints);
1089:   PetscSFGatherBegin(sf, MPIU_2INT, leafPoints, rootPoints);
1090:   PetscSFGatherEnd(sf, MPIU_2INT, leafPoints, rootPoints);
1091:   PetscSFCreate(comm,& sfLabel);
1092:   PetscSFSetGraph(sfLabel, nroots, nmultiroots, NULL, PETSC_OWN_POINTER, rootPoints, PETSC_OWN_POINTER);
1093:   /* Migrate label over inverted SF to pull stratum values at leaves into roots. */
1094:   DMLabelDistribute_Internal(label, sfLabel, &rootSection, &rootStrata);
1095:   /* Rebuild the point strata on the receiver */
1096:   for (p = 0, idx = 0; p < nroots; p++) {
1097:     for (d = 0; d < rootDegree[p]; d++) {
1098:       PetscSectionGetDof(rootSection, idx+d, &dof);
1099:       PetscSectionGetOffset(rootSection, idx+d, &offset);
1100:       for (s = 0; s < dof; s++) {DMLabelSetValue(*labelNew, p, rootStrata[offset+s]);}
1101:     }
1102:     idx += rootDegree[p];
1103:   }
1104:   PetscFree(leafPoints);
1105:   PetscFree(rootStrata);
1106:   PetscSectionDestroy(&rootSection);
1107:   PetscSFDestroy(&sfLabel);
1108:   return(0);
1109: }

1113: PetscErrorCode DMLabelConvertToSection(DMLabel label, PetscSection *section, IS *is)
1114: {
1115:   IS              vIS;
1116:   const PetscInt *values;
1117:   PetscInt       *points;
1118:   PetscInt        nV, vS = 0, vE = 0, v, N;
1119:   PetscErrorCode  ierr;

1122:   DMLabelGetNumValues(label, &nV);
1123:   DMLabelGetValueIS(label, &vIS);
1124:   ISGetIndices(vIS, &values);
1125:   if (nV) {vS = values[0]; vE = values[0]+1;}
1126:   for (v = 1; v < nV; ++v) {
1127:     vS = PetscMin(vS, values[v]);
1128:     vE = PetscMax(vE, values[v]+1);
1129:   }
1130:   PetscSectionCreate(PETSC_COMM_SELF, section);
1131:   PetscSectionSetChart(*section, vS, vE);
1132:   for (v = 0; v < nV; ++v) {
1133:     PetscInt n;

1135:     DMLabelGetStratumSize(label, values[v], &n);
1136:     PetscSectionSetDof(*section, values[v], n);
1137:   }
1138:   PetscSectionSetUp(*section);
1139:   PetscSectionGetStorageSize(*section, &N);
1140:   PetscMalloc1(N, &points);
1141:   for (v = 0; v < nV; ++v) {
1142:     IS              is;
1143:     const PetscInt *spoints;
1144:     PetscInt        dof, off, p;

1146:     PetscSectionGetDof(*section, values[v], &dof);
1147:     PetscSectionGetOffset(*section, values[v], &off);
1148:     DMLabelGetStratumIS(label, values[v], &is);
1149:     ISGetIndices(is, &spoints);
1150:     for (p = 0; p < dof; ++p) points[off+p] = spoints[p];
1151:     ISRestoreIndices(is, &spoints);
1152:     ISDestroy(&is);
1153:   }
1154:   ISRestoreIndices(vIS, &values);
1155:   ISDestroy(&vIS);
1156:   ISCreateGeneral(PETSC_COMM_SELF, N, points, PETSC_OWN_POINTER, is);
1157:   return(0);
1158: }

1162: /*@C
1163:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
1164:   the local section and an SF describing the section point overlap.

1166:   Input Parameters:
1167:   + s - The PetscSection for the local field layout
1168:   . sf - The SF describing parallel layout of the section points
1169:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1170:   . label - The label specifying the points
1171:   - labelValue - The label stratum specifying the points

1173:   Output Parameter:
1174:   . gsection - The PetscSection for the global field layout

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

1178:   Level: developer

1180: .seealso: PetscSectionCreate()
1181: @*/
1182: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
1183: {
1184:   PetscInt      *neg = NULL, *tmpOff = NULL;
1185:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

1189:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1190:   PetscSectionGetChart(s, &pStart, &pEnd);
1191:   PetscSectionSetChart(*gsection, pStart, pEnd);
1192:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1193:   if (nroots >= 0) {
1194:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1195:     PetscCalloc1(nroots, &neg);
1196:     if (nroots > pEnd-pStart) {
1197:       PetscCalloc1(nroots, &tmpOff);
1198:     } else {
1199:       tmpOff = &(*gsection)->atlasDof[-pStart];
1200:     }
1201:   }
1202:   /* Mark ghost points with negative dof */
1203:   for (p = pStart; p < pEnd; ++p) {
1204:     PetscInt value;

1206:     DMLabelGetValue(label, p, &value);
1207:     if (value != labelValue) continue;
1208:     PetscSectionGetDof(s, p, &dof);
1209:     PetscSectionSetDof(*gsection, p, dof);
1210:     PetscSectionGetConstraintDof(s, p, &cdof);
1211:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1212:     if (neg) neg[p] = -(dof+1);
1213:   }
1214:   PetscSectionSetUpBC(*gsection);
1215:   if (nroots >= 0) {
1216:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1217:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1218:     if (nroots > pEnd-pStart) {
1219:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1220:     }
1221:   }
1222:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1223:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1224:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
1225:     (*gsection)->atlasOff[p] = off;
1226:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
1227:   }
1228:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1229:   globalOff -= off;
1230:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1231:     (*gsection)->atlasOff[p] += globalOff;
1232:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
1233:   }
1234:   /* Put in negative offsets for ghost points */
1235:   if (nroots >= 0) {
1236:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1237:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1238:     if (nroots > pEnd-pStart) {
1239:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1240:     }
1241:   }
1242:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1243:   PetscFree(neg);
1244:   return(0);
1245: }