Actual source code: vsectionis.c

petsc-master 2014-12-24
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,_p_PetscSection,int,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:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1129:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1130:     if (nroots > pEnd - pStart) {
1131:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1132:     }
1133:   }
1134:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1135:   PetscFree(neg);
1136:   return(0);
1137: }

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

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

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

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

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

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

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

1189:   Not collective

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

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

1198:   Level: intermediate

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

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

1217:   Not collective

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

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

1226:   Level: intermediate

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

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

1243:   Not collective

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

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

1253:   Level: intermediate

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

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

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

1272:   Not collective

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

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

1282:   Level: intermediate

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

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

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

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

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

1317:   Not collective

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

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

1326:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1538:   Collective on PetscSection

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

1544:   Level: developer

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

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

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

1578:   Not collective

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

1583:   Level: developer

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

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

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

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

1621:   Not collective

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

1626:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1762:   Level: advanced

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

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

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

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

1787:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1789:   Level: advanced

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

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

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

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

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

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

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

1833:   Collective on PetscSection

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

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

1842:   Level: intermediate

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

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

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

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

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

1912: /*
1913:   I need extract and merge routines for section based on fields. This should be trivial except for updating the
1914:   constraint indices.

1916:   Then I need a new interface for DMCreateDecomposition that takes groups of fields and returns a real DMPlex
1917:   that shares the mesh parts, and has the extracted section
1918: */

1922: PetscErrorCode PetscSFConvertPartition(PetscSF sfPart, PetscSection partSection, IS partition, ISLocalToGlobalMapping *renumbering, PetscSF *sf)
1923: {
1924:   MPI_Comm       comm;
1925:   PetscSF        sfPoints;
1926:   PetscInt       *partSizes,*partOffsets,p,i,numParts,numMyPoints,numPoints,count;
1927:   const PetscInt *partArray;
1928:   PetscSFNode    *sendPoints;
1929:   PetscMPIInt    rank;

1933:   PetscObjectGetComm((PetscObject)sfPart,&comm);
1934:   MPI_Comm_rank(comm, &rank);

1936:   /* Get the number of parts and sizes that I have to distribute */
1937:   PetscSFGetGraph(sfPart,NULL,&numParts,NULL,NULL);
1938:   PetscMalloc2(numParts,&partSizes,numParts,&partOffsets);
1939:   for (p=0,numPoints=0; p<numParts; p++) {
1940:     PetscSectionGetDof(partSection, p, &partSizes[p]);
1941:     numPoints += partSizes[p];
1942:   }
1943:   numMyPoints = 0;
1944:   PetscSFFetchAndOpBegin(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1945:   PetscSFFetchAndOpEnd(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1946:   /* I will receive numMyPoints. I will send a total of numPoints, to be placed on remote procs at partOffsets */

1948:   /* Create SF mapping locations (addressed through partition, as indexed leaves) to new owners (roots) */
1949:   PetscMalloc1(numPoints,&sendPoints);
1950:   for (p=0,count=0; p<numParts; p++) {
1951:     for (i=0; i<partSizes[p]; i++) {
1952:       sendPoints[count].rank = p;
1953:       sendPoints[count].index = partOffsets[p]+i;
1954:       count++;
1955:     }
1956:   }
1957:   if (count != numPoints) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count %D should equal numPoints=%D",count,numPoints);
1958:   PetscFree2(partSizes,partOffsets);
1959:   ISGetIndices(partition,&partArray);
1960:   PetscSFCreate(comm,&sfPoints);
1961:   PetscSFSetFromOptions(sfPoints);
1962:   PetscSFSetGraph(sfPoints,numMyPoints,numPoints,partArray,PETSC_USE_POINTER,sendPoints,PETSC_OWN_POINTER);

1964:   /* Invert SF so that the new owners are leaves and the locations indexed through partition are the roots */
1965:   PetscSFCreateInverseSF(sfPoints,sf);
1966:   PetscSFDestroy(&sfPoints);
1967:   ISRestoreIndices(partition,&partArray);

1969:   /* Create the new local-to-global mapping */
1970:   ISLocalToGlobalMappingCreateSF(*sf,0,renumbering);
1971:   return(0);
1972: }

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

1979:   Collective

1981:   Input Parameters:
1982: + sf - The SF
1983: - rootSection - Section defined on root space

1985:   Output Parameters:
1986: + remoteOffsets - root offsets in leaf storage, or NULL
1987: - leafSection - Section defined on the leaf space

1989:   Level: intermediate

1991: .seealso: PetscSFCreate()
1992: @*/
1993: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1994: {
1995:   PetscSF        embedSF;
1996:   const PetscInt *ilocal, *indices;
1997:   IS             selected;
1998:   PetscInt       numFields, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

2002:   PetscSectionGetNumFields(rootSection, &numFields);
2003:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2004:   for (f = 0; f < numFields; ++f) {
2005:     PetscInt numComp = 0;
2006:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2007:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2008:   }
2009:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2010:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2011:   ISGetIndices(selected, &indices);
2012:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2013:   ISRestoreIndices(selected, &indices);
2014:   ISDestroy(&selected);
2015:   PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2016:   if (nleaves && ilocal) {
2017:     for (i = 0; i < nleaves; ++i) {
2018:       lpStart = PetscMin(lpStart, ilocal[i]);
2019:       lpEnd   = PetscMax(lpEnd,   ilocal[i]);
2020:     }
2021:     ++lpEnd;
2022:   } else {
2023:     lpStart = 0;
2024:     lpEnd   = nleaves;
2025:   }
2026:   PetscSectionSetChart(leafSection, lpStart, lpEnd);
2027:   /* Could fuse these at the cost of a copy and extra allocation */
2028:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2029:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2030:   for (f = 0; f < numFields; ++f) {
2031:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2032:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2033:   }
2034:   if (remoteOffsets) {
2035:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2036:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2037:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2038:   }
2039:   PetscSFDestroy(&embedSF);
2040:   PetscSectionSetUp(leafSection);
2041:   return(0);
2042: }

2046: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2047: {
2048:   PetscSF         embedSF;
2049:   const PetscInt *indices;
2050:   IS              selected;
2051:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2052:   PetscErrorCode  ierr;

2055:   *remoteOffsets = NULL;
2056:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2057:   if (numRoots < 0) return(0);
2058:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2059:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2060:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2061:   ISGetIndices(selected, &indices);
2062:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2063:   ISRestoreIndices(selected, &indices);
2064:   ISDestroy(&selected);
2065:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2066:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2067:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2068:   PetscSFDestroy(&embedSF);
2069:   return(0);
2070: }

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

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

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

2086:   Note: Either rootSection or remoteOffsets can be specified

2088:   Level: intermediate

2090: .seealso: PetscSFCreate()
2091: @*/
2092: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2093: {
2094:   MPI_Comm          comm;
2095:   const PetscInt    *localPoints;
2096:   const PetscSFNode *remotePoints;
2097:   PetscInt          lpStart, lpEnd;
2098:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2099:   PetscInt          *localIndices;
2100:   PetscSFNode       *remoteIndices;
2101:   PetscInt          i, ind;
2102:   PetscErrorCode    ierr;

2110:   PetscObjectGetComm((PetscObject)sf,&comm);
2111:   PetscSFCreate(comm, sectionSF);
2112:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2113:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2114:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2115:   if (numRoots < 0) return(0);
2116:   for (i = 0; i < numPoints; ++i) {
2117:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2118:     PetscInt dof;

2120:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2121:       PetscSectionGetDof(leafSection, localPoint, &dof);
2122:       numIndices += dof;
2123:     }
2124:   }
2125:   PetscMalloc1(numIndices, &localIndices);
2126:   PetscMalloc1(numIndices, &remoteIndices);
2127:   /* Create new index graph */
2128:   for (i = 0, ind = 0; i < numPoints; ++i) {
2129:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2130:     PetscInt rank       = remotePoints[i].rank;

2132:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2133:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2134:       PetscInt loff, dof, d;

2136:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2137:       PetscSectionGetDof(leafSection, localPoint, &dof);
2138:       for (d = 0; d < dof; ++d, ++ind) {
2139:         localIndices[ind]        = loff+d;
2140:         remoteIndices[ind].rank  = rank;
2141:         remoteIndices[ind].index = remoteOffset+d;
2142:       }
2143:     }
2144:   }
2145:   PetscFree(remoteOffsets);
2146:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2147:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2148:   return(0);
2149: }

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

2156:   Input Parameters:
2157: + section   - The PetscSection
2158: . obj       - A PetscObject which serves as the key for this index
2159: . clSection - Section giving the size of the closure of each point
2160: - clPoints  - IS giving the points in each closure

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

2164:   Level: intermediate

2166: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2167: @*/
2168: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2169: {

2173:   section->clObj     = obj;
2174:   PetscSectionDestroy(&section->clSection);
2175:   ISDestroy(&section->clPoints);
2176:   section->clSection = clSection;
2177:   section->clPoints  = clPoints;
2178:   return(0);
2179: }

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

2186:   Input Parameters:
2187: + section   - The PetscSection
2188: - obj       - A PetscObject which serves as the key for this index

2190:   Output Parameters:
2191: + clSection - Section giving the size of the closure of each point
2192: - clPoints  - IS giving the points in each closure

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

2196:   Level: intermediate

2198: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2199: @*/
2200: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2201: {
2203:   if (section->clObj == obj) {
2204:     if (clSection) *clSection = section->clSection;
2205:     if (clPoints)  *clPoints  = section->clPoints;
2206:   } else {
2207:     if (clSection) *clSection = NULL;
2208:     if (clPoints)  *clPoints  = NULL;
2209:   }
2210:   return(0);
2211: }

2215: /*@
2216:   PetscSectionGetField - Get the subsection associated with a single field

2218:   Input Parameters:
2219: + s     - The PetscSection
2220: - field - The field number

2222:   Output Parameter:
2223: . subs  - The subsection for the given field

2225:   Level: intermediate

2227: .seealso: PetscSectionSetNumFields()
2228: @*/
2229: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2230: {

2236:   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);
2237:   PetscObjectReference((PetscObject) s->field[field]);
2238:   *subs = s->field[field];
2239:   return(0);
2240: }