Actual source code: vsectionis.c

petsc-3.7.7 2017-09-25
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: - size - the allocated size

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
964:   - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points

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

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

971:   Level: developer

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

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

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

1049: /*@
1050:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1051:   the local section and an SF describing the section point overlap.

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

1060:   Output Parameter:
1061:   . gsection - The PetscSection for the global field layout

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

1065:   Level: developer

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

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

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

1145: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1146: {
1147:   PetscInt       pStart, pEnd, p, localSize = 0;

1151:   PetscSectionGetChart(s, &pStart, &pEnd);
1152:   for (p = pStart; p < pEnd; ++p) {
1153:     PetscInt dof;

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

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

1170:   Input Parameters:
1171: + comm - The MPI_Comm
1172: - s    - The PetscSection

1174:   Output Parameter:
1175: . layout - The layout for the section

1177:   Level: developer

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

1187:   PetscSectionGetChart(s, &pStart, &pEnd);
1188:   for (p = pStart; p < pEnd; ++p) {
1189:     PetscInt dof,cdof;

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

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

1207:   Not collective

1209:   Input Parameters:
1210: + s - the PetscSection
1211: - point - the point

1213:   Output Parameter:
1214: . offset - the offset

1216:   Level: intermediate

1218: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1219: @*/
1220: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1221: {
1223: #if defined(PETSC_USE_DEBUG)
1224:   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);
1225: #endif
1226:   *offset = s->atlasOff[point - s->pStart];
1227:   return(0);
1228: }

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

1235:   Not collective

1237:   Input Parameters:
1238: + s - the PetscSection
1239: . point - the point
1240: - offset - the offset

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

1244:   Level: intermediate

1246: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1247: @*/
1248: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1249: {
1251:   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);
1252:   s->atlasOff[point - s->pStart] = offset;
1253:   return(0);
1254: }

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

1261:   Not collective

1263:   Input Parameters:
1264: + s - the PetscSection
1265: . point - the point
1266: - field - the field

1268:   Output Parameter:
1269: . offset - the offset

1271:   Level: intermediate

1273: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1274: @*/
1275: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1276: {

1280:   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);
1281:   PetscSectionGetOffset(s->field[field], point, offset);
1282:   return(0);
1283: }

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

1290:   Not collective

1292:   Input Parameters:
1293: + s - the PetscSection
1294: . point - the point
1295: . field - the field
1296: - offset - the offset

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

1300:   Level: intermediate

1302: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1303: @*/
1304: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1305: {

1309:   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);
1310:   PetscSectionSetOffset(s->field[field], point, offset);
1311:   return(0);
1312: }

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

1323:   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);
1324:   PetscSectionGetOffset(s, point, &off);
1325:   PetscSectionGetOffset(s->field[field], point, &foff);
1326:   *offset = foff - off;
1327:   return(0);
1328: }

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

1335:   Not collective

1337:   Input Parameter:
1338: . s - the PetscSection

1340:   Output Parameters:
1341: + start - the minimum offset
1342: - end   - one more than the maximum offset

1344:   Level: intermediate

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

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

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

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

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

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

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

1412:     PetscMalloc1(maxCdof, &indices);
1413:     for (p = pStart; p < pEnd; ++p) {
1414:       PetscInt cdof;

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

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

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

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

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

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

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

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

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

1524: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1525: {
1526:   PetscInt       p;
1527:   PetscMPIInt    rank;

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

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

1554: /*@C
1555:   PetscSectionView - Views a PetscSection

1557:   Collective on PetscSection

1559:   Input Parameters:
1560: + s - the PetscSection object to view
1561: - v - the viewer

1563:   Level: developer

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

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

1595: /*@
1596:   PetscSectionReset - Frees all section data.

1598:   Not collective

1600:   Input Parameters:
1601: . s - the PetscSection

1603:   Level: developer

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

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

1627:   s->pStart    = -1;
1628:   s->pEnd      = -1;
1629:   s->maxDof    = 0;
1630:   s->setup     = PETSC_FALSE;
1631:   s->numFields = 0;
1632:   s->clObj     = NULL;
1633:   return(0);
1634: }

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

1641:   Not collective

1643:   Input Parameters:
1644: . s - the PetscSection

1646:   Level: developer

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

1651: .seealso: PetscSection, PetscSectionCreate()
1652: @*/
1653: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1654: {

1658:   if (!*s) return(0);
1660:   if (--((PetscObject)(*s))->refct > 0) {
1661:     *s = NULL;
1662:     return(0);
1663:   }
1664:   PetscSectionReset(*s);
1665:   PetscHeaderDestroy(s);
1666:   return(0);
1667: }

1671: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1672: {
1673:   const PetscInt p = point - s->pStart;

1676:   *values = &baseArray[s->atlasOff[p]];
1677:   return(0);
1678: }

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

1691:   PetscSectionGetConstraintDof(s, p, &cDim);
1692:   array = &baseArray[s->atlasOff[p]];
1693:   if (!cDim) {
1694:     if (orientation >= 0) {
1695:       const PetscInt dim = s->atlasDof[p];
1696:       PetscInt       i;

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

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

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

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

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

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

1759: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1760: {
1764:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1765:   return(0);
1766: }

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

1773:   Input Parameters:
1774: + s     - The PetscSection
1775: - point - The point

1777:   Output Parameter:
1778: . indices - The constrained dofs

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

1782:   Level: advanced

1784: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1785: @*/
1786: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1787: {

1791:   if (s->bc) {
1792:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1793:   } else *indices = NULL;
1794:   return(0);
1795: }

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

1802:   Input Parameters:
1803: + s     - The PetscSection
1804: . point - The point
1805: - indices - The constrained dofs

1807:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1809:   Level: advanced

1811: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1812: @*/
1813: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1814: {

1818:   if (s->bc) {
1819:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1820:   }
1821:   return(0);
1822: }

1826: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1827: {

1831:   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);
1832:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1833:   return(0);
1834: }

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

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

1850: /*@
1851:   PetscSectionPermute - Reorder the section according to the input point permutation

1853:   Collective on PetscSection

1855:   Input Parameter:
1856: + section - The PetscSection object
1857: - perm - The point permutation, old point p becomes new point perm[p]

1859:   Output Parameter:
1860: . sectionNew - The permuted PetscSection

1862:   Level: intermediate

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

1878:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1879:   PetscSectionGetNumFields(s, &numFields);
1880:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1881:   for (f = 0; f < numFields; ++f) {
1882:     const char *name;
1883:     PetscInt    numComp;

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

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

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

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

1937:   Collective

1939:   Input Parameters:
1940: + sf - The SF
1941: - rootSection - Section defined on root space

1943:   Output Parameters:
1944: + remoteOffsets - root offsets in leaf storage, or NULL
1945: - leafSection - Section defined on the leaf space

1947:   Level: intermediate

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

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

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

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

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

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

2044:   Output Parameters:
2045: - sectionSF - The new SF

2047:   Note: Either rootSection or remoteOffsets can be specified

2049:   Level: intermediate

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

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

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

2093:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2094:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2095:       PetscInt loff, dof, d;

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

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

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

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

2124:   Level: intermediate

2126: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2127: @*/
2128: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2129: {

2133:   section->clObj     = obj;
2134:   PetscSectionDestroy(&section->clSection);
2135:   ISDestroy(&section->clPoints);
2136:   section->clSection = clSection;
2137:   section->clPoints  = clPoints;
2138:   return(0);
2139: }

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

2146:   Input Parameters:
2147: + section   - The PetscSection
2148: - obj       - A PetscObject which serves as the key for this index

2150:   Output Parameters:
2151: + clSection - Section giving the size of the closure of each point
2152: - clPoints  - IS giving the points in each closure

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

2156:   Level: intermediate

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

2175: /*@
2176:   PetscSectionGetField - Get the subsection associated with a single field

2178:   Input Parameters:
2179: + s     - The PetscSection
2180: - field - The field number

2182:   Output Parameter:
2183: . subs  - The subsection for the given field

2185:   Level: intermediate

2187: .seealso: PetscSectionSetNumFields()
2188: @*/
2189: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2190: {

2196:   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);
2197:   PetscObjectReference((PetscObject) s->field[field]);
2198:   *subs = s->field[field];
2199:   return(0);
2200: }