Actual source code: vsectionis.c

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

 5:  #include <petsc/private/isimpl.h>
 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:   if (s->sym) {
1553:     PetscViewerASCIIPushTab(viewer);
1554:     PetscSectionSymView(s->sym,viewer);
1555:     PetscViewerASCIIPopTab(viewer);
1556:   }
1557:   return(0);
1558: }

1562: /*@C
1563:   PetscSectionView - Views a PetscSection

1565:   Collective on PetscSection

1567:   Input Parameters:
1568: + s - the PetscSection object to view
1569: - v - the viewer

1571:   Level: developer

1573: .seealso PetscSectionCreate(), PetscSectionDestroy()
1574: @*/
1575: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1576: {
1577:   PetscBool      isascii;
1578:   PetscInt       f;

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

1602: /*@
1603:   PetscSectionReset - Frees all section data.

1605:   Not collective

1607:   Input Parameters:
1608: . s - the PetscSection

1610:   Level: developer

1612: .seealso: PetscSection, PetscSectionCreate()
1613: @*/
1614: PetscErrorCode PetscSectionReset(PetscSection s)
1615: {
1616:   PetscInt       f;

1620:   PetscFree(s->numFieldComponents);
1621:   for (f = 0; f < s->numFields; ++f) {
1622:     PetscSectionDestroy(&s->field[f]);
1623:     PetscFree(s->fieldNames[f]);
1624:   }
1625:   PetscFree(s->fieldNames);
1626:   PetscFree(s->field);
1627:   PetscSectionDestroy(&s->bc);
1628:   PetscFree(s->bcIndices);
1629:   PetscFree2(s->atlasDof, s->atlasOff);
1630:   PetscSectionDestroy(&s->clSection);
1631:   ISDestroy(&s->clPoints);
1632:   ISDestroy(&s->perm);
1633:   PetscFree(s->clPerm);
1634:   PetscFree(s->clInvPerm);
1635:   PetscSectionSymDestroy(&s->sym);

1637:   s->pStart    = -1;
1638:   s->pEnd      = -1;
1639:   s->maxDof    = 0;
1640:   s->setup     = PETSC_FALSE;
1641:   s->numFields = 0;
1642:   s->clObj     = NULL;
1643:   return(0);
1644: }

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

1651:   Not collective

1653:   Input Parameters:
1654: . s - the PetscSection

1656:   Level: developer

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

1661: .seealso: PetscSection, PetscSectionCreate()
1662: @*/
1663: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1664: {

1668:   if (!*s) return(0);
1670:   if (--((PetscObject)(*s))->refct > 0) {
1671:     *s = NULL;
1672:     return(0);
1673:   }
1674:   PetscSectionReset(*s);
1675:   PetscHeaderDestroy(s);
1676:   return(0);
1677: }

1681: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1682: {
1683:   const PetscInt p = point - s->pStart;

1686:   *values = &baseArray[s->atlasOff[p]];
1687:   return(0);
1688: }

1692: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1693: {
1694:   PetscInt       *array;
1695:   const PetscInt p           = point - s->pStart;
1696:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1697:   PetscInt       cDim        = 0;

1701:   PetscSectionGetConstraintDof(s, p, &cDim);
1702:   array = &baseArray[s->atlasOff[p]];
1703:   if (!cDim) {
1704:     if (orientation >= 0) {
1705:       const PetscInt dim = s->atlasDof[p];
1706:       PetscInt       i;

1708:       if (mode == INSERT_VALUES) {
1709:         for (i = 0; i < dim; ++i) array[i] = values[i];
1710:       } else {
1711:         for (i = 0; i < dim; ++i) array[i] += values[i];
1712:       }
1713:     } else {
1714:       PetscInt offset = 0;
1715:       PetscInt j      = -1, field, i;

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

1720:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1721:         offset += dim;
1722:       }
1723:     }
1724:   } else {
1725:     if (orientation >= 0) {
1726:       const PetscInt dim  = s->atlasDof[p];
1727:       PetscInt       cInd = 0, i;
1728:       const PetscInt *cDof;

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

1748:       PetscSectionGetConstraintIndices(s, point, &cDof);
1749:       for (field = 0; field < s->numFields; ++field) {
1750:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1751:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1752:         const PetscInt sDim = dim - tDim;
1753:         PetscInt       cInd = 0, i,k;

1755:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1756:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1757:           array[j] = values[k];
1758:         }
1759:         offset  += dim;
1760:         cOffset += dim - tDim;
1761:       }
1762:     }
1763:   }
1764:   return(0);
1765: }

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

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

1783:   Input Parameters:
1784: + s     - The PetscSection
1785: - point - The point

1787:   Output Parameter:
1788: . indices - The constrained dofs

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

1792:   Level: advanced

1794: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1795: @*/
1796: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1797: {

1801:   if (s->bc) {
1802:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1803:   } else *indices = NULL;
1804:   return(0);
1805: }

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

1812:   Input Parameters:
1813: + s     - The PetscSection
1814: . point - The point
1815: - indices - The constrained dofs

1817:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1819:   Level: advanced

1821: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1822: @*/
1823: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1824: {

1828:   if (s->bc) {
1829:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1830:   }
1831:   return(0);
1832: }

1836: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1837: {

1841:   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);
1842:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1843:   return(0);
1844: }

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

1853:   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);
1854:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1855:   return(0);
1856: }

1860: /*@
1861:   PetscSectionPermute - Reorder the section according to the input point permutation

1863:   Collective on PetscSection

1865:   Input Parameter:
1866: + section - The PetscSection object
1867: - perm - The point permutation, old point p becomes new point perm[p]

1869:   Output Parameter:
1870: . sectionNew - The permuted PetscSection

1872:   Level: intermediate

1874: .keywords: mesh
1875: .seealso: MatPermute()
1876: @*/
1877: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1878: {
1879:   PetscSection    s = section, sNew;
1880:   const PetscInt *perm;
1881:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1882:   PetscErrorCode  ierr;

1888:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1889:   PetscSectionGetNumFields(s, &numFields);
1890:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1891:   for (f = 0; f < numFields; ++f) {
1892:     const char *name;
1893:     PetscInt    numComp;

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

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

1924:     PetscSectionGetConstraintDof(s, p, &cdof);
1925:     if (cdof) {
1926:       PetscSectionGetConstraintIndices(s, p, &cind);
1927:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1928:     }
1929:     for (f = 0; f < numFields; ++f) {
1930:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1931:       if (cdof) {
1932:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1933:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1934:       }
1935:     }
1936:   }
1937:   ISRestoreIndices(permutation, &perm);
1938:   *sectionNew = sNew;
1939:   return(0);
1940: }

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

1947:   Collective

1949:   Input Parameters:
1950: + sf - The SF
1951: - rootSection - Section defined on root space

1953:   Output Parameters:
1954: + remoteOffsets - root offsets in leaf storage, or NULL
1955: - leafSection - Section defined on the leaf space

1957:   Level: intermediate

1959: .seealso: PetscSFCreate()
1960: @*/
1961: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1962: {
1963:   PetscSF        embedSF;
1964:   const PetscInt *ilocal, *indices;
1965:   IS             selected;
1966:   PetscInt       numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

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

2017: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2018: {
2019:   PetscSF         embedSF;
2020:   const PetscInt *indices;
2021:   IS              selected;
2022:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2023:   PetscErrorCode  ierr;

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

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

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

2054:   Output Parameters:
2055: - sectionSF - The new SF

2057:   Note: Either rootSection or remoteOffsets can be specified

2059:   Level: intermediate

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

2081:   PetscObjectGetComm((PetscObject)sf,&comm);
2082:   PetscSFCreate(comm, sectionSF);
2083:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2084:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2085:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2086:   if (numRoots < 0) return(0);
2087:   for (i = 0; i < numPoints; ++i) {
2088:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2089:     PetscInt dof;

2091:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2092:       PetscSectionGetDof(leafSection, localPoint, &dof);
2093:       numIndices += dof;
2094:     }
2095:   }
2096:   PetscMalloc1(numIndices, &localIndices);
2097:   PetscMalloc1(numIndices, &remoteIndices);
2098:   /* Create new index graph */
2099:   for (i = 0, ind = 0; i < numPoints; ++i) {
2100:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2101:     PetscInt rank       = remotePoints[i].rank;

2103:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2104:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2105:       PetscInt loff, dof, d;

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

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

2126:   Input Parameters:
2127: + section   - The PetscSection
2128: . obj       - A PetscObject which serves as the key for this index
2129: . clSection - Section giving the size of the closure of each point
2130: - clPoints  - IS giving the points in each closure

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

2134:   Level: intermediate

2136: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2137: @*/
2138: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2139: {

2143:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2144:   section->clObj     = obj;
2145:   PetscSectionDestroy(&section->clSection);
2146:   ISDestroy(&section->clPoints);
2147:   section->clSection = clSection;
2148:   section->clPoints  = clPoints;
2149:   return(0);
2150: }

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

2157:   Input Parameters:
2158: + section   - The PetscSection
2159: - obj       - A PetscObject which serves as the key for this index

2161:   Output Parameters:
2162: + clSection - Section giving the size of the closure of each point
2163: - clPoints  - IS giving the points in each closure

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

2167:   Level: intermediate

2169: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2170: @*/
2171: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2172: {
2174:   if (section->clObj == obj) {
2175:     if (clSection) *clSection = section->clSection;
2176:     if (clPoints)  *clPoints  = section->clPoints;
2177:   } else {
2178:     if (clSection) *clSection = NULL;
2179:     if (clPoints)  *clPoints  = NULL;
2180:   }
2181:   return(0);
2182: }

2186: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2187: {
2188:   PetscInt       i;

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

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

2217:   Input Parameters:
2218: + section - The PetscSection
2219: . obj     - A PetscObject which serves as the key for this index
2220: - perm    - Permutation of the cell dof closure

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

2225:   Level: intermediate

2227: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2228: @*/
2229: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2230: {
2231:   const PetscInt *clPerm = NULL;
2232:   PetscInt        clSize = 0;
2233:   PetscErrorCode  ierr;

2236:   if (perm) {
2237:     ISGetLocalSize(perm, &clSize);
2238:     ISGetIndices(perm, &clPerm);
2239:   }
2240:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2241:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2242:   return(0);
2243: }

2247: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2248: {
2250:   if (section->clObj == obj) {
2251:     if (size) *size = section->clSize;
2252:     if (perm) *perm = section->clPerm;
2253:   } else {
2254:     if (size) *size = 0;
2255:     if (perm) *perm = NULL;
2256:   }
2257:   return(0);
2258: }

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

2265:   Input Parameters:
2266: + section   - The PetscSection
2267: - obj       - A PetscObject which serves as the key for this index

2269:   Output Parameter:
2270: . perm - The dof closure permutation

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

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

2277:   Level: intermediate

2279: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2280: @*/
2281: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2282: {
2283:   const PetscInt *clPerm;
2284:   PetscInt        clSize;
2285:   PetscErrorCode  ierr;

2288:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2289:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2290:   return(0);
2291: }

2295: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2296: {
2298:   if (section->clObj == obj) {
2299:     if (size) *size = section->clSize;
2300:     if (perm) *perm = section->clInvPerm;
2301:   } else {
2302:     if (size) *size = 0;
2303:     if (perm) *perm = NULL;
2304:   }
2305:   return(0);
2306: }

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

2313:   Input Parameters:
2314: + section   - The PetscSection
2315: - obj       - A PetscObject which serves as the key for this index

2317:   Output Parameters:
2318: + size - The dof closure size
2319: - perm - The dof closure permutation

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

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

2326:   Level: intermediate

2328: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2329: @*/
2330: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2331: {
2332:   const PetscInt *clPerm;
2333:   PetscInt        clSize;
2334:   PetscErrorCode  ierr;

2337:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2338:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2339:   return(0);
2340: }

2344: /*@
2345:   PetscSectionGetField - Get the subsection associated with a single field

2347:   Input Parameters:
2348: + s     - The PetscSection
2349: - field - The field number

2351:   Output Parameter:
2352: . subs  - The subsection for the given field

2354:   Level: intermediate

2356: .seealso: PetscSectionSetNumFields()
2357: @*/
2358: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2359: {

2365:   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);
2366:   PetscObjectReference((PetscObject) s->field[field]);
2367:   *subs = s->field[field];
2368:   return(0);
2369: }

2371: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2372: PetscFunctionList PetscSectionSymList = NULL;

2376: /*@
2377:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2379:   Collective on MPI_Comm

2381:   Input Parameter:
2382: . comm - the MPI communicator

2384:   Output Parameter:
2385: . sym - pointer to the new set of symmetries

2387:   Level: developer

2389: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2390: @*/
2391: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2392: {

2397:   ISInitializePackage();
2398:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2399:   return(0);
2400: }

2404: /*@C
2405:   PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.

2407:   Collective on PetscSectionSym

2409:   Input Parameters:
2410: + sym    - The section symmetry object
2411: - method - The name of the section symmetry type

2413:   Level: developer

2415: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2416: @*/
2417: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2418: {
2419:   PetscErrorCode (*r)(PetscSectionSym);
2420:   PetscBool      match;

2425:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2426:   if (match) return(0);

2428:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2429:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2430:   if (sym->ops->destroy) {
2431:     (*sym->ops->destroy)(sym);
2432:     sym->ops->destroy = NULL;
2433:   }
2434:   (*r)(sym);
2435:   PetscObjectChangeTypeName((PetscObject)sym,method);
2436:   return(0);
2437: }


2442: /*@C
2443:   PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.

2445:   Not Collective

2447:   Input Parameter:
2448: . sym  - The section symmetry

2450:   Output Parameter:
2451: . type - The index set type name

2453:   Level: developer

2455: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2456: @*/
2457: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2458: {
2462:   *type = ((PetscObject)sym)->type_name;
2463:   return(0);
2464: }

2468: /*@C
2469:   PetscSectionSymRegister - Adds a new section symmetry implementation

2471:   Not Collective

2473:   Input Parameters:
2474: + name        - The name of a new user-defined creation routine
2475: - create_func - The creation routine itself

2477:   Notes:
2478:   PetscSectionSymRegister() may be called multiple times to add several user-defined vectors

2480:   Level: developer

2482: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2483: @*/
2484: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2485: {

2489:   ISInitializePackage();
2490:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2491:   return(0);
2492: }

2496: /*@
2497:    PetscSectionSymDestroy - Destroys a section symmetry.

2499:    Collective on PetscSectionSym

2501:    Input Parameters:
2502: .  sym - the section symmetry

2504:    Level: developer

2506: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2507: @*/
2508: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2509: {
2510:   SymWorkLink    link,next;

2514:   if (!*sym) return(0);
2516:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2517:   if ((*sym)->ops->destroy) {
2518:     (*(*sym)->ops->destroy)(*sym);
2519:   }
2520:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2521:   for (link=(*sym)->workin; link; link=next) {
2522:     next = link->next;
2523:     PetscFree2(link->perms,link->rots);
2524:     PetscFree(link);
2525:   }
2526:   (*sym)->workin = NULL;
2527:   PetscHeaderDestroy(sym);
2528:   return(0);
2529: }

2533: /*@C
2534:    PetscSectionSymView - Displays a section symmetry

2536:    Collective on PetscSectionSym

2538:    Input Parameters:
2539: +  sym - the index set
2540: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2542:    Level: developer

2544: .seealso: PetscViewerASCIIOpen()
2545: @*/
2546: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2547: {

2552:   if (!viewer) {
2553:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2554:   }
2557:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2558:   if (sym->ops->view) {
2559:     (*sym->ops->view)(sym,viewer);
2560:   }
2561:   return(0);
2562: }

2566: /*@
2567:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2569:   Collective on PetscSection

2571:   Input Parameters:
2572: + section - the section describing data layout
2573: - sym - the symmetry describing the affect of orientation on the access of the data

2575:   Level: developer

2577: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2578: @*/
2579: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2580: {

2587:   PetscObjectReference((PetscObject)sym);
2588:   PetscSectionSymDestroy(&(section->sym));
2589:   section->sym = sym;
2590:   return(0);
2591: }

2595: /*@
2596:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

2598:   Not collective

2600:   Input Parameters:
2601: . section - the section describing data layout

2603:   Output Parameters:
2604: . sym - the symmetry describing the affect of orientation on the access of the data

2606:   Level: developer

2608: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2609: @*/
2610: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2611: {
2614:   *sym = section->sym;
2615:   return(0);
2616: }

2620: /*@
2621:   PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section

2623:   Collective on PetscSection

2625:   Input Parameters:
2626: + section - the section describing data layout
2627: . field - the field number
2628: - sym - the symmetry describing the affect of orientation on the access of the data

2630:   Level: developer

2632: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2633: @*/
2634: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2635: {

2640:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
2641:   PetscSectionSetSym(section->field[field],sym);
2642:   return(0);
2643: }

2647: /*@
2648:   PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section

2650:   Collective on PetscSection

2652:   Input Parameters:
2653: + section - the section describing data layout
2654: - field - the field number

2656:   Output Parameters:
2657: . sym - the symmetry describing the affect of orientation on the access of the data

2659:   Level: developer

2661: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2662: @*/
2663: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2664: {
2667:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
2668:   *sym = section->field[field]->sym;
2669:   return(0);
2670: }

2674: /*@C
2675:   PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.

2677:   Not collective

2679:   Input Parameters:
2680: + section - the section
2681: . numPoints - the number of points
2682: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2683:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2684:     context, see DMPlexGetConeOrientation()).

2686:   Output Parameter:
2687: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2688: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2689:     identity).

2691:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2692: .vb
2693:      const PetscInt    **perms;
2694:      const PetscScalar **rots;
2695:      PetscInt            lOffset;

2697:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2698:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2699:        PetscInt           point = points[2*i], dof, sOffset;
2700:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2701:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2703:        PetscSectionGetDof(section,point,&dof);
2704:        PetscSectionGetOffset(section,point,&sOffset);

2706:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
2707:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
2708:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
2709:        lOffset += dof;
2710:      }
2711:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2712: .ve

2714:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2715: .vb
2716:      const PetscInt    **perms;
2717:      const PetscScalar **rots;
2718:      PetscInt            lOffset;

2720:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2721:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2722:        PetscInt           point = points[2*i], dof, sOffset;
2723:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2724:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2726:        PetscSectionGetDof(section,point,&dof);
2727:        PetscSectionGetOffset(section,point,&sOff);

2729:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2730:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
2731:        offset += dof;
2732:      }
2733:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2734: .ve

2736:   Level: developer

2738: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2739: @*/
2740: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2741: {
2742:   PetscSectionSym sym;
2743:   PetscErrorCode  ierr;

2748:   if (perms) *perms = NULL;
2749:   if (rots)  *rots  = NULL;
2750:   sym = section->sym;
2751:   if (sym && (perms || rots)) {
2752:     SymWorkLink link;

2754:     if (sym->workin) {
2755:       link        = sym->workin;
2756:       sym->workin = sym->workin->next;
2757:     } else {
2758:       PetscNewLog(sym,&link);
2759:     }
2760:     if (numPoints > link->numPoints) {
2761:       PetscFree2(link->perms,link->rots);
2762:       PetscMalloc2(numPoints,&link->perms,numPoints,&link->rots);
2763:       link->numPoints = numPoints;
2764:     }
2765:     link->next   = sym->workout;
2766:     sym->workout = link;
2767:     PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2768:     PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2769:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2770:     if (perms) *perms = link->perms;
2771:     if (rots)  *rots  = link->rots;
2772:   }
2773:   return(0);
2774: }

2778: /*@C
2779:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

2781:   Not collective

2783:   Input Parameters:
2784: + section - the section
2785: . numPoints - the number of points
2786: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2787:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2788:     context, see DMPlexGetConeOrientation()).

2790:   Output Parameter:
2791: + perms - The permutations for the given orientations: set to NULL at conclusion
2792: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2794:   Level: developer

2796: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2797: @*/
2798: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2799: {
2800:   PetscSectionSym sym;

2804:   sym = section->sym;
2805:   if (sym && (perms || rots)) {
2806:     SymWorkLink *p,link;

2808:     for (p=&sym->workout; (link=*p); p=&link->next) {
2809:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2810:         *p          = link->next;
2811:         link->next  = sym->workin;
2812:         sym->workin = link;
2813:         if (perms) *perms = NULL;
2814:         if (rots)  *rots  = NULL;
2815:         return(0);
2816:       }
2817:     }
2818:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2819:   }
2820:   return(0);
2821: }

2825: /*@C
2826:   PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.

2828:   Not collective

2830:   Input Parameters:
2831: + section - the section
2832: . field - the field of the section
2833: . numPoints - the number of points
2834: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2835:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2836:     context, see DMPlexGetConeOrientation()).

2838:   Output Parameter:
2839: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2840: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2841:     identity).

2843:   Level: developer

2845: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2846: @*/
2847: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2848: {

2853:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
2854:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2855:   return(0);
2856: }

2860: /*@C
2861:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

2863:   Not collective

2865:   Input Parameters:
2866: + section - the section
2867: . field - the field number
2868: . numPoints - the number of points
2869: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2870:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2871:     context, see DMPlexGetConeOrientation()).

2873:   Output Parameter:
2874: + perms - The permutations for the given orientations: set to NULL at conclusion
2875: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2877:   Level: developer

2879: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2880: @*/
2881: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2882: {

2887:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
2888:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
2889:   return(0);
2890: }