Actual source code: vsectionis.c

petsc-master 2017-01-16
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;

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

 14:   Collective on MPI_Comm

 16:   Input Parameters:
 17: + comm - the MPI communicator
 18: - s    - pointer to the section

 20:   Level: developer

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

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

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

 42:   ISInitializePackage();

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

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

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

 70:   Collective on MPI_Comm

 72:   Input Parameter:
 73: . section - the PetscSection

 75:   Output Parameter:
 76: . newSection - the copy

 78:   Level: developer

 80: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
 81: @*/
 82: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 83: {
 84:   IS             perm;
 85:   PetscInt       numFields, f, pStart, pEnd, p;

 89:   PetscSectionGetNumFields(section, &numFields);
 90:   if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
 91:   for (f = 0; f < numFields; ++f) {
 92:     const char *name   = NULL;
 93:     PetscInt   numComp = 0;

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

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

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

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

147:   Collective on MPI_Comm

149:   Input Parameter:
150: . section - the PetscSection

152:   Output Parameter:
153: . newSection - the copy

155:   Level: developer

157: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
158: @*/
159: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
160: {

164:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
165:   PetscSectionCopy(section, *newSection);
166:   return(0);
167: }

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

172:   Not collective

174:   Input Parameter:
175: . s - the PetscSection

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

180:   Level: intermediate

182: .seealso: PetscSectionSetNumFields()
183: @*/
184: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
185: {
188:   *numFields = s->numFields;
189:   return(0);
190: }

192: /*@
193:   PetscSectionSetNumFields - Sets the number of fields.

195:   Not collective

197:   Input Parameters:
198: + s - the PetscSection
199: - numFields - the number of fields

201:   Level: intermediate

203: .seealso: PetscSectionGetNumFields()
204: @*/
205: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
206: {
207:   PetscInt       f;

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

214:   s->numFields = numFields;
215:   PetscMalloc1(s->numFields, &s->numFieldComponents);
216:   PetscMalloc1(s->numFields, &s->fieldNames);
217:   PetscMalloc1(s->numFields, &s->field);
218:   for (f = 0; f < s->numFields; ++f) {
219:     char name[64];

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

223:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
224:     PetscSNPrintf(name, 64, "Field_%D", f);
225:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
226:   }
227:   return(0);
228: }

230: /*@C
231:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

233:   Not Collective

235:   Input Parameters:
236: + s     - the PetscSection
237: - field - the field number

239:   Output Parameter:
240: . fieldName - the field name

242:   Level: developer

244: .seealso: PetscSectionSetFieldName()
245: @*/
246: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
247: {
250:   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);
251:   *fieldName = s->fieldNames[field];
252:   return(0);
253: }

255: /*@C
256:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

258:   Not Collective

260:   Input Parameters:
261: + s     - the PetscSection
262: . field - the field number
263: - fieldName - the field name

265:   Level: developer

267: .seealso: PetscSectionGetFieldName()
268: @*/
269: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
270: {

275:   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);
276:   PetscFree(s->fieldNames[field]);
277:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
278:   return(0);
279: }

281: /*@
282:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

284:   Not collective

286:   Input Parameters:
287: + s - the PetscSection
288: - field - the field number

290:   Output Parameter:
291: . numComp - the number of field components

293:   Level: intermediate

295: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
296: @*/
297: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
298: {
301:   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);
302:   *numComp = s->numFieldComponents[field];
303:   return(0);
304: }

306: /*@
307:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

309:   Not collective

311:   Input Parameters:
312: + s - the PetscSection
313: . field - the field number
314: - numComp - the number of field components

316:   Level: intermediate

318: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
319: @*/
320: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
321: {
323:   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);
324:   s->numFieldComponents[field] = numComp;
325:   return(0);
326: }

328: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
329: {

333:   if (!s->bc) {
334:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
335:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
336:   }
337:   return(0);
338: }

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

343:   Not collective

345:   Input Parameter:
346: . s - the PetscSection

348:   Output Parameters:
349: + pStart - the first point
350: - pEnd - one past the last point

352:   Level: intermediate

354: .seealso: PetscSectionSetChart(), PetscSectionCreate()
355: @*/
356: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
357: {
360:   if (pStart) *pStart = s->pStart;
361:   if (pEnd)   *pEnd   = s->pEnd;
362:   return(0);
363: }

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

368:   Not collective

370:   Input Parameters:
371: + s - the PetscSection
372: . pStart - the first point
373: - pEnd - one past the last point

375:   Level: intermediate

377: .seealso: PetscSectionGetChart(), PetscSectionCreate()
378: @*/
379: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
380: {
381:   PetscInt       f;

386:   /* Cannot Reset() because it destroys field information */
387:   s->setup = PETSC_FALSE;
388:   PetscSectionDestroy(&s->bc);
389:   PetscFree(s->bcIndices);
390:   PetscFree2(s->atlasDof, s->atlasOff);

392:   s->pStart = pStart;
393:   s->pEnd   = pEnd;
394:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
395:   PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
396:   for (f = 0; f < s->numFields; ++f) {
397:     PetscSectionSetChart(s->field[f], pStart, pEnd);
398:   }
399:   return(0);
400: }

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

405:   Not collective

407:   Input Parameter:
408: . s - the PetscSection

410:   Output Parameters:
411: . perm - The permutation as an IS

413:   Level: intermediate

415: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
416: @*/
417: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
418: {
422:   return(0);
423: }

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

428:   Not collective

430:   Input Parameters:
431: + s - the PetscSection
432: - perm - the permutation of points

434:   Level: intermediate

436: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
437: @*/
438: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
439: {

443:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
444:   if (s->perm != perm) {
445:     ISDestroy(&s->perm);
446:     s->perm = perm;
447:     PetscObjectReference((PetscObject) s->perm);
448:   }
449:   return(0);
450: }

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

455:   Not collective

457:   Input Parameters:
458: + s - the PetscSection
459: - point - the point

461:   Output Parameter:
462: . numDof - the number of dof

464:   Level: intermediate

466: .seealso: PetscSectionSetDof(), PetscSectionCreate()
467: @*/
468: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
469: {
471: #if defined(PETSC_USE_DEBUG)
472:   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);
473: #endif
474:   *numDof = s->atlasDof[point - s->pStart];
475:   return(0);
476: }

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

481:   Not collective

483:   Input Parameters:
484: + s - the PetscSection
485: . point - the point
486: - numDof - the number of dof

488:   Level: intermediate

490: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
491: @*/
492: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
493: {
495:   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);
496:   s->atlasDof[point - s->pStart] = numDof;
497:   return(0);
498: }

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

503:   Not collective

505:   Input Parameters:
506: + s - the PetscSection
507: . point - the point
508: - numDof - the number of additional dof

510:   Level: intermediate

512: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
513: @*/
514: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
515: {
517:   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);
518:   s->atlasDof[point - s->pStart] += numDof;
519:   return(0);
520: }

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

525:   Not collective

527:   Input Parameters:
528: + s - the PetscSection
529: . point - the point
530: - field - the field

532:   Output Parameter:
533: . numDof - the number of dof

535:   Level: intermediate

537: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
538: @*/
539: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
540: {

544:   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);
545:   PetscSectionGetDof(s->field[field], point, numDof);
546:   return(0);
547: }

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

552:   Not collective

554:   Input Parameters:
555: + s - the PetscSection
556: . point - the point
557: . field - the field
558: - numDof - the number of dof

560:   Level: intermediate

562: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
563: @*/
564: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
565: {

569:   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);
570:   PetscSectionSetDof(s->field[field], point, numDof);
571:   return(0);
572: }

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

577:   Not collective

579:   Input Parameters:
580: + s - the PetscSection
581: . point - the point
582: . field - the field
583: - numDof - the number of dof

585:   Level: intermediate

587: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
588: @*/
589: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
590: {

594:   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);
595:   PetscSectionAddDof(s->field[field], point, numDof);
596:   return(0);
597: }

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

602:   Not collective

604:   Input Parameters:
605: + s - the PetscSection
606: - point - the point

608:   Output Parameter:
609: . numDof - the number of dof which are fixed by constraints

611:   Level: intermediate

613: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
614: @*/
615: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
616: {

620:   if (s->bc) {
621:     PetscSectionGetDof(s->bc, point, numDof);
622:   } else *numDof = 0;
623:   return(0);
624: }

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

629:   Not collective

631:   Input Parameters:
632: + s - the PetscSection
633: . point - the point
634: - numDof - the number of dof which are fixed by constraints

636:   Level: intermediate

638: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
639: @*/
640: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
641: {

645:   if (numDof) {
646:     PetscSectionCheckConstraints_Static(s);
647:     PetscSectionSetDof(s->bc, point, numDof);
648:   }
649:   return(0);
650: }

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

655:   Not collective

657:   Input Parameters:
658: + s - the PetscSection
659: . point - the point
660: - numDof - the number of additional dof which are fixed by constraints

662:   Level: intermediate

664: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
665: @*/
666: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
667: {

671:   if (numDof) {
672:     PetscSectionCheckConstraints_Static(s);
673:     PetscSectionAddDof(s->bc, point, numDof);
674:   }
675:   return(0);
676: }

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

681:   Not collective

683:   Input Parameters:
684: + s - the PetscSection
685: . point - the point
686: - field - the field

688:   Output Parameter:
689: . numDof - the number of dof which are fixed by constraints

691:   Level: intermediate

693: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
694: @*/
695: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
696: {

700:   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);
701:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
702:   return(0);
703: }

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

708:   Not collective

710:   Input Parameters:
711: + s - the PetscSection
712: . point - the point
713: . field - the field
714: - numDof - the number of dof which are fixed by constraints

716:   Level: intermediate

718: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
719: @*/
720: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
721: {

725:   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);
726:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
727:   return(0);
728: }

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

733:   Not collective

735:   Input Parameters:
736: + s - the PetscSection
737: . point - the point
738: . field - the field
739: - numDof - the number of additional dof which are fixed by constraints

741:   Level: intermediate

743: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
744: @*/
745: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
746: {

750:   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);
751:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
752:   return(0);
753: }

755: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
756: {

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

763:     PetscSectionSetUp(s->bc);
764:     PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
765:   }
766:   return(0);
767: }

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

772:   Not collective

774:   Input Parameter:
775: . s - the PetscSection

777:   Level: intermediate

779: .seealso: PetscSectionCreate()
780: @*/
781: PetscErrorCode PetscSectionSetUp(PetscSection s)
782: {
783:   const PetscInt *pind   = NULL;
784:   PetscInt        offset = 0, p, f;
785:   PetscErrorCode  ierr;

788:   if (s->setup) return(0);
789:   s->setup = PETSC_TRUE;
790:   if (s->perm) {ISGetIndices(s->perm, &pind);}
791:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
792:     const PetscInt q = pind ? pind[p] : p;

794:     s->atlasOff[q] = offset;
795:     offset        += s->atlasDof[q];
796:     s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
797:   }
798:   PetscSectionSetUpBC(s);
799:   /* Assume that all fields have the same chart */
800:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
801:     const PetscInt q   = pind ? pind[p] : p;
802:     PetscInt       off = s->atlasOff[q];

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

807:       sf->atlasOff[q] = off;
808:       off += sf->atlasDof[q];
809:     }
810:   }
811:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
812:   for (f = 0; f < s->numFields; ++f) {
813:     PetscSectionSetUpBC(s->field[f]);
814:   }
815:   return(0);
816: }

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

821:   Not collective

823:   Input Parameters:
824: . s - the PetscSection

826:   Output Parameter:
827: . maxDof - the maximum dof

829:   Level: intermediate

831: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
832: @*/
833: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
834: {
836:   *maxDof = s->maxDof;
837:   return(0);
838: }

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

843:   Not collective

845:   Input Parameters:
846: + s - the PetscSection
847: - size - the allocated size

849:   Output Parameter:
850: . size - the size of an array which can hold all the dofs

852:   Level: intermediate

854: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
855: @*/
856: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
857: {
858:   PetscInt p, n = 0;

861:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
862:   *size = n;
863:   return(0);
864: }

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

869:   Not collective

871:   Input Parameters:
872: + s - the PetscSection
873: - point - the point

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

878:   Level: intermediate

880: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
881: @*/
882: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
883: {
884:   PetscInt p, n = 0;

887:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
888:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
889:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
890:   }
891:   *size = n;
892:   return(0);
893: }

895: /*@
896:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
897:   the local section and an SF describing the section point overlap.

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

905:   Output Parameter:
906:   . gsection - The PetscSection for the global field layout

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

910:   Level: developer

912: .seealso: PetscSectionCreate()
913: @*/
914: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
915: {
916:   const PetscInt *pind = NULL;
917:   PetscInt       *recv = NULL, *neg = NULL;
918:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
919:   PetscErrorCode  ierr;

922:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
923:   PetscSectionGetChart(s, &pStart, &pEnd);
924:   PetscSectionSetChart(*gsection, pStart, pEnd);
925:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
926:   nlocal = nroots;              /* The local/leaf space matches global/root space */
927:   /* Must allocate for all points visible to SF, which may be more than this section */
928:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
929:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
930:     PetscMalloc2(nroots,&neg,nlocal,&recv);
931:     PetscMemzero(neg,nroots*sizeof(PetscInt));
932:   }
933:   /* Mark all local points with negative dof */
934:   for (p = pStart; p < pEnd; ++p) {
935:     PetscSectionGetDof(s, p, &dof);
936:     PetscSectionSetDof(*gsection, p, dof);
937:     PetscSectionGetConstraintDof(s, p, &cdof);
938:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
939:     if (neg) neg[p] = -(dof+1);
940:   }
941:   PetscSectionSetUpBC(*gsection);
942:   if (nroots >= 0) {
943:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
944:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
945:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
946:     for (p = pStart; p < pEnd; ++p) {
947:       if (recv[p] < 0) {
948:         (*gsection)->atlasDof[p-pStart] = recv[p];
949:         PetscSectionGetDof(s, p, &dof);
950:         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);
951:       }
952:     }
953:   }
954:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
955:   if (s->perm) {ISGetIndices(s->perm, &pind);}
956:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
957:     const PetscInt q = pind ? pind[p] : p;

959:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
960:     (*gsection)->atlasOff[q] = off;
961:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
962:   }
963:   if (!localOffsets) {
964:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
965:     globalOff -= off;
966:   }
967:   for (p = pStart, off = 0; p < pEnd; ++p) {
968:     (*gsection)->atlasOff[p-pStart] += globalOff;
969:     if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
970:   }
971:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
972:   /* Put in negative offsets for ghost points */
973:   if (nroots >= 0) {
974:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
975:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
976:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
977:     for (p = pStart; p < pEnd; ++p) {
978:       if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
979:     }
980:   }
981:   PetscFree2(neg,recv);
982:   PetscSectionViewFromOptions(*gsection,NULL,"-global_section_view");
983:   return(0);
984: }

986: /*@
987:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
988:   the local section and an SF describing the section point overlap.

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

997:   Output Parameter:
998:   . gsection - The PetscSection for the global field layout

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

1002:   Level: developer

1004: .seealso: PetscSectionCreate()
1005: @*/
1006: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1007: {
1008:   const PetscInt *pind = NULL;
1009:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1010:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1011:   PetscErrorCode  ierr;

1014:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1015:   PetscSectionGetChart(s, &pStart, &pEnd);
1016:   PetscSectionSetChart(*gsection, pStart, pEnd);
1017:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1018:   if (nroots >= 0) {
1019:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1020:     PetscCalloc1(nroots, &neg);
1021:     if (nroots > pEnd-pStart) {
1022:       PetscCalloc1(nroots, &tmpOff);
1023:     } else {
1024:       tmpOff = &(*gsection)->atlasDof[-pStart];
1025:     }
1026:   }
1027:   /* Mark ghost points with negative dof */
1028:   for (p = pStart; p < pEnd; ++p) {
1029:     for (e = 0; e < numExcludes; ++e) {
1030:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1031:         PetscSectionSetDof(*gsection, p, 0);
1032:         break;
1033:       }
1034:     }
1035:     if (e < numExcludes) continue;
1036:     PetscSectionGetDof(s, p, &dof);
1037:     PetscSectionSetDof(*gsection, p, dof);
1038:     PetscSectionGetConstraintDof(s, p, &cdof);
1039:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1040:     if (neg) neg[p] = -(dof+1);
1041:   }
1042:   PetscSectionSetUpBC(*gsection);
1043:   if (nroots >= 0) {
1044:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1045:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1046:     if (nroots > pEnd - pStart) {
1047:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1048:     }
1049:   }
1050:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1051:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1052:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1053:     const PetscInt q = pind ? pind[p] : p;

1055:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1056:     (*gsection)->atlasOff[q] = off;
1057:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1058:   }
1059:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1060:   globalOff -= off;
1061:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1062:     (*gsection)->atlasOff[p] += globalOff;
1063:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1064:   }
1065:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1066:   /* Put in negative offsets for ghost points */
1067:   if (nroots >= 0) {
1068:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1069:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1070:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1071:     if (nroots > pEnd - pStart) {
1072:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1073:     }
1074:   }
1075:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1076:   PetscFree(neg);
1077:   return(0);
1078: }

1080: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1081: {
1082:   PetscInt       pStart, pEnd, p, localSize = 0;

1086:   PetscSectionGetChart(s, &pStart, &pEnd);
1087:   for (p = pStart; p < pEnd; ++p) {
1088:     PetscInt dof;

1090:     PetscSectionGetDof(s, p, &dof);
1091:     if (dof > 0) ++localSize;
1092:   }
1093:   PetscLayoutCreate(comm, layout);
1094:   PetscLayoutSetLocalSize(*layout, localSize);
1095:   PetscLayoutSetBlockSize(*layout, 1);
1096:   PetscLayoutSetUp(*layout);
1097:   return(0);
1098: }

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

1103:   Input Parameters:
1104: + comm - The MPI_Comm
1105: - s    - The PetscSection

1107:   Output Parameter:
1108: . layout - The layout for the section

1110:   Level: developer

1112: .seealso: PetscSectionCreate()
1113: @*/
1114: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1115: {
1116:   PetscInt       pStart, pEnd, p, localSize = 0;

1120:   PetscSectionGetChart(s, &pStart, &pEnd);
1121:   for (p = pStart; p < pEnd; ++p) {
1122:     PetscInt dof,cdof;

1124:     PetscSectionGetDof(s, p, &dof);
1125:     PetscSectionGetConstraintDof(s, p, &cdof);
1126:     if (dof-cdof > 0) localSize += dof-cdof;
1127:   }
1128:   PetscLayoutCreate(comm, layout);
1129:   PetscLayoutSetLocalSize(*layout, localSize);
1130:   PetscLayoutSetBlockSize(*layout, 1);
1131:   PetscLayoutSetUp(*layout);
1132:   return(0);
1133: }

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

1138:   Not collective

1140:   Input Parameters:
1141: + s - the PetscSection
1142: - point - the point

1144:   Output Parameter:
1145: . offset - the offset

1147:   Level: intermediate

1149: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1150: @*/
1151: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1152: {
1154: #if defined(PETSC_USE_DEBUG)
1155:   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);
1156: #endif
1157:   *offset = s->atlasOff[point - s->pStart];
1158:   return(0);
1159: }

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

1164:   Not collective

1166:   Input Parameters:
1167: + s - the PetscSection
1168: . point - the point
1169: - offset - the offset

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

1173:   Level: intermediate

1175: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1176: @*/
1177: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1178: {
1180:   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);
1181:   s->atlasOff[point - s->pStart] = offset;
1182:   return(0);
1183: }

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

1188:   Not collective

1190:   Input Parameters:
1191: + s - the PetscSection
1192: . point - the point
1193: - field - the field

1195:   Output Parameter:
1196: . offset - the offset

1198:   Level: intermediate

1200: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1201: @*/
1202: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1203: {

1207:   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);
1208:   PetscSectionGetOffset(s->field[field], point, offset);
1209:   return(0);
1210: }

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

1215:   Not collective

1217:   Input Parameters:
1218: + s - the PetscSection
1219: . point - the point
1220: . field - the field
1221: - offset - the offset

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

1225:   Level: intermediate

1227: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1228: @*/
1229: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1230: {

1234:   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);
1235:   PetscSectionSetOffset(s->field[field], point, offset);
1236:   return(0);
1237: }

1239: /* This gives the offset on a point of the field, ignoring constraints */
1240: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1241: {
1242:   PetscInt       off, foff;

1246:   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);
1247:   PetscSectionGetOffset(s, point, &off);
1248:   PetscSectionGetOffset(s->field[field], point, &foff);
1249:   *offset = foff - off;
1250:   return(0);
1251: }

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

1256:   Not collective

1258:   Input Parameter:
1259: . s - the PetscSection

1261:   Output Parameters:
1262: + start - the minimum offset
1263: - end   - one more than the maximum offset

1265:   Level: intermediate

1267: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1268: @*/
1269: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1270: {
1271:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1275:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1276:   PetscSectionGetChart(s, &pStart, &pEnd);
1277:   for (p = 0; p < pEnd-pStart; ++p) {
1278:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1280:     if (off >= 0) {
1281:       os = PetscMin(os, off);
1282:       oe = PetscMax(oe, off+dof);
1283:     }
1284:   }
1285:   if (start) *start = os;
1286:   if (end)   *end   = oe;
1287:   return(0);
1288: }

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

1296:   if (!numFields) return(0);
1297:   PetscSectionGetNumFields(s, &nF);
1298:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1299:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1300:   PetscSectionSetNumFields(*subs, numFields);
1301:   for (f = 0; f < numFields; ++f) {
1302:     const char *name   = NULL;
1303:     PetscInt   numComp = 0;

1305:     PetscSectionGetFieldName(s, fields[f], &name);
1306:     PetscSectionSetFieldName(*subs, f, name);
1307:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1308:     PetscSectionSetFieldComponents(*subs, f, numComp);
1309:   }
1310:   PetscSectionGetChart(s, &pStart, &pEnd);
1311:   PetscSectionSetChart(*subs, pStart, pEnd);
1312:   for (p = pStart; p < pEnd; ++p) {
1313:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1315:     for (f = 0; f < numFields; ++f) {
1316:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1317:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1318:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1319:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1320:       dof  += fdof;
1321:       cdof += cfdof;
1322:     }
1323:     PetscSectionSetDof(*subs, p, dof);
1324:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1325:     maxCdof = PetscMax(cdof, maxCdof);
1326:   }
1327:   PetscSectionSetUp(*subs);
1328:   if (maxCdof) {
1329:     PetscInt *indices;

1331:     PetscMalloc1(maxCdof, &indices);
1332:     for (p = pStart; p < pEnd; ++p) {
1333:       PetscInt cdof;

1335:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1336:       if (cdof) {
1337:         const PetscInt *oldIndices = NULL;
1338:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1343:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1344:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1345:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1346:           /* This can be sped up if we assume sorted fields */
1347:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1348:             PetscInt oldfdof = 0;
1349:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1350:             oldFoff += oldfdof;
1351:           }
1352:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1353:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1354:           numConst += cfdof;
1355:           fOff     += fdof;
1356:         }
1357:         PetscSectionSetConstraintIndices(*subs, p, indices);
1358:       }
1359:     }
1360:     PetscFree(indices);
1361:   }
1362:   return(0);
1363: }

1365: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1366: {
1367:   const PetscInt *points = NULL, *indices = NULL;
1368:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1372:   PetscSectionGetNumFields(s, &numFields);
1373:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1374:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1375:   for (f = 0; f < numFields; ++f) {
1376:     const char *name   = NULL;
1377:     PetscInt   numComp = 0;

1379:     PetscSectionGetFieldName(s, f, &name);
1380:     PetscSectionSetFieldName(*subs, f, name);
1381:     PetscSectionGetFieldComponents(s, f, &numComp);
1382:     PetscSectionSetFieldComponents(*subs, f, numComp);
1383:   }
1384:   /* For right now, we do not try to squeeze the subchart */
1385:   if (subpointMap) {
1386:     ISGetSize(subpointMap, &numSubpoints);
1387:     ISGetIndices(subpointMap, &points);
1388:   }
1389:   PetscSectionGetChart(s, &pStart, &pEnd);
1390:   PetscSectionSetChart(*subs, 0, numSubpoints);
1391:   for (p = pStart; p < pEnd; ++p) {
1392:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1394:     PetscFindInt(p, numSubpoints, points, &subp);
1395:     if (subp < 0) continue;
1396:     for (f = 0; f < numFields; ++f) {
1397:       PetscSectionGetFieldDof(s, p, f, &fdof);
1398:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1399:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1400:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1401:     }
1402:     PetscSectionGetDof(s, p, &dof);
1403:     PetscSectionSetDof(*subs, subp, dof);
1404:     PetscSectionGetConstraintDof(s, p, &cdof);
1405:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1406:   }
1407:   PetscSectionSetUp(*subs);
1408:   /* Change offsets to original offsets */
1409:   for (p = pStart; p < pEnd; ++p) {
1410:     PetscInt off, foff = 0;

1412:     PetscFindInt(p, numSubpoints, points, &subp);
1413:     if (subp < 0) continue;
1414:     for (f = 0; f < numFields; ++f) {
1415:       PetscSectionGetFieldOffset(s, p, f, &foff);
1416:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1417:     }
1418:     PetscSectionGetOffset(s, p, &off);
1419:     PetscSectionSetOffset(*subs, subp, off);
1420:   }
1421:   /* Copy constraint indices */
1422:   for (subp = 0; subp < numSubpoints; ++subp) {
1423:     PetscInt cdof;

1425:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1426:     if (cdof) {
1427:       for (f = 0; f < numFields; ++f) {
1428:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1429:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1430:       }
1431:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1432:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1433:     }
1434:   }
1435:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1436:   return(0);
1437: }

1439: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1440: {
1441:   PetscInt       p;
1442:   PetscMPIInt    rank;

1446:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1447:   PetscViewerASCIIPushSynchronized(viewer);
1448:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1449:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1450:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1451:       PetscInt b;

1453:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1454:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1455:         PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1456:       }
1457:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1458:     } else {
1459:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1460:     }
1461:   }
1462:   PetscViewerFlush(viewer);
1463:   PetscViewerASCIIPopSynchronized(viewer);
1464:   if (s->sym) {
1465:     PetscViewerASCIIPushTab(viewer);
1466:     PetscSectionSymView(s->sym,viewer);
1467:     PetscViewerASCIIPopTab(viewer);
1468:   }
1469:   return(0);
1470: }

1472: /*@C
1473:   PetscSectionView - Views a PetscSection

1475:   Collective on PetscSection

1477:   Input Parameters:
1478: + s - the PetscSection object to view
1479: - v - the viewer

1481:   Level: developer

1483: .seealso PetscSectionCreate(), PetscSectionDestroy()
1484: @*/
1485: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1486: {
1487:   PetscBool      isascii;
1488:   PetscInt       f;

1493:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1495:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1496:   if (isascii) {
1497:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1498:     if (s->numFields) {
1499:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1500:       for (f = 0; f < s->numFields; ++f) {
1501:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1502:         PetscSectionView_ASCII(s->field[f], viewer);
1503:       }
1504:     } else {
1505:       PetscSectionView_ASCII(s, viewer);
1506:     }
1507:   }
1508:   return(0);
1509: }

1511: /*@
1512:   PetscSectionReset - Frees all section data.

1514:   Not collective

1516:   Input Parameters:
1517: . s - the PetscSection

1519:   Level: developer

1521: .seealso: PetscSection, PetscSectionCreate()
1522: @*/
1523: PetscErrorCode PetscSectionReset(PetscSection s)
1524: {
1525:   PetscInt       f;

1529:   PetscFree(s->numFieldComponents);
1530:   for (f = 0; f < s->numFields; ++f) {
1531:     PetscSectionDestroy(&s->field[f]);
1532:     PetscFree(s->fieldNames[f]);
1533:   }
1534:   PetscFree(s->fieldNames);
1535:   PetscFree(s->field);
1536:   PetscSectionDestroy(&s->bc);
1537:   PetscFree(s->bcIndices);
1538:   PetscFree2(s->atlasDof, s->atlasOff);
1539:   PetscSectionDestroy(&s->clSection);
1540:   ISDestroy(&s->clPoints);
1541:   ISDestroy(&s->perm);
1542:   PetscFree(s->clPerm);
1543:   PetscFree(s->clInvPerm);
1544:   PetscSectionSymDestroy(&s->sym);

1546:   s->pStart    = -1;
1547:   s->pEnd      = -1;
1548:   s->maxDof    = 0;
1549:   s->setup     = PETSC_FALSE;
1550:   s->numFields = 0;
1551:   s->clObj     = NULL;
1552:   return(0);
1553: }

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

1558:   Not collective

1560:   Input Parameters:
1561: . s - the PetscSection

1563:   Level: developer

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

1568: .seealso: PetscSection, PetscSectionCreate()
1569: @*/
1570: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1571: {

1575:   if (!*s) return(0);
1577:   if (--((PetscObject)(*s))->refct > 0) {
1578:     *s = NULL;
1579:     return(0);
1580:   }
1581:   PetscSectionReset(*s);
1582:   PetscHeaderDestroy(s);
1583:   return(0);
1584: }

1586: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1587: {
1588:   const PetscInt p = point - s->pStart;

1591:   *values = &baseArray[s->atlasOff[p]];
1592:   return(0);
1593: }

1595: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1596: {
1597:   PetscInt       *array;
1598:   const PetscInt p           = point - s->pStart;
1599:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1600:   PetscInt       cDim        = 0;

1604:   PetscSectionGetConstraintDof(s, p, &cDim);
1605:   array = &baseArray[s->atlasOff[p]];
1606:   if (!cDim) {
1607:     if (orientation >= 0) {
1608:       const PetscInt dim = s->atlasDof[p];
1609:       PetscInt       i;

1611:       if (mode == INSERT_VALUES) {
1612:         for (i = 0; i < dim; ++i) array[i] = values[i];
1613:       } else {
1614:         for (i = 0; i < dim; ++i) array[i] += values[i];
1615:       }
1616:     } else {
1617:       PetscInt offset = 0;
1618:       PetscInt j      = -1, field, i;

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

1623:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1624:         offset += dim;
1625:       }
1626:     }
1627:   } else {
1628:     if (orientation >= 0) {
1629:       const PetscInt dim  = s->atlasDof[p];
1630:       PetscInt       cInd = 0, i;
1631:       const PetscInt *cDof;

1633:       PetscSectionGetConstraintIndices(s, point, &cDof);
1634:       if (mode == INSERT_VALUES) {
1635:         for (i = 0; i < dim; ++i) {
1636:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1637:           array[i] = values[i];
1638:         }
1639:       } else {
1640:         for (i = 0; i < dim; ++i) {
1641:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1642:           array[i] += values[i];
1643:         }
1644:       }
1645:     } else {
1646:       const PetscInt *cDof;
1647:       PetscInt       offset  = 0;
1648:       PetscInt       cOffset = 0;
1649:       PetscInt       j       = 0, field;

1651:       PetscSectionGetConstraintIndices(s, point, &cDof);
1652:       for (field = 0; field < s->numFields; ++field) {
1653:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1654:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1655:         const PetscInt sDim = dim - tDim;
1656:         PetscInt       cInd = 0, i,k;

1658:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1659:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1660:           array[j] = values[k];
1661:         }
1662:         offset  += dim;
1663:         cOffset += dim - tDim;
1664:       }
1665:     }
1666:   }
1667:   return(0);
1668: }

1670: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1671: {
1675:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1676:   return(0);
1677: }

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

1682:   Input Parameters:
1683: + s     - The PetscSection
1684: - point - The point

1686:   Output Parameter:
1687: . indices - The constrained dofs

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

1691:   Level: advanced

1693: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1694: @*/
1695: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1696: {

1700:   if (s->bc) {
1701:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1702:   } else *indices = NULL;
1703:   return(0);
1704: }

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

1709:   Input Parameters:
1710: + s     - The PetscSection
1711: . point - The point
1712: - indices - The constrained dofs

1714:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1716:   Level: advanced

1718: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1719: @*/
1720: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1721: {

1725:   if (s->bc) {
1726:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1727:   }
1728:   return(0);
1729: }

1731: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1732: {

1736:   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);
1737:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1738:   return(0);
1739: }

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

1746:   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);
1747:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1748:   return(0);
1749: }

1751: /*@
1752:   PetscSectionPermute - Reorder the section according to the input point permutation

1754:   Collective on PetscSection

1756:   Input Parameter:
1757: + section - The PetscSection object
1758: - perm - The point permutation, old point p becomes new point perm[p]

1760:   Output Parameter:
1761: . sectionNew - The permuted PetscSection

1763:   Level: intermediate

1765: .keywords: mesh
1766: .seealso: MatPermute()
1767: @*/
1768: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1769: {
1770:   PetscSection    s = section, sNew;
1771:   const PetscInt *perm;
1772:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1773:   PetscErrorCode  ierr;

1779:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1780:   PetscSectionGetNumFields(s, &numFields);
1781:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1782:   for (f = 0; f < numFields; ++f) {
1783:     const char *name;
1784:     PetscInt    numComp;

1786:     PetscSectionGetFieldName(s, f, &name);
1787:     PetscSectionSetFieldName(sNew, f, name);
1788:     PetscSectionGetFieldComponents(s, f, &numComp);
1789:     PetscSectionSetFieldComponents(sNew, f, numComp);
1790:   }
1791:   ISGetLocalSize(permutation, &numPoints);
1792:   ISGetIndices(permutation, &perm);
1793:   PetscSectionGetChart(s, &pStart, &pEnd);
1794:   PetscSectionSetChart(sNew, pStart, pEnd);
1795:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1796:   for (p = pStart; p < pEnd; ++p) {
1797:     PetscInt dof, cdof;

1799:     PetscSectionGetDof(s, p, &dof);
1800:     PetscSectionSetDof(sNew, perm[p], dof);
1801:     PetscSectionGetConstraintDof(s, p, &cdof);
1802:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1803:     for (f = 0; f < numFields; ++f) {
1804:       PetscSectionGetFieldDof(s, p, f, &dof);
1805:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1806:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1807:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1808:     }
1809:   }
1810:   PetscSectionSetUp(sNew);
1811:   for (p = pStart; p < pEnd; ++p) {
1812:     const PetscInt *cind;
1813:     PetscInt        cdof;

1815:     PetscSectionGetConstraintDof(s, p, &cdof);
1816:     if (cdof) {
1817:       PetscSectionGetConstraintIndices(s, p, &cind);
1818:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1819:     }
1820:     for (f = 0; f < numFields; ++f) {
1821:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1822:       if (cdof) {
1823:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1824:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1825:       }
1826:     }
1827:   }
1828:   ISRestoreIndices(permutation, &perm);
1829:   *sectionNew = sNew;
1830:   return(0);
1831: }

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

1836:   Collective

1838:   Input Parameters:
1839: + sf - The SF
1840: - rootSection - Section defined on root space

1842:   Output Parameters:
1843: + remoteOffsets - root offsets in leaf storage, or NULL
1844: - leafSection - Section defined on the leaf space

1846:   Level: intermediate

1848: .seealso: PetscSFCreate()
1849: @*/
1850: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1851: {
1852:   PetscSF        embedSF;
1853:   const PetscInt *ilocal, *indices;
1854:   IS             selected;
1855:   PetscInt       numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

1859:   PetscSectionGetNumFields(rootSection, &numFields);
1860:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1861:   for (f = 0; f < numFields; ++f) {
1862:     PetscInt numComp = 0;
1863:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
1864:     PetscSectionSetFieldComponents(leafSection, f, numComp);
1865:   }
1866:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1867:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
1868:   rpEnd = PetscMin(rpEnd,nroots);
1869:   rpEnd = PetscMax(rpStart,rpEnd);
1870:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1871:   ISGetIndices(selected, &indices);
1872:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1873:   ISRestoreIndices(selected, &indices);
1874:   ISDestroy(&selected);
1875:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1876:   if (nleaves && ilocal) {
1877:     for (i = 0; i < nleaves; ++i) {
1878:       lpStart = PetscMin(lpStart, ilocal[i]);
1879:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
1880:     }
1881:     ++lpEnd;
1882:   } else {
1883:     lpStart = 0;
1884:     lpEnd   = nleaves;
1885:   }
1886:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
1887:   /* Could fuse these at the cost of a copy and extra allocation */
1888:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1889:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1890:   for (f = 0; f < numFields; ++f) {
1891:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1892:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1893:   }
1894:   if (remoteOffsets) {
1895:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
1896:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1897:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1898:   }
1899:   PetscSFDestroy(&embedSF);
1900:   PetscSectionSetUp(leafSection);
1901:   return(0);
1902: }

1904: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
1905: {
1906:   PetscSF         embedSF;
1907:   const PetscInt *indices;
1908:   IS              selected;
1909:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
1910:   PetscErrorCode  ierr;

1913:   *remoteOffsets = NULL;
1914:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
1915:   if (numRoots < 0) return(0);
1916:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1917:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1918:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1919:   ISGetIndices(selected, &indices);
1920:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1921:   ISRestoreIndices(selected, &indices);
1922:   ISDestroy(&selected);
1923:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
1924:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1925:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1926:   PetscSFDestroy(&embedSF);
1927:   return(0);
1928: }

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

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

1939:   Output Parameters:
1940: - sectionSF - The new SF

1942:   Note: Either rootSection or remoteOffsets can be specified

1944:   Level: intermediate

1946: .seealso: PetscSFCreate()
1947: @*/
1948: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1949: {
1950:   MPI_Comm          comm;
1951:   const PetscInt    *localPoints;
1952:   const PetscSFNode *remotePoints;
1953:   PetscInt          lpStart, lpEnd;
1954:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
1955:   PetscInt          *localIndices;
1956:   PetscSFNode       *remoteIndices;
1957:   PetscInt          i, ind;
1958:   PetscErrorCode    ierr;

1966:   PetscObjectGetComm((PetscObject)sf,&comm);
1967:   PetscSFCreate(comm, sectionSF);
1968:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1969:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
1970:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
1971:   if (numRoots < 0) return(0);
1972:   for (i = 0; i < numPoints; ++i) {
1973:     PetscInt localPoint = localPoints ? localPoints[i] : i;
1974:     PetscInt dof;

1976:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1977:       PetscSectionGetDof(leafSection, localPoint, &dof);
1978:       numIndices += dof;
1979:     }
1980:   }
1981:   PetscMalloc1(numIndices, &localIndices);
1982:   PetscMalloc1(numIndices, &remoteIndices);
1983:   /* Create new index graph */
1984:   for (i = 0, ind = 0; i < numPoints; ++i) {
1985:     PetscInt localPoint = localPoints ? localPoints[i] : i;
1986:     PetscInt rank       = remotePoints[i].rank;

1988:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1989:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
1990:       PetscInt loff, dof, d;

1992:       PetscSectionGetOffset(leafSection, localPoint, &loff);
1993:       PetscSectionGetDof(leafSection, localPoint, &dof);
1994:       for (d = 0; d < dof; ++d, ++ind) {
1995:         localIndices[ind]        = loff+d;
1996:         remoteIndices[ind].rank  = rank;
1997:         remoteIndices[ind].index = remoteOffset+d;
1998:       }
1999:     }
2000:   }
2001:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2002:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2003:   return(0);
2004: }

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

2009:   Input Parameters:
2010: + section   - The PetscSection
2011: . obj       - A PetscObject which serves as the key for this index
2012: . clSection - Section giving the size of the closure of each point
2013: - clPoints  - IS giving the points in each closure

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

2017:   Level: intermediate

2019: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2020: @*/
2021: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2022: {

2026:   if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2027:   section->clObj     = obj;
2028:   PetscSectionDestroy(&section->clSection);
2029:   ISDestroy(&section->clPoints);
2030:   section->clSection = clSection;
2031:   section->clPoints  = clPoints;
2032:   return(0);
2033: }

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

2038:   Input Parameters:
2039: + section   - The PetscSection
2040: - obj       - A PetscObject which serves as the key for this index

2042:   Output Parameters:
2043: + clSection - Section giving the size of the closure of each point
2044: - clPoints  - IS giving the points in each closure

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

2048:   Level: intermediate

2050: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2051: @*/
2052: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2053: {
2055:   if (section->clObj == obj) {
2056:     if (clSection) *clSection = section->clSection;
2057:     if (clPoints)  *clPoints  = section->clPoints;
2058:   } else {
2059:     if (clSection) *clSection = NULL;
2060:     if (clPoints)  *clPoints  = NULL;
2061:   }
2062:   return(0);
2063: }

2065: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2066: {
2067:   PetscInt       i;

2071:   if (section->clObj != obj) {
2072:     PetscSectionDestroy(&section->clSection);
2073:     ISDestroy(&section->clPoints);
2074:   }
2075:   section->clObj  = obj;
2076:   PetscFree(section->clPerm);
2077:   PetscFree(section->clInvPerm);
2078:   section->clSize = clSize;
2079:   if (mode == PETSC_COPY_VALUES) {
2080:     PetscMalloc1(clSize, &section->clPerm);
2081:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2082:     PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2083:   } else if (mode == PETSC_OWN_POINTER) {
2084:     section->clPerm = clPerm;
2085:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2086:   PetscMalloc1(clSize, &section->clInvPerm);
2087:   for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2088:   return(0);
2089: }

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

2094:   Input Parameters:
2095: + section - The PetscSection
2096: . obj     - A PetscObject which serves as the key for this index
2097: - perm    - Permutation of the cell dof closure

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

2102:   Level: intermediate

2104: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2105: @*/
2106: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2107: {
2108:   const PetscInt *clPerm = NULL;
2109:   PetscInt        clSize = 0;
2110:   PetscErrorCode  ierr;

2113:   if (perm) {
2114:     ISGetLocalSize(perm, &clSize);
2115:     ISGetIndices(perm, &clPerm);
2116:   }
2117:   PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2118:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2119:   return(0);
2120: }

2122: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2123: {
2125:   if (section->clObj == obj) {
2126:     if (size) *size = section->clSize;
2127:     if (perm) *perm = section->clPerm;
2128:   } else {
2129:     if (size) *size = 0;
2130:     if (perm) *perm = NULL;
2131:   }
2132:   return(0);
2133: }

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

2138:   Input Parameters:
2139: + section   - The PetscSection
2140: - obj       - A PetscObject which serves as the key for this index

2142:   Output Parameter:
2143: . perm - The dof closure permutation

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

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

2150:   Level: intermediate

2152: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2153: @*/
2154: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2155: {
2156:   const PetscInt *clPerm;
2157:   PetscInt        clSize;
2158:   PetscErrorCode  ierr;

2161:   PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2162:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2163:   return(0);
2164: }

2166: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2167: {
2169:   if (section->clObj == obj) {
2170:     if (size) *size = section->clSize;
2171:     if (perm) *perm = section->clInvPerm;
2172:   } else {
2173:     if (size) *size = 0;
2174:     if (perm) *perm = NULL;
2175:   }
2176:   return(0);
2177: }

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

2182:   Input Parameters:
2183: + section   - The PetscSection
2184: - obj       - A PetscObject which serves as the key for this index

2186:   Output Parameters:
2187: + size - The dof closure size
2188: - perm - The dof closure permutation

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

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

2195:   Level: intermediate

2197: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2198: @*/
2199: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2200: {
2201:   const PetscInt *clPerm;
2202:   PetscInt        clSize;
2203:   PetscErrorCode  ierr;

2206:   PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2207:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2208:   return(0);
2209: }

2211: /*@
2212:   PetscSectionGetField - Get the subsection associated with a single field

2214:   Input Parameters:
2215: + s     - The PetscSection
2216: - field - The field number

2218:   Output Parameter:
2219: . subs  - The subsection for the given field

2221:   Level: intermediate

2223: .seealso: PetscSectionSetNumFields()
2224: @*/
2225: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2226: {

2232:   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);
2233:   PetscObjectReference((PetscObject) s->field[field]);
2234:   *subs = s->field[field];
2235:   return(0);
2236: }

2238: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2239: PetscFunctionList PetscSectionSymList = NULL;

2241: /*@
2242:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2244:   Collective on MPI_Comm

2246:   Input Parameter:
2247: . comm - the MPI communicator

2249:   Output Parameter:
2250: . sym - pointer to the new set of symmetries

2252:   Level: developer

2254: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2255: @*/
2256: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2257: {

2262:   ISInitializePackage();
2263:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2264:   return(0);
2265: }

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

2270:   Collective on PetscSectionSym

2272:   Input Parameters:
2273: + sym    - The section symmetry object
2274: - method - The name of the section symmetry type

2276:   Level: developer

2278: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2279: @*/
2280: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2281: {
2282:   PetscErrorCode (*r)(PetscSectionSym);
2283:   PetscBool      match;

2288:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
2289:   if (match) return(0);

2291:   PetscFunctionListFind(PetscSectionSymList,method,&r);
2292:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2293:   if (sym->ops->destroy) {
2294:     (*sym->ops->destroy)(sym);
2295:     sym->ops->destroy = NULL;
2296:   }
2297:   (*r)(sym);
2298:   PetscObjectChangeTypeName((PetscObject)sym,method);
2299:   return(0);
2300: }


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

2306:   Not Collective

2308:   Input Parameter:
2309: . sym  - The section symmetry

2311:   Output Parameter:
2312: . type - The index set type name

2314:   Level: developer

2316: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2317: @*/
2318: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2319: {
2323:   *type = ((PetscObject)sym)->type_name;
2324:   return(0);
2325: }

2327: /*@C
2328:   PetscSectionSymRegister - Adds a new section symmetry implementation

2330:   Not Collective

2332:   Input Parameters:
2333: + name        - The name of a new user-defined creation routine
2334: - create_func - The creation routine itself

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

2339:   Level: developer

2341: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2342: @*/
2343: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2344: {

2348:   ISInitializePackage();
2349:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2350:   return(0);
2351: }

2353: /*@
2354:    PetscSectionSymDestroy - Destroys a section symmetry.

2356:    Collective on PetscSectionSym

2358:    Input Parameters:
2359: .  sym - the section symmetry

2361:    Level: developer

2363: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2364: @*/
2365: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2366: {
2367:   SymWorkLink    link,next;

2371:   if (!*sym) return(0);
2373:   if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2374:   if ((*sym)->ops->destroy) {
2375:     (*(*sym)->ops->destroy)(*sym);
2376:   }
2377:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2378:   for (link=(*sym)->workin; link; link=next) {
2379:     next = link->next;
2380:     PetscFree2(link->perms,link->rots);
2381:     PetscFree(link);
2382:   }
2383:   (*sym)->workin = NULL;
2384:   PetscHeaderDestroy(sym);
2385:   return(0);
2386: }

2388: /*@C
2389:    PetscSectionSymView - Displays a section symmetry

2391:    Collective on PetscSectionSym

2393:    Input Parameters:
2394: +  sym - the index set
2395: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2397:    Level: developer

2399: .seealso: PetscViewerASCIIOpen()
2400: @*/
2401: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2402: {

2407:   if (!viewer) {
2408:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2409:   }
2412:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2413:   if (sym->ops->view) {
2414:     (*sym->ops->view)(sym,viewer);
2415:   }
2416:   return(0);
2417: }

2419: /*@
2420:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2422:   Collective on PetscSection

2424:   Input Parameters:
2425: + section - the section describing data layout
2426: - sym - the symmetry describing the affect of orientation on the access of the data

2428:   Level: developer

2430: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2431: @*/
2432: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2433: {

2440:   PetscObjectReference((PetscObject)sym);
2441:   PetscSectionSymDestroy(&(section->sym));
2442:   section->sym = sym;
2443:   return(0);
2444: }

2446: /*@
2447:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

2449:   Not collective

2451:   Input Parameters:
2452: . section - the section describing data layout

2454:   Output Parameters:
2455: . sym - the symmetry describing the affect of orientation on the access of the data

2457:   Level: developer

2459: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2460: @*/
2461: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2462: {
2465:   *sym = section->sym;
2466:   return(0);
2467: }

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

2472:   Collective on PetscSection

2474:   Input Parameters:
2475: + section - the section describing data layout
2476: . field - the field number
2477: - sym - the symmetry describing the affect of orientation on the access of the data

2479:   Level: developer

2481: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2482: @*/
2483: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2484: {

2489:   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);
2490:   PetscSectionSetSym(section->field[field],sym);
2491:   return(0);
2492: }

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

2497:   Collective on PetscSection

2499:   Input Parameters:
2500: + section - the section describing data layout
2501: - field - the field number

2503:   Output Parameters:
2504: . sym - the symmetry describing the affect of orientation on the access of the data

2506:   Level: developer

2508: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2509: @*/
2510: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2511: {
2514:   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);
2515:   *sym = section->field[field]->sym;
2516:   return(0);
2517: }

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

2522:   Not collective

2524:   Input Parameters:
2525: + section - the section
2526: . numPoints - the number of points
2527: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2528:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2529:     context, see DMPlexGetConeOrientation()).

2531:   Output Parameter:
2532: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2533: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2534:     identity).

2536:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2537: .vb
2538:      const PetscInt    **perms;
2539:      const PetscScalar **rots;
2540:      PetscInt            lOffset;

2542:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2543:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2544:        PetscInt           point = points[2*i], dof, sOffset;
2545:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2546:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2548:        PetscSectionGetDof(section,point,&dof);
2549:        PetscSectionGetOffset(section,point,&sOffset);

2551:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
2552:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
2553:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
2554:        lOffset += dof;
2555:      }
2556:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2557: .ve

2559:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2560: .vb
2561:      const PetscInt    **perms;
2562:      const PetscScalar **rots;
2563:      PetscInt            lOffset;

2565:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2566:      for (i = 0, lOffset = 0; i < numPoints; i++) {
2567:        PetscInt           point = points[2*i], dof, sOffset;
2568:        const PetscInt    *perm  = perms ? perms[i] : NULL;
2569:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

2571:        PetscSectionGetDof(section,point,&dof);
2572:        PetscSectionGetOffset(section,point,&sOff);

2574:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2575:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
2576:        offset += dof;
2577:      }
2578:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2579: .ve

2581:   Level: developer

2583: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2584: @*/
2585: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2586: {
2587:   PetscSectionSym sym;
2588:   PetscErrorCode  ierr;

2593:   if (perms) *perms = NULL;
2594:   if (rots)  *rots  = NULL;
2595:   sym = section->sym;
2596:   if (sym && (perms || rots)) {
2597:     SymWorkLink link;

2599:     if (sym->workin) {
2600:       link        = sym->workin;
2601:       sym->workin = sym->workin->next;
2602:     } else {
2603:       PetscNewLog(sym,&link);
2604:     }
2605:     if (numPoints > link->numPoints) {
2606:       PetscFree2(link->perms,link->rots);
2607:       PetscMalloc2(numPoints,&link->perms,numPoints,&link->rots);
2608:       link->numPoints = numPoints;
2609:     }
2610:     link->next   = sym->workout;
2611:     sym->workout = link;
2612:     PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2613:     PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2614:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2615:     if (perms) *perms = link->perms;
2616:     if (rots)  *rots  = link->rots;
2617:   }
2618:   return(0);
2619: }

2621: /*@C
2622:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

2624:   Not collective

2626:   Input Parameters:
2627: + section - the section
2628: . numPoints - the number of points
2629: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2630:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2631:     context, see DMPlexGetConeOrientation()).

2633:   Output Parameter:
2634: + perms - The permutations for the given orientations: set to NULL at conclusion
2635: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2637:   Level: developer

2639: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2640: @*/
2641: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2642: {
2643:   PetscSectionSym sym;

2647:   sym = section->sym;
2648:   if (sym && (perms || rots)) {
2649:     SymWorkLink *p,link;

2651:     for (p=&sym->workout; (link=*p); p=&link->next) {
2652:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2653:         *p          = link->next;
2654:         link->next  = sym->workin;
2655:         sym->workin = link;
2656:         if (perms) *perms = NULL;
2657:         if (rots)  *rots  = NULL;
2658:         return(0);
2659:       }
2660:     }
2661:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2662:   }
2663:   return(0);
2664: }

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

2669:   Not collective

2671:   Input Parameters:
2672: + section - the section
2673: . field - the field of the section
2674: . numPoints - the number of points
2675: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2676:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2677:     context, see DMPlexGetConeOrientation()).

2679:   Output Parameter:
2680: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2681: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2682:     identity).

2684:   Level: developer

2686: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2687: @*/
2688: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2689: {

2694:   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);
2695:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2696:   return(0);
2697: }

2699: /*@C
2700:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

2702:   Not collective

2704:   Input Parameters:
2705: + section - the section
2706: . field - the field number
2707: . numPoints - the number of points
2708: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2709:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
2710:     context, see DMPlexGetConeOrientation()).

2712:   Output Parameter:
2713: + perms - The permutations for the given orientations: set to NULL at conclusion
2714: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

2716:   Level: developer

2718: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2719: @*/
2720: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2721: {

2726:   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);
2727:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
2728:   return(0);
2729: }