Actual source code: vsectionis.c

petsc-3.5.4 2015-05-23
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:     PetscMalloc2(nroots,&neg,nlocal,&recv);
964:     PetscMemzero(neg,nroots*sizeof(PetscInt));
965:   }
966:   /* Mark all local points with negative dof */
967:   for (p = pStart; p < pEnd; ++p) {
968:     PetscSectionGetDof(s, p, &dof);
969:     PetscSectionSetDof(*gsection, p, dof);
970:     PetscSectionGetConstraintDof(s, p, &cdof);
971:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
972:     if (neg) neg[p] = -(dof+1);
973:   }
974:   PetscSectionSetUpBC(*gsection);
975:   if (nroots >= 0) {
976:     PetscMemzero(recv,nlocal*sizeof(PetscInt));
977:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
978:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
979:     for (p = pStart; p < pEnd; ++p) {
980:       if (recv[p] < 0) {
981:         (*gsection)->atlasDof[p-pStart] = recv[p];
982:         PetscSectionGetDof(s, p, &dof);
983:         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);
984:       }
985:     }
986:   }
987:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
988:   if (s->perm) {ISGetIndices(s->perm, &pind);}
989:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
990:     const PetscInt q = pind ? pind[p] : p;

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

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

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

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

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

1035:   Level: developer

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

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

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

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

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

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

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

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

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

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

1162:   Not collective

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

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

1171:   Level: intermediate

1173: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1174: @*/
1175: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1176: {
1178: #if defined(PETSC_USE_DEBUG)
1179:   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);
1180: #endif
1181:   *offset = s->atlasOff[point - s->pStart];
1182:   return(0);
1183: }

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

1190:   Not collective

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

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

1199:   Level: intermediate

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

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

1216:   Not collective

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

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

1226:   Level: intermediate

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

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

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

1245:   Not collective

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

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

1255:   Level: intermediate

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

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

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

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

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

1290:   Not collective

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

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

1299:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1511:   Collective on PetscSection

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

1517:   Level: developer

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

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

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

1551:   Not collective

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

1556:   Level: developer

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

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

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

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

1593:   Not collective

1595:   Input Parameters:
1596: . s - the PetscSection

1598:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

1711: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1712: {
1716:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1717:   return(0);
1718: }

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

1725:   Input Parameters:
1726: + s     - The PetscSection
1727: - point - The point

1729:   Output Parameter:
1730: . indices - The constrained dofs

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

1734:   Level: advanced

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

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

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

1754:   Input Parameters:
1755: + s     - The PetscSection
1756: . point - The point
1757: - indices - The constrained dofs

1759:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

1761:   Level: advanced

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

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

1778: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1779: {

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

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

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

1802: /*@
1803:   PetscSectionPermute - Reorder the section according to the input point permutation

1805:   Collective on PetscSection

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

1811:   Output Parameter:
1812: . sectionNew - The permuted PetscSection

1814:   Level: intermediate

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

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

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

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

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

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

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

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

1905:   PetscObjectGetComm((PetscObject)sfPart,&comm);
1906:   MPI_Comm_rank(comm, &rank);

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

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

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

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

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

1951:   Collective

1953:   Input Parameters:
1954: + sf - The SF
1955: - rootSection - Section defined on root space

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

1961:   Level: intermediate

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

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

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

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

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

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

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

2058:   Note: Either rootSection or remoteOffsets can be specified

2060:   Level: intermediate

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

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

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

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

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

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

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

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

2136:   Level: intermediate

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

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

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

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

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

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

2168:   Level: intermediate

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

2187: /*@
2188:   PetscSectionGetField - Get the subsection associated with a single field

2190:   Input Parameters:
2191: + s     - The PetscSection
2192: - field - The field number

2194:   Output Parameter:
2195: . subs  - The subsection for the given field

2197:   Level: intermediate

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

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