Actual source code: vsectionis.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

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

  9: PetscClassId PETSC_SECTION_CLASSID;

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

 16:   Collective on MPI_Comm

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

 22:   Level: developer

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

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

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

 44:   ISInitializePackage();

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

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

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

 71:   Collective on MPI_Comm

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

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

 79:   Level: developer

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

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

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

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

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

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

150:   Collective on MPI_Comm

152:   Input Parameter:
153: . section - the PetscSection

155:   Output Parameter:
156: . newSection - the copy

158:   Level: developer

160: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
161: @*/
162: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
163: {

167:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
168:   PetscSectionCopy(section, *newSection);
169:   return(0);
170: }

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

177:   Not collective

179:   Input Parameter:
180: . s - the PetscSection

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

185:   Level: intermediate

187: .seealso: PetscSectionSetNumFields()
188: @*/
189: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
190: {
193:   *numFields = s->numFields;
194:   return(0);
195: }

199: /*@
200:   PetscSectionSetNumFields - Sets the number of fields.

202:   Not collective

204:   Input Parameters:
205: + s - the PetscSection
206: - numFields - the number of fields

208:   Level: intermediate

210: .seealso: PetscSectionGetNumFields()
211: @*/
212: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
213: {
214:   PetscInt       f;

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

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

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

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

239: /*@C
240:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

242:   Not Collective

244:   Input Parameters:
245: + s     - the PetscSection
246: - field - the field number

248:   Output Parameter:
249: . fieldName - the field name

251:   Level: developer

253: .seealso: PetscSectionSetFieldName()
254: @*/
255: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
256: {
259:   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);
260:   *fieldName = s->fieldNames[field];
261:   return(0);
262: }

266: /*@C
267:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

269:   Not Collective

271:   Input Parameters:
272: + s     - the PetscSection
273: . field - the field number
274: - fieldName - the field name

276:   Level: developer

278: .seealso: PetscSectionGetFieldName()
279: @*/
280: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
281: {

286:   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);
287:   PetscFree(s->fieldNames[field]);
288:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
289:   return(0);
290: }

294: /*@
295:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

297:   Not collective

299:   Input Parameters:
300: + s - the PetscSection
301: - field - the field number

303:   Output Parameter:
304: . numComp - the number of field components

306:   Level: intermediate

308: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
309: @*/
310: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
311: {
314:   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);
315:   *numComp = s->numFieldComponents[field];
316:   return(0);
317: }

321: /*@
322:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

324:   Not collective

326:   Input Parameters:
327: + s - the PetscSection
328: . field - the field number
329: - numComp - the number of field components

331:   Level: intermediate

333: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
334: @*/
335: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
336: {
338:   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);
339:   s->numFieldComponents[field] = numComp;
340:   return(0);
341: }

345: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
346: {

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

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

362:   Not collective

364:   Input Parameter:
365: . s - the PetscSection

367:   Output Parameters:
368: + pStart - the first point
369: - pEnd - one past the last point

371:   Level: intermediate

373: .seealso: PetscSectionSetChart(), PetscSectionCreate()
374: @*/
375: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
376: {
379:   if (pStart) *pStart = s->pStart;
380:   if (pEnd)   *pEnd   = s->pEnd;
381:   return(0);
382: }

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

389:   Not collective

391:   Input Parameters:
392: + s - the PetscSection
393: . pStart - the first point
394: - pEnd - one past the last point

396:   Level: intermediate

398: .seealso: PetscSectionGetChart(), PetscSectionCreate()
399: @*/
400: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
401: {
402:   PetscInt       f;

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

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

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

428:   Not collective

430:   Input Parameter:
431: . s - the PetscSection

433:   Output Parameters:
434: . perm - The permutation as an IS

436:   Level: intermediate

438: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
439: @*/
440: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
441: {
445:   return(0);
446: }

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

453:   Not collective

455:   Input Parameters:
456: + s - the PetscSection
457: - perm - the permutation of points

459:   Level: intermediate

461: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
462: @*/
463: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
464: {

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

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

482:   Not collective

484:   Input Parameters:
485: + s - the PetscSection
486: - point - the point

488:   Output Parameter:
489: . numDof - the number of dof

491:   Level: intermediate

493: .seealso: PetscSectionSetDof(), PetscSectionCreate()
494: @*/
495: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
496: {
498: #if defined(PETSC_USE_DEBUG)
499:   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);
500: #endif
501:   *numDof = s->atlasDof[point - s->pStart];
502:   return(0);
503: }

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

510:   Not collective

512:   Input Parameters:
513: + s - the PetscSection
514: . point - the point
515: - numDof - the number of dof

517:   Level: intermediate

519: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
520: @*/
521: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
522: {
524:   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);
525:   s->atlasDof[point - s->pStart] = numDof;
526:   return(0);
527: }

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

534:   Not collective

536:   Input Parameters:
537: + s - the PetscSection
538: . point - the point
539: - numDof - the number of additional dof

541:   Level: intermediate

543: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
544: @*/
545: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
546: {
548:   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);
549:   s->atlasDof[point - s->pStart] += numDof;
550:   return(0);
551: }

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

558:   Not collective

560:   Input Parameters:
561: + s - the PetscSection
562: . point - the point
563: - field - the field

565:   Output Parameter:
566: . numDof - the number of dof

568:   Level: intermediate

570: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
571: @*/
572: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
573: {

577:   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);
578:   PetscSectionGetDof(s->field[field], point, numDof);
579:   return(0);
580: }

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

587:   Not collective

589:   Input Parameters:
590: + s - the PetscSection
591: . point - the point
592: . field - the field
593: - numDof - the number of dof

595:   Level: intermediate

597: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
598: @*/
599: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
600: {

604:   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);
605:   PetscSectionSetDof(s->field[field], point, numDof);
606:   return(0);
607: }

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

614:   Not collective

616:   Input Parameters:
617: + s - the PetscSection
618: . point - the point
619: . field - the field
620: - numDof - the number of dof

622:   Level: intermediate

624: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
625: @*/
626: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
627: {

631:   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);
632:   PetscSectionAddDof(s->field[field], point, numDof);
633:   return(0);
634: }

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

641:   Not collective

643:   Input Parameters:
644: + s - the PetscSection
645: - point - the point

647:   Output Parameter:
648: . numDof - the number of dof which are fixed by constraints

650:   Level: intermediate

652: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
653: @*/
654: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
655: {

659:   if (s->bc) {
660:     PetscSectionGetDof(s->bc, point, numDof);
661:   } else *numDof = 0;
662:   return(0);
663: }

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

670:   Not collective

672:   Input Parameters:
673: + s - the PetscSection
674: . point - the point
675: - numDof - the number of dof which are fixed by constraints

677:   Level: intermediate

679: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
680: @*/
681: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
682: {

686:   if (numDof) {
687:     PetscSectionCheckConstraints_Static(s);
688:     PetscSectionSetDof(s->bc, point, numDof);
689:   }
690:   return(0);
691: }

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

698:   Not collective

700:   Input Parameters:
701: + s - the PetscSection
702: . point - the point
703: - numDof - the number of additional dof which are fixed by constraints

705:   Level: intermediate

707: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
708: @*/
709: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
710: {

714:   if (numDof) {
715:     PetscSectionCheckConstraints_Static(s);
716:     PetscSectionAddDof(s->bc, point, numDof);
717:   }
718:   return(0);
719: }

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

726:   Not collective

728:   Input Parameters:
729: + s - the PetscSection
730: . point - the point
731: - field - the field

733:   Output Parameter:
734: . numDof - the number of dof which are fixed by constraints

736:   Level: intermediate

738: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
739: @*/
740: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
741: {

745:   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);
746:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
747:   return(0);
748: }

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

755:   Not collective

757:   Input Parameters:
758: + s - the PetscSection
759: . point - the point
760: . field - the field
761: - numDof - the number of dof which are fixed by constraints

763:   Level: intermediate

765: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
766: @*/
767: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
768: {

772:   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);
773:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
774:   return(0);
775: }

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

782:   Not collective

784:   Input Parameters:
785: + s - the PetscSection
786: . point - the point
787: . field - the field
788: - numDof - the number of additional dof which are fixed by constraints

790:   Level: intermediate

792: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
793: @*/
794: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
795: {

799:   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);
800:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
801:   return(0);
802: }

806: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
807: {

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

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

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

825:   Not collective

827:   Input Parameter:
828: . s - the PetscSection

830:   Level: intermediate

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

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

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

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

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

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

876:   Not collective

878:   Input Parameters:
879: . s - the PetscSection

881:   Output Parameter:
882: . maxDof - the maximum dof

884:   Level: intermediate

886: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
887: @*/
888: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
889: {
891:   *maxDof = s->maxDof;
892:   return(0);
893: }

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

900:   Not collective

902:   Input Parameters:
903: + s - the PetscSection
904: - point - the point

906:   Output Parameter:
907: . size - the size of an array which can hold all the dofs

909:   Level: intermediate

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

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

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

928:   Not collective

930:   Input Parameters:
931: + s - the PetscSection
932: - point - the point

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

937:   Level: intermediate

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

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

956: /*@
957:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
958:   the local section and an SF describing the section point overlap.

960:   Input Parameters:
961:   + s - The PetscSection for the local field layout
962:   . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
963:   - includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs

965:   Output Parameter:
966:   . gsection - The PetscSection for the global field layout

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

970:   Level: developer

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

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

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

1046: /*@
1047:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1048:   the local section and an SF describing the section point overlap.

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

1057:   Output Parameter:
1058:   . gsection - The PetscSection for the global field layout

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

1062:   Level: developer

1064: .seealso: PetscSectionCreate()
1065: @*/
1066: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1067: {
1068:   const PetscInt *pind = NULL;
1069:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1070:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1071:   PetscErrorCode  ierr;

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

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

1142: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1143: {
1144:   PetscInt       pStart, pEnd, p, localSize = 0;

1148:   PetscSectionGetChart(s, &pStart, &pEnd);
1149:   for (p = pStart; p < pEnd; ++p) {
1150:     PetscInt dof;

1152:     PetscSectionGetDof(s, p, &dof);
1153:     if (dof > 0) ++localSize;
1154:   }
1155:   PetscLayoutCreate(comm, layout);
1156:   PetscLayoutSetLocalSize(*layout, localSize);
1157:   PetscLayoutSetBlockSize(*layout, 1);
1158:   PetscLayoutSetUp(*layout);
1159:   return(0);
1160: }

1164: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1165: {
1166:   PetscInt       pStart, pEnd, p, localSize = 0;

1170:   PetscSectionGetChart(s, &pStart, &pEnd);
1171:   for (p = pStart; p < pEnd; ++p) {
1172:     PetscInt dof,cdof;

1174:     PetscSectionGetDof(s, p, &dof);
1175:     PetscSectionGetConstraintDof(s, p, &cdof);
1176:     if (dof-cdof > 0) localSize += dof-cdof;
1177:   }
1178:   PetscLayoutCreate(comm, layout);
1179:   PetscLayoutSetLocalSize(*layout, localSize);
1180:   PetscLayoutSetBlockSize(*layout, 1);
1181:   PetscLayoutSetUp(*layout);
1182:   return(0);
1183: }

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

1190:   Not collective

1192:   Input Parameters:
1193: + s - the PetscSection
1194: - point - the point

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

1199:   Level: intermediate

1201: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1202: @*/
1203: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1204: {
1206: #if defined(PETSC_USE_DEBUG)
1207:   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);
1208: #endif
1209:   *offset = s->atlasOff[point - s->pStart];
1210:   return(0);
1211: }

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

1218:   Not collective

1220:   Input Parameters:
1221: + s - the PetscSection
1222: . point - the point
1223: - offset - the offset

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

1227:   Level: intermediate

1229: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1230: @*/
1231: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1232: {
1234:   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);
1235:   s->atlasOff[point - s->pStart] = offset;
1236:   return(0);
1237: }

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

1244:   Not collective

1246:   Input Parameters:
1247: + s - the PetscSection
1248: . point - the point
1249: - field - the field

1251:   Output Parameter:
1252: . offset - the offset

1254:   Level: intermediate

1256: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1257: @*/
1258: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1259: {

1263:   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);
1264:   PetscSectionGetOffset(s->field[field], point, offset);
1265:   return(0);
1266: }

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

1273:   Not collective

1275:   Input Parameters:
1276: + s - the PetscSection
1277: . point - the point
1278: . field - the field
1279: - offset - the offset

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

1283:   Level: intermediate

1285: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1286: @*/
1287: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1288: {

1292:   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);
1293:   PetscSectionSetOffset(s->field[field], point, offset);
1294:   return(0);
1295: }

1299: /* This gives the offset on a point of the field, ignoring constraints */
1300: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1301: {
1302:   PetscInt       off, foff;

1306:   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);
1307:   PetscSectionGetOffset(s, point, &off);
1308:   PetscSectionGetOffset(s->field[field], point, &foff);
1309:   *offset = foff - off;
1310:   return(0);
1311: }

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

1318:   Not collective

1320:   Input Parameter:
1321: . s - the PetscSection

1323:   Output Parameters:
1324: + start - the minimum offset
1325: - end   - one more than the maximum offset

1327:   Level: intermediate

1329: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1330: @*/
1331: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1332: {
1333:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1337:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1338:   PetscSectionGetChart(s, &pStart, &pEnd);
1339:   for (p = 0; p < pEnd-pStart; ++p) {
1340:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1342:     if (off >= 0) {
1343:       os = PetscMin(os, off);
1344:       oe = PetscMax(oe, off+dof);
1345:     }
1346:   }
1347:   if (start) *start = os;
1348:   if (end)   *end   = oe;
1349:   return(0);
1350: }

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

1360:   if (!numFields) return(0);
1361:   PetscSectionGetNumFields(s, &nF);
1362:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1363:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1364:   PetscSectionSetNumFields(*subs, numFields);
1365:   for (f = 0; f < numFields; ++f) {
1366:     const char *name   = NULL;
1367:     PetscInt   numComp = 0;

1369:     PetscSectionGetFieldName(s, fields[f], &name);
1370:     PetscSectionSetFieldName(*subs, f, name);
1371:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1372:     PetscSectionSetFieldComponents(*subs, f, numComp);
1373:   }
1374:   PetscSectionGetChart(s, &pStart, &pEnd);
1375:   PetscSectionSetChart(*subs, pStart, pEnd);
1376:   for (p = pStart; p < pEnd; ++p) {
1377:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1379:     for (f = 0; f < numFields; ++f) {
1380:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1381:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1382:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1383:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1384:       dof  += fdof;
1385:       cdof += cfdof;
1386:     }
1387:     PetscSectionSetDof(*subs, p, dof);
1388:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1389:     maxCdof = PetscMax(cdof, maxCdof);
1390:   }
1391:   PetscSectionSetUp(*subs);
1392:   if (maxCdof) {
1393:     PetscInt *indices;

1395:     PetscMalloc1(maxCdof, &indices);
1396:     for (p = pStart; p < pEnd; ++p) {
1397:       PetscInt cdof;

1399:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1400:       if (cdof) {
1401:         const PetscInt *oldIndices = NULL;
1402:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1407:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1408:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1409:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1410:           /* This can be sped up if we assume sorted fields */
1411:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1412:             PetscInt oldfdof = 0;
1413:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1414:             oldFoff += oldfdof;
1415:           }
1416:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1417:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1418:           numConst += cfdof;
1419:           fOff     += fdof;
1420:         }
1421:         PetscSectionSetConstraintIndices(*subs, p, indices);
1422:       }
1423:     }
1424:     PetscFree(indices);
1425:   }
1426:   return(0);
1427: }

1431: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1432: {
1433:   const PetscInt *points = NULL, *indices = NULL;
1434:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1438:   PetscSectionGetNumFields(s, &numFields);
1439:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1440:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1441:   for (f = 0; f < numFields; ++f) {
1442:     const char *name   = NULL;
1443:     PetscInt   numComp = 0;

1445:     PetscSectionGetFieldName(s, f, &name);
1446:     PetscSectionSetFieldName(*subs, f, name);
1447:     PetscSectionGetFieldComponents(s, f, &numComp);
1448:     PetscSectionSetFieldComponents(*subs, f, numComp);
1449:   }
1450:   /* For right now, we do not try to squeeze the subchart */
1451:   if (subpointMap) {
1452:     ISGetSize(subpointMap, &numSubpoints);
1453:     ISGetIndices(subpointMap, &points);
1454:   }
1455:   PetscSectionGetChart(s, &pStart, &pEnd);
1456:   PetscSectionSetChart(*subs, 0, numSubpoints);
1457:   for (p = pStart; p < pEnd; ++p) {
1458:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1460:     PetscFindInt(p, numSubpoints, points, &subp);
1461:     if (subp < 0) continue;
1462:     for (f = 0; f < numFields; ++f) {
1463:       PetscSectionGetFieldDof(s, p, f, &fdof);
1464:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1465:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1466:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1467:     }
1468:     PetscSectionGetDof(s, p, &dof);
1469:     PetscSectionSetDof(*subs, subp, dof);
1470:     PetscSectionGetConstraintDof(s, p, &cdof);
1471:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1472:   }
1473:   PetscSectionSetUp(*subs);
1474:   /* Change offsets to original offsets */
1475:   for (p = pStart; p < pEnd; ++p) {
1476:     PetscInt off, foff = 0;

1478:     PetscFindInt(p, numSubpoints, points, &subp);
1479:     if (subp < 0) continue;
1480:     for (f = 0; f < numFields; ++f) {
1481:       PetscSectionGetFieldOffset(s, p, f, &foff);
1482:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1483:     }
1484:     PetscSectionGetOffset(s, p, &off);
1485:     PetscSectionSetOffset(*subs, subp, off);
1486:   }
1487:   /* Copy constraint indices */
1488:   for (subp = 0; subp < numSubpoints; ++subp) {
1489:     PetscInt cdof;

1491:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1492:     if (cdof) {
1493:       for (f = 0; f < numFields; ++f) {
1494:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1495:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1496:       }
1497:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1498:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1499:     }
1500:   }
1501:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1502:   return(0);
1503: }

1507: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1508: {
1509:   PetscInt       p;
1510:   PetscMPIInt    rank;

1514:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1515:   PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
1516:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1517:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1518:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1519:       PetscInt b;

1521:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1522:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1523:         PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1524:       }
1525:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1526:     } else {
1527:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1528:     }
1529:   }
1530:   PetscViewerFlush(viewer);
1531:   return(0);
1532: }

1536: /*@C
1537:   PetscSectionView - Views a PetscSection

1539:   Collective on PetscSection

1541:   Input Parameters:
1542: + s - the PetscSection object to view
1543: - v - the viewer

1545:   Level: developer

1547: .seealso PetscSectionCreate(), PetscSectionDestroy()
1548: @*/
1549: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1550: {
1551:   PetscBool      isascii;
1552:   PetscInt       f;

1556:   if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
1558:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1559:   if (isascii) {
1560:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1561:     if (s->numFields) {
1562:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1563:       for (f = 0; f < s->numFields; ++f) {
1564:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1565:         PetscSectionView_ASCII(s->field[f], viewer);
1566:       }
1567:     } else {
1568:       PetscSectionView_ASCII(s, viewer);
1569:     }
1570:   }
1571:   return(0);
1572: }

1576: /*@
1577:   PetscSectionReset - Frees all section data.

1579:   Not collective

1581:   Input Parameters:
1582: . s - the PetscSection

1584:   Level: developer

1586: .seealso: PetscSection, PetscSectionCreate()
1587: @*/
1588: PetscErrorCode PetscSectionReset(PetscSection s)
1589: {
1590:   PetscInt       f;

1594:   PetscFree(s->numFieldComponents);
1595:   for (f = 0; f < s->numFields; ++f) {
1596:     PetscSectionDestroy(&s->field[f]);
1597:     PetscFree(s->fieldNames[f]);
1598:   }
1599:   PetscFree(s->fieldNames);
1600:   PetscFree(s->field);
1601:   PetscSectionDestroy(&s->bc);
1602:   PetscFree(s->bcIndices);
1603:   PetscFree2(s->atlasDof, s->atlasOff);
1604:   PetscSectionDestroy(&s->clSection);
1605:   ISDestroy(&s->clPoints);
1606:   ISDestroy(&s->perm);

1608:   s->pStart    = -1;
1609:   s->pEnd      = -1;
1610:   s->maxDof    = 0;
1611:   s->setup     = PETSC_FALSE;
1612:   s->numFields = 0;
1613:   s->clObj     = NULL;
1614:   return(0);
1615: }

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

1622:   Not collective

1624:   Input Parameters:
1625: . s - the PetscSection

1627:   Level: developer

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

1632: .seealso: PetscSection, PetscSectionCreate()
1633: @*/
1634: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1635: {

1639:   if (!*s) return(0);
1641:   if (--((PetscObject)(*s))->refct > 0) {
1642:     *s = NULL;
1643:     return(0);
1644:   }
1645:   PetscSectionReset(*s);
1646:   PetscHeaderDestroy(s);
1647:   return(0);
1648: }

1652: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1653: {
1654:   const PetscInt p = point - s->pStart;

1657:   *values = &baseArray[s->atlasOff[p]];
1658:   return(0);
1659: }

1663: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1664: {
1665:   PetscInt       *array;
1666:   const PetscInt p           = point - s->pStart;
1667:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1668:   PetscInt       cDim        = 0;

1672:   PetscSectionGetConstraintDof(s, p, &cDim);
1673:   array = &baseArray[s->atlasOff[p]];
1674:   if (!cDim) {
1675:     if (orientation >= 0) {
1676:       const PetscInt dim = s->atlasDof[p];
1677:       PetscInt       i;

1679:       if (mode == INSERT_VALUES) {
1680:         for (i = 0; i < dim; ++i) array[i] = values[i];
1681:       } else {
1682:         for (i = 0; i < dim; ++i) array[i] += values[i];
1683:       }
1684:     } else {
1685:       PetscInt offset = 0;
1686:       PetscInt j      = -1, field, i;

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

1691:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1692:         offset += dim;
1693:       }
1694:     }
1695:   } else {
1696:     if (orientation >= 0) {
1697:       const PetscInt dim  = s->atlasDof[p];
1698:       PetscInt       cInd = 0, i;
1699:       const PetscInt *cDof;

1701:       PetscSectionGetConstraintIndices(s, point, &cDof);
1702:       if (mode == INSERT_VALUES) {
1703:         for (i = 0; i < dim; ++i) {
1704:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1705:           array[i] = values[i];
1706:         }
1707:       } else {
1708:         for (i = 0; i < dim; ++i) {
1709:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1710:           array[i] += values[i];
1711:         }
1712:       }
1713:     } else {
1714:       const PetscInt *cDof;
1715:       PetscInt       offset  = 0;
1716:       PetscInt       cOffset = 0;
1717:       PetscInt       j       = 0, field;

1719:       PetscSectionGetConstraintIndices(s, point, &cDof);
1720:       for (field = 0; field < s->numFields; ++field) {
1721:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1722:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1723:         const PetscInt sDim = dim - tDim;
1724:         PetscInt       cInd = 0, i,k;

1726:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1727:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1728:           array[j] = values[k];
1729:         }
1730:         offset  += dim;
1731:         cOffset += dim - tDim;
1732:       }
1733:     }
1734:   }
1735:   return(0);
1736: }

1740: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1741: {
1745:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1746:   return(0);
1747: }

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

1754:   Input Parameters:
1755: + s     - The PetscSection
1756: - point - The point

1758:   Output Parameter:
1759: . indices - The constrained dofs

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

1763:   Level: advanced

1765: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1766: @*/
1767: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1768: {

1772:   if (s->bc) {
1773:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1774:   } else *indices = NULL;
1775:   return(0);
1776: }

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

1783:   Input Parameters:
1784: + s     - The PetscSection
1785: . point - The point
1786: - indices - The constrained dofs

1788:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1790:   Level: advanced

1792: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1793: @*/
1794: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1795: {

1799:   if (s->bc) {
1800:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1801:   }
1802:   return(0);
1803: }

1807: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1808: {

1812:   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);
1813:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1814:   return(0);
1815: }

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

1824:   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);
1825:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1826:   return(0);
1827: }

1831: /*@
1832:   PetscSectionPermute - Reorder the section according to the input point permutation

1834:   Collective on PetscSection

1836:   Input Parameter:
1837: + section - The PetscSection object
1838: - perm - The point permutation, old point p becomes new point perm[p]

1840:   Output Parameter:
1841: . sectionNew - The permuted PetscSection

1843:   Level: intermediate

1845: .keywords: mesh
1846: .seealso: MatPermute()
1847: @*/
1848: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1849: {
1850:   PetscSection    s = section, sNew;
1851:   const PetscInt *perm;
1852:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1853:   PetscErrorCode  ierr;

1859:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1860:   PetscSectionGetNumFields(s, &numFields);
1861:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1862:   for (f = 0; f < numFields; ++f) {
1863:     const char *name;
1864:     PetscInt    numComp;

1866:     PetscSectionGetFieldName(s, f, &name);
1867:     PetscSectionSetFieldName(sNew, f, name);
1868:     PetscSectionGetFieldComponents(s, f, &numComp);
1869:     PetscSectionSetFieldComponents(sNew, f, numComp);
1870:   }
1871:   ISGetLocalSize(permutation, &numPoints);
1872:   ISGetIndices(permutation, &perm);
1873:   PetscSectionGetChart(s, &pStart, &pEnd);
1874:   PetscSectionSetChart(sNew, pStart, pEnd);
1875:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1876:   for (p = pStart; p < pEnd; ++p) {
1877:     PetscInt dof, cdof;

1879:     PetscSectionGetDof(s, p, &dof);
1880:     PetscSectionSetDof(sNew, perm[p], dof);
1881:     PetscSectionGetConstraintDof(s, p, &cdof);
1882:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1883:     for (f = 0; f < numFields; ++f) {
1884:       PetscSectionGetFieldDof(s, p, f, &dof);
1885:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1886:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1887:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1888:     }
1889:   }
1890:   PetscSectionSetUp(sNew);
1891:   for (p = pStart; p < pEnd; ++p) {
1892:     const PetscInt *cind;
1893:     PetscInt        cdof;

1895:     PetscSectionGetConstraintDof(s, p, &cdof);
1896:     if (cdof) {
1897:       PetscSectionGetConstraintIndices(s, p, &cind);
1898:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1899:     }
1900:     for (f = 0; f < numFields; ++f) {
1901:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1902:       if (cdof) {
1903:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1904:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1905:       }
1906:     }
1907:   }
1908:   ISRestoreIndices(permutation, &perm);
1909:   *sectionNew = sNew;
1910:   return(0);
1911: }

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

1918:   Collective

1920:   Input Parameters:
1921: + sf - The SF
1922: - rootSection - Section defined on root space

1924:   Output Parameters:
1925: + remoteOffsets - root offsets in leaf storage, or NULL
1926: - leafSection - Section defined on the leaf space

1928:   Level: intermediate

1930: .seealso: PetscSFCreate()
1931: @*/
1932: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1933: {
1934:   PetscSF        embedSF;
1935:   const PetscInt *ilocal, *indices;
1936:   IS             selected;
1937:   PetscInt       numFields, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

1941:   PetscSectionGetNumFields(rootSection, &numFields);
1942:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1943:   for (f = 0; f < numFields; ++f) {
1944:     PetscInt numComp = 0;
1945:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
1946:     PetscSectionSetFieldComponents(leafSection, f, numComp);
1947:   }
1948:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1949:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1950:   ISGetIndices(selected, &indices);
1951:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1952:   ISRestoreIndices(selected, &indices);
1953:   ISDestroy(&selected);
1954:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1955:   if (nleaves && ilocal) {
1956:     for (i = 0; i < nleaves; ++i) {
1957:       lpStart = PetscMin(lpStart, ilocal[i]);
1958:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
1959:     }
1960:     ++lpEnd;
1961:   } else {
1962:     lpStart = 0;
1963:     lpEnd   = nleaves;
1964:   }
1965:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
1966:   /* Could fuse these at the cost of a copy and extra allocation */
1967:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1968:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1969:   for (f = 0; f < numFields; ++f) {
1970:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1971:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1972:   }
1973:   if (remoteOffsets) {
1974:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
1975:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1976:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1977:   }
1978:   PetscSFDestroy(&embedSF);
1979:   PetscSectionSetUp(leafSection);
1980:   return(0);
1981: }

1985: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
1986: {
1987:   PetscSF         embedSF;
1988:   const PetscInt *indices;
1989:   IS              selected;
1990:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
1991:   PetscErrorCode  ierr;

1994:   *remoteOffsets = NULL;
1995:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
1996:   if (numRoots < 0) return(0);
1997:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1998:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1999:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2000:   ISGetIndices(selected, &indices);
2001:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2002:   ISRestoreIndices(selected, &indices);
2003:   ISDestroy(&selected);
2004:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2005:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2006:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2007:   PetscSFDestroy(&embedSF);
2008:   return(0);
2009: }

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

2016:   Input Parameters:
2017: + sf - The SF
2018: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section), or NULL
2019: - remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL

2021:   Output Parameters:
2022: + leafSection - Data layout of local points for incoming data  (this is the distributed section)
2023: - sectionSF - The new SF

2025:   Note: Either rootSection or remoteOffsets can be specified

2027:   Level: intermediate

2029: .seealso: PetscSFCreate()
2030: @*/
2031: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2032: {
2033:   MPI_Comm          comm;
2034:   const PetscInt    *localPoints;
2035:   const PetscSFNode *remotePoints;
2036:   PetscInt          lpStart, lpEnd;
2037:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2038:   PetscInt          *localIndices;
2039:   PetscSFNode       *remoteIndices;
2040:   PetscInt          i, ind;
2041:   PetscErrorCode    ierr;

2049:   PetscObjectGetComm((PetscObject)sf,&comm);
2050:   PetscSFCreate(comm, sectionSF);
2051:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2052:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2053:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2054:   if (numRoots < 0) return(0);
2055:   for (i = 0; i < numPoints; ++i) {
2056:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2057:     PetscInt dof;

2059:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2060:       PetscSectionGetDof(leafSection, localPoint, &dof);
2061:       numIndices += dof;
2062:     }
2063:   }
2064:   PetscMalloc1(numIndices, &localIndices);
2065:   PetscMalloc1(numIndices, &remoteIndices);
2066:   /* Create new index graph */
2067:   for (i = 0, ind = 0; i < numPoints; ++i) {
2068:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2069:     PetscInt rank       = remotePoints[i].rank;

2071:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2072:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2073:       PetscInt loff, dof, d;

2075:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2076:       PetscSectionGetDof(leafSection, localPoint, &dof);
2077:       for (d = 0; d < dof; ++d, ++ind) {
2078:         localIndices[ind]        = loff+d;
2079:         remoteIndices[ind].rank  = rank;
2080:         remoteIndices[ind].index = remoteOffset+d;
2081:       }
2082:     }
2083:   }
2084:   PetscFree(remoteOffsets);
2085:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2086:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2087:   return(0);
2088: }

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

2095:   Input Parameters:
2096: + section   - The PetscSection
2097: . obj       - A PetscObject which serves as the key for this index
2098: . clSection - Section giving the size of the closure of each point
2099: - clPoints  - IS giving the points in each closure

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

2103:   Level: intermediate

2105: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2106: @*/
2107: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2108: {

2112:   section->clObj     = obj;
2113:   PetscSectionDestroy(&section->clSection);
2114:   ISDestroy(&section->clPoints);
2115:   section->clSection = clSection;
2116:   section->clPoints  = clPoints;
2117:   return(0);
2118: }

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

2125:   Input Parameters:
2126: + section   - The PetscSection
2127: - obj       - A PetscObject which serves as the key for this index

2129:   Output Parameters:
2130: + clSection - Section giving the size of the closure of each point
2131: - clPoints  - IS giving the points in each closure

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

2135:   Level: intermediate

2137: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2138: @*/
2139: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2140: {
2142:   if (section->clObj == obj) {
2143:     if (clSection) *clSection = section->clSection;
2144:     if (clPoints)  *clPoints  = section->clPoints;
2145:   } else {
2146:     if (clSection) *clSection = NULL;
2147:     if (clPoints)  *clPoints  = NULL;
2148:   }
2149:   return(0);
2150: }

2154: /*@
2155:   PetscSectionGetField - Get the subsection associated with a single field

2157:   Input Parameters:
2158: + s     - The PetscSection
2159: - field - The field number

2161:   Output Parameter:
2162: . subs  - The subsection for the given field

2164:   Level: intermediate

2166: .seealso: PetscSectionSetNumFields()
2167: @*/
2168: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2169: {

2175:   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);
2176:   PetscObjectReference((PetscObject) s->field[field]);
2177:   *subs = s->field[field];
2178:   return(0);
2179: }