Actual source code: vsectionis.c

petsc-dev 2014-08-27
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:   PetscSectionClone - 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 PetscSectionClone(PetscSection section, PetscSection *newSection)
 84: {
 85:   IS             perm;
 86:   PetscInt       numFields, f, pStart, pEnd, p;

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

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

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

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

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

151:   Not collective

153:   Input Parameter:
154: . s - the PetscSection

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

159:   Level: intermediate

161: .seealso: PetscSectionSetNumFields()
162: @*/
163: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
164: {
167:   *numFields = s->numFields;
168:   return(0);
169: }

173: /*@
174:   PetscSectionSetNumFields - Sets the number of fields.

176:   Not collective

178:   Input Parameters:
179: + s - the PetscSection
180: - numFields - the number of fields

182:   Level: intermediate

184: .seealso: PetscSectionGetNumFields()
185: @*/
186: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
187: {
188:   PetscInt       f;

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

195:   s->numFields = numFields;
196:   PetscMalloc1(s->numFields, &s->numFieldComponents);
197:   PetscMalloc1(s->numFields, &s->fieldNames);
198:   PetscMalloc1(s->numFields, &s->field);
199:   for (f = 0; f < s->numFields; ++f) {
200:     char name[64];

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

204:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
205:     PetscSNPrintf(name, 64, "Field_%D", f);
206:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
207:   }
208:   return(0);
209: }

213: /*@C
214:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

216:   Not Collective

218:   Input Parameters:
219: + s     - the PetscSection
220: - field - the field number

222:   Output Parameter:
223: . fieldName - the field name

225:   Level: developer

227: .seealso: PetscSectionSetFieldName()
228: @*/
229: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
230: {
233:   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);
234:   *fieldName = s->fieldNames[field];
235:   return(0);
236: }

240: /*@C
241:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

243:   Not Collective

245:   Input Parameters:
246: + s     - the PetscSection
247: . field - the field number
248: - fieldName - the field name

250:   Level: developer

252: .seealso: PetscSectionGetFieldName()
253: @*/
254: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
255: {

260:   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);
261:   PetscFree(s->fieldNames[field]);
262:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
263:   return(0);
264: }

268: /*@
269:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

271:   Not collective

273:   Input Parameters:
274: + s - the PetscSection
275: - field - the field number

277:   Output Parameter:
278: . numComp - the number of field components

280:   Level: intermediate

282: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
283: @*/
284: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
285: {
288:   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);
289:   *numComp = s->numFieldComponents[field];
290:   return(0);
291: }

295: /*@
296:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

298:   Not collective

300:   Input Parameters:
301: + s - the PetscSection
302: . field - the field number
303: - numComp - the number of field components

305:   Level: intermediate

307: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
308: @*/
309: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
310: {
312:   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);
313:   s->numFieldComponents[field] = numComp;
314:   return(0);
315: }

319: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
320: {

324:   if (!s->bc) {
325:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
326:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
327:   }
328:   return(0);
329: }

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

336:   Not collective

338:   Input Parameter:
339: . s - the PetscSection

341:   Output Parameters:
342: + pStart - the first point
343: - pEnd - one past the last point

345:   Level: intermediate

347: .seealso: PetscSectionSetChart(), PetscSectionCreate()
348: @*/
349: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
350: {
353:   if (pStart) *pStart = s->pStart;
354:   if (pEnd)   *pEnd   = s->pEnd;
355:   return(0);
356: }

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

363:   Not collective

365:   Input Parameters:
366: + s - the PetscSection
367: . pStart - the first point
368: - pEnd - one past the last point

370:   Level: intermediate

372: .seealso: PetscSectionGetChart(), PetscSectionCreate()
373: @*/
374: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
375: {
376:   PetscInt       f;

381:   /* Cannot Reset() because it destroys field information */
382:   s->setup = PETSC_FALSE;
383:   PetscSectionDestroy(&s->bc);
384:   PetscFree(s->bcIndices);
385:   PetscFree2(s->atlasDof, s->atlasOff);

387:   s->pStart = pStart;
388:   s->pEnd   = pEnd;
389:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
390:   PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
391:   for (f = 0; f < s->numFields; ++f) {
392:     PetscSectionSetChart(s->field[f], pStart, pEnd);
393:   }
394:   return(0);
395: }

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

402:   Not collective

404:   Input Parameter:
405: . s - the PetscSection

407:   Output Parameters:
408: . perm - The permutation as an IS

410:   Level: intermediate

412: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
413: @*/
414: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
415: {
419:   return(0);
420: }

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

427:   Not collective

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

433:   Level: intermediate

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

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

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

456:   Not collective

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

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

465:   Level: intermediate

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

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

484:   Not collective

486:   Input Parameters:
487: + s - the PetscSection
488: . point - the point
489: - numDof - the number of dof

491:   Level: intermediate

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

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

508:   Not collective

510:   Input Parameters:
511: + s - the PetscSection
512: . point - the point
513: - numDof - the number of additional dof

515:   Level: intermediate

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

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

532:   Not collective

534:   Input Parameters:
535: + s - the PetscSection
536: . point - the point
537: - field - the field

539:   Output Parameter:
540: . numDof - the number of dof

542:   Level: intermediate

544: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
545: @*/
546: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
547: {

551:   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);
552:   PetscSectionGetDof(s->field[field], point, numDof);
553:   return(0);
554: }

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

561:   Not collective

563:   Input Parameters:
564: + s - the PetscSection
565: . point - the point
566: . field - the field
567: - numDof - the number of dof

569:   Level: intermediate

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

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

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

588:   Not collective

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

596:   Level: intermediate

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

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

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

615:   Not collective

617:   Input Parameters:
618: + s - the PetscSection
619: - point - the point

621:   Output Parameter:
622: . numDof - the number of dof which are fixed by constraints

624:   Level: intermediate

626: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
627: @*/
628: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
629: {

633:   if (s->bc) {
634:     PetscSectionGetDof(s->bc, point, numDof);
635:   } else *numDof = 0;
636:   return(0);
637: }

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

644:   Not collective

646:   Input Parameters:
647: + s - the PetscSection
648: . point - the point
649: - numDof - the number of dof which are fixed by constraints

651:   Level: intermediate

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

660:   if (numDof) {
661:     PetscSectionCheckConstraints_Static(s);
662:     PetscSectionSetDof(s->bc, point, numDof);
663:   }
664:   return(0);
665: }

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

672:   Not collective

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

679:   Level: intermediate

681: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
682: @*/
683: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
684: {

688:   if (numDof) {
689:     PetscSectionCheckConstraints_Static(s);
690:     PetscSectionAddDof(s->bc, point, numDof);
691:   }
692:   return(0);
693: }

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

700:   Not collective

702:   Input Parameters:
703: + s - the PetscSection
704: . point - the point
705: - field - the field

707:   Output Parameter:
708: . numDof - the number of dof which are fixed by constraints

710:   Level: intermediate

712: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
713: @*/
714: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
715: {

719:   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);
720:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
721:   return(0);
722: }

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

729:   Not collective

731:   Input Parameters:
732: + s - the PetscSection
733: . point - the point
734: . field - the field
735: - numDof - the number of dof which are fixed by constraints

737:   Level: intermediate

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

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

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

756:   Not collective

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

764:   Level: intermediate

766: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
767: @*/
768: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
769: {

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

780: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
781: {

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

788:     PetscSectionSetUp(s->bc);
789:     PetscMalloc1((s->bc->atlasOff[last] + s->bc->atlasDof[last]), &s->bcIndices);
790:   }
791:   return(0);
792: }

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

799:   Not collective

801:   Input Parameter:
802: . s - the PetscSection

804:   Level: intermediate

806: .seealso: PetscSectionCreate()
807: @*/
808: PetscErrorCode PetscSectionSetUp(PetscSection s)
809: {
810:   const PetscInt *pind   = NULL;
811:   PetscInt        offset = 0, p, f;
812:   PetscErrorCode  ierr;

815:   if (s->setup) return(0);
816:   s->setup = PETSC_TRUE;
817:   if (s->perm) {ISGetIndices(s->perm, &pind);}
818:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
819:     const PetscInt q = pind ? pind[p] : p;

821:     s->atlasOff[q] = offset;
822:     offset        += s->atlasDof[q];
823:     s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
824:   }
825:   PetscSectionSetUpBC(s);
826:   /* Assume that all fields have the same chart */
827:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
828:     const PetscInt q   = pind ? pind[p] : p;
829:     PetscInt       off = s->atlasOff[q];

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

834:       sf->atlasOff[q] = off;
835:       off += sf->atlasDof[q];
836:     }
837:   }
838:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
839:   for (f = 0; f < s->numFields; ++f) {
840:     PetscSectionSetUpBC(s->field[f]);
841:   }
842:   return(0);
843: }

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

850:   Not collective

852:   Input Parameters:
853: . s - the PetscSection

855:   Output Parameter:
856: . maxDof - the maximum dof

858:   Level: intermediate

860: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
861: @*/
862: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
863: {
865:   *maxDof = s->maxDof;
866:   return(0);
867: }

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

874:   Not collective

876:   Input Parameters:
877: + s - the PetscSection
878: - point - the point

880:   Output Parameter:
881: . size - the size of an array which can hold all the dofs

883:   Level: intermediate

885: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
886: @*/
887: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
888: {
889:   PetscInt p, n = 0;

892:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
893:   *size = n;
894:   return(0);
895: }

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

902:   Not collective

904:   Input Parameters:
905: + s - the PetscSection
906: - point - the point

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

911:   Level: intermediate

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

920:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
921:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
922:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
923:   }
924:   *size = n;
925:   return(0);
926: }

930: /*@
931:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
932:   the local section and an SF describing the section point overlap.

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

939:   Output Parameter:
940:   . gsection - The PetscSection for the global field layout

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

944:   Level: developer

946: .seealso: PetscSectionCreate()
947: @*/
948: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscSection *gsection)
949: {
950:   const PetscInt *pind = NULL;
951:   PetscInt       *recv = NULL, *neg = NULL;
952:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
953:   PetscErrorCode  ierr;

956:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
957:   PetscSectionGetChart(s, &pStart, &pEnd);
958:   PetscSectionSetChart(*gsection, pStart, pEnd);
959:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
960:   nlocal = nroots;              /* The local/leaf space matches global/root space */
961:   /* Must allocate for all points visible to SF, which may be more than this section */
962:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
963:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
964:     PetscMalloc2(nroots,&neg,nlocal,&recv);
965:     PetscMemzero(neg,nroots*sizeof(PetscInt));
966:   }
967:   /* Mark all local points with negative dof */
968:   for (p = pStart; p < pEnd; ++p) {
969:     PetscSectionGetDof(s, p, &dof);
970:     PetscSectionSetDof(*gsection, p, dof);
971:     PetscSectionGetConstraintDof(s, p, &cdof);
972:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
973:     if (neg) neg[p] = -(dof+1);
974:   }
975:   PetscSectionSetUpBC(*gsection);
976:   if (nroots >= 0) {
977:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
978:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
979:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
980:     for (p = pStart; p < pEnd; ++p) {
981:       if (recv[p] < 0) {
982:         (*gsection)->atlasDof[p-pStart] = recv[p];
983:         PetscSectionGetDof(s, p, &dof);
984:         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);
985:       }
986:     }
987:   }
988:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
989:   if (s->perm) {ISGetIndices(s->perm, &pind);}
990:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
991:     const PetscInt q = pind ? pind[p] : p;

993:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
994:     (*gsection)->atlasOff[q] = off;
995:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
996:   }
997:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
998:   globalOff -= off;
999:   for (p = pStart, off = 0; p < pEnd; ++p) {
1000:     (*gsection)->atlasOff[p-pStart] += globalOff;
1001:     if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
1002:   }
1003:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1004:   /* Put in negative offsets for ghost points */
1005:   if (nroots >= 0) {
1006:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
1007:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1008:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1009:     for (p = pStart; p < pEnd; ++p) {
1010:       if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
1011:     }
1012:   }
1013:   PetscFree2(neg,recv);
1014:   PetscSectionViewFromOptions(*gsection,NULL,"-section_global_view");
1015:   return(0);
1016: }

1020: /*@
1021:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1022:   the local section and an SF describing the section point overlap.

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

1031:   Output Parameter:
1032:   . gsection - The PetscSection for the global field layout

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

1036:   Level: developer

1038: .seealso: PetscSectionCreate()
1039: @*/
1040: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1041: {
1042:   const PetscInt *pind = NULL;
1043:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1044:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1045:   PetscErrorCode  ierr;

1048:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1049:   PetscSectionGetChart(s, &pStart, &pEnd);
1050:   PetscSectionSetChart(*gsection, pStart, pEnd);
1051:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1052:   if (nroots >= 0) {
1053:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1054:     PetscCalloc1(nroots, &neg);
1055:     if (nroots > pEnd-pStart) {
1056:       PetscCalloc1(nroots, &tmpOff);
1057:     } else {
1058:       tmpOff = &(*gsection)->atlasDof[-pStart];
1059:     }
1060:   }
1061:   /* Mark ghost points with negative dof */
1062:   for (p = pStart; p < pEnd; ++p) {
1063:     for (e = 0; e < numExcludes; ++e) {
1064:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1065:         PetscSectionSetDof(*gsection, p, 0);
1066:         break;
1067:       }
1068:     }
1069:     if (e < numExcludes) continue;
1070:     PetscSectionGetDof(s, p, &dof);
1071:     PetscSectionSetDof(*gsection, p, dof);
1072:     PetscSectionGetConstraintDof(s, p, &cdof);
1073:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1074:     if (neg) neg[p] = -(dof+1);
1075:   }
1076:   PetscSectionSetUpBC(*gsection);
1077:   if (nroots >= 0) {
1078:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1079:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1080:     if (nroots > pEnd - pStart) {
1081:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1082:     }
1083:   }
1084:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1085:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1086:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1087:     const PetscInt q = pind ? pind[p] : p;

1089:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1090:     (*gsection)->atlasOff[q] = off;
1091:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1092:   }
1093:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1094:   globalOff -= off;
1095:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1096:     (*gsection)->atlasOff[p] += globalOff;
1097:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1098:   }
1099:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1100:   /* Put in negative offsets for ghost points */
1101:   if (nroots >= 0) {
1102:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1103:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1104:     if (nroots > pEnd - pStart) {
1105:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1106:     }
1107:   }
1108:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1109:   PetscFree(neg);
1110:   return(0);
1111: }

1115: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1116: {
1117:   PetscInt       pStart, pEnd, p, localSize = 0;

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

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

1137: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1138: {
1139:   PetscInt       pStart, pEnd, p, localSize = 0;

1143:   PetscSectionGetChart(s, &pStart, &pEnd);
1144:   for (p = pStart; p < pEnd; ++p) {
1145:     PetscInt dof,cdof;

1147:     PetscSectionGetDof(s, p, &dof);
1148:     PetscSectionGetConstraintDof(s, p, &cdof);
1149:     if (dof-cdof > 0) localSize += dof-cdof;
1150:   }
1151:   PetscLayoutCreate(comm, layout);
1152:   PetscLayoutSetLocalSize(*layout, localSize);
1153:   PetscLayoutSetBlockSize(*layout, 1);
1154:   PetscLayoutSetUp(*layout);
1155:   return(0);
1156: }

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

1163:   Not collective

1165:   Input Parameters:
1166: + s - the PetscSection
1167: - point - the point

1169:   Output Parameter:
1170: . offset - the offset

1172:   Level: intermediate

1174: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1175: @*/
1176: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1177: {
1179: #if defined(PETSC_USE_DEBUG)
1180:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1181: #endif
1182:   *offset = s->atlasOff[point - s->pStart];
1183:   return(0);
1184: }

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

1191:   Not collective

1193:   Input Parameters:
1194: + s - the PetscSection
1195: . point - the point
1196: - offset - the offset

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

1200:   Level: intermediate

1202: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1203: @*/
1204: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1205: {
1207:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1208:   s->atlasOff[point - s->pStart] = offset;
1209:   return(0);
1210: }

1214: /*@
1215:   PetscSectionGetFieldOffset - Return 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: - field - the field

1224:   Output Parameter:
1225: . offset - the offset

1227:   Level: intermediate

1229: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1230: @*/
1231: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1232: {

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

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

1246:   Not collective

1248:   Input Parameters:
1249: + s - the PetscSection
1250: . point - the point
1251: . field - the field
1252: - offset - the offset

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

1256:   Level: intermediate

1258: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1259: @*/
1260: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1261: {

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

1272: /* This gives the offset on a point of the field, ignoring constraints */
1273: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1274: {
1275:   PetscInt       off, foff;

1279:   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);
1280:   PetscSectionGetOffset(s, point, &off);
1281:   PetscSectionGetOffset(s->field[field], point, &foff);
1282:   *offset = foff - off;
1283:   return(0);
1284: }

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

1291:   Not collective

1293:   Input Parameter:
1294: . s - the PetscSection

1296:   Output Parameters:
1297: + start - the minimum offset
1298: - end   - one more than the maximum offset

1300:   Level: intermediate

1302: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1303: @*/
1304: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1305: {
1306:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1310:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1311:   PetscSectionGetChart(s, &pStart, &pEnd);
1312:   for (p = 0; p < pEnd-pStart; ++p) {
1313:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1315:     if (off >= 0) {
1316:       os = PetscMin(os, off);
1317:       oe = PetscMax(oe, off+dof);
1318:     }
1319:   }
1320:   if (start) *start = os;
1321:   if (end)   *end   = oe;
1322:   return(0);
1323: }

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

1333:   if (!numFields) return(0);
1334:   PetscSectionGetNumFields(s, &nF);
1335:   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1336:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1337:   PetscSectionSetNumFields(*subs, numFields);
1338:   for (f = 0; f < numFields; ++f) {
1339:     const char *name   = NULL;
1340:     PetscInt   numComp = 0;

1342:     PetscSectionGetFieldName(s, fields[f], &name);
1343:     PetscSectionSetFieldName(*subs, f, name);
1344:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1345:     PetscSectionSetFieldComponents(*subs, f, numComp);
1346:   }
1347:   PetscSectionGetChart(s, &pStart, &pEnd);
1348:   PetscSectionSetChart(*subs, pStart, pEnd);
1349:   for (p = pStart; p < pEnd; ++p) {
1350:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1352:     for (f = 0; f < numFields; ++f) {
1353:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1354:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1355:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1356:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1357:       dof  += fdof;
1358:       cdof += cfdof;
1359:     }
1360:     PetscSectionSetDof(*subs, p, dof);
1361:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1362:     maxCdof = PetscMax(cdof, maxCdof);
1363:   }
1364:   PetscSectionSetUp(*subs);
1365:   if (maxCdof) {
1366:     PetscInt *indices;

1368:     PetscMalloc1(maxCdof, &indices);
1369:     for (p = pStart; p < pEnd; ++p) {
1370:       PetscInt cdof;

1372:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1373:       if (cdof) {
1374:         const PetscInt *oldIndices = NULL;
1375:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

1380:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1381:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1382:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1383:           /* This can be sped up if we assume sorted fields */
1384:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1385:             PetscInt oldfdof = 0;
1386:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1387:             oldFoff += oldfdof;
1388:           }
1389:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1390:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1391:           numConst += cfdof;
1392:           fOff     += fdof;
1393:         }
1394:         PetscSectionSetConstraintIndices(*subs, p, indices);
1395:       }
1396:     }
1397:     PetscFree(indices);
1398:   }
1399:   return(0);
1400: }

1404: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1405: {
1406:   const PetscInt *points = NULL, *indices = NULL;
1407:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

1411:   PetscSectionGetNumFields(s, &numFields);
1412:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1413:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1414:   for (f = 0; f < numFields; ++f) {
1415:     const char *name   = NULL;
1416:     PetscInt   numComp = 0;

1418:     PetscSectionGetFieldName(s, f, &name);
1419:     PetscSectionSetFieldName(*subs, f, name);
1420:     PetscSectionGetFieldComponents(s, f, &numComp);
1421:     PetscSectionSetFieldComponents(*subs, f, numComp);
1422:   }
1423:   /* For right now, we do not try to squeeze the subchart */
1424:   if (subpointMap) {
1425:     ISGetSize(subpointMap, &numSubpoints);
1426:     ISGetIndices(subpointMap, &points);
1427:   }
1428:   PetscSectionGetChart(s, &pStart, &pEnd);
1429:   PetscSectionSetChart(*subs, 0, numSubpoints);
1430:   for (p = pStart; p < pEnd; ++p) {
1431:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1433:     PetscFindInt(p, numSubpoints, points, &subp);
1434:     if (subp < 0) continue;
1435:     for (f = 0; f < numFields; ++f) {
1436:       PetscSectionGetFieldDof(s, p, f, &fdof);
1437:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1438:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1439:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1440:     }
1441:     PetscSectionGetDof(s, p, &dof);
1442:     PetscSectionSetDof(*subs, subp, dof);
1443:     PetscSectionGetConstraintDof(s, p, &cdof);
1444:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1445:   }
1446:   PetscSectionSetUp(*subs);
1447:   /* Change offsets to original offsets */
1448:   for (p = pStart; p < pEnd; ++p) {
1449:     PetscInt off, foff = 0;

1451:     PetscFindInt(p, numSubpoints, points, &subp);
1452:     if (subp < 0) continue;
1453:     for (f = 0; f < numFields; ++f) {
1454:       PetscSectionGetFieldOffset(s, p, f, &foff);
1455:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1456:     }
1457:     PetscSectionGetOffset(s, p, &off);
1458:     PetscSectionSetOffset(*subs, subp, off);
1459:   }
1460:   /* Copy constraint indices */
1461:   for (subp = 0; subp < numSubpoints; ++subp) {
1462:     PetscInt cdof;

1464:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1465:     if (cdof) {
1466:       for (f = 0; f < numFields; ++f) {
1467:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1468:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1469:       }
1470:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
1471:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1472:     }
1473:   }
1474:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1475:   return(0);
1476: }

1480: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1481: {
1482:   PetscInt       p;
1483:   PetscMPIInt    rank;

1487:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1488:   PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
1489:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1490:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1491:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1492:       PetscInt b;

1494:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1495:       for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1496:         PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1497:       }
1498:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1499:     } else {
1500:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1501:     }
1502:   }
1503:   PetscViewerFlush(viewer);
1504:   return(0);
1505: }

1509: /*@C
1510:   PetscSectionView - Views a PetscSection

1512:   Collective on PetscSection

1514:   Input Parameters:
1515: + s - the PetscSection object to view
1516: - v - the viewer

1518:   Level: developer

1520: .seealso PetscSectionCreate(), PetscSectionDestroy()
1521: @*/
1522: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1523: {
1524:   PetscBool      isascii;
1525:   PetscInt       f;

1529:   if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
1531:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1532:   if (isascii) {
1533:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1534:     if (s->numFields) {
1535:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1536:       for (f = 0; f < s->numFields; ++f) {
1537:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
1538:         PetscSectionView_ASCII(s->field[f], viewer);
1539:       }
1540:     } else {
1541:       PetscSectionView_ASCII(s, viewer);
1542:     }
1543:   }
1544:   return(0);
1545: }

1549: /*@
1550:   PetscSectionReset - Frees all section data.

1552:   Not collective

1554:   Input Parameters:
1555: . s - the PetscSection

1557:   Level: developer

1559: .seealso: PetscSection, PetscSectionCreate()
1560: @*/
1561: PetscErrorCode PetscSectionReset(PetscSection s)
1562: {
1563:   PetscInt       f;

1567:   PetscFree(s->numFieldComponents);
1568:   for (f = 0; f < s->numFields; ++f) {
1569:     PetscSectionDestroy(&s->field[f]);
1570:     PetscFree(s->fieldNames[f]);
1571:   }
1572:   PetscFree(s->fieldNames);
1573:   PetscFree(s->field);
1574:   PetscSectionDestroy(&s->bc);
1575:   PetscFree(s->bcIndices);
1576:   PetscFree2(s->atlasDof, s->atlasOff);
1577:   PetscSectionDestroy(&s->clSection);
1578:   ISDestroy(&s->clPoints);
1579:   ISDestroy(&s->perm);

1581:   s->pStart    = -1;
1582:   s->pEnd      = -1;
1583:   s->maxDof    = 0;
1584:   s->setup     = PETSC_FALSE;
1585:   s->numFields = 0;
1586:   s->clObj     = NULL;
1587:   return(0);
1588: }

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

1595:   Not collective

1597:   Input Parameters:
1598: . s - the PetscSection

1600:   Level: developer

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

1605: .seealso: PetscSection, PetscSectionCreate()
1606: @*/
1607: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1608: {

1612:   if (!*s) return(0);
1614:   if (--((PetscObject)(*s))->refct > 0) {
1615:     *s = NULL;
1616:     return(0);
1617:   }
1618:   PetscSectionReset(*s);
1619:   PetscHeaderDestroy(s);
1620:   return(0);
1621: }

1625: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1626: {
1627:   const PetscInt p = point - s->pStart;

1630:   *values = &baseArray[s->atlasOff[p]];
1631:   return(0);
1632: }

1636: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1637: {
1638:   PetscInt       *array;
1639:   const PetscInt p           = point - s->pStart;
1640:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1641:   PetscInt       cDim        = 0;

1645:   PetscSectionGetConstraintDof(s, p, &cDim);
1646:   array = &baseArray[s->atlasOff[p]];
1647:   if (!cDim) {
1648:     if (orientation >= 0) {
1649:       const PetscInt dim = s->atlasDof[p];
1650:       PetscInt       i;

1652:       if (mode == INSERT_VALUES) {
1653:         for (i = 0; i < dim; ++i) array[i] = values[i];
1654:       } else {
1655:         for (i = 0; i < dim; ++i) array[i] += values[i];
1656:       }
1657:     } else {
1658:       PetscInt offset = 0;
1659:       PetscInt j      = -1, field, i;

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

1664:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1665:         offset += dim;
1666:       }
1667:     }
1668:   } else {
1669:     if (orientation >= 0) {
1670:       const PetscInt dim  = s->atlasDof[p];
1671:       PetscInt       cInd = 0, i;
1672:       const PetscInt *cDof;

1674:       PetscSectionGetConstraintIndices(s, point, &cDof);
1675:       if (mode == INSERT_VALUES) {
1676:         for (i = 0; i < dim; ++i) {
1677:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1678:           array[i] = values[i];
1679:         }
1680:       } else {
1681:         for (i = 0; i < dim; ++i) {
1682:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1683:           array[i] += values[i];
1684:         }
1685:       }
1686:     } else {
1687:       const PetscInt *cDof;
1688:       PetscInt       offset  = 0;
1689:       PetscInt       cOffset = 0;
1690:       PetscInt       j       = 0, field;

1692:       PetscSectionGetConstraintIndices(s, point, &cDof);
1693:       for (field = 0; field < s->numFields; ++field) {
1694:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
1695:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1696:         const PetscInt sDim = dim - tDim;
1697:         PetscInt       cInd = 0, i,k;

1699:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1700:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1701:           array[j] = values[k];
1702:         }
1703:         offset  += dim;
1704:         cOffset += dim - tDim;
1705:       }
1706:     }
1707:   }
1708:   return(0);
1709: }

1713: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1714: {
1718:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1719:   return(0);
1720: }

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

1727:   Input Parameters:
1728: + s     - The PetscSection
1729: - point - The point

1731:   Output Parameter:
1732: . indices - The constrained dofs

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

1736:   Level: advanced

1738: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1739: @*/
1740: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1741: {

1745:   if (s->bc) {
1746:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1747:   } else *indices = NULL;
1748:   return(0);
1749: }

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

1756:   Input Parameters:
1757: + s     - The PetscSection
1758: . point - The point
1759: - indices - The constrained dofs

1761:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1763:   Level: advanced

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

1772:   if (s->bc) {
1773:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1774:   }
1775:   return(0);
1776: }

1780: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1781: {

1785:   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);
1786:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1787:   return(0);
1788: }

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

1797:   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);
1798:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1799:   return(0);
1800: }

1804: /*@
1805:   PetscSectionPermute - Reorder the section according to the input point permutation

1807:   Collective on PetscSection

1809:   Input Parameter:
1810: + section - The PetscSection object
1811: - perm - The point permutation, old point p becomes new point perm[p]

1813:   Output Parameter:
1814: . sectionNew - The permuted PetscSection

1816:   Level: intermediate

1818: .keywords: mesh
1819: .seealso: MatPermute()
1820: @*/
1821: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1822: {
1823:   PetscSection    s = section, sNew;
1824:   const PetscInt *perm;
1825:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1826:   PetscErrorCode  ierr;

1832:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1833:   PetscSectionGetNumFields(s, &numFields);
1834:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1835:   for (f = 0; f < numFields; ++f) {
1836:     const char *name;
1837:     PetscInt    numComp;

1839:     PetscSectionGetFieldName(s, f, &name);
1840:     PetscSectionSetFieldName(sNew, f, name);
1841:     PetscSectionGetFieldComponents(s, f, &numComp);
1842:     PetscSectionSetFieldComponents(sNew, f, numComp);
1843:   }
1844:   ISGetLocalSize(permutation, &numPoints);
1845:   ISGetIndices(permutation, &perm);
1846:   PetscSectionGetChart(s, &pStart, &pEnd);
1847:   PetscSectionSetChart(sNew, pStart, pEnd);
1848:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1849:   for (p = pStart; p < pEnd; ++p) {
1850:     PetscInt dof, cdof;

1852:     PetscSectionGetDof(s, p, &dof);
1853:     PetscSectionSetDof(sNew, perm[p], dof);
1854:     PetscSectionGetConstraintDof(s, p, &cdof);
1855:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1856:     for (f = 0; f < numFields; ++f) {
1857:       PetscSectionGetFieldDof(s, p, f, &dof);
1858:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1859:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1860:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1861:     }
1862:   }
1863:   PetscSectionSetUp(sNew);
1864:   for (p = pStart; p < pEnd; ++p) {
1865:     const PetscInt *cind;
1866:     PetscInt        cdof;

1868:     PetscSectionGetConstraintDof(s, p, &cdof);
1869:     if (cdof) {
1870:       PetscSectionGetConstraintIndices(s, p, &cind);
1871:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1872:     }
1873:     for (f = 0; f < numFields; ++f) {
1874:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1875:       if (cdof) {
1876:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1877:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1878:       }
1879:     }
1880:   }
1881:   ISRestoreIndices(permutation, &perm);
1882:   *sectionNew = sNew;
1883:   return(0);
1884: }

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

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

1896: PetscErrorCode PetscSFConvertPartition(PetscSF sfPart, PetscSection partSection, IS partition, ISLocalToGlobalMapping *renumbering, PetscSF *sf)
1897: {
1898:   MPI_Comm       comm;
1899:   PetscSF        sfPoints;
1900:   PetscInt       *partSizes,*partOffsets,p,i,numParts,numMyPoints,numPoints,count;
1901:   const PetscInt *partArray;
1902:   PetscSFNode    *sendPoints;
1903:   PetscMPIInt    rank;

1907:   PetscObjectGetComm((PetscObject)sfPart,&comm);
1908:   MPI_Comm_rank(comm, &rank);

1910:   /* Get the number of parts and sizes that I have to distribute */
1911:   PetscSFGetGraph(sfPart,NULL,&numParts,NULL,NULL);
1912:   PetscMalloc2(numParts,&partSizes,numParts,&partOffsets);
1913:   for (p=0,numPoints=0; p<numParts; p++) {
1914:     PetscSectionGetDof(partSection, p, &partSizes[p]);
1915:     numPoints += partSizes[p];
1916:   }
1917:   numMyPoints = 0;
1918:   PetscSFFetchAndOpBegin(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1919:   PetscSFFetchAndOpEnd(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1920:   /* I will receive numMyPoints. I will send a total of numPoints, to be placed on remote procs at partOffsets */

1922:   /* Create SF mapping locations (addressed through partition, as indexed leaves) to new owners (roots) */
1923:   PetscMalloc1(numPoints,&sendPoints);
1924:   for (p=0,count=0; p<numParts; p++) {
1925:     for (i=0; i<partSizes[p]; i++) {
1926:       sendPoints[count].rank = p;
1927:       sendPoints[count].index = partOffsets[p]+i;
1928:       count++;
1929:     }
1930:   }
1931:   if (count != numPoints) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count %D should equal numPoints=%D",count,numPoints);
1932:   PetscFree2(partSizes,partOffsets);
1933:   ISGetIndices(partition,&partArray);
1934:   PetscSFCreate(comm,&sfPoints);
1935:   PetscSFSetFromOptions(sfPoints);
1936:   PetscSFSetGraph(sfPoints,numMyPoints,numPoints,partArray,PETSC_USE_POINTER,sendPoints,PETSC_OWN_POINTER);

1938:   /* Invert SF so that the new owners are leaves and the locations indexed through partition are the roots */
1939:   PetscSFCreateInverseSF(sfPoints,sf);
1940:   PetscSFDestroy(&sfPoints);
1941:   ISRestoreIndices(partition,&partArray);

1943:   /* Create the new local-to-global mapping */
1944:   ISLocalToGlobalMappingCreateSF(*sf,0,renumbering);
1945:   return(0);
1946: }

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

1953:   Collective

1955:   Input Parameters:
1956: + sf - The SF
1957: - rootSection - Section defined on root space

1959:   Output Parameters:
1960: + remoteOffsets - root offsets in leaf storage, or NULL
1961: - leafSection - Section defined on the leaf space

1963:   Level: intermediate

1965: .seealso: PetscSFCreate()
1966: @*/
1967: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1968: {
1969:   PetscSF        embedSF;
1970:   const PetscInt *ilocal, *indices;
1971:   IS             selected;
1972:   PetscInt       numFields, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;

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

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

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

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

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

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

2060:   Note: Either rootSection or remoteOffsets can be specified

2062:   Level: intermediate

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

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

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

2106:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2107:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2108:       PetscInt loff, dof, d;

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

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

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

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

2138:   Level: intermediate

2140: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2141: @*/
2142: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2143: {

2147:   section->clObj     = obj;
2148:   PetscSectionDestroy(&section->clSection);
2149:   ISDestroy(&section->clPoints);
2150:   section->clSection = clSection;
2151:   section->clPoints  = clPoints;
2152:   return(0);
2153: }

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

2160:   Input Parameters:
2161: + section   - The PetscSection
2162: - obj       - A PetscObject which serves as the key for this index

2164:   Output Parameters:
2165: + clSection - Section giving the size of the closure of each point
2166: - clPoints  - IS giving the points in each closure

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

2170:   Level: intermediate

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

2189: /*@
2190:   PetscSectionGetField - Get the subsection associated with a single field

2192:   Input Parameters:
2193: + s     - The PetscSection
2194: - field - The field number

2196:   Output Parameter:
2197: . subs  - The subsection for the given field

2199:   Level: intermediate

2201: .seealso: PetscSectionSetNumFields()
2202: @*/
2203: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2204: {

2210:   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);
2211:   PetscObjectReference((PetscObject) s->field[field]);
2212:   *subs = s->field[field];
2213:   return(0);
2214: }