Actual source code: vsectionis.c

petsc-master 2016-07-29
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

  5: #include <petsc/private/isimpl.h>   /*I  "petscvec.h"   I*/
  6: #include <petscsf.h>
  7: #include <petscviewer.h>

  9: PetscClassId PETSC_SECTION_CLASSID;

 13: /*@
 14:   PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.

 16:   Collective on MPI_Comm

 18:   Input Parameters:
 19: + comm - the MPI communicator
 20: - s    - pointer to the section

 22:   Level: developer

 24:   Notes: Typical calling sequence
 25: $       PetscSectionCreate(MPI_Comm,PetscSection *);
 26: $       PetscSectionSetNumFields(PetscSection, numFields);
 27: $       PetscSectionSetChart(PetscSection,low,high);
 28: $       PetscSectionSetDof(PetscSection,point,numdof);
 29: $       PetscSectionSetUp(PetscSection);
 30: $       PetscSectionGetOffset(PetscSection,point,PetscInt *);
 31: $       PetscSectionDestroy(PetscSection);

 33:        The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
 34:        recommended they not be used in user codes unless you really gain something in their use.

 36: .seealso: PetscSection, PetscSectionDestroy()
 37: @*/
 38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
 39: {

 44:   ISInitializePackage();

 46:   PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);

 48:   (*s)->pStart             = -1;
 49:   (*s)->pEnd               = -1;
 50:   (*s)->perm               = NULL;
 51:   (*s)->maxDof             = 0;
 52:   (*s)->atlasDof           = NULL;
 53:   (*s)->atlasOff           = NULL;
 54:   (*s)->bc                 = NULL;
 55:   (*s)->bcIndices          = NULL;
 56:   (*s)->setup              = PETSC_FALSE;
 57:   (*s)->numFields          = 0;
 58:   (*s)->fieldNames         = NULL;
 59:   (*s)->field              = NULL;
 60:   (*s)->clObj              = NULL;
 61:   (*s)->clSection          = NULL;
 62:   (*s)->clPoints           = NULL;
 63:   (*s)->clSize             = 0;
 64:   (*s)->clPerm             = NULL;
 65:   (*s)->clInvPerm          = NULL;
 66:   return(0);
 67: }

 71: /*@
 72:   PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection

 74:   Collective on MPI_Comm

 76:   Input Parameter:
 77: . section - the PetscSection

 79:   Output Parameter:
 80: . newSection - the copy

 82:   Level: developer

 84: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
 85: @*/
 86: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 87: {
 88:   IS             perm;
 89:   PetscInt       numFields, f, pStart, pEnd, p;

 93:   PetscSectionGetNumFields(section, &numFields);
 94:   if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
 95:   for (f = 0; f < numFields; ++f) {
 96:     const char *name   = NULL;
 97:     PetscInt   numComp = 0;

 99:     PetscSectionGetFieldName(section, f, &name);
100:     PetscSectionSetFieldName(newSection, f, name);
101:     PetscSectionGetFieldComponents(section, f, &numComp);
102:     PetscSectionSetFieldComponents(newSection, f, numComp);
103:   }
104:   PetscSectionGetChart(section, &pStart, &pEnd);
105:   PetscSectionSetChart(newSection, pStart, pEnd);
106:   PetscSectionGetPermutation(section, &perm);
107:   PetscSectionSetPermutation(newSection, perm);
108:   for (p = pStart; p < pEnd; ++p) {
109:     PetscInt dof, cdof, fcdof = 0;

111:     PetscSectionGetDof(section, p, &dof);
112:     PetscSectionSetDof(newSection, p, dof);
113:     PetscSectionGetConstraintDof(section, p, &cdof);
114:     if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
115:     for (f = 0; f < numFields; ++f) {
116:       PetscSectionGetFieldDof(section, p, f, &dof);
117:       PetscSectionSetFieldDof(newSection, p, f, dof);
118:       if (cdof) {
119:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
120:         if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
121:       }
122:     }
123:   }
124:   PetscSectionSetUp(newSection);
125:   for (p = pStart; p < pEnd; ++p) {
126:     PetscInt       off, cdof, fcdof = 0;
127:     const PetscInt *cInd;

129:     /* Must set offsets in case they do not agree with the prefix sums */
130:     PetscSectionGetOffset(section, p, &off);
131:     PetscSectionSetOffset(newSection, p, off);
132:     PetscSectionGetConstraintDof(section, p, &cdof);
133:     if (cdof) {
134:       PetscSectionGetConstraintIndices(section, p, &cInd);
135:       PetscSectionSetConstraintIndices(newSection, p, cInd);
136:       for (f = 0; f < numFields; ++f) {
137:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
138:         if (fcdof) {
139:           PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
140:           PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
141:         }
142:       }
143:     }
144:   }
145:   return(0);
146: }

150: /*@
151:   PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection

153:   Collective on MPI_Comm

155:   Input Parameter:
156: . section - the PetscSection

158:   Output Parameter:
159: . newSection - the copy

161:   Level: developer

163: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
164: @*/
165: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
166: {

170:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
171:   PetscSectionCopy(section, *newSection);
172:   return(0);
173: }

177: /*@
178:   PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.

180:   Not collective

182:   Input Parameter:
183: . s - the PetscSection

185:   Output Parameter:
186: . numFields - the number of fields defined, or 0 if none were defined

188:   Level: intermediate

190: .seealso: PetscSectionSetNumFields()
191: @*/
192: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
193: {
196:   *numFields = s->numFields;
197:   return(0);
198: }

202: /*@
203:   PetscSectionSetNumFields - Sets the number of fields.

205:   Not collective

207:   Input Parameters:
208: + s - the PetscSection
209: - numFields - the number of fields

211:   Level: intermediate

213: .seealso: PetscSectionGetNumFields()
214: @*/
215: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
216: {
217:   PetscInt       f;

221:   if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
222:   PetscSectionReset(s);

224:   s->numFields = numFields;
225:   PetscMalloc1(s->numFields, &s->numFieldComponents);
226:   PetscMalloc1(s->numFields, &s->fieldNames);
227:   PetscMalloc1(s->numFields, &s->field);
228:   for (f = 0; f < s->numFields; ++f) {
229:     char name[64];

231:     s->numFieldComponents[f] = 1;

233:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
234:     PetscSNPrintf(name, 64, "Field_%D", f);
235:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
236:   }
237:   return(0);
238: }

242: /*@C
243:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

245:   Not Collective

247:   Input Parameters:
248: + s     - the PetscSection
249: - field - the field number

251:   Output Parameter:
252: . fieldName - the field name

254:   Level: developer

256: .seealso: PetscSectionSetFieldName()
257: @*/
258: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
259: {
262:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
263:   *fieldName = s->fieldNames[field];
264:   return(0);
265: }

269: /*@C
270:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

272:   Not Collective

274:   Input Parameters:
275: + s     - the PetscSection
276: . field - the field number
277: - fieldName - the field name

279:   Level: developer

281: .seealso: PetscSectionGetFieldName()
282: @*/
283: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
284: {

289:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
290:   PetscFree(s->fieldNames[field]);
291:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
292:   return(0);
293: }

297: /*@
298:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

300:   Not collective

302:   Input Parameters:
303: + s - the PetscSection
304: - field - the field number

306:   Output Parameter:
307: . numComp - the number of field components

309:   Level: intermediate

311: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
312: @*/
313: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
314: {
317:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
318:   *numComp = s->numFieldComponents[field];
319:   return(0);
320: }

324: /*@
325:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

327:   Not collective

329:   Input Parameters:
330: + s - the PetscSection
331: . field - the field number
332: - numComp - the number of field components

334:   Level: intermediate

336: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
337: @*/
338: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
339: {
341:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
342:   s->numFieldComponents[field] = numComp;
343:   return(0);
344: }

348: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
349: {

353:   if (!s->bc) {
354:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
355:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
356:   }
357:   return(0);
358: }

362: /*@
363:   PetscSectionGetChart - Returns the range [pStart, pEnd) in which points in the lie.

365:   Not collective

367:   Input Parameter:
368: . s - the PetscSection

370:   Output Parameters:
371: + pStart - the first point
372: - pEnd - one past the last point

374:   Level: intermediate

376: .seealso: PetscSectionSetChart(), PetscSectionCreate()
377: @*/
378: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
379: {
382:   if (pStart) *pStart = s->pStart;
383:   if (pEnd)   *pEnd   = s->pEnd;
384:   return(0);
385: }

389: /*@
390:   PetscSectionSetChart - Sets the range [pStart, pEnd) in which points in the lie.

392:   Not collective

394:   Input Parameters:
395: + s - the PetscSection
396: . pStart - the first point
397: - pEnd - one past the last point

399:   Level: intermediate

401: .seealso: PetscSectionGetChart(), PetscSectionCreate()
402: @*/
403: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
404: {
405:   PetscInt       f;

410:   /* Cannot Reset() because it destroys field information */
411:   s->setup = PETSC_FALSE;
412:   PetscSectionDestroy(&s->bc);
413:   PetscFree(s->bcIndices);
414:   PetscFree2(s->atlasDof, s->atlasOff);

416:   s->pStart = pStart;
417:   s->pEnd   = pEnd;
418:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
419:   PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
420:   for (f = 0; f < s->numFields; ++f) {
421:     PetscSectionSetChart(s->field[f], pStart, pEnd);
422:   }
423:   return(0);
424: }

428: /*@
429:   PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL

431:   Not collective

433:   Input Parameter:
434: . s - the PetscSection

436:   Output Parameters:
437: . perm - The permutation as an IS

439:   Level: intermediate

441: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
442: @*/
443: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
444: {
448:   return(0);
449: }

453: /*@
454:   PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)

456:   Not collective

458:   Input Parameters:
459: + s - the PetscSection
460: - perm - the permutation of points

462:   Level: intermediate

464: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
465: @*/
466: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
467: {

471:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
472:   if (s->perm != perm) {
473:     ISDestroy(&s->perm);
474:     s->perm = perm;
475:     PetscObjectReference((PetscObject) s->perm);
476:   }
477:   return(0);
478: }

482: /*@
483:   PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.

485:   Not collective

487:   Input Parameters:
488: + s - the PetscSection
489: - point - the point

491:   Output Parameter:
492: . numDof - the number of dof

494:   Level: intermediate

496: .seealso: PetscSectionSetDof(), PetscSectionCreate()
497: @*/
498: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
499: {
501: #if defined(PETSC_USE_DEBUG)
502:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
503: #endif
504:   *numDof = s->atlasDof[point - s->pStart];
505:   return(0);
506: }

510: /*@
511:   PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.

513:   Not collective

515:   Input Parameters:
516: + s - the PetscSection
517: . point - the point
518: - numDof - the number of dof

520:   Level: intermediate

522: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
523: @*/
524: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
525: {
527:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
528:   s->atlasDof[point - s->pStart] = numDof;
529:   return(0);
530: }

534: /*@
535:   PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.

537:   Not collective

539:   Input Parameters:
540: + s - the PetscSection
541: . point - the point
542: - numDof - the number of additional dof

544:   Level: intermediate

546: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
547: @*/
548: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
549: {
551:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
552:   s->atlasDof[point - s->pStart] += numDof;
553:   return(0);
554: }

558: /*@
559:   PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.

561:   Not collective

563:   Input Parameters:
564: + s - the PetscSection
565: . point - the point
566: - field - the field

568:   Output Parameter:
569: . numDof - the number of dof

571:   Level: intermediate

573: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
574: @*/
575: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
576: {

580:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
581:   PetscSectionGetDof(s->field[field], point, numDof);
582:   return(0);
583: }

587: /*@
588:   PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.

590:   Not collective

592:   Input Parameters:
593: + s - the PetscSection
594: . point - the point
595: . field - the field
596: - numDof - the number of dof

598:   Level: intermediate

600: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
601: @*/
602: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
603: {

607:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
608:   PetscSectionSetDof(s->field[field], point, numDof);
609:   return(0);
610: }

614: /*@
615:   PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.

617:   Not collective

619:   Input Parameters:
620: + s - the PetscSection
621: . point - the point
622: . field - the field
623: - numDof - the number of dof

625:   Level: intermediate

627: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
628: @*/
629: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
630: {

634:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
635:   PetscSectionAddDof(s->field[field], point, numDof);
636:   return(0);
637: }

641: /*@
642:   PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.

644:   Not collective

646:   Input Parameters:
647: + s - the PetscSection
648: - point - the point

650:   Output Parameter:
651: . numDof - the number of dof which are fixed by constraints

653:   Level: intermediate

655: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
656: @*/
657: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
658: {

662:   if (s->bc) {
663:     PetscSectionGetDof(s->bc, point, numDof);
664:   } else *numDof = 0;
665:   return(0);
666: }

670: /*@
671:   PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.

673:   Not collective

675:   Input Parameters:
676: + s - the PetscSection
677: . point - the point
678: - numDof - the number of dof which are fixed by constraints

680:   Level: intermediate

682: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
683: @*/
684: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
685: {

689:   if (numDof) {
690:     PetscSectionCheckConstraints_Static(s);
691:     PetscSectionSetDof(s->bc, point, numDof);
692:   }
693:   return(0);
694: }

698: /*@
699:   PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.

701:   Not collective

703:   Input Parameters:
704: + s - the PetscSection
705: . point - the point
706: - numDof - the number of additional dof which are fixed by constraints

708:   Level: intermediate

710: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
711: @*/
712: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
713: {

717:   if (numDof) {
718:     PetscSectionCheckConstraints_Static(s);
719:     PetscSectionAddDof(s->bc, point, numDof);
720:   }
721:   return(0);
722: }

726: /*@
727:   PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.

729:   Not collective

731:   Input Parameters:
732: + s - the PetscSection
733: . point - the point
734: - field - the field

736:   Output Parameter:
737: . numDof - the number of dof which are fixed by constraints

739:   Level: intermediate

741: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
742: @*/
743: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
744: {

748:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
749:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
750:   return(0);
751: }

755: /*@
756:   PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.

758:   Not collective

760:   Input Parameters:
761: + s - the PetscSection
762: . point - the point
763: . field - the field
764: - numDof - the number of dof which are fixed by constraints

766:   Level: intermediate

768: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
769: @*/
770: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
771: {

775:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
776:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
777:   return(0);
778: }

782: /*@
783:   PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.

785:   Not collective

787:   Input Parameters:
788: + s - the PetscSection
789: . point - the point
790: . field - the field
791: - numDof - the number of additional dof which are fixed by constraints

793:   Level: intermediate

795: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
796: @*/
797: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
798: {

802:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
803:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
804:   return(0);
805: }

809: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
810: {

814:   if (s->bc) {
815:     const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;

817:     PetscSectionSetUp(s->bc);
818:     PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
819:   }
820:   return(0);
821: }

825: /*@
826:   PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.

828:   Not collective

830:   Input Parameter:
831: . s - the PetscSection

833:   Level: intermediate

835: .seealso: PetscSectionCreate()
836: @*/
837: PetscErrorCode PetscSectionSetUp(PetscSection s)
838: {
839:   const PetscInt *pind   = NULL;
840:   PetscInt        offset = 0, p, f;
841:   PetscErrorCode  ierr;

844:   if (s->setup) return(0);
845:   s->setup = PETSC_TRUE;
846:   if (s->perm) {ISGetIndices(s->perm, &pind);}
847:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
848:     const PetscInt q = pind ? pind[p] : p;

850:     s->atlasOff[q] = offset;
851:     offset        += s->atlasDof[q];
852:     s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
853:   }
854:   PetscSectionSetUpBC(s);
855:   /* Assume that all fields have the same chart */
856:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
857:     const PetscInt q   = pind ? pind[p] : p;
858:     PetscInt       off = s->atlasOff[q];

860:     for (f = 0; f < s->numFields; ++f) {
861:       PetscSection sf = s->field[f];

863:       sf->atlasOff[q] = off;
864:       off += sf->atlasDof[q];
865:     }
866:   }
867:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
868:   for (f = 0; f < s->numFields; ++f) {
869:     PetscSectionSetUpBC(s->field[f]);
870:   }
871:   return(0);
872: }

876: /*@
877:   PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart

879:   Not collective

881:   Input Parameters:
882: . s - the PetscSection

884:   Output Parameter:
885: . maxDof - the maximum dof

887:   Level: intermediate

889: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
890: @*/
891: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
892: {
894:   *maxDof = s->maxDof;
895:   return(0);
896: }

900: /*@
901:   PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.

903:   Not collective

905:   Input Parameters:
906: + s - the PetscSection
907: - size - the allocated size

909:   Output Parameter:
910: . size - the size of an array which can hold all the dofs

912:   Level: intermediate

914: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
915: @*/
916: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
917: {
918:   PetscInt p, n = 0;

921:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
922:   *size = n;
923:   return(0);
924: }

928: /*@
929:   PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.

931:   Not collective

933:   Input Parameters:
934: + s - the PetscSection
935: - point - the point

937:   Output Parameter:
938: . size - the size of an array which can hold all unconstrained dofs

940:   Level: intermediate

942: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
943: @*/
944: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
945: {
946:   PetscInt p, n = 0;

949:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
950:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
951:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
952:   }
953:   *size = n;
954:   return(0);
955: }

959: /*@
960:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
961:   the local section and an SF describing the section point overlap.

963:   Input Parameters:
964:   + s - The PetscSection for the local field layout
965:   . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
966:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
967:   - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points

969:   Output Parameter:
970:   . gsection - The PetscSection for the global field layout

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

974:   Level: developer

976: .seealso: PetscSectionCreate()
977: @*/
978: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
979: {
980:   const PetscInt *pind = NULL;
981:   PetscInt       *recv = NULL, *neg = NULL;
982:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
983:   PetscErrorCode  ierr;

986:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
987:   PetscSectionGetChart(s, &pStart, &pEnd);
988:   PetscSectionSetChart(*gsection, pStart, pEnd);
989:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
990:   nlocal = nroots;              /* The local/leaf space matches global/root space */
991:   /* Must allocate for all points visible to SF, which may be more than this section */
992:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
993:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
994:     PetscMalloc2(nroots,&neg,nlocal,&recv);
995:     PetscMemzero(neg,nroots*sizeof(PetscInt));
996:   }
997:   /* Mark all local points with negative dof */
998:   for (p = pStart; p < pEnd; ++p) {
999:     PetscSectionGetDof(s, p, &dof);
1000:     PetscSectionSetDof(*gsection, p, dof);
1001:     PetscSectionGetConstraintDof(s, p, &cdof);
1002:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1003:     if (neg) neg[p] = -(dof+1);
1004:   }
1005:   PetscSectionSetUpBC(*gsection);
1006:   if (nroots >= 0) {
1007:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1008:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1009:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1010:     for (p = pStart; p < pEnd; ++p) {
1011:       if (recv[p] < 0) {
1012:         (*gsection)->atlasDof[p-pStart] = recv[p];
1013:         PetscSectionGetDof(s, p, &dof);
1014:         if (-(recv[p]+1) != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %d for point %d is not the unconstrained %d", -(recv[p]+1), p, dof);
1015:       }
1016:     }
1017:   }
1018:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1019:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1020:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1021:     const PetscInt q = pind ? pind[p] : p;

1023:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1024:     (*gsection)->atlasOff[q] = off;
1025:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1026:   }
1027:   if (!localOffsets) {
1028:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1029:     globalOff -= off;
1030:   }
1031:   for (p = pStart, off = 0; p < pEnd; ++p) {
1032:     (*gsection)->atlasOff[p-pStart] += globalOff;
1033:     if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
1034:   }
1035:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1036:   /* Put in negative offsets for ghost points */
1037:   if (nroots >= 0) {
1038:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1039:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1040:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1041:     for (p = pStart; p < pEnd; ++p) {
1042:       if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
1043:     }
1044:   }
1045:   PetscFree2(neg,recv);
1046:   PetscSectionViewFromOptions(*gsection,NULL,"-global_section_view");
1047:   return(0);
1048: }

1052: /*@
1053:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1054:   the local section and an SF describing the section point overlap.

1056:   Input Parameters:
1057:   + s - The PetscSection for the local field layout
1058:   . sf - The SF describing parallel layout of the section points
1059:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1060:   . numExcludes - The number of exclusion ranges
1061:   - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs

1063:   Output Parameter:
1064:   . gsection - The PetscSection for the global field layout

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

1068:   Level: developer

1070: .seealso: PetscSectionCreate()
1071: @*/
1072: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1073: {
1074:   const PetscInt *pind = NULL;
1075:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1076:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1077:   PetscErrorCode  ierr;

1080:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1081:   PetscSectionGetChart(s, &pStart, &pEnd);
1082:   PetscSectionSetChart(*gsection, pStart, pEnd);
1083:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1084:   if (nroots >= 0) {
1085:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1086:     PetscCalloc1(nroots, &neg);
1087:     if (nroots > pEnd-pStart) {
1088:       PetscCalloc1(nroots, &tmpOff);
1089:     } else {
1090:       tmpOff = &(*gsection)->atlasDof[-pStart];
1091:     }
1092:   }
1093:   /* Mark ghost points with negative dof */
1094:   for (p = pStart; p < pEnd; ++p) {
1095:     for (e = 0; e < numExcludes; ++e) {
1096:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1097:         PetscSectionSetDof(*gsection, p, 0);
1098:         break;
1099:       }
1100:     }
1101:     if (e < numExcludes) continue;
1102:     PetscSectionGetDof(s, p, &dof);
1103:     PetscSectionSetDof(*gsection, p, dof);
1104:     PetscSectionGetConstraintDof(s, p, &cdof);
1105:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1106:     if (neg) neg[p] = -(dof+1);
1107:   }
1108:   PetscSectionSetUpBC(*gsection);
1109:   if (nroots >= 0) {
1110:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1111:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1112:     if (nroots > pEnd - pStart) {
1113:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1114:     }
1115:   }
1116:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1117:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1118:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1119:     const PetscInt q = pind ? pind[p] : p;

1121:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1122:     (*gsection)->atlasOff[q] = off;
1123:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1124:   }
1125:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1126:   globalOff -= off;
1127:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1128:     (*gsection)->atlasOff[p] += globalOff;
1129:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1130:   }
1131:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1132:   /* Put in negative offsets for ghost points */
1133:   if (nroots >= 0) {
1134:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1135:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1136:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1137:     if (nroots > pEnd - pStart) {
1138:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1139:     }
1140:   }
1141:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1142:   PetscFree(neg);
1143:   return(0);
1144: }

1148: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1149: {
1150:   PetscInt       pStart, pEnd, p, localSize = 0;

1154:   PetscSectionGetChart(s, &pStart, &pEnd);
1155:   for (p = pStart; p < pEnd; ++p) {
1156:     PetscInt dof;

1158:     PetscSectionGetDof(s, p, &dof);
1159:     if (dof > 0) ++localSize;
1160:   }
1161:   PetscLayoutCreate(comm, layout);
1162:   PetscLayoutSetLocalSize(*layout, localSize);
1163:   PetscLayoutSetBlockSize(*layout, 1);
1164:   PetscLayoutSetUp(*layout);
1165:   return(0);
1166: }

1170: /*@
1171:   PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.

1173:   Input Parameters:
1174: + comm - The MPI_Comm
1175: - s    - The PetscSection

1177:   Output Parameter:
1178: . layout - The layout for the section

1180:   Level: developer

1182: .seealso: PetscSectionCreate()
1183: @*/
1184: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1185: {
1186:   PetscInt       pStart, pEnd, p, localSize = 0;

1190:   PetscSectionGetChart(s, &pStart, &pEnd);
1191:   for (p = pStart; p < pEnd; ++p) {
1192:     PetscInt dof,cdof;

1194:     PetscSectionGetDof(s, p, &dof);
1195:     PetscSectionGetConstraintDof(s, p, &cdof);
1196:     if (dof-cdof > 0) localSize += dof-cdof;
1197:   }
1198:   PetscLayoutCreate(comm, layout);
1199:   PetscLayoutSetLocalSize(*layout, localSize);
1200:   PetscLayoutSetBlockSize(*layout, 1);
1201:   PetscLayoutSetUp(*layout);
1202:   return(0);
1203: }

1207: /*@
1208:   PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.

1210:   Not collective

1212:   Input Parameters:
1213: + s - the PetscSection
1214: - point - the point

1216:   Output Parameter:
1217: . offset - the offset

1219:   Level: intermediate

1221: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1222: @*/
1223: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1224: {
1226: #if defined(PETSC_USE_DEBUG)
1227:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1228: #endif
1229:   *offset = s->atlasOff[point - s->pStart];
1230:   return(0);
1231: }

1235: /*@
1236:   PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.

1238:   Not collective

1240:   Input Parameters:
1241: + s - the PetscSection
1242: . point - the point
1243: - offset - the offset

1245:   Note: The user usually does not call this function, but uses PetscSectionSetUp()

1247:   Level: intermediate

1249: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1250: @*/
1251: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1252: {
1254:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1255:   s->atlasOff[point - s->pStart] = offset;
1256:   return(0);
1257: }

1261: /*@
1262:   PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.

1264:   Not collective

1266:   Input Parameters:
1267: + s - the PetscSection
1268: . point - the point
1269: - field - the field

1271:   Output Parameter:
1272: . offset - the offset

1274:   Level: intermediate

1276: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1277: @*/
1278: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1279: {

1283:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1284:   PetscSectionGetOffset(s->field[field], point, offset);
1285:   return(0);
1286: }

1290: /*@
1291:   PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.

1293:   Not collective

1295:   Input Parameters:
1296: + s - the PetscSection
1297: . point - the point
1298: . field - the field
1299: - offset - the offset

1301:   Note: The user usually does not call this function, but uses PetscSectionSetUp()

1303:   Level: intermediate

1305: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1306: @*/
1307: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1308: {

1312:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1313:   PetscSectionSetOffset(s->field[field], point, offset);
1314:   return(0);
1315: }

1319: /* This gives the offset on a point of the field, ignoring constraints */
1320: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1321: {
1322:   PetscInt       off, foff;

1326:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1327:   PetscSectionGetOffset(s, point, &off);
1328:   PetscSectionGetOffset(s->field[field], point, &foff);
1329:   *offset = foff - off;
1330:   return(0);
1331: }

1335: /*@
1336:   PetscSectionGetOffsetRange - Return the full range of offsets [start, end)

1338:   Not collective

1340:   Input Parameter:
1341: . s - the PetscSection

1343:   Output Parameters:
1344: + start - the minimum offset
1345: - end   - one more than the maximum offset

1347:   Level: intermediate

1349: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1350: @*/
1351: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1352: {
1353:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1357:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1358:   PetscSectionGetChart(s, &pStart, &pEnd);
1359:   for (p = 0; p < pEnd-pStart; ++p) {
1360:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1362:     if (off >= 0) {
1363:       os = PetscMin(os, off);
1364:       oe = PetscMax(oe, off+dof);
1365:     }
1366:   }
1367:   if (start) *start = os;
1368:   if (end)   *end   = oe;
1369:   return(0);
1370: }

1374: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, PetscInt fields[], PetscSection *subs)
1375: {
1376:   PetscInt       nF, f, pStart, pEnd, p, maxCdof = 0;

1380:   if (!numFields) return(0);
1381:   PetscSectionGetNumFields(s, &nF);
1382:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1383:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1384:   PetscSectionSetNumFields(*subs, numFields);
1385:   for (f = 0; f < numFields; ++f) {
1386:     const char *name   = NULL;
1387:     PetscInt   numComp = 0;

1389:     PetscSectionGetFieldName(s, fields[f], &name);
1390:     PetscSectionSetFieldName(*subs, f, name);
1391:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1392:     PetscSectionSetFieldComponents(*subs, f, numComp);
1393:   }
1394:   PetscSectionGetChart(s, &pStart, &pEnd);
1395:   PetscSectionSetChart(*subs, pStart, pEnd);
1396:   for (p = pStart; p < pEnd; ++p) {
1397:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1399:     for (f = 0; f < numFields; ++f) {
1400:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1401:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1402:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1403:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1404:       dof  += fdof;
1405:       cdof += cfdof;
1406:     }
1407:     PetscSectionSetDof(*subs, p, dof);
1408:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1409:     maxCdof = PetscMax(cdof, maxCdof);
1410:   }
1411:   PetscSectionSetUp(*subs);
1412:   if (maxCdof) {
1413:     PetscInt *indices;

1415:     PetscMalloc1(maxCdof, &indices);
1416:     for (p = pStart; p < pEnd; ++p) {
1417:       PetscInt cdof;

1419:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1420:       if (cdof) {
1421:         const PetscInt *oldIndices = NULL;
1422:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1424:         for (f = 0; f < numFields; ++f) {
1425:           PetscInt oldFoff = 0, oldf;

1427:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1428:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1429:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1430:           /* This can be sped up if we assume sorted fields */
1431:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1432:             PetscInt oldfdof = 0;
1433:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1434:             oldFoff += oldfdof;
1435:           }
1436:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1437:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1438:           numConst += cfdof;
1439:           fOff     += fdof;
1440:         }
1441:         PetscSectionSetConstraintIndices(*subs, p, indices);
1442:       }
1443:     }
1444:     PetscFree(indices);
1445:   }
1446:   return(0);
1447: }

1451: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1452: {
1453:   const PetscInt *points = NULL, *indices = NULL;
1454:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1458:   PetscSectionGetNumFields(s, &numFields);
1459:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1460:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1461:   for (f = 0; f < numFields; ++f) {
1462:     const char *name   = NULL;
1463:     PetscInt   numComp = 0;

1465:     PetscSectionGetFieldName(s, f, &name);
1466:     PetscSectionSetFieldName(*subs, f, name);
1467:     PetscSectionGetFieldComponents(s, f, &numComp);
1468:     PetscSectionSetFieldComponents(*subs, f, numComp);
1469:   }
1470:   /* For right now, we do not try to squeeze the subchart */
1471:   if (subpointMap) {
1472:     ISGetSize(subpointMap, &numSubpoints);
1473:     ISGetIndices(subpointMap, &points);
1474:   }
1475:   PetscSectionGetChart(s, &pStart, &pEnd);
1476:   PetscSectionSetChart(*subs, 0, numSubpoints);
1477:   for (p = pStart; p < pEnd; ++p) {
1478:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1480:     PetscFindInt(p, numSubpoints, points, &subp);
1481:     if (subp < 0) continue;
1482:     for (f = 0; f < numFields; ++f) {
1483:       PetscSectionGetFieldDof(s, p, f, &fdof);
1484:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1485:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1486:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1487:     }
1488:     PetscSectionGetDof(s, p, &dof);
1489:     PetscSectionSetDof(*subs, subp, dof);
1490:     PetscSectionGetConstraintDof(s, p, &cdof);
1491:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1492:   }
1493:   PetscSectionSetUp(*subs);
1494:   /* Change offsets to original offsets */
1495:   for (p = pStart; p < pEnd; ++p) {
1496:     PetscInt off, foff = 0;

1498:     PetscFindInt(p, numSubpoints, points, &subp);
1499:     if (subp < 0) continue;
1500:     for (f = 0; f < numFields; ++f) {
1501:       PetscSectionGetFieldOffset(s, p, f, &foff);
1502:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1503:     }
1504:     PetscSectionGetOffset(s, p, &off);
1505:     PetscSectionSetOffset(*subs, subp, off);
1506:   }
1507:   /* Copy constraint indices */
1508:   for (subp = 0; subp < numSubpoints; ++subp) {
1509:     PetscInt cdof;

1511:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1512:     if (cdof) {
1513:       for (f = 0; f < numFields; ++f) {
1514:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1515:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1516:       }
1517:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1518:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1519:     }
1520:   }
1521:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1522:   return(0);
1523: }

1527: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1528: {
1529:   PetscInt       p;
1530:   PetscMPIInt    rank;

1534:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1535:   PetscViewerASCIIPushSynchronized(viewer);
1536:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1537:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1538:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1539:       PetscInt b;

1541:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1542:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1543:         PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1544:       }
1545:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1546:     } else {
1547:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1548:     }
1549:   }
1550:   PetscViewerFlush(viewer);
1551:   PetscViewerASCIIPopSynchronized(viewer);
1552:   return(0);
1553: }

1557: /*@C
1558:   PetscSectionView - Views a PetscSection

1560:   Collective on PetscSection

1562:   Input Parameters:
1563: + s - the PetscSection object to view
1564: - v - the viewer

1566:   Level: developer

1568: .seealso PetscSectionCreate(), PetscSectionDestroy()
1569: @*/
1570: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1571: {
1572:   PetscBool      isascii;
1573:   PetscInt       f;

1577:   if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
1579:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1580:   if (isascii) {
1581:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1582:     if (s->numFields) {
1583:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1584:       for (f = 0; f < s->numFields; ++f) {
1585:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1586:         PetscSectionView_ASCII(s->field[f], viewer);
1587:       }
1588:     } else {
1589:       PetscSectionView_ASCII(s, viewer);
1590:     }
1591:   }
1592:   return(0);
1593: }

1597: /*@
1598:   PetscSectionReset - Frees all section data.

1600:   Not collective

1602:   Input Parameters:
1603: . s - the PetscSection

1605:   Level: developer

1607: .seealso: PetscSection, PetscSectionCreate()
1608: @*/
1609: PetscErrorCode PetscSectionReset(PetscSection s)
1610: {
1611:   PetscInt       f;

1615:   PetscFree(s->numFieldComponents);
1616:   for (f = 0; f < s->numFields; ++f) {
1617:     PetscSectionDestroy(&s->field[f]);
1618:     PetscFree(s->fieldNames[f]);
1619:   }
1620:   PetscFree(s->fieldNames);
1621:   PetscFree(s->field);
1622:   PetscSectionDestroy(&s->bc);
1623:   PetscFree(s->bcIndices);
1624:   PetscFree2(s->atlasDof, s->atlasOff);
1625:   PetscSectionDestroy(&s->clSection);
1626:   ISDestroy(&s->clPoints);
1627:   ISDestroy(&s->perm);
1628:   PetscFree(s->clPerm);
1629:   PetscFree(s->clInvPerm);

1631:   s->pStart    = -1;
1632:   s->pEnd      = -1;
1633:   s->maxDof    = 0;
1634:   s->setup     = PETSC_FALSE;
1635:   s->numFields = 0;
1636:   s->clObj     = NULL;
1637:   return(0);
1638: }

1642: /*@
1643:   PetscSectionDestroy - Frees a section object and frees its range if that exists.

1645:   Not collective

1647:   Input Parameters:
1648: . s - the PetscSection

1650:   Level: developer

1652:     The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1653:     recommended they not be used in user codes unless you really gain something in their use.

1655: .seealso: PetscSection, PetscSectionCreate()
1656: @*/
1657: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1658: {

1662:   if (!*s) return(0);
1664:   if (--((PetscObject)(*s))->refct > 0) {
1665:     *s = NULL;
1666:     return(0);
1667:   }
1668:   PetscSectionReset(*s);
1669:   PetscHeaderDestroy(s);
1670:   return(0);
1671: }

1675: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1676: {
1677:   const PetscInt p = point - s->pStart;

1680:   *values = &baseArray[s->atlasOff[p]];
1681:   return(0);
1682: }

1686: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1687: {
1688:   PetscInt       *array;
1689:   const PetscInt p           = point - s->pStart;
1690:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1691:   PetscInt       cDim        = 0;

1695:   PetscSectionGetConstraintDof(s, p, &cDim);
1696:   array = &baseArray[s->atlasOff[p]];
1697:   if (!cDim) {
1698:     if (orientation >= 0) {
1699:       const PetscInt dim = s->atlasDof[p];
1700:       PetscInt       i;

1702:       if (mode == INSERT_VALUES) {
1703:         for (i = 0; i < dim; ++i) array[i] = values[i];
1704:       } else {
1705:         for (i = 0; i < dim; ++i) array[i] += values[i];
1706:       }
1707:     } else {
1708:       PetscInt offset = 0;
1709:       PetscInt j      = -1, field, i;

1711:       for (field = 0; field < s->numFields; ++field) {
1712:         const PetscInt dim = s->field[field]->atlasDof[p];

1714:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1715:         offset += dim;
1716:       }
1717:     }
1718:   } else {
1719:     if (orientation >= 0) {
1720:       const PetscInt dim  = s->atlasDof[p];
1721:       PetscInt       cInd = 0, i;
1722:       const PetscInt *cDof;

1724:       PetscSectionGetConstraintIndices(s, point, &cDof);
1725:       if (mode == INSERT_VALUES) {
1726:         for (i = 0; i < dim; ++i) {
1727:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1728:           array[i] = values[i];
1729:         }
1730:       } else {
1731:         for (i = 0; i < dim; ++i) {
1732:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1733:           array[i] += values[i];
1734:         }
1735:       }
1736:     } else {
1737:       const PetscInt *cDof;
1738:       PetscInt       offset  = 0;
1739:       PetscInt       cOffset = 0;
1740:       PetscInt       j       = 0, field;

1742:       PetscSectionGetConstraintIndices(s, point, &cDof);
1743:       for (field = 0; field < s->numFields; ++field) {
1744:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1745:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1746:         const PetscInt sDim = dim - tDim;
1747:         PetscInt       cInd = 0, i,k;

1749:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1750:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1751:           array[j] = values[k];
1752:         }
1753:         offset  += dim;
1754:         cOffset += dim - tDim;
1755:       }
1756:     }
1757:   }
1758:   return(0);
1759: }

1763: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1764: {
1768:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1769:   return(0);
1770: }

1774: /*@C
1775:   PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained

1777:   Input Parameters:
1778: + s     - The PetscSection
1779: - point - The point

1781:   Output Parameter:
1782: . indices - The constrained dofs

1784:   Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()

1786:   Level: advanced

1788: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1789: @*/
1790: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1791: {

1795:   if (s->bc) {
1796:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1797:   } else *indices = NULL;
1798:   return(0);
1799: }

1803: /*@C
1804:   PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained

1806:   Input Parameters:
1807: + s     - The PetscSection
1808: . point - The point
1809: - indices - The constrained dofs

1811:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1813:   Level: advanced

1815: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1816: @*/
1817: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1818: {

1822:   if (s->bc) {
1823:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1824:   }
1825:   return(0);
1826: }

1830: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1831: {

1835:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1836:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1837:   return(0);
1838: }

1842: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
1843: {

1847:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1848:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1849:   return(0);
1850: }

1854: /*@
1855:   PetscSectionPermute - Reorder the section according to the input point permutation

1857:   Collective on PetscSection

1859:   Input Parameter:
1860: + section - The PetscSection object
1861: - perm - The point permutation, old point p becomes new point perm[p]

1863:   Output Parameter:
1864: . sectionNew - The permuted PetscSection

1866:   Level: intermediate

1868: .keywords: mesh
1869: .seealso: MatPermute()
1870: @*/
1871: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1872: {
1873:   PetscSection    s = section, sNew;
1874:   const PetscInt *perm;
1875:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1876:   PetscErrorCode  ierr;

1882:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1883:   PetscSectionGetNumFields(s, &numFields);
1884:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1885:   for (f = 0; f < numFields; ++f) {
1886:     const char *name;
1887:     PetscInt    numComp;

1889:     PetscSectionGetFieldName(s, f, &name);
1890:     PetscSectionSetFieldName(sNew, f, name);
1891:     PetscSectionGetFieldComponents(s, f, &numComp);
1892:     PetscSectionSetFieldComponents(sNew, f, numComp);
1893:   }
1894:   ISGetLocalSize(permutation, &numPoints);
1895:   ISGetIndices(permutation, &perm);
1896:   PetscSectionGetChart(s, &pStart, &pEnd);
1897:   PetscSectionSetChart(sNew, pStart, pEnd);
1898:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1899:   for (p = pStart; p < pEnd; ++p) {
1900:     PetscInt dof, cdof;

1902:     PetscSectionGetDof(s, p, &dof);
1903:     PetscSectionSetDof(sNew, perm[p], dof);
1904:     PetscSectionGetConstraintDof(s, p, &cdof);
1905:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1906:     for (f = 0; f < numFields; ++f) {
1907:       PetscSectionGetFieldDof(s, p, f, &dof);
1908:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1909:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1910:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1911:     }
1912:   }
1913:   PetscSectionSetUp(sNew);
1914:   for (p = pStart; p < pEnd; ++p) {
1915:     const PetscInt *cind;
1916:     PetscInt        cdof;

1918:     PetscSectionGetConstraintDof(s, p, &cdof);
1919:     if (cdof) {
1920:       PetscSectionGetConstraintIndices(s, p, &cind);
1921:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1922:     }
1923:     for (f = 0; f < numFields; ++f) {
1924:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1925:       if (cdof) {
1926:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1927:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1928:       }
1929:     }
1930:   }
1931:   ISRestoreIndices(permutation, &perm);
1932:   *sectionNew = sNew;
1933:   return(0);
1934: }

1938: /*@C
1939:   PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF

1941:   Collective

1943:   Input Parameters:
1944: + sf - The SF
1945: - rootSection - Section defined on root space

1947:   Output Parameters:
1948: + remoteOffsets - root offsets in leaf storage, or NULL
1949: - leafSection - Section defined on the leaf space

1951:   Level: intermediate

1953: .seealso: PetscSFCreate()
1954: @*/
1955: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1956: {
1957:   PetscSF        embedSF;
1958:   const PetscInt *ilocal, *indices;
1959:   IS             selected;
1960:   PetscInt       numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

1964:   PetscSectionGetNumFields(rootSection, &numFields);
1965:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1966:   for (f = 0; f < numFields; ++f) {
1967:     PetscInt numComp = 0;
1968:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
1969:     PetscSectionSetFieldComponents(leafSection, f, numComp);
1970:   }
1971:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1972:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
1973:   rpEnd = PetscMin(rpEnd,nroots);
1974:   rpEnd = PetscMax(rpStart,rpEnd);
1975:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1976:   ISGetIndices(selected, &indices);
1977:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1978:   ISRestoreIndices(selected, &indices);
1979:   ISDestroy(&selected);
1980:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1981:   if (nleaves && ilocal) {
1982:     for (i = 0; i < nleaves; ++i) {
1983:       lpStart = PetscMin(lpStart, ilocal[i]);
1984:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
1985:     }
1986:     ++lpEnd;
1987:   } else {
1988:     lpStart = 0;
1989:     lpEnd   = nleaves;
1990:   }
1991:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
1992:   /* Could fuse these at the cost of a copy and extra allocation */
1993:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1994:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1995:   for (f = 0; f < numFields; ++f) {
1996:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1997:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1998:   }
1999:   if (remoteOffsets) {
2000:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2001:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2002:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2003:   }
2004:   PetscSFDestroy(&embedSF);
2005:   PetscSectionSetUp(leafSection);
2006:   return(0);
2007: }

2011: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2012: {
2013:   PetscSF         embedSF;
2014:   const PetscInt *indices;
2015:   IS              selected;
2016:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2017:   PetscErrorCode  ierr;

2020:   *remoteOffsets = NULL;
2021:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2022:   if (numRoots < 0) return(0);
2023:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2024:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2025:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2026:   ISGetIndices(selected, &indices);
2027:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2028:   ISRestoreIndices(selected, &indices);
2029:   ISDestroy(&selected);
2030:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2031:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2032:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2033:   PetscSFDestroy(&embedSF);
2034:   return(0);
2035: }

2039: /*@C
2040:   PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points

2042:   Input Parameters:
2043: + sf - The SF
2044: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2045: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2046: - leafSection - Data layout of local points for incoming data  (this is the distributed section)

2048:   Output Parameters:
2049: - sectionSF - The new SF

2051:   Note: Either rootSection or remoteOffsets can be specified

2053:   Level: intermediate

2055: .seealso: PetscSFCreate()
2056: @*/
2057: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2058: {
2059:   MPI_Comm          comm;
2060:   const PetscInt    *localPoints;
2061:   const PetscSFNode *remotePoints;
2062:   PetscInt          lpStart, lpEnd;
2063:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2064:   PetscInt          *localIndices;
2065:   PetscSFNode       *remoteIndices;
2066:   PetscInt          i, ind;
2067:   PetscErrorCode    ierr;

2075:   PetscObjectGetComm((PetscObject)sf,&comm);
2076:   PetscSFCreate(comm, sectionSF);
2077:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2078:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2079:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2080:   if (numRoots < 0) return(0);
2081:   for (i = 0; i < numPoints; ++i) {
2082:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2083:     PetscInt dof;

2085:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2086:       PetscSectionGetDof(leafSection, localPoint, &dof);
2087:       numIndices += dof;
2088:     }
2089:   }
2090:   PetscMalloc1(numIndices, &localIndices);
2091:   PetscMalloc1(numIndices, &remoteIndices);
2092:   /* Create new index graph */
2093:   for (i = 0, ind = 0; i < numPoints; ++i) {
2094:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2095:     PetscInt rank       = remotePoints[i].rank;

2097:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2098:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2099:       PetscInt loff, dof, d;

2101:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2102:       PetscSectionGetDof(leafSection, localPoint, &dof);
2103:       for (d = 0; d < dof; ++d, ++ind) {
2104:         localIndices[ind]        = loff+d;
2105:         remoteIndices[ind].rank  = rank;
2106:         remoteIndices[ind].index = remoteOffset+d;
2107:       }
2108:     }
2109:   }
2110:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2111:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2112:   return(0);
2113: }

2117: /*@
2118:   PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section

2120:   Input Parameters:
2121: + section   - The PetscSection
2122: . obj       - A PetscObject which serves as the key for this index
2123: . clSection - Section giving the size of the closure of each point
2124: - clPoints  - IS giving the points in each closure

2126:   Note: We compress out closure points with no dofs in this section

2128:   Level: intermediate

2130: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2131: @*/
2132: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2133: {

2137:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2138:   section->clObj     = obj;
2139:   PetscSectionDestroy(&section->clSection);
2140:   ISDestroy(&section->clPoints);
2141:   section->clSection = clSection;
2142:   section->clPoints  = clPoints;
2143:   return(0);
2144: }

2148: /*@
2149:   PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section

2151:   Input Parameters:
2152: + section   - The PetscSection
2153: - obj       - A PetscObject which serves as the key for this index

2155:   Output Parameters:
2156: + clSection - Section giving the size of the closure of each point
2157: - clPoints  - IS giving the points in each closure

2159:   Note: We compress out closure points with no dofs in this section

2161:   Level: intermediate

2163: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2164: @*/
2165: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2166: {
2168:   if (section->clObj == obj) {
2169:     if (clSection) *clSection = section->clSection;
2170:     if (clPoints)  *clPoints  = section->clPoints;
2171:   } else {
2172:     if (clSection) *clSection = NULL;
2173:     if (clPoints)  *clPoints  = NULL;
2174:   }
2175:   return(0);
2176: }

2180: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2181: {
2182:   PetscInt       i;

2186:   if (section->clObj != obj) {
2187:     PetscSectionDestroy(&section->clSection);
2188:     ISDestroy(&section->clPoints);
2189:   }
2190:   section->clObj  = obj;
2191:   PetscFree(section->clPerm);
2192:   PetscFree(section->clInvPerm);
2193:   section->clSize = clSize;
2194:   if (mode == PETSC_COPY_VALUES) {
2195:     PetscMalloc1(clSize, &section->clPerm);
2196:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2197:     PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2198:   } else if (mode == PETSC_OWN_POINTER) {
2199:     section->clPerm = clPerm;
2200:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2201:   PetscMalloc1(clSize, &section->clInvPerm);
2202:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2203:   return(0);
2204: }

2208: /*@
2209:   PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.

2211:   Input Parameters:
2212: + section - The PetscSection
2213: . obj     - A PetscObject which serves as the key for this index
2214: - perm    - Permutation of the cell dof closure

2216:   Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2217:   other points (like faces).

2219:   Level: intermediate

2221: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2222: @*/
2223: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2224: {
2225:   const PetscInt *clPerm = NULL;
2226:   PetscInt        clSize = 0;
2227:   PetscErrorCode  ierr;

2230:   if (perm) {
2231:     ISGetLocalSize(perm, &clSize);
2232:     ISGetIndices(perm, &clPerm);
2233:   }
2234:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2235:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2236:   return(0);
2237: }

2241: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2242: {
2244:   if (section->clObj == obj) {
2245:     if (size) *size = section->clSize;
2246:     if (perm) *perm = section->clPerm;
2247:   } else {
2248:     if (size) *size = 0;
2249:     if (perm) *perm = NULL;
2250:   }
2251:   return(0);
2252: }

2256: /*@
2257:   PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.

2259:   Input Parameters:
2260: + section   - The PetscSection
2261: - obj       - A PetscObject which serves as the key for this index

2263:   Output Parameter:
2264: . perm - The dof closure permutation

2266:   Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2267:   other points (like faces).

2269:   The user must destroy the IS that is returned.

2271:   Level: intermediate

2273: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2274: @*/
2275: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2276: {
2277:   const PetscInt *clPerm;
2278:   PetscInt        clSize;
2279:   PetscErrorCode  ierr;

2282:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2283:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2284:   return(0);
2285: }

2289: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2290: {
2292:   if (section->clObj == obj) {
2293:     if (size) *size = section->clSize;
2294:     if (perm) *perm = section->clInvPerm;
2295:   } else {
2296:     if (size) *size = 0;
2297:     if (perm) *perm = NULL;
2298:   }
2299:   return(0);
2300: }

2304: /*@
2305:   PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.

2307:   Input Parameters:
2308: + section   - The PetscSection
2309: - obj       - A PetscObject which serves as the key for this index

2311:   Output Parameters:
2312: + size - The dof closure size
2313: - perm - The dof closure permutation

2315:   Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2316:   other points (like faces).

2318:   The user must destroy the IS that is returned.

2320:   Level: intermediate

2322: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2323: @*/
2324: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2325: {
2326:   const PetscInt *clPerm;
2327:   PetscInt        clSize;
2328:   PetscErrorCode  ierr;

2331:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2332:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2333:   return(0);
2334: }

2338: /*@
2339:   PetscSectionGetField - Get the subsection associated with a single field

2341:   Input Parameters:
2342: + s     - The PetscSection
2343: - field - The field number

2345:   Output Parameter:
2346: . subs  - The subsection for the given field

2348:   Level: intermediate

2350: .seealso: PetscSectionSetNumFields()
2351: @*/
2352: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2353: {

2359:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
2360:   PetscObjectReference((PetscObject) s->field[field]);
2361:   *subs = s->field[field];
2362:   return(0);
2363: }