Actual source code: section.c

petsc-3.14.0 2020-09-29
Report Typos and Errors
  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

  5: #include <petsc/private/sectionimpl.h>
  6: #include <petscsection.h>
  7: #include <petscsf.h>
  8: #include <petscviewer.h>

 10: PetscClassId PETSC_SECTION_CLASSID;

 12: /*@
 13:   PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.

 15:   Collective

 17:   Input Parameters:
 18: + comm - the MPI communicator
 19: - s    - pointer to the section

 21:   Level: beginner

 23:   Notes:
 24:   Typical calling sequence
 25: $       PetscSectionCreate(MPI_Comm,PetscSection *);
 26: $       PetscSectionSetNumFields(PetscSection, numFields);
 27: $       PetscSectionSetChart(PetscSection,low,high);
 28: $       PetscSectionSetDof(PetscSection,point,numdof);
 29: $       PetscSectionSetUp(PetscSection);
 30: $       PetscSectionGetOffset(PetscSection,point,PetscInt *);
 31: $       PetscSectionDestroy(PetscSection);

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

 36: .seealso: PetscSection, PetscSectionDestroy()
 37: @*/
 38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
 39: {

 44:   ISInitializePackage();

 46:   PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);

 48:   (*s)->pStart             = -1;
 49:   (*s)->pEnd               = -1;
 50:   (*s)->perm               = NULL;
 51:   (*s)->pointMajor         = PETSC_TRUE;
 52:   (*s)->maxDof             = 0;
 53:   (*s)->atlasDof           = NULL;
 54:   (*s)->atlasOff           = NULL;
 55:   (*s)->bc                 = NULL;
 56:   (*s)->bcIndices          = NULL;
 57:   (*s)->setup              = PETSC_FALSE;
 58:   (*s)->numFields          = 0;
 59:   (*s)->fieldNames         = NULL;
 60:   (*s)->field              = NULL;
 61:   (*s)->useFieldOff        = PETSC_FALSE;
 62:   (*s)->compNames          = NULL;
 63:   (*s)->clObj              = NULL;
 64:   (*s)->clHash             = NULL;
 65:   (*s)->clSection          = NULL;
 66:   (*s)->clPoints           = NULL;
 67:   return(0);
 68: }

 70: /*@
 71:   PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection

 73:   Collective

 75:   Input Parameter:
 76: . section - the PetscSection

 78:   Output Parameter:
 79: . newSection - the copy

 81:   Level: intermediate

 83: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
 84: @*/
 85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 86: {
 87:   PetscSectionSym sym;
 88:   IS              perm;
 89:   PetscInt        numFields, f, c, pStart, pEnd, p;
 90:   PetscErrorCode  ierr;

 95:   PetscSectionReset(newSection);
 96:   PetscSectionGetNumFields(section, &numFields);
 97:   if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
 98:   for (f = 0; f < numFields; ++f) {
 99:     const char *fieldName = NULL, *compName = NULL;
100:     PetscInt   numComp = 0;

102:     PetscSectionGetFieldName(section, f, &fieldName);
103:     PetscSectionSetFieldName(newSection, f, fieldName);
104:     PetscSectionGetFieldComponents(section, f, &numComp);
105:     PetscSectionSetFieldComponents(newSection, f, numComp);
106:     for (c = 0; c < numComp; ++c) {
107:       PetscSectionGetComponentName(section, f, c, &compName);
108:       PetscSectionSetComponentName(newSection, f, c, compName);
109:     }
110:     PetscSectionGetFieldSym(section, f, &sym);
111:     PetscSectionSetFieldSym(newSection, f, sym);
112:   }
113:   PetscSectionGetChart(section, &pStart, &pEnd);
114:   PetscSectionSetChart(newSection, pStart, pEnd);
115:   PetscSectionGetPermutation(section, &perm);
116:   PetscSectionSetPermutation(newSection, perm);
117:   PetscSectionGetSym(section, &sym);
118:   PetscSectionSetSym(newSection, sym);
119:   for (p = pStart; p < pEnd; ++p) {
120:     PetscInt dof, cdof, fcdof = 0;

122:     PetscSectionGetDof(section, p, &dof);
123:     PetscSectionSetDof(newSection, p, dof);
124:     PetscSectionGetConstraintDof(section, p, &cdof);
125:     if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
126:     for (f = 0; f < numFields; ++f) {
127:       PetscSectionGetFieldDof(section, p, f, &dof);
128:       PetscSectionSetFieldDof(newSection, p, f, dof);
129:       if (cdof) {
130:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
131:         if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
132:       }
133:     }
134:   }
135:   PetscSectionSetUp(newSection);
136:   for (p = pStart; p < pEnd; ++p) {
137:     PetscInt       off, cdof, fcdof = 0;
138:     const PetscInt *cInd;

140:     /* Must set offsets in case they do not agree with the prefix sums */
141:     PetscSectionGetOffset(section, p, &off);
142:     PetscSectionSetOffset(newSection, p, off);
143:     PetscSectionGetConstraintDof(section, p, &cdof);
144:     if (cdof) {
145:       PetscSectionGetConstraintIndices(section, p, &cInd);
146:       PetscSectionSetConstraintIndices(newSection, p, cInd);
147:       for (f = 0; f < numFields; ++f) {
148:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
149:         if (fcdof) {
150:           PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
151:           PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
152:         }
153:       }
154:     }
155:   }
156:   return(0);
157: }

159: /*@
160:   PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection

162:   Collective

164:   Input Parameter:
165: . section - the PetscSection

167:   Output Parameter:
168: . newSection - the copy

170:   Level: beginner

172: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
173: @*/
174: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
175: {

181:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
182:   PetscSectionCopy(section, *newSection);
183:   return(0);
184: }

186: /*@
187:   PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database

189:   Collective on PetscSection

191:   Input Parameter:
192: . section - the PetscSection

194:   Options Database:
195: . -petscsection_point_major the dof order

197:   Level: intermediate

199: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
200: @*/
201: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
202: {

207:   PetscObjectOptionsBegin((PetscObject) s);
208:   PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
209:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
210:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
211:   PetscOptionsEnd();
212:   PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
213:   return(0);
214: }

216: /*@
217:   PetscSectionCompare - Compares two sections

219:   Collective on PetscSection

221:   Input Parameters:
222: + s1 - the first PetscSection
223: - s2 - the second PetscSection

225:   Output Parameter:
226: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise

228:   Level: intermediate

230:   Notes:
231:   Field names are disregarded.

233: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
234: @*/
235: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
236: {
237:   PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
238:   const PetscInt *idx1, *idx2;
239:   IS perm1, perm2;
240:   PetscBool flg;
241:   PetscMPIInt mflg;

248:   flg = PETSC_FALSE;

250:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
251:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
252:     *congruent = PETSC_FALSE;
253:     return(0);
254:   }

256:   PetscSectionGetChart(s1, &pStart, &pEnd);
257:   PetscSectionGetChart(s2, &n1, &n2);
258:   if (pStart != n1 || pEnd != n2) goto not_congruent;

260:   PetscSectionGetPermutation(s1, &perm1);
261:   PetscSectionGetPermutation(s2, &perm2);
262:   if (perm1 && perm2) {
263:     ISEqual(perm1, perm2, congruent);
264:     if (!(*congruent)) goto not_congruent;
265:   } else if (perm1 != perm2) goto not_congruent;

267:   for (p = pStart; p < pEnd; ++p) {
268:     PetscSectionGetOffset(s1, p, &n1);
269:     PetscSectionGetOffset(s2, p, &n2);
270:     if (n1 != n2) goto not_congruent;

272:     PetscSectionGetDof(s1, p, &n1);
273:     PetscSectionGetDof(s2, p, &n2);
274:     if (n1 != n2) goto not_congruent;

276:     PetscSectionGetConstraintDof(s1, p, &ncdof);
277:     PetscSectionGetConstraintDof(s2, p, &n2);
278:     if (ncdof != n2) goto not_congruent;

280:     PetscSectionGetConstraintIndices(s1, p, &idx1);
281:     PetscSectionGetConstraintIndices(s2, p, &idx2);
282:     PetscArraycmp(idx1, idx2, ncdof, congruent);
283:     if (!(*congruent)) goto not_congruent;
284:   }

286:   PetscSectionGetNumFields(s1, &nfields);
287:   PetscSectionGetNumFields(s2, &n2);
288:   if (nfields != n2) goto not_congruent;

290:   for (f = 0; f < nfields; ++f) {
291:     PetscSectionGetFieldComponents(s1, f, &n1);
292:     PetscSectionGetFieldComponents(s2, f, &n2);
293:     if (n1 != n2) goto not_congruent;

295:     for (p = pStart; p < pEnd; ++p) {
296:       PetscSectionGetFieldOffset(s1, p, f, &n1);
297:       PetscSectionGetFieldOffset(s2, p, f, &n2);
298:       if (n1 != n2) goto not_congruent;

300:       PetscSectionGetFieldDof(s1, p, f, &n1);
301:       PetscSectionGetFieldDof(s2, p, f, &n2);
302:       if (n1 != n2) goto not_congruent;

304:       PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
305:       PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
306:       if (nfcdof != n2) goto not_congruent;

308:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
309:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
310:       PetscArraycmp(idx1, idx2, nfcdof, congruent);
311:       if (!(*congruent)) goto not_congruent;
312:     }
313:   }

315:   flg = PETSC_TRUE;
316: not_congruent:
317:   MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
318:   return(0);
319: }

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

324:   Not collective

326:   Input Parameter:
327: . s - the PetscSection

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

332:   Level: intermediate

334: .seealso: PetscSectionSetNumFields()
335: @*/
336: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
337: {
341:   *numFields = s->numFields;
342:   return(0);
343: }

345: /*@
346:   PetscSectionSetNumFields - Sets the number of fields.

348:   Not collective

350:   Input Parameters:
351: + s - the PetscSection
352: - numFields - the number of fields

354:   Level: intermediate

356: .seealso: PetscSectionGetNumFields()
357: @*/
358: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
359: {
360:   PetscInt       f;

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

368:   s->numFields = numFields;
369:   PetscMalloc1(s->numFields, &s->numFieldComponents);
370:   PetscMalloc1(s->numFields, &s->fieldNames);
371:   PetscMalloc1(s->numFields, &s->compNames);
372:   PetscMalloc1(s->numFields, &s->field);
373:   for (f = 0; f < s->numFields; ++f) {
374:     char name[64];

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

378:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
379:     PetscSNPrintf(name, 64, "Field_%D", f);
380:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
381:     PetscSNPrintf(name, 64, "Component_0");
382:     PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
383:     PetscStrallocpy(name, (char **) &s->compNames[f][0]);
384:   }
385:   return(0);
386: }

388: /*@C
389:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

391:   Not Collective

393:   Input Parameters:
394: + s     - the PetscSection
395: - field - the field number

397:   Output Parameter:
398: . fieldName - the field name

400:   Level: intermediate

402: .seealso: PetscSectionSetFieldName()
403: @*/
404: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
405: {
409:   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);
410:   *fieldName = s->fieldNames[field];
411:   return(0);
412: }

414: /*@C
415:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

417:   Not Collective

419:   Input Parameters:
420: + s     - the PetscSection
421: . field - the field number
422: - fieldName - the field name

424:   Level: intermediate

426: .seealso: PetscSectionGetFieldName()
427: @*/
428: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
429: {

435:   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);
436:   PetscFree(s->fieldNames[field]);
437:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
438:   return(0);
439: }

441: /*@C
442:   PetscSectionGetComponentName - Gets the name of a field component in the PetscSection

444:   Not Collective

446:   Input Parameters:
447: + s     - the PetscSection
448: . field - the field number
449: . comp  - the component number
450: - compName - the component name

452:   Level: intermediate

454: .seealso: PetscSectionSetComponentName()
455: @*/
456: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
457: {
461:   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);
462:   if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
463:   *compName = s->compNames[field][comp];
464:   return(0);
465: }

467: /*@C
468:   PetscSectionSetComponentName - Sets the name of a field component in the PetscSection

470:   Not Collective

472:   Input Parameters:
473: + s     - the PetscSection
474: . field - the field number
475: . comp  - the component number
476: - compName - the component name

478:   Level: intermediate

480: .seealso: PetscSectionGetComponentName()
481: @*/
482: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
483: {

489:   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);
490:   if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
491:   PetscFree(s->compNames[field][comp]);
492:   PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
493:   return(0);
494: }

496: /*@
497:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

499:   Not collective

501:   Input Parameters:
502: + s - the PetscSection
503: - field - the field number

505:   Output Parameter:
506: . numComp - the number of field components

508:   Level: intermediate

510: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
511: @*/
512: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
513: {
517:   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);
518:   *numComp = s->numFieldComponents[field];
519:   return(0);
520: }

522: /*@
523:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

525:   Not collective

527:   Input Parameters:
528: + s - the PetscSection
529: . field - the field number
530: - numComp - the number of field components

532:   Level: intermediate

534: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
535: @*/
536: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
537: {
539:   PetscInt c;
540:   char name[64];

544:   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);
545:   if (s->compNames) {
546:     for (c = 0; c < s->numFieldComponents[field]; ++c) {
547:       PetscFree(s->compNames[field][c]);
548:     }
549:     PetscFree(s->compNames[field]);
550:   }

552:   s->numFieldComponents[field] = numComp;
553:   if (numComp) {
554:     PetscMalloc1(numComp, (char ***) &s->compNames[field]);
555:     for (c = 0; c < numComp; ++c) {
556:       PetscSNPrintf(name, 64, "%D", c);
557:       PetscStrallocpy(name, (char **) &s->compNames[field][c]);
558:     }
559:   }
560:   return(0);
561: }

563: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
564: {

568:   if (!s->bc) {
569:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
570:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
571:   }
572:   return(0);
573: }

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

578:   Not collective

580:   Input Parameter:
581: . s - the PetscSection

583:   Output Parameters:
584: + pStart - the first point
585: - pEnd - one past the last point

587:   Level: intermediate

589: .seealso: PetscSectionSetChart(), PetscSectionCreate()
590: @*/
591: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
592: {
595:   if (pStart) *pStart = s->pStart;
596:   if (pEnd)   *pEnd   = s->pEnd;
597:   return(0);
598: }

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

603:   Not collective

605:   Input Parameters:
606: + s - the PetscSection
607: . pStart - the first point
608: - pEnd - one past the last point

610:   Level: intermediate

612: .seealso: PetscSectionGetChart(), PetscSectionCreate()
613: @*/
614: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
615: {
616:   PetscInt       f;

621:   /* Cannot Reset() because it destroys field information */
622:   s->setup = PETSC_FALSE;
623:   PetscSectionDestroy(&s->bc);
624:   PetscFree(s->bcIndices);
625:   PetscFree2(s->atlasDof, s->atlasOff);

627:   s->pStart = pStart;
628:   s->pEnd   = pEnd;
629:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
630:   PetscArrayzero(s->atlasDof, pEnd - pStart);
631:   for (f = 0; f < s->numFields; ++f) {
632:     PetscSectionSetChart(s->field[f], pStart, pEnd);
633:   }
634:   return(0);
635: }

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

640:   Not collective

642:   Input Parameter:
643: . s - the PetscSection

645:   Output Parameters:
646: . perm - The permutation as an IS

648:   Level: intermediate

650: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
651: @*/
652: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
653: {
657:   return(0);
658: }

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

663:   Not collective

665:   Input Parameters:
666: + s - the PetscSection
667: - perm - the permutation of points

669:   Level: intermediate

671: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
672: @*/
673: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
674: {

680:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
681:   if (s->perm != perm) {
682:     ISDestroy(&s->perm);
683:     if (perm) {
684:       s->perm = perm;
685:       PetscObjectReference((PetscObject) s->perm);
686:     }
687:   }
688:   return(0);
689: }

691: /*@
692:   PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major

694:   Not collective

696:   Input Parameter:
697: . s - the PetscSection

699:   Output Parameter:
700: . pm - the flag for point major ordering

702:   Level: intermediate

704: .seealso: PetscSectionSetPointMajor()
705: @*/
706: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
707: {
711:   *pm = s->pointMajor;
712:   return(0);
713: }

715: /*@
716:   PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major

718:   Not collective

720:   Input Parameters:
721: + s  - the PetscSection
722: - pm - the flag for point major ordering

724:   Not collective

726:   Level: intermediate

728: .seealso: PetscSectionGetPointMajor()
729: @*/
730: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
731: {
734:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
735:   s->pointMajor = pm;
736:   return(0);
737: }

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

742:   Not collective

744:   Input Parameters:
745: + s - the PetscSection
746: - point - the point

748:   Output Parameter:
749: . numDof - the number of dof

751:   Level: intermediate

753: .seealso: PetscSectionSetDof(), PetscSectionCreate()
754: @*/
755: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
756: {
760:   if (PetscDefined(USE_DEBUG)) {
761:     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);
762:   }
763:   *numDof = s->atlasDof[point - s->pStart];
764:   return(0);
765: }

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

770:   Not collective

772:   Input Parameters:
773: + s - the PetscSection
774: . point - the point
775: - numDof - the number of dof

777:   Level: intermediate

779: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
780: @*/
781: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
782: {
785:   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);
786:   s->atlasDof[point - s->pStart] = numDof;
787:   return(0);
788: }

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

793:   Not collective

795:   Input Parameters:
796: + s - the PetscSection
797: . point - the point
798: - numDof - the number of additional dof

800:   Level: intermediate

802: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
803: @*/
804: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
805: {
808:   if (PetscDefined(USE_DEBUG)) {
809:     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);
810:   }
811:   s->atlasDof[point - s->pStart] += numDof;
812:   return(0);
813: }

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

818:   Not collective

820:   Input Parameters:
821: + s - the PetscSection
822: . point - the point
823: - field - the field

825:   Output Parameter:
826: . numDof - the number of dof

828:   Level: intermediate

830: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
831: @*/
832: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
833: {

838:   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);
839:   PetscSectionGetDof(s->field[field], point, numDof);
840:   return(0);
841: }

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

846:   Not collective

848:   Input Parameters:
849: + s - the PetscSection
850: . point - the point
851: . field - the field
852: - numDof - the number of dof

854:   Level: intermediate

856: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
857: @*/
858: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
859: {

864:   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);
865:   PetscSectionSetDof(s->field[field], point, numDof);
866:   return(0);
867: }

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

872:   Not collective

874:   Input Parameters:
875: + s - the PetscSection
876: . point - the point
877: . field - the field
878: - numDof - the number of dof

880:   Level: intermediate

882: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
883: @*/
884: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
885: {

890:   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);
891:   PetscSectionAddDof(s->field[field], point, numDof);
892:   return(0);
893: }

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

898:   Not collective

900:   Input Parameters:
901: + s - the PetscSection
902: - point - the point

904:   Output Parameter:
905: . numDof - the number of dof which are fixed by constraints

907:   Level: intermediate

909: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
910: @*/
911: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
912: {

918:   if (s->bc) {
919:     PetscSectionGetDof(s->bc, point, numDof);
920:   } else *numDof = 0;
921:   return(0);
922: }

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

927:   Not collective

929:   Input Parameters:
930: + s - the PetscSection
931: . point - the point
932: - numDof - the number of dof which are fixed by constraints

934:   Level: intermediate

936: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
937: @*/
938: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
939: {

944:   if (numDof) {
945:     PetscSectionCheckConstraints_Static(s);
946:     PetscSectionSetDof(s->bc, point, numDof);
947:   }
948:   return(0);
949: }

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

954:   Not collective

956:   Input Parameters:
957: + s - the PetscSection
958: . point - the point
959: - numDof - the number of additional dof which are fixed by constraints

961:   Level: intermediate

963: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
964: @*/
965: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
966: {

971:   if (numDof) {
972:     PetscSectionCheckConstraints_Static(s);
973:     PetscSectionAddDof(s->bc, point, numDof);
974:   }
975:   return(0);
976: }

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

981:   Not collective

983:   Input Parameters:
984: + s - the PetscSection
985: . point - the point
986: - field - the field

988:   Output Parameter:
989: . numDof - the number of dof which are fixed by constraints

991:   Level: intermediate

993: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
994: @*/
995: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
996: {

1002:   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);
1003:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
1004:   return(0);
1005: }

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

1010:   Not collective

1012:   Input Parameters:
1013: + s - the PetscSection
1014: . point - the point
1015: . field - the field
1016: - numDof - the number of dof which are fixed by constraints

1018:   Level: intermediate

1020: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1021: @*/
1022: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1023: {

1028:   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);
1029:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
1030:   return(0);
1031: }

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

1036:   Not collective

1038:   Input Parameters:
1039: + s - the PetscSection
1040: . point - the point
1041: . field - the field
1042: - numDof - the number of additional dof which are fixed by constraints

1044:   Level: intermediate

1046: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1047: @*/
1048: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1049: {

1054:   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);
1055:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
1056:   return(0);
1057: }

1059: /*@
1060:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1062:   Not collective

1064:   Input Parameter:
1065: . s - the PetscSection

1067:   Level: advanced

1069: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1070: @*/
1071: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1072: {

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

1080:     PetscSectionSetUp(s->bc);
1081:     PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1082:   }
1083:   return(0);
1084: }

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

1089:   Not collective

1091:   Input Parameter:
1092: . s - the PetscSection

1094:   Level: intermediate

1096: .seealso: PetscSectionCreate()
1097: @*/
1098: PetscErrorCode PetscSectionSetUp(PetscSection s)
1099: {
1100:   const PetscInt *pind   = NULL;
1101:   PetscInt        offset = 0, foff, p, f;
1102:   PetscErrorCode  ierr;

1106:   if (s->setup) return(0);
1107:   s->setup = PETSC_TRUE;
1108:   /* Set offsets and field offsets for all points */
1109:   /*   Assume that all fields have the same chart */
1110:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1111:   if (s->pointMajor) {
1112:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1113:       const PetscInt q = pind ? pind[p] : p;

1115:       /* Set point offset */
1116:       s->atlasOff[q] = offset;
1117:       offset        += s->atlasDof[q];
1118:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
1119:       /* Set field offset */
1120:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1121:         PetscSection sf = s->field[f];

1123:         sf->atlasOff[q] = foff;
1124:         foff += sf->atlasDof[q];
1125:       }
1126:     }
1127:   } else {
1128:     /* Set field offsets for all points */
1129:     for (f = 0; f < s->numFields; ++f) {
1130:       PetscSection sf = s->field[f];

1132:       for (p = 0; p < s->pEnd - s->pStart; ++p) {
1133:         const PetscInt q = pind ? pind[p] : p;

1135:         sf->atlasOff[q] = offset;
1136:         offset += sf->atlasDof[q];
1137:       }
1138:     }
1139:     /* Disable point offsets since these are unused */
1140:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1141:       s->atlasOff[p] = -1;
1142:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[p]);
1143:     }
1144:   }
1145:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1146:   /* Setup BC sections */
1147:   PetscSectionSetUpBC(s);
1148:   for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1149:   return(0);
1150: }

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

1155:   Not collective

1157:   Input Parameters:
1158: . s - the PetscSection

1160:   Output Parameter:
1161: . maxDof - the maximum dof

1163:   Level: intermediate

1165: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1166: @*/
1167: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1168: {
1172:   *maxDof = s->maxDof;
1173:   return(0);
1174: }

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

1179:   Not collective

1181:   Input Parameter:
1182: . s - the PetscSection

1184:   Output Parameter:
1185: . size - the size of an array which can hold all the dofs

1187:   Level: intermediate

1189: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1190: @*/
1191: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1192: {
1193:   PetscInt p, n = 0;

1198:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1199:   *size = n;
1200:   return(0);
1201: }

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

1206:   Not collective

1208:   Input Parameter:
1209: . s - the PetscSection

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

1214:   Level: intermediate

1216: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1217: @*/
1218: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1219: {
1220:   PetscInt p, n = 0;

1225:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1226:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1227:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1228:   }
1229:   *size = n;
1230:   return(0);
1231: }

1233: /*@
1234:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1235:   the local section and an SF describing the section point overlap.

1237:   Input Parameters:
1238: + s - The PetscSection for the local field layout
1239: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1240: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1241: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points

1243:   Output Parameter:
1244: . gsection - The PetscSection for the global field layout

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

1248:   Level: intermediate

1250: .seealso: PetscSectionCreate()
1251: @*/
1252: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1253: {
1254:   PetscSection    gs;
1255:   const PetscInt *pind = NULL;
1256:   PetscInt       *recv = NULL, *neg = NULL;
1257:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1258:   PetscErrorCode  ierr;

1266:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1267:   PetscSectionGetChart(s, &pStart, &pEnd);
1268:   PetscSectionSetChart(gs, pStart, pEnd);
1269:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1270:   nlocal = nroots;              /* The local/leaf space matches global/root space */
1271:   /* Must allocate for all points visible to SF, which may be more than this section */
1272:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
1273:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
1274:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1275:     if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1276:     PetscMalloc2(nroots,&neg,nlocal,&recv);
1277:     PetscArrayzero(neg,nroots);
1278:   }
1279:   /* Mark all local points with negative dof */
1280:   for (p = pStart; p < pEnd; ++p) {
1281:     PetscSectionGetDof(s, p, &dof);
1282:     PetscSectionSetDof(gs, p, dof);
1283:     PetscSectionGetConstraintDof(s, p, &cdof);
1284:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1285:     if (neg) neg[p] = -(dof+1);
1286:   }
1287:   PetscSectionSetUpBC(gs);
1288:   if (gs->bcIndices) {PetscArraycpy(gs->bcIndices, s->bcIndices,gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]);}
1289:   if (nroots >= 0) {
1290:     PetscArrayzero(recv,nlocal);
1291:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1292:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1293:     for (p = pStart; p < pEnd; ++p) {
1294:       if (recv[p] < 0) {
1295:         gs->atlasDof[p-pStart] = recv[p];
1296:         PetscSectionGetDof(s, p, &dof);
1297:         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);
1298:       }
1299:     }
1300:   }
1301:   /* Calculate new sizes, get process offset, and calculate point offsets */
1302:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1303:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1304:     const PetscInt q = pind ? pind[p] : p;

1306:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1307:     gs->atlasOff[q] = off;
1308:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1309:   }
1310:   if (!localOffsets) {
1311:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1312:     globalOff -= off;
1313:   }
1314:   for (p = pStart, off = 0; p < pEnd; ++p) {
1315:     gs->atlasOff[p-pStart] += globalOff;
1316:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1317:   }
1318:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1319:   /* Put in negative offsets for ghost points */
1320:   if (nroots >= 0) {
1321:     PetscArrayzero(recv,nlocal);
1322:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1323:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1324:     for (p = pStart; p < pEnd; ++p) {
1325:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1326:     }
1327:   }
1328:   PetscFree2(neg,recv);
1329:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1330:   *gsection = gs;
1331:   return(0);
1332: }

1334: /*@
1335:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1336:   the local section and an SF describing the section point overlap.

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

1345:   Output Parameter:
1346:   . gsection - The PetscSection for the global field layout

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

1350:   Level: advanced

1352: .seealso: PetscSectionCreate()
1353: @*/
1354: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1355: {
1356:   const PetscInt *pind = NULL;
1357:   PetscInt       *neg  = NULL, *tmpOff = NULL;
1358:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1359:   PetscErrorCode  ierr;

1365:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1366:   PetscSectionGetChart(s, &pStart, &pEnd);
1367:   PetscSectionSetChart(*gsection, pStart, pEnd);
1368:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1369:   if (nroots >= 0) {
1370:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1371:     PetscCalloc1(nroots, &neg);
1372:     if (nroots > pEnd-pStart) {
1373:       PetscCalloc1(nroots, &tmpOff);
1374:     } else {
1375:       tmpOff = &(*gsection)->atlasDof[-pStart];
1376:     }
1377:   }
1378:   /* Mark ghost points with negative dof */
1379:   for (p = pStart; p < pEnd; ++p) {
1380:     for (e = 0; e < numExcludes; ++e) {
1381:       if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1382:         PetscSectionSetDof(*gsection, p, 0);
1383:         break;
1384:       }
1385:     }
1386:     if (e < numExcludes) continue;
1387:     PetscSectionGetDof(s, p, &dof);
1388:     PetscSectionSetDof(*gsection, p, dof);
1389:     PetscSectionGetConstraintDof(s, p, &cdof);
1390:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1391:     if (neg) neg[p] = -(dof+1);
1392:   }
1393:   PetscSectionSetUpBC(*gsection);
1394:   if (nroots >= 0) {
1395:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1396:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1397:     if (nroots > pEnd - pStart) {
1398:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1399:     }
1400:   }
1401:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
1402:   if (s->perm) {ISGetIndices(s->perm, &pind);}
1403:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1404:     const PetscInt q = pind ? pind[p] : p;

1406:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1407:     (*gsection)->atlasOff[q] = off;
1408:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1409:   }
1410:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1411:   globalOff -= off;
1412:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1413:     (*gsection)->atlasOff[p] += globalOff;
1414:     if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1415:   }
1416:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1417:   /* Put in negative offsets for ghost points */
1418:   if (nroots >= 0) {
1419:     if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1420:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1421:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1422:     if (nroots > pEnd - pStart) {
1423:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1424:     }
1425:   }
1426:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1427:   PetscFree(neg);
1428:   return(0);
1429: }

1431: /*@
1432:   PetscSectionGetPointLayout - Get the PetscLayout associated with the section points

1434:   Collective on comm

1436:   Input Parameters:
1437: + comm - The MPI_Comm
1438: - s    - The PetscSection

1440:   Output Parameter:
1441: . layout - The point layout for the section

1443:   Note: This is usually caleld for the default global section.

1445:   Level: advanced

1447: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1448: @*/
1449: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1450: {
1451:   PetscInt       pStart, pEnd, p, localSize = 0;

1455:   PetscSectionGetChart(s, &pStart, &pEnd);
1456:   for (p = pStart; p < pEnd; ++p) {
1457:     PetscInt dof;

1459:     PetscSectionGetDof(s, p, &dof);
1460:     if (dof > 0) ++localSize;
1461:   }
1462:   PetscLayoutCreate(comm, layout);
1463:   PetscLayoutSetLocalSize(*layout, localSize);
1464:   PetscLayoutSetBlockSize(*layout, 1);
1465:   PetscLayoutSetUp(*layout);
1466:   return(0);
1467: }

1469: /*@
1470:   PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.

1472:   Collective on comm

1474:   Input Parameters:
1475: + comm - The MPI_Comm
1476: - s    - The PetscSection

1478:   Output Parameter:
1479: . layout - The dof layout for the section

1481:   Note: This is usually called for the default global section.

1483:   Level: advanced

1485: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1486: @*/
1487: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1488: {
1489:   PetscInt       pStart, pEnd, p, localSize = 0;

1495:   PetscSectionGetChart(s, &pStart, &pEnd);
1496:   for (p = pStart; p < pEnd; ++p) {
1497:     PetscInt dof,cdof;

1499:     PetscSectionGetDof(s, p, &dof);
1500:     PetscSectionGetConstraintDof(s, p, &cdof);
1501:     if (dof-cdof > 0) localSize += dof-cdof;
1502:   }
1503:   PetscLayoutCreate(comm, layout);
1504:   PetscLayoutSetLocalSize(*layout, localSize);
1505:   PetscLayoutSetBlockSize(*layout, 1);
1506:   PetscLayoutSetUp(*layout);
1507:   return(0);
1508: }

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

1513:   Not collective

1515:   Input Parameters:
1516: + s - the PetscSection
1517: - point - the point

1519:   Output Parameter:
1520: . offset - the offset

1522:   Level: intermediate

1524: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1525: @*/
1526: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1527: {
1531:   if (PetscDefined(USE_DEBUG)) {
1532:     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);
1533:   }
1534:   *offset = s->atlasOff[point - s->pStart];
1535:   return(0);
1536: }

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

1541:   Not collective

1543:   Input Parameters:
1544: + s - the PetscSection
1545: . point - the point
1546: - offset - the offset

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

1550:   Level: intermediate

1552: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1553: @*/
1554: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1555: {
1558:   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);
1559:   s->atlasOff[point - s->pStart] = offset;
1560:   return(0);
1561: }

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

1566:   Not collective

1568:   Input Parameters:
1569: + s - the PetscSection
1570: . point - the point
1571: - field - the field

1573:   Output Parameter:
1574: . offset - the offset

1576:   Level: intermediate

1578: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1579: @*/
1580: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1581: {

1587:   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);
1588:   PetscSectionGetOffset(s->field[field], point, offset);
1589:   return(0);
1590: }

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

1595:   Not collective

1597:   Input Parameters:
1598: + s - the PetscSection
1599: . point - the point
1600: . field - the field
1601: - offset - the offset

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

1605:   Level: intermediate

1607: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1608: @*/
1609: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1610: {

1615:   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);
1616:   PetscSectionSetOffset(s->field[field], point, offset);
1617:   return(0);
1618: }

1620: /*@
1621:   PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.

1623:   Not collective

1625:   Input Parameters:
1626: + s - the PetscSection
1627: . point - the point
1628: - field - the field

1630:   Output Parameter:
1631: . offset - the offset

1633:   Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1634:         this point, what number is the first dof with this field.

1636:   Level: advanced

1638: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1639: @*/
1640: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1641: {
1642:   PetscInt       off, foff;

1648:   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);
1649:   PetscSectionGetOffset(s, point, &off);
1650:   PetscSectionGetOffset(s->field[field], point, &foff);
1651:   *offset = foff - off;
1652:   return(0);
1653: }

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

1658:   Not collective

1660:   Input Parameter:
1661: . s - the PetscSection

1663:   Output Parameters:
1664: + start - the minimum offset
1665: - end   - one more than the maximum offset

1667:   Level: intermediate

1669: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1670: @*/
1671: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1672: {
1673:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

1678:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1679:   PetscSectionGetChart(s, &pStart, &pEnd);
1680:   for (p = 0; p < pEnd-pStart; ++p) {
1681:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1683:     if (off >= 0) {
1684:       os = PetscMin(os, off);
1685:       oe = PetscMax(oe, off+dof);
1686:     }
1687:   }
1688:   if (start) *start = os;
1689:   if (end)   *end   = oe;
1690:   return(0);
1691: }

1693: /*@
1694:   PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields

1696:   Collective on s

1698:   Input Parameter:
1699: + s      - the PetscSection
1700: . len    - the number of subfields
1701: - fields - the subfield numbers

1703:   Output Parameter:
1704: . subs   - the subsection

1706:   Note: The section offsets now refer to a new, smaller vector.

1708:   Level: advanced

1710: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1711: @*/
1712: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1713: {
1714:   PetscInt       nF, f, c, pStart, pEnd, p, maxCdof = 0;

1718:   if (!len) return(0);
1722:   PetscSectionGetNumFields(s, &nF);
1723:   if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1724:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1725:   PetscSectionSetNumFields(*subs, len);
1726:   for (f = 0; f < len; ++f) {
1727:     const char *name   = NULL;
1728:     PetscInt   numComp = 0;

1730:     PetscSectionGetFieldName(s, fields[f], &name);
1731:     PetscSectionSetFieldName(*subs, f, name);
1732:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1733:     PetscSectionSetFieldComponents(*subs, f, numComp);
1734:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1735:       PetscSectionGetComponentName(s, fields[f], c, &name);
1736:       PetscSectionSetComponentName(*subs, f, c, name);
1737:     }
1738:   }
1739:   PetscSectionGetChart(s, &pStart, &pEnd);
1740:   PetscSectionSetChart(*subs, pStart, pEnd);
1741:   for (p = pStart; p < pEnd; ++p) {
1742:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1744:     for (f = 0; f < len; ++f) {
1745:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1746:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1747:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1748:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1749:       dof  += fdof;
1750:       cdof += cfdof;
1751:     }
1752:     PetscSectionSetDof(*subs, p, dof);
1753:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1754:     maxCdof = PetscMax(cdof, maxCdof);
1755:   }
1756:   PetscSectionSetUp(*subs);
1757:   if (maxCdof) {
1758:     PetscInt *indices;

1760:     PetscMalloc1(maxCdof, &indices);
1761:     for (p = pStart; p < pEnd; ++p) {
1762:       PetscInt cdof;

1764:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1765:       if (cdof) {
1766:         const PetscInt *oldIndices = NULL;
1767:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1769:         for (f = 0; f < len; ++f) {
1770:           PetscInt oldFoff = 0, oldf;

1772:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1773:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1774:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1775:           /* This can be sped up if we assume sorted fields */
1776:           for (oldf = 0; oldf < fields[f]; ++oldf) {
1777:             PetscInt oldfdof = 0;
1778:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1779:             oldFoff += oldfdof;
1780:           }
1781:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1782:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1783:           numConst += cfdof;
1784:           fOff     += fdof;
1785:         }
1786:         PetscSectionSetConstraintIndices(*subs, p, indices);
1787:       }
1788:     }
1789:     PetscFree(indices);
1790:   }
1791:   return(0);
1792: }

1794: /*@
1795:   PetscSectionCreateSupersection - Create a new, larger section composed of the input sections

1797:   Collective on s

1799:   Input Parameters:
1800: + s     - the input sections
1801: - len   - the number of input sections

1803:   Output Parameter:
1804: . supers - the supersection

1806:   Note: The section offsets now refer to a new, larger vector.

1808:   Level: advanced

1810: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1811: @*/
1812: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1813: {
1814:   PetscInt       Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1818:   if (!len) return(0);
1819:   for (i = 0; i < len; ++i) {
1820:     PetscInt nf, pStarti, pEndi;

1822:     PetscSectionGetNumFields(s[i], &nf);
1823:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1824:     pStart = PetscMin(pStart, pStarti);
1825:     pEnd   = PetscMax(pEnd,   pEndi);
1826:     Nf += nf;
1827:   }
1828:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1829:   PetscSectionSetNumFields(*supers, Nf);
1830:   for (i = 0, f = 0; i < len; ++i) {
1831:     PetscInt nf, fi, ci;

1833:     PetscSectionGetNumFields(s[i], &nf);
1834:     for (fi = 0; fi < nf; ++fi, ++f) {
1835:       const char *name   = NULL;
1836:       PetscInt   numComp = 0;

1838:       PetscSectionGetFieldName(s[i], fi, &name);
1839:       PetscSectionSetFieldName(*supers, f, name);
1840:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1841:       PetscSectionSetFieldComponents(*supers, f, numComp);
1842:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1843:         PetscSectionGetComponentName(s[i], fi, ci, &name);
1844:         PetscSectionSetComponentName(*supers, f, ci, name);
1845:       }
1846:     }
1847:   }
1848:   PetscSectionSetChart(*supers, pStart, pEnd);
1849:   for (p = pStart; p < pEnd; ++p) {
1850:     PetscInt dof = 0, cdof = 0;

1852:     for (i = 0, f = 0; i < len; ++i) {
1853:       PetscInt nf, fi, pStarti, pEndi;
1854:       PetscInt fdof = 0, cfdof = 0;

1856:       PetscSectionGetNumFields(s[i], &nf);
1857:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1858:       if ((p < pStarti) || (p >= pEndi)) continue;
1859:       for (fi = 0; fi < nf; ++fi, ++f) {
1860:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1861:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1862:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1863:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1864:         dof  += fdof;
1865:         cdof += cfdof;
1866:       }
1867:     }
1868:     PetscSectionSetDof(*supers, p, dof);
1869:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1870:     maxCdof = PetscMax(cdof, maxCdof);
1871:   }
1872:   PetscSectionSetUp(*supers);
1873:   if (maxCdof) {
1874:     PetscInt *indices;

1876:     PetscMalloc1(maxCdof, &indices);
1877:     for (p = pStart; p < pEnd; ++p) {
1878:       PetscInt cdof;

1880:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1881:       if (cdof) {
1882:         PetscInt dof, numConst = 0, fOff = 0;

1884:         for (i = 0, f = 0; i < len; ++i) {
1885:           const PetscInt *oldIndices = NULL;
1886:           PetscInt        nf, fi, pStarti, pEndi, fdof, cfdof, fc;

1888:           PetscSectionGetNumFields(s[i], &nf);
1889:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1890:           if ((p < pStarti) || (p >= pEndi)) continue;
1891:           for (fi = 0; fi < nf; ++fi, ++f) {
1892:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1893:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1894:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1895:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1896:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1897:             numConst += cfdof;
1898:           }
1899:           PetscSectionGetDof(s[i], p, &dof);
1900:           fOff += dof;
1901:         }
1902:         PetscSectionSetConstraintIndices(*supers, p, indices);
1903:       }
1904:     }
1905:     PetscFree(indices);
1906:   }
1907:   return(0);
1908: }

1910: /*@
1911:   PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh

1913:   Collective on s

1915:   Input Parameters:
1916: + s           - the PetscSection
1917: - subpointMap - a sorted list of points in the original mesh which are in the submesh

1919:   Output Parameter:
1920: . subs - the subsection

1922:   Note: The section offsets now refer to a new, smaller vector.

1924:   Level: advanced

1926: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1927: @*/
1928: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1929: {
1930:   const PetscInt *points = NULL, *indices = NULL;
1931:   PetscInt       numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;

1938:   PetscSectionGetNumFields(s, &numFields);
1939:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1940:   if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1941:   for (f = 0; f < numFields; ++f) {
1942:     const char *name   = NULL;
1943:     PetscInt   numComp = 0;

1945:     PetscSectionGetFieldName(s, f, &name);
1946:     PetscSectionSetFieldName(*subs, f, name);
1947:     PetscSectionGetFieldComponents(s, f, &numComp);
1948:     PetscSectionSetFieldComponents(*subs, f, numComp);
1949:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
1950:       PetscSectionGetComponentName(s, f, c, &name);
1951:       PetscSectionSetComponentName(*subs, f, c, name);
1952:     }
1953:   }
1954:   /* For right now, we do not try to squeeze the subchart */
1955:   if (subpointMap) {
1956:     ISGetSize(subpointMap, &numSubpoints);
1957:     ISGetIndices(subpointMap, &points);
1958:   }
1959:   PetscSectionGetChart(s, &pStart, &pEnd);
1960:   PetscSectionSetChart(*subs, 0, numSubpoints);
1961:   for (p = pStart; p < pEnd; ++p) {
1962:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1964:     PetscFindInt(p, numSubpoints, points, &subp);
1965:     if (subp < 0) continue;
1966:     for (f = 0; f < numFields; ++f) {
1967:       PetscSectionGetFieldDof(s, p, f, &fdof);
1968:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1969:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1970:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1971:     }
1972:     PetscSectionGetDof(s, p, &dof);
1973:     PetscSectionSetDof(*subs, subp, dof);
1974:     PetscSectionGetConstraintDof(s, p, &cdof);
1975:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1976:   }
1977:   PetscSectionSetUp(*subs);
1978:   /* Change offsets to original offsets */
1979:   for (p = pStart; p < pEnd; ++p) {
1980:     PetscInt off, foff = 0;

1982:     PetscFindInt(p, numSubpoints, points, &subp);
1983:     if (subp < 0) continue;
1984:     for (f = 0; f < numFields; ++f) {
1985:       PetscSectionGetFieldOffset(s, p, f, &foff);
1986:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1987:     }
1988:     PetscSectionGetOffset(s, p, &off);
1989:     PetscSectionSetOffset(*subs, subp, off);
1990:   }
1991:   /* Copy constraint indices */
1992:   for (subp = 0; subp < numSubpoints; ++subp) {
1993:     PetscInt cdof;

1995:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1996:     if (cdof) {
1997:       for (f = 0; f < numFields; ++f) {
1998:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1999:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2000:       }
2001:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
2002:       PetscSectionSetConstraintIndices(*subs, subp, indices);
2003:     }
2004:   }
2005:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
2006:   return(0);
2007: }

2009: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2010: {
2011:   PetscInt       p;
2012:   PetscMPIInt    rank;

2016:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2017:   PetscViewerASCIIPushSynchronized(viewer);
2018:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2019:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2020:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2021:       PetscInt b;

2023:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2024:       if (s->bcIndices) {
2025:         for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2026:           PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
2027:         }
2028:       }
2029:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2030:     } else {
2031:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2032:     }
2033:   }
2034:   PetscViewerFlush(viewer);
2035:   PetscViewerASCIIPopSynchronized(viewer);
2036:   if (s->sym) {
2037:     PetscViewerASCIIPushTab(viewer);
2038:     PetscSectionSymView(s->sym,viewer);
2039:     PetscViewerASCIIPopTab(viewer);
2040:   }
2041:   return(0);
2042: }

2044: /*@C
2045:    PetscSectionViewFromOptions - View from Options

2047:    Collective on PetscSection

2049:    Input Parameters:
2050: +  A - the PetscSection object to view
2051: .  obj - Optional object
2052: -  name - command line option

2054:    Level: intermediate
2055: .seealso:  PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2056: @*/
2057: PetscErrorCode  PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2058: {

2063:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
2064:   return(0);
2065: }

2067: /*@C
2068:   PetscSectionView - Views a PetscSection

2070:   Collective on PetscSection

2072:   Input Parameters:
2073: + s - the PetscSection object to view
2074: - v - the viewer

2076:   Level: beginner

2078: .seealso PetscSectionCreate(), PetscSectionDestroy()
2079: @*/
2080: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2081: {
2082:   PetscBool      isascii;
2083:   PetscInt       f;

2088:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
2090:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2091:   if (isascii) {
2092:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2093:     if (s->numFields) {
2094:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2095:       for (f = 0; f < s->numFields; ++f) {
2096:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
2097:         PetscSectionView_ASCII(s->field[f], viewer);
2098:       }
2099:     } else {
2100:       PetscSectionView_ASCII(s, viewer);
2101:     }
2102:   }
2103:   return(0);
2104: }

2106: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2107: {
2109:   PetscSectionClosurePermVal clVal;

2112:   if (!section->clHash) return(0);
2113:   kh_foreach_value(section->clHash, clVal, {
2114:       PetscFree(clVal.perm);
2115:       PetscFree(clVal.invPerm);
2116:     });
2117:   kh_destroy(ClPerm, section->clHash);
2118:   section->clHash = NULL;
2119:   return(0);
2120: }

2122: /*@
2123:   PetscSectionReset - Frees all section data.

2125:   Not collective

2127:   Input Parameters:
2128: . s - the PetscSection

2130:   Level: beginner

2132: .seealso: PetscSection, PetscSectionCreate()
2133: @*/
2134: PetscErrorCode PetscSectionReset(PetscSection s)
2135: {
2136:   PetscInt       f, c;

2141:   for (f = 0; f < s->numFields; ++f) {
2142:     PetscSectionDestroy(&s->field[f]);
2143:     PetscFree(s->fieldNames[f]);
2144:     for (c = 0; c < s->numFieldComponents[f]; ++c)
2145:       PetscFree(s->compNames[f][c]);
2146:     PetscFree(s->compNames[f]);
2147:   }
2148:   PetscFree(s->numFieldComponents);
2149:   PetscFree(s->fieldNames);
2150:   PetscFree(s->compNames);
2151:   PetscFree(s->field);
2152:   PetscSectionDestroy(&s->bc);
2153:   PetscFree(s->bcIndices);
2154:   PetscFree2(s->atlasDof, s->atlasOff);
2155:   PetscSectionDestroy(&s->clSection);
2156:   ISDestroy(&s->clPoints);
2157:   ISDestroy(&s->perm);
2158:   PetscSectionResetClosurePermutation(s);
2159:   PetscSectionSymDestroy(&s->sym);
2160:   PetscSectionDestroy(&s->clSection);
2161:   ISDestroy(&s->clPoints);

2163:   s->pStart    = -1;
2164:   s->pEnd      = -1;
2165:   s->maxDof    = 0;
2166:   s->setup     = PETSC_FALSE;
2167:   s->numFields = 0;
2168:   s->clObj     = NULL;
2169:   return(0);
2170: }

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

2175:   Not collective

2177:   Input Parameters:
2178: . s - the PetscSection

2180:   Level: beginner

2182: .seealso: PetscSection, PetscSectionCreate()
2183: @*/
2184: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2185: {

2189:   if (!*s) return(0);
2191:   if (--((PetscObject)(*s))->refct > 0) {
2192:     *s = NULL;
2193:     return(0);
2194:   }
2195:   PetscSectionReset(*s);
2196:   PetscHeaderDestroy(s);
2197:   return(0);
2198: }

2200: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2201: {
2202:   const PetscInt p = point - s->pStart;

2206:   *values = &baseArray[s->atlasOff[p]];
2207:   return(0);
2208: }

2210: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2211: {
2212:   PetscInt       *array;
2213:   const PetscInt p           = point - s->pStart;
2214:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2215:   PetscInt       cDim        = 0;

2220:   PetscSectionGetConstraintDof(s, p, &cDim);
2221:   array = &baseArray[s->atlasOff[p]];
2222:   if (!cDim) {
2223:     if (orientation >= 0) {
2224:       const PetscInt dim = s->atlasDof[p];
2225:       PetscInt       i;

2227:       if (mode == INSERT_VALUES) {
2228:         for (i = 0; i < dim; ++i) array[i] = values[i];
2229:       } else {
2230:         for (i = 0; i < dim; ++i) array[i] += values[i];
2231:       }
2232:     } else {
2233:       PetscInt offset = 0;
2234:       PetscInt j      = -1, field, i;

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

2239:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2240:         offset += dim;
2241:       }
2242:     }
2243:   } else {
2244:     if (orientation >= 0) {
2245:       const PetscInt dim  = s->atlasDof[p];
2246:       PetscInt       cInd = 0, i;
2247:       const PetscInt *cDof;

2249:       PetscSectionGetConstraintIndices(s, point, &cDof);
2250:       if (mode == INSERT_VALUES) {
2251:         for (i = 0; i < dim; ++i) {
2252:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2253:           array[i] = values[i];
2254:         }
2255:       } else {
2256:         for (i = 0; i < dim; ++i) {
2257:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2258:           array[i] += values[i];
2259:         }
2260:       }
2261:     } else {
2262:       const PetscInt *cDof;
2263:       PetscInt       offset  = 0;
2264:       PetscInt       cOffset = 0;
2265:       PetscInt       j       = 0, field;

2267:       PetscSectionGetConstraintIndices(s, point, &cDof);
2268:       for (field = 0; field < s->numFields; ++field) {
2269:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2270:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2271:         const PetscInt sDim = dim - tDim;
2272:         PetscInt       cInd = 0, i,k;

2274:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2275:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2276:           array[j] = values[k];
2277:         }
2278:         offset  += dim;
2279:         cOffset += dim - tDim;
2280:       }
2281:     }
2282:   }
2283:   return(0);
2284: }

2286: /*@C
2287:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2289:   Not collective

2291:   Input Parameter:
2292: . s - The PetscSection

2294:   Output Parameter:
2295: . hasConstraints - flag indicating that the section has constrained dofs

2297:   Level: intermediate

2299: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2300: @*/
2301: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2302: {
2306:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2307:   return(0);
2308: }

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

2313:   Not collective

2315:   Input Parameters:
2316: + s     - The PetscSection
2317: - point - The point

2319:   Output Parameter:
2320: . indices - The constrained dofs

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

2324:   Level: intermediate

2326: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2327: @*/
2328: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2329: {

2334:   if (s->bc) {
2335:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2336:   } else *indices = NULL;
2337:   return(0);
2338: }

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

2343:   Not collective

2345:   Input Parameters:
2346: + s     - The PetscSection
2347: . point - The point
2348: - indices - The constrained dofs

2350:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2352:   Level: intermediate

2354: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2355: @*/
2356: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2357: {

2362:   if (s->bc) {
2363:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2364:   }
2365:   return(0);
2366: }

2368: /*@C
2369:   PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained

2371:   Not collective

2373:   Input Parameters:
2374: + s     - The PetscSection
2375: . field  - The field number
2376: - point - The point

2378:   Output Parameter:
2379: . indices - The constrained dofs sorted in ascending order

2381:   Notes:
2382:   The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof().

2384:   Fortran Note:
2385:   In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()

2387:   Level: intermediate

2389: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2390: @*/
2391: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2392: {

2397:   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);
2398:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2399:   return(0);
2400: }

2402: /*@C
2403:   PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained

2405:   Not collective

2407:   Input Parameters:
2408: + s       - The PetscSection
2409: . point   - The point
2410: . field   - The field number
2411: - indices - The constrained dofs

2413:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

2415:   Level: intermediate

2417: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2418: @*/
2419: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2420: {

2425:   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);
2426:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2427:   return(0);
2428: }

2430: /*@
2431:   PetscSectionPermute - Reorder the section according to the input point permutation

2433:   Collective on PetscSection

2435:   Input Parameter:
2436: + section - The PetscSection object
2437: - perm - The point permutation, old point p becomes new point perm[p]

2439:   Output Parameter:
2440: . sectionNew - The permuted PetscSection

2442:   Level: intermediate

2444: .seealso: MatPermute()
2445: @*/
2446: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2447: {
2448:   PetscSection    s = section, sNew;
2449:   const PetscInt *perm;
2450:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;
2451:   PetscErrorCode  ierr;

2457:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2458:   PetscSectionGetNumFields(s, &numFields);
2459:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2460:   for (f = 0; f < numFields; ++f) {
2461:     const char *name;
2462:     PetscInt    numComp;

2464:     PetscSectionGetFieldName(s, f, &name);
2465:     PetscSectionSetFieldName(sNew, f, name);
2466:     PetscSectionGetFieldComponents(s, f, &numComp);
2467:     PetscSectionSetFieldComponents(sNew, f, numComp);
2468:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2469:       PetscSectionGetComponentName(s, f, c, &name);
2470:       PetscSectionSetComponentName(sNew, f, c, name);
2471:     }
2472:   }
2473:   ISGetLocalSize(permutation, &numPoints);
2474:   ISGetIndices(permutation, &perm);
2475:   PetscSectionGetChart(s, &pStart, &pEnd);
2476:   PetscSectionSetChart(sNew, pStart, pEnd);
2477:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2478:   for (p = pStart; p < pEnd; ++p) {
2479:     PetscInt dof, cdof;

2481:     PetscSectionGetDof(s, p, &dof);
2482:     PetscSectionSetDof(sNew, perm[p], dof);
2483:     PetscSectionGetConstraintDof(s, p, &cdof);
2484:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2485:     for (f = 0; f < numFields; ++f) {
2486:       PetscSectionGetFieldDof(s, p, f, &dof);
2487:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2488:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2489:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2490:     }
2491:   }
2492:   PetscSectionSetUp(sNew);
2493:   for (p = pStart; p < pEnd; ++p) {
2494:     const PetscInt *cind;
2495:     PetscInt        cdof;

2497:     PetscSectionGetConstraintDof(s, p, &cdof);
2498:     if (cdof) {
2499:       PetscSectionGetConstraintIndices(s, p, &cind);
2500:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2501:     }
2502:     for (f = 0; f < numFields; ++f) {
2503:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2504:       if (cdof) {
2505:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2506:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2507:       }
2508:     }
2509:   }
2510:   ISRestoreIndices(permutation, &perm);
2511:   *sectionNew = sNew;
2512:   return(0);
2513: }

2515: /* TODO: the next three functions should be moved to sf/utils */
2516: #include <petsc/private/sfimpl.h>

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

2521:   Collective on sf

2523:   Input Parameters:
2524: + sf - The SF
2525: - rootSection - Section defined on root space

2527:   Output Parameters:
2528: + remoteOffsets - root offsets in leaf storage, or NULL
2529: - leafSection - Section defined on the leaf space

2531:   Level: advanced

2533: .seealso: PetscSFCreate()
2534: @*/
2535: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2536: {
2537:   PetscSF        embedSF;
2538:   const PetscInt *indices;
2539:   IS             selected;
2540:   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
2541:   PetscBool      *sub, hasc;

2545:   PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2546:   PetscSectionGetNumFields(rootSection, &numFields);
2547:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2548:   PetscMalloc1(numFields+2, &sub);
2549:   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2550:   for (f = 0; f < numFields; ++f) {
2551:     PetscSectionSym sym;
2552:     const char      *name   = NULL;
2553:     PetscInt        numComp = 0;

2555:     sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2556:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
2557:     PetscSectionGetFieldName(rootSection, f, &name);
2558:     PetscSectionGetFieldSym(rootSection, f, &sym);
2559:     PetscSectionSetFieldComponents(leafSection, f, numComp);
2560:     PetscSectionSetFieldName(leafSection, f, name);
2561:     PetscSectionSetFieldSym(leafSection, f, sym);
2562:     for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
2563:       PetscSectionGetComponentName(rootSection, f, c, &name);
2564:       PetscSectionSetComponentName(leafSection, f, c, name);
2565:     }
2566:   }
2567:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2568:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2569:   rpEnd = PetscMin(rpEnd,nroots);
2570:   rpEnd = PetscMax(rpStart,rpEnd);
2571:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2572:   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2573:   MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
2574:   if (sub[0]) {
2575:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2576:     ISGetIndices(selected, &indices);
2577:     PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2578:     ISRestoreIndices(selected, &indices);
2579:     ISDestroy(&selected);
2580:   } else {
2581:     PetscObjectReference((PetscObject)sf);
2582:     embedSF = sf;
2583:   }
2584:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2585:   lpEnd++;

2587:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

2589:   /* Constrained dof section */
2590:   hasc = sub[1];
2591:   for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);

2593:   /* Could fuse these at the cost of copies and extra allocation */
2594:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2595:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2596:   if (sub[1]) {
2597:     PetscSectionCheckConstraints_Static(rootSection);
2598:     PetscSectionCheckConstraints_Static(leafSection);
2599:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2600:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2601:   }
2602:   for (f = 0; f < numFields; ++f) {
2603:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2604:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2605:     if (sub[2+f]) {
2606:       PetscSectionCheckConstraints_Static(rootSection->field[f]);
2607:       PetscSectionCheckConstraints_Static(leafSection->field[f]);
2608:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2609:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2610:     }
2611:   }
2612:   if (remoteOffsets) {
2613:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2614:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2615:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2616:   }
2617:   PetscSectionSetUp(leafSection);
2618:   if (hasc) { /* need to communicate bcIndices */
2619:     PetscSF  bcSF;
2620:     PetscInt *rOffBc;

2622:     PetscMalloc1(lpEnd - lpStart, &rOffBc);
2623:     if (sub[1]) {
2624:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2625:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2626:       PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
2627:       PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2628:       PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2629:       PetscSFDestroy(&bcSF);
2630:     }
2631:     for (f = 0; f < numFields; ++f) {
2632:       if (sub[2+f]) {
2633:         PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2634:         PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2635:         PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
2636:         PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2637:         PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2638:         PetscSFDestroy(&bcSF);
2639:       }
2640:     }
2641:     PetscFree(rOffBc);
2642:   }
2643:   PetscSFDestroy(&embedSF);
2644:   PetscFree(sub);
2645:   PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2646:   return(0);
2647: }

2649: /*@C
2650:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

2652:   Collective on sf

2654:   Input Parameters:
2655: + sf          - The SF
2656: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2657: - leafSection - Data layout of local points for incoming data  (this is layout for SF leaves)

2659:   Output Parameter:
2660: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL

2662:   Level: developer

2664: .seealso: PetscSFCreate()
2665: @*/
2666: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2667: {
2668:   PetscSF         embedSF;
2669:   const PetscInt *indices;
2670:   IS              selected;
2671:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2672:   PetscErrorCode  ierr;

2675:   *remoteOffsets = NULL;
2676:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2677:   if (numRoots < 0) return(0);
2678:   PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2679:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2680:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2681:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2682:   ISGetIndices(selected, &indices);
2683:   PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2684:   ISRestoreIndices(selected, &indices);
2685:   ISDestroy(&selected);
2686:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2687:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2688:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2689:   PetscSFDestroy(&embedSF);
2690:   PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2691:   return(0);
2692: }

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

2697:   Collective on sf

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

2705:   Output Parameters:
2706: - sectionSF - The new SF

2708:   Note: Either rootSection or remoteOffsets can be specified

2710:   Level: advanced

2712: .seealso: PetscSFCreate()
2713: @*/
2714: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2715: {
2716:   MPI_Comm          comm;
2717:   const PetscInt    *localPoints;
2718:   const PetscSFNode *remotePoints;
2719:   PetscInt          lpStart, lpEnd;
2720:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
2721:   PetscInt          *localIndices;
2722:   PetscSFNode       *remoteIndices;
2723:   PetscInt          i, ind;
2724:   PetscErrorCode    ierr;

2732:   PetscObjectGetComm((PetscObject)sf,&comm);
2733:   PetscSFCreate(comm, sectionSF);
2734:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2735:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2736:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2737:   if (numRoots < 0) return(0);
2738:   PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2739:   for (i = 0; i < numPoints; ++i) {
2740:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2741:     PetscInt dof;

2743:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2744:       PetscSectionGetDof(leafSection, localPoint, &dof);
2745:       numIndices += dof;
2746:     }
2747:   }
2748:   PetscMalloc1(numIndices, &localIndices);
2749:   PetscMalloc1(numIndices, &remoteIndices);
2750:   /* Create new index graph */
2751:   for (i = 0, ind = 0; i < numPoints; ++i) {
2752:     PetscInt localPoint = localPoints ? localPoints[i] : i;
2753:     PetscInt rank       = remotePoints[i].rank;

2755:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2756:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2757:       PetscInt loff, dof, d;

2759:       PetscSectionGetOffset(leafSection, localPoint, &loff);
2760:       PetscSectionGetDof(leafSection, localPoint, &dof);
2761:       for (d = 0; d < dof; ++d, ++ind) {
2762:         localIndices[ind]        = loff+d;
2763:         remoteIndices[ind].rank  = rank;
2764:         remoteIndices[ind].index = remoteOffset+d;
2765:       }
2766:     }
2767:   }
2768:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2769:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2770:   PetscSFSetUp(*sectionSF);
2771:   PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2772:   return(0);
2773: }

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

2778:   Collective on section

2780:   Input Parameters:
2781: + section   - The PetscSection
2782: . obj       - A PetscObject which serves as the key for this index
2783: . clSection - Section giving the size of the closure of each point
2784: - clPoints  - IS giving the points in each closure

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

2788:   Level: advanced

2790: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2791: @*/
2792: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2793: {

2800:   if (section->clObj != obj) {PetscSectionResetClosurePermutation(section);}
2801:   section->clObj     = obj;
2802:   PetscObjectReference((PetscObject)clSection);
2803:   PetscObjectReference((PetscObject)clPoints);
2804:   PetscSectionDestroy(&section->clSection);
2805:   ISDestroy(&section->clPoints);
2806:   section->clSection = clSection;
2807:   section->clPoints  = clPoints;
2808:   return(0);
2809: }

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

2814:   Collective on section

2816:   Input Parameters:
2817: + section   - The PetscSection
2818: - obj       - A PetscObject which serves as the key for this index

2820:   Output Parameters:
2821: + clSection - Section giving the size of the closure of each point
2822: - clPoints  - IS giving the points in each closure

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

2826:   Level: advanced

2828: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2829: @*/
2830: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2831: {
2833:   if (section->clObj == obj) {
2834:     if (clSection) *clSection = section->clSection;
2835:     if (clPoints)  *clPoints  = section->clPoints;
2836:   } else {
2837:     if (clSection) *clSection = NULL;
2838:     if (clPoints)  *clPoints  = NULL;
2839:   }
2840:   return(0);
2841: }

2843: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2844: {
2845:   PetscInt       i;
2846:   khiter_t iter;
2847:   int new_entry;
2848:   PetscSectionClosurePermKey key = {depth, clSize};
2849:   PetscSectionClosurePermVal *val;

2853:   if (section->clObj != obj) {
2854:     PetscSectionDestroy(&section->clSection);
2855:     ISDestroy(&section->clPoints);
2856:   }
2857:   section->clObj = obj;
2858:   if (!section->clHash) {PetscClPermCreate(&section->clHash);}
2859:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2860:   val = &kh_val(section->clHash, iter);
2861:   if (!new_entry) {
2862:     PetscFree(val->perm);
2863:     PetscFree(val->invPerm);
2864:   }
2865:   if (mode == PETSC_COPY_VALUES) {
2866:     PetscMalloc1(clSize, &val->perm);
2867:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2868:     PetscArraycpy(val->perm, clPerm, clSize);
2869:   } else if (mode == PETSC_OWN_POINTER) {
2870:     val->perm = clPerm;
2871:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2872:   PetscMalloc1(clSize, &val->invPerm);
2873:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2874:   return(0);
2875: }

2877: /*@
2878:   PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.

2880:   Not Collective

2882:   Input Parameters:
2883: + section - The PetscSection
2884: . obj     - A PetscObject which serves as the key for this index (usually a DM)
2885: . depth   - Depth of points on which to apply the given permutation
2886: - perm    - Permutation of the cell dof closure

2888:   Note:
2889:   The specified permutation will only be applied to points at depth whose closure size matches the length of perm.  In a
2890:   mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2891:   each topology and degree.

2893:   This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2894:   exotic/enriched spaces on mixed topology meshes.

2896:   Level: intermediate

2898: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2899: @*/
2900: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2901: {
2902:   const PetscInt *clPerm = NULL;
2903:   PetscInt        clSize = 0;
2904:   PetscErrorCode  ierr;

2907:   if (perm) {
2908:     ISGetLocalSize(perm, &clSize);
2909:     ISGetIndices(perm, &clPerm);
2910:   }
2911:   PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2912:   if (perm) {ISRestoreIndices(perm, &clPerm);}
2913:   return(0);
2914: }

2916: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2917: {

2921:   if (section->clObj == obj) {
2922:     PetscSectionClosurePermKey k = {depth, size};
2923:     PetscSectionClosurePermVal v;
2924:     PetscClPermGet(section->clHash, k, &v);
2925:     if (perm) *perm = v.perm;
2926:   } else {
2927:     if (perm) *perm = NULL;
2928:   }
2929:   return(0);
2930: }

2932: /*@
2933:   PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.

2935:   Not collective

2937:   Input Parameters:
2938: + section   - The PetscSection
2939: . obj       - A PetscObject which serves as the key for this index (usually a DM)
2940: . depth     - Depth stratum on which to obtain closure permutation
2941: - clSize    - Closure size to be permuted (e.g., may vary with element topology and degree)

2943:   Output Parameter:
2944: . perm - The dof closure permutation

2946:   Note:
2947:   The user must destroy the IS that is returned.

2949:   Level: intermediate

2951: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2952: @*/
2953: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2954: {
2955:   const PetscInt *clPerm;
2956:   PetscErrorCode  ierr;

2959:   PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2960:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2961:   return(0);
2962: }

2964: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2965: {

2969:   if (section->clObj == obj && section->clHash) {
2970:     PetscSectionClosurePermKey k = {depth, size};
2971:     PetscSectionClosurePermVal v;
2972:     PetscClPermGet(section->clHash, k, &v);
2973:     if (perm) *perm = v.invPerm;
2974:   } else {
2975:     if (perm) *perm = NULL;
2976:   }
2977:   return(0);
2978: }

2980: /*@
2981:   PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.

2983:   Not collective

2985:   Input Parameters:
2986: + section   - The PetscSection
2987: . obj       - A PetscObject which serves as the key for this index (usually a DM)
2988: . depth     - Depth stratum on which to obtain closure permutation
2989: - clSize    - Closure size to be permuted (e.g., may vary with element topology and degree)

2991:   Output Parameters:
2992: . perm - The dof closure permutation

2994:   Note:
2995:   The user must destroy the IS that is returned.

2997:   Level: intermediate

2999: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
3000: @*/
3001: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3002: {
3003:   const PetscInt *clPerm;
3004:   PetscErrorCode  ierr;

3007:   PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
3008:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
3009:   return(0);
3010: }

3012: /*@
3013:   PetscSectionGetField - Get the subsection associated with a single field

3015:   Input Parameters:
3016: + s     - The PetscSection
3017: - field - The field number

3019:   Output Parameter:
3020: . subs  - The subsection for the given field

3022:   Level: intermediate

3024: .seealso: PetscSectionSetNumFields()
3025: @*/
3026: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3027: {
3031:   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);
3032:   *subs = s->field[field];
3033:   return(0);
3034: }

3036: PetscClassId      PETSC_SECTION_SYM_CLASSID;
3037: PetscFunctionList PetscSectionSymList = NULL;

3039: /*@
3040:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

3042:   Collective

3044:   Input Parameter:
3045: . comm - the MPI communicator

3047:   Output Parameter:
3048: . sym - pointer to the new set of symmetries

3050:   Level: developer

3052: .seealso: PetscSectionSym, PetscSectionSymDestroy()
3053: @*/
3054: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3055: {

3060:   ISInitializePackage();
3061:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
3062:   return(0);
3063: }

3065: /*@C
3066:   PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.

3068:   Collective on PetscSectionSym

3070:   Input Parameters:
3071: + sym    - The section symmetry object
3072: - method - The name of the section symmetry type

3074:   Level: developer

3076: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
3077: @*/
3078: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3079: {
3080:   PetscErrorCode (*r)(PetscSectionSym);
3081:   PetscBool      match;

3086:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
3087:   if (match) return(0);

3089:   PetscFunctionListFind(PetscSectionSymList,method,&r);
3090:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3091:   if (sym->ops->destroy) {
3092:     (*sym->ops->destroy)(sym);
3093:     sym->ops->destroy = NULL;
3094:   }
3095:   (*r)(sym);
3096:   PetscObjectChangeTypeName((PetscObject)sym,method);
3097:   return(0);
3098: }


3101: /*@C
3102:   PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.

3104:   Not Collective

3106:   Input Parameter:
3107: . sym  - The section symmetry

3109:   Output Parameter:
3110: . type - The index set type name

3112:   Level: developer

3114: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
3115: @*/
3116: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3117: {
3121:   *type = ((PetscObject)sym)->type_name;
3122:   return(0);
3123: }

3125: /*@C
3126:   PetscSectionSymRegister - Adds a new section symmetry implementation

3128:   Not Collective

3130:   Input Parameters:
3131: + name        - The name of a new user-defined creation routine
3132: - create_func - The creation routine itself

3134:   Notes:
3135:   PetscSectionSymRegister() may be called multiple times to add several user-defined vectors

3137:   Level: developer

3139: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
3140: @*/
3141: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3142: {

3146:   ISInitializePackage();
3147:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
3148:   return(0);
3149: }

3151: /*@
3152:    PetscSectionSymDestroy - Destroys a section symmetry.

3154:    Collective on PetscSectionSym

3156:    Input Parameters:
3157: .  sym - the section symmetry

3159:    Level: developer

3161: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3162: @*/
3163: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3164: {
3165:   SymWorkLink    link,next;

3169:   if (!*sym) return(0);
3171:   if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return(0);}
3172:   if ((*sym)->ops->destroy) {
3173:     (*(*sym)->ops->destroy)(*sym);
3174:   }
3175:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3176:   for (link=(*sym)->workin; link; link=next) {
3177:     next = link->next;
3178:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3179:     PetscFree(link);
3180:   }
3181:   (*sym)->workin = NULL;
3182:   PetscHeaderDestroy(sym);
3183:   return(0);
3184: }

3186: /*@C
3187:    PetscSectionSymView - Displays a section symmetry

3189:    Collective on PetscSectionSym

3191:    Input Parameters:
3192: +  sym - the index set
3193: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

3195:    Level: developer

3197: .seealso: PetscViewerASCIIOpen()
3198: @*/
3199: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3200: {

3205:   if (!viewer) {
3206:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3207:   }
3210:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3211:   if (sym->ops->view) {
3212:     (*sym->ops->view)(sym,viewer);
3213:   }
3214:   return(0);
3215: }

3217: /*@
3218:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3220:   Collective on PetscSection

3222:   Input Parameters:
3223: + section - the section describing data layout
3224: - sym - the symmetry describing the affect of orientation on the access of the data

3226:   Level: developer

3228: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3229: @*/
3230: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3231: {

3236:   PetscSectionSymDestroy(&(section->sym));
3237:   if (sym) {
3240:     PetscObjectReference((PetscObject) sym);
3241:   }
3242:   section->sym = sym;
3243:   return(0);
3244: }

3246: /*@
3247:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3249:   Not collective

3251:   Input Parameters:
3252: . section - the section describing data layout

3254:   Output Parameters:
3255: . sym - the symmetry describing the affect of orientation on the access of the data

3257:   Level: developer

3259: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3260: @*/
3261: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3262: {
3265:   *sym = section->sym;
3266:   return(0);
3267: }

3269: /*@
3270:   PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section

3272:   Collective on PetscSection

3274:   Input Parameters:
3275: + section - the section describing data layout
3276: . field - the field number
3277: - sym - the symmetry describing the affect of orientation on the access of the data

3279:   Level: developer

3281: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3282: @*/
3283: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3284: {

3289:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
3290:   PetscSectionSetSym(section->field[field],sym);
3291:   return(0);
3292: }

3294: /*@
3295:   PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section

3297:   Collective on PetscSection

3299:   Input Parameters:
3300: + section - the section describing data layout
3301: - field - the field number

3303:   Output Parameters:
3304: . sym - the symmetry describing the affect of orientation on the access of the data

3306:   Level: developer

3308: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3309: @*/
3310: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3311: {
3314:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
3315:   *sym = section->field[field]->sym;
3316:   return(0);
3317: }

3319: /*@C
3320:   PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.

3322:   Not collective

3324:   Input Parameters:
3325: + section - the section
3326: . numPoints - the number of points
3327: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3328:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3329:     context, see DMPlexGetConeOrientation()).

3331:   Output Parameter:
3332: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3333: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3334:     identity).

3336:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3337: .vb
3338:      const PetscInt    **perms;
3339:      const PetscScalar **rots;
3340:      PetscInt            lOffset;

3342:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3343:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3344:        PetscInt           point = points[2*i], dof, sOffset;
3345:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3346:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3348:        PetscSectionGetDof(section,point,&dof);
3349:        PetscSectionGetOffset(section,point,&sOffset);

3351:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3352:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3353:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3354:        lOffset += dof;
3355:      }
3356:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3357: .ve

3359:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3360: .vb
3361:      const PetscInt    **perms;
3362:      const PetscScalar **rots;
3363:      PetscInt            lOffset;

3365:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3366:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3367:        PetscInt           point = points[2*i], dof, sOffset;
3368:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3369:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3371:        PetscSectionGetDof(section,point,&dof);
3372:        PetscSectionGetOffset(section,point,&sOff);

3374:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3375:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3376:        offset += dof;
3377:      }
3378:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3379: .ve

3381:   Level: developer

3383: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3384: @*/
3385: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3386: {
3387:   PetscSectionSym sym;
3388:   PetscErrorCode  ierr;

3393:   if (perms) *perms = NULL;
3394:   if (rots)  *rots  = NULL;
3395:   sym = section->sym;
3396:   if (sym && (perms || rots)) {
3397:     SymWorkLink link;

3399:     if (sym->workin) {
3400:       link        = sym->workin;
3401:       sym->workin = sym->workin->next;
3402:     } else {
3403:       PetscNewLog(sym,&link);
3404:     }
3405:     if (numPoints > link->numPoints) {
3406:       PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3407:       PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3408:       link->numPoints = numPoints;
3409:     }
3410:     link->next   = sym->workout;
3411:     sym->workout = link;
3412:     PetscArrayzero((PetscInt**)link->perms,numPoints);
3413:     PetscArrayzero((PetscInt**)link->rots,numPoints);
3414:     (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3415:     if (perms) *perms = link->perms;
3416:     if (rots)  *rots  = link->rots;
3417:   }
3418:   return(0);
3419: }

3421: /*@C
3422:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3424:   Not collective

3426:   Input Parameters:
3427: + section - the section
3428: . numPoints - the number of points
3429: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3430:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3431:     context, see DMPlexGetConeOrientation()).

3433:   Output Parameter:
3434: + perms - The permutations for the given orientations: set to NULL at conclusion
3435: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3437:   Level: developer

3439: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3440: @*/
3441: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3442: {
3443:   PetscSectionSym sym;

3447:   sym = section->sym;
3448:   if (sym && (perms || rots)) {
3449:     SymWorkLink *p,link;

3451:     for (p=&sym->workout; (link=*p); p=&link->next) {
3452:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3453:         *p          = link->next;
3454:         link->next  = sym->workin;
3455:         sym->workin = link;
3456:         if (perms) *perms = NULL;
3457:         if (rots)  *rots  = NULL;
3458:         return(0);
3459:       }
3460:     }
3461:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3462:   }
3463:   return(0);
3464: }

3466: /*@C
3467:   PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.

3469:   Not collective

3471:   Input Parameters:
3472: + section - the section
3473: . field - the field of the section
3474: . numPoints - the number of points
3475: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3476:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3477:     context, see DMPlexGetConeOrientation()).

3479:   Output Parameter:
3480: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3481: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3482:     identity).

3484:   Level: developer

3486: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3487: @*/
3488: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3489: {

3494:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
3495:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3496:   return(0);
3497: }

3499: /*@C
3500:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3502:   Not collective

3504:   Input Parameters:
3505: + section - the section
3506: . field - the field number
3507: . numPoints - the number of points
3508: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3509:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3510:     context, see DMPlexGetConeOrientation()).

3512:   Output Parameter:
3513: + perms - The permutations for the given orientations: set to NULL at conclusion
3514: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3516:   Level: developer

3518: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3519: @*/
3520: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3521: {

3526:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
3527:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3528:   return(0);
3529: }

3531: /*@
3532:   PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset

3534:   Not collective

3536:   Input Parameter:
3537: . s - the global PetscSection

3539:   Output Parameters:
3540: . flg - the flag

3542:   Level: developer

3544: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3545: @*/
3546: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3547: {
3550:   *flg = s->useFieldOff;
3551:   return(0);
3552: }

3554: /*@
3555:   PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset

3557:   Not collective

3559:   Input Parameters:
3560: + s   - the global PetscSection
3561: - flg - the flag

3563:   Level: developer

3565: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3566: @*/
3567: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3568: {
3571:   s->useFieldOff = flg;
3572:   return(0);
3573: }

3575: #define PetscSectionExpandPoints_Loop(TYPE) \
3576: { \
3577:   PetscInt i, n, o0, o1, size; \
3578:   TYPE *a0 = (TYPE*)origArray, *a1; \
3579:   PetscSectionGetStorageSize(s, &size); \
3580:   PetscMalloc1(size, &a1); \
3581:   for (i=0; i<npoints; i++) { \
3582:     PetscSectionGetOffset(origSection, points_[i], &o0); \
3583:     PetscSectionGetOffset(s, i, &o1); \
3584:     PetscSectionGetDof(s, i, &n); \
3585:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3586:   } \
3587:   *newArray = (void*)a1; \
3588: }

3590: /*@
3591:   PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.

3593:   Not collective

3595:   Input Parameters:
3596: + origSection - the PetscSection describing the layout of the array
3597: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3598: . origArray - the array; its size must be equal to the storage size of origSection
3599: - points - IS with points to extract; its indices must lie in the chart of origSection

3601:   Output Parameters:
3602: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3603: - newArray - the array of the extracted DOFs; its size is the storage size of newSection

3605:   Level: developer

3607: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3608: @*/
3609: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3610: {
3611:   PetscSection        s;
3612:   const PetscInt      *points_;
3613:   PetscInt            i, n, npoints, pStart, pEnd;
3614:   PetscMPIInt         unitsize;
3615:   PetscErrorCode      ierr;

3623:   MPI_Type_size(dataType, &unitsize);
3624:   ISGetLocalSize(points, &npoints);
3625:   ISGetIndices(points, &points_);
3626:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3627:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3628:   PetscSectionSetChart(s, 0, npoints);
3629:   for (i=0; i<npoints; i++) {
3630:     if (PetscUnlikely(points_[i] < pStart || points_[i] >= pEnd)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %d (index %d) in input IS out of input section's chart", points_[i], i);
3631:     PetscSectionGetDof(origSection, points_[i], &n);
3632:     PetscSectionSetDof(s, i, n);
3633:   }
3634:   PetscSectionSetUp(s);
3635:   if (newArray) {
3636:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
3637:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
3638:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
3639:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3640:   }
3641:   if (newSection) {
3642:     *newSection = s;
3643:   } else {
3644:     PetscSectionDestroy(&s);
3645:   }
3646:   return(0);
3647: }