Actual source code: vsectionis.c

petsc-dev 2014-04-16
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) (*gsection)->atlasDof[p-pStart] = recv[p];
981:     }
982:   }
983:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
984:   if (s->perm) {ISGetIndices(s->perm, &pind);}
985:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
986:     const PetscInt q = pind ? pind[p] : p;

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

1015: /*@
1016:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1017:   the local section and an SF describing the section point overlap.

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

1026:   Output Parameter:
1027:   . gsection - The PetscSection for the global field layout

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

1031:   Level: developer

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

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

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

1110: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1111: {
1112:   PetscInt       pStart, pEnd, p, localSize = 0;

1116:   PetscSectionGetChart(s, &pStart, &pEnd);
1117:   for (p = pStart; p < pEnd; ++p) {
1118:     PetscInt dof;

1120:     PetscSectionGetDof(s, p, &dof);
1121:     if (dof > 0) ++localSize;
1122:   }
1123:   PetscLayoutCreate(comm, layout);
1124:   PetscLayoutSetLocalSize(*layout, localSize);
1125:   PetscLayoutSetBlockSize(*layout, 1);
1126:   PetscLayoutSetUp(*layout);
1127:   return(0);
1128: }

1132: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1133: {
1134:   PetscInt       pStart, pEnd, p, localSize = 0;

1138:   PetscSectionGetChart(s, &pStart, &pEnd);
1139:   for (p = pStart; p < pEnd; ++p) {
1140:     PetscInt dof,cdof;

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

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

1158:   Not collective

1160:   Input Parameters:
1161: + s - the PetscSection
1162: - point - the point

1164:   Output Parameter:
1165: . offset - the offset

1167:   Level: intermediate

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

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

1186:   Not collective

1188:   Input Parameters:
1189: + s - the PetscSection
1190: . point - the point
1191: - offset - the offset

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

1195:   Level: intermediate

1197: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1198: @*/
1199: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1200: {
1202:   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);
1203:   s->atlasOff[point - s->pStart] = offset;
1204:   return(0);
1205: }

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

1212:   Not collective

1214:   Input Parameters:
1215: + s - the PetscSection
1216: . point - the point
1217: - field - the field

1219:   Output Parameter:
1220: . offset - the offset

1222:   Level: intermediate

1224: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1225: @*/
1226: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1227: {

1231:   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);
1232:   PetscSectionGetOffset(s->field[field], point, offset);
1233:   return(0);
1234: }

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

1241:   Not collective

1243:   Input Parameters:
1244: + s - the PetscSection
1245: . point - the point
1246: . field - the field
1247: - offset - the offset

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

1251:   Level: intermediate

1253: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1254: @*/
1255: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1256: {

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

1267: /* This gives the offset on a point of the field, ignoring constraints */
1268: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1269: {
1270:   PetscInt       off, foff;

1274:   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);
1275:   PetscSectionGetOffset(s, point, &off);
1276:   PetscSectionGetOffset(s->field[field], point, &foff);
1277:   *offset = foff - off;
1278:   return(0);
1279: }

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

1286:   Not collective

1288:   Input Parameter:
1289: . s - the PetscSection

1291:   Output Parameters:
1292: + start - the minimum offset
1293: - end   - one more than the maximum offset

1295:   Level: intermediate

1297: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1298: @*/
1299: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1300: {
1301:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1305:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1306:   PetscSectionGetChart(s, &pStart, &pEnd);
1307:   for (p = 0; p < pEnd-pStart; ++p) {
1308:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1310:     if (off >= 0) {
1311:       os = PetscMin(os, off);
1312:       oe = PetscMax(oe, off+dof);
1313:     }
1314:   }
1315:   if (start) *start = os;
1316:   if (end)   *end   = oe;
1317:   return(0);
1318: }

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

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

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

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

1363:     PetscMalloc1(maxCdof, &indices);
1364:     for (p = pStart; p < pEnd; ++p) {
1365:       PetscInt cdof;

1367:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1368:       if (cdof) {
1369:         const PetscInt *oldIndices = NULL;
1370:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

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

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

1399: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1400: {
1401:   const PetscInt *points = NULL, *indices = NULL;
1402:   PetscInt       numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;

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

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

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

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

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

1475: PetscErrorCode  PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1476: {
1477:   PetscInt       p;
1478:   PetscMPIInt    rank;

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

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

1504: /*@C
1505:   PetscSectionView - Views a PetscSection

1507:   Collective on PetscSection

1509:   Input Parameters:
1510: + s - the PetscSection object to view
1511: - v - the viewer

1513:   Level: developer

1515: .seealso PetscSectionCreate(), PetscSectionDestroy()
1516: @*/
1517: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1518: {
1519:   PetscBool      isascii;
1520:   PetscInt       f;

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

1544: /*@
1545:   PetscSectionReset - Frees all section data.

1547:   Not collective

1549:   Input Parameters:
1550: . s - the PetscSection

1552:   Level: developer

1554: .seealso: PetscSection, PetscSectionCreate()
1555: @*/
1556: PetscErrorCode PetscSectionReset(PetscSection s)
1557: {
1558:   PetscInt       f;

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

1575:   s->pStart    = -1;
1576:   s->pEnd      = -1;
1577:   s->maxDof    = 0;
1578:   s->setup     = PETSC_FALSE;
1579:   s->numFields = 0;
1580:   s->clObj     = NULL;
1581:   return(0);
1582: }

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

1589:   Not collective

1591:   Input Parameters:
1592: . s - the PetscSection

1594:   Level: developer

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

1599: .seealso: PetscSection, PetscSectionCreate()
1600: @*/
1601: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1602: {

1606:   if (!*s) return(0);
1608:   if (--((PetscObject)(*s))->refct > 0) {
1609:     *s = NULL;
1610:     return(0);
1611:   }
1612:   PetscSectionReset(*s);
1613:   PetscHeaderDestroy(s);
1614:   return(0);
1615: }

1619: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1620: {
1621:   const PetscInt p = point - s->pStart;

1624:   *values = &baseArray[s->atlasOff[p]];
1625:   return(0);
1626: }

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

1639:   PetscSectionGetConstraintDof(s, p, &cDim);
1640:   array = &baseArray[s->atlasOff[p]];
1641:   if (!cDim) {
1642:     if (orientation >= 0) {
1643:       const PetscInt dim = s->atlasDof[p];
1644:       PetscInt       i;

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

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

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

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

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

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

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

1718: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1719: {

1723:   if (s->bc) {
1724:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1725:   } else *indices = NULL;
1726:   return(0);
1727: }

1731: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1732: {

1736:   if (s->bc) {
1737:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1738:   }
1739:   return(0);
1740: }

1744: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1745: {

1749:   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);
1750:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
1751:   return(0);
1752: }

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

1761:   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);
1762:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
1763:   return(0);
1764: }

1768: /*@
1769:   PetscSectionPermute - Reorder the section according to the input point permutation

1771:   Collective on PetscSection

1773:   Input Parameter:
1774: + section - The PetscSection object
1775: - perm - The point permutation, old point p becomes new point perm[p]

1777:   Output Parameter:
1778: . sectionNew - The permuted PetscSection

1780:   Level: intermediate

1782: .keywords: mesh
1783: .seealso: MatPermute()
1784: @*/
1785: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1786: {
1787:   PetscSection    s = section, sNew;
1788:   const PetscInt *perm;
1789:   PetscInt        numFields, f, numPoints, pStart, pEnd, p;
1790:   PetscErrorCode  ierr;

1796:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1797:   PetscSectionGetNumFields(s, &numFields);
1798:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1799:   for (f = 0; f < numFields; ++f) {
1800:     const char *name;
1801:     PetscInt    numComp;

1803:     PetscSectionGetFieldName(s, f, &name);
1804:     PetscSectionSetFieldName(sNew, f, name);
1805:     PetscSectionGetFieldComponents(s, f, &numComp);
1806:     PetscSectionSetFieldComponents(sNew, f, numComp);
1807:   }
1808:   ISGetLocalSize(permutation, &numPoints);
1809:   ISGetIndices(permutation, &perm);
1810:   PetscSectionGetChart(s, &pStart, &pEnd);
1811:   PetscSectionSetChart(sNew, pStart, pEnd);
1812:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1813:   for (p = pStart; p < pEnd; ++p) {
1814:     PetscInt dof, cdof;

1816:     PetscSectionGetDof(s, p, &dof);
1817:     PetscSectionSetDof(sNew, perm[p], dof);
1818:     PetscSectionGetConstraintDof(s, p, &cdof);
1819:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1820:     for (f = 0; f < numFields; ++f) {
1821:       PetscSectionGetFieldDof(s, p, f, &dof);
1822:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1823:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1824:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1825:     }
1826:   }
1827:   PetscSectionSetUp(sNew);
1828:   for (p = pStart; p < pEnd; ++p) {
1829:     const PetscInt *cind;
1830:     PetscInt        cdof;

1832:     PetscSectionGetConstraintDof(s, p, &cdof);
1833:     if (cdof) {
1834:       PetscSectionGetConstraintIndices(s, p, &cind);
1835:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1836:     }
1837:     for (f = 0; f < numFields; ++f) {
1838:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1839:       if (cdof) {
1840:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1841:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1842:       }
1843:     }
1844:   }
1845:   ISRestoreIndices(permutation, &perm);
1846:   *sectionNew = sNew;
1847:   return(0);
1848: }

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

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

1860: PetscErrorCode PetscSFConvertPartition(PetscSF sfPart, PetscSection partSection, IS partition, ISLocalToGlobalMapping *renumbering, PetscSF *sf)
1861: {
1862:   MPI_Comm       comm;
1863:   PetscSF        sfPoints;
1864:   PetscInt       *partSizes,*partOffsets,p,i,numParts,numMyPoints,numPoints,count;
1865:   const PetscInt *partArray;
1866:   PetscSFNode    *sendPoints;
1867:   PetscMPIInt    rank;

1871:   PetscObjectGetComm((PetscObject)sfPart,&comm);
1872:   MPI_Comm_rank(comm, &rank);

1874:   /* Get the number of parts and sizes that I have to distribute */
1875:   PetscSFGetGraph(sfPart,NULL,&numParts,NULL,NULL);
1876:   PetscMalloc2(numParts,&partSizes,numParts,&partOffsets);
1877:   for (p=0,numPoints=0; p<numParts; p++) {
1878:     PetscSectionGetDof(partSection, p, &partSizes[p]);
1879:     numPoints += partSizes[p];
1880:   }
1881:   numMyPoints = 0;
1882:   PetscSFFetchAndOpBegin(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1883:   PetscSFFetchAndOpEnd(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1884:   /* I will receive numMyPoints. I will send a total of numPoints, to be placed on remote procs at partOffsets */

1886:   /* Create SF mapping locations (addressed through partition, as indexed leaves) to new owners (roots) */
1887:   PetscMalloc1(numPoints,&sendPoints);
1888:   for (p=0,count=0; p<numParts; p++) {
1889:     for (i=0; i<partSizes[p]; i++) {
1890:       sendPoints[count].rank = p;
1891:       sendPoints[count].index = partOffsets[p]+i;
1892:       count++;
1893:     }
1894:   }
1895:   if (count != numPoints) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count %D should equal numPoints=%D",count,numPoints);
1896:   PetscFree2(partSizes,partOffsets);
1897:   ISGetIndices(partition,&partArray);
1898:   PetscSFCreate(comm,&sfPoints);
1899:   PetscSFSetFromOptions(sfPoints);
1900:   PetscSFSetGraph(sfPoints,numMyPoints,numPoints,partArray,PETSC_USE_POINTER,sendPoints,PETSC_OWN_POINTER);

1902:   /* Invert SF so that the new owners are leaves and the locations indexed through partition are the roots */
1903:   PetscSFCreateInverseSF(sfPoints,sf);
1904:   PetscSFDestroy(&sfPoints);
1905:   ISRestoreIndices(partition,&partArray);

1907:   /* Create the new local-to-global mapping */
1908:   ISLocalToGlobalMappingCreateSF(*sf,0,renumbering);
1909:   return(0);
1910: }

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

1917:   Collective

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

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

1927:   Level: intermediate

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

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

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

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

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

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

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

2024:   Note: Either rootSection or remoteOffsets can be specified

2026:   Level: intermediate

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

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

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

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

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

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

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

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

2102:   Level: intermediate

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

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

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

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

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

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

2134:   Level: intermediate

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

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

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

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

2163:   Level: intermediate

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

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