Actual source code: section.c

  1: /*
  2:    This file contains routines for basic section object implementation.
  3: */

  5: #include <petsc/private/sectionimpl.h>
  6: #include <petscsf.h>

  8: PetscClassId PETSC_SECTION_CLASSID;

 10: /*@
 11:   PetscSectionCreate - Allocates `PetscSection` and sets the map contents to the default.

 13:   Collective

 15:   Input Parameters:
 16: + comm - the MPI communicator
 17: - s    - pointer to the section

 19:   Level: beginner

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

 33:   The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.

 35: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
 36: @*/
 37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
 38: {
 40:   ISInitializePackage();

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

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

 67: /*@
 68:   PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`

 70:   Collective

 72:   Input Parameter:
 73: . section - the `PetscSection`

 75:   Output Parameter:
 76: . newSection - the copy

 78:   Level: intermediate

 80:   Developer Note:
 81:   What exactly does shallow mean in this context?

 83: .seealso: [PetscSection](sec_petscsection), `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;

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

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

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

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

157: /*@
158:   PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`

160:   Collective

162:   Input Parameter:
163: . section - the `PetscSection`

165:   Output Parameter:
166: . newSection - the copy

168:   Level: beginner

170:   Developer Note:
171:   With standard PETSc terminology this should be called `PetscSectionDuplicate()`

173: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
174: @*/
175: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
176: {
179:   PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection);
180:   PetscSectionCopy(section, *newSection);
181:   return 0;
182: }

184: /*@
185:   PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database

187:   Collective

189:   Input Parameter:
190: . section - the `PetscSection`

192:   Options Database:
193: . -petscsection_point_major - `PETSC_TRUE` for point-major order

195:   Level: intermediate

197: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
198: @*/
199: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
200: {
202:   PetscObjectOptionsBegin((PetscObject)s);
203:   PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
204:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
205:   PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject);
206:   PetscOptionsEnd();
207:   PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view");
208:   return 0;
209: }

211: /*@
212:   PetscSectionCompare - Compares two sections

214:   Collective

216:   Input Parameters:
217: + s1 - the first `PetscSection`
218: - s2 - the second `PetscSection`

220:   Output Parameter:
221: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise

223:   Level: intermediate

225:   Note:
226:   Field names are disregarded.

228: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
229: @*/
230: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
231: {
232:   PetscInt        pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
233:   const PetscInt *idx1, *idx2;
234:   IS              perm1, perm2;
235:   PetscBool       flg;
236:   PetscMPIInt     mflg;

241:   flg = PETSC_FALSE;

243:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg);
244:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
245:     *congruent = PETSC_FALSE;
246:     return 0;
247:   }

249:   PetscSectionGetChart(s1, &pStart, &pEnd);
250:   PetscSectionGetChart(s2, &n1, &n2);
251:   if (pStart != n1 || pEnd != n2) goto not_congruent;

253:   PetscSectionGetPermutation(s1, &perm1);
254:   PetscSectionGetPermutation(s2, &perm2);
255:   if (perm1 && perm2) {
256:     ISEqual(perm1, perm2, congruent);
257:     if (!(*congruent)) goto not_congruent;
258:   } else if (perm1 != perm2) goto not_congruent;

260:   for (p = pStart; p < pEnd; ++p) {
261:     PetscSectionGetOffset(s1, p, &n1);
262:     PetscSectionGetOffset(s2, p, &n2);
263:     if (n1 != n2) goto not_congruent;

265:     PetscSectionGetDof(s1, p, &n1);
266:     PetscSectionGetDof(s2, p, &n2);
267:     if (n1 != n2) goto not_congruent;

269:     PetscSectionGetConstraintDof(s1, p, &ncdof);
270:     PetscSectionGetConstraintDof(s2, p, &n2);
271:     if (ncdof != n2) goto not_congruent;

273:     PetscSectionGetConstraintIndices(s1, p, &idx1);
274:     PetscSectionGetConstraintIndices(s2, p, &idx2);
275:     PetscArraycmp(idx1, idx2, ncdof, congruent);
276:     if (!(*congruent)) goto not_congruent;
277:   }

279:   PetscSectionGetNumFields(s1, &nfields);
280:   PetscSectionGetNumFields(s2, &n2);
281:   if (nfields != n2) goto not_congruent;

283:   for (f = 0; f < nfields; ++f) {
284:     PetscSectionGetFieldComponents(s1, f, &n1);
285:     PetscSectionGetFieldComponents(s2, f, &n2);
286:     if (n1 != n2) goto not_congruent;

288:     for (p = pStart; p < pEnd; ++p) {
289:       PetscSectionGetFieldOffset(s1, p, f, &n1);
290:       PetscSectionGetFieldOffset(s2, p, f, &n2);
291:       if (n1 != n2) goto not_congruent;

293:       PetscSectionGetFieldDof(s1, p, f, &n1);
294:       PetscSectionGetFieldDof(s2, p, f, &n2);
295:       if (n1 != n2) goto not_congruent;

297:       PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
298:       PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
299:       if (nfcdof != n2) goto not_congruent;

301:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
302:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
303:       PetscArraycmp(idx1, idx2, nfcdof, congruent);
304:       if (!(*congruent)) goto not_congruent;
305:     }
306:   }

308:   flg = PETSC_TRUE;
309: not_congruent:
310:   MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1));
311:   return 0;
312: }

314: /*@
315:   PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.

317:   Not Collective

319:   Input Parameter:
320: . s - the `PetscSection`

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

325:   Level: intermediate

327: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
328: @*/
329: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
330: {
333:   *numFields = s->numFields;
334:   return 0;
335: }

337: /*@
338:   PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`

340:   Not Collective

342:   Input Parameters:
343: + s - the PetscSection
344: - numFields - the number of fields

346:   Level: intermediate

348: .seealso: [PetscSection](sec_petscsection), `PetscSection()`, `PetscSectionGetNumFields()`
349: @*/
350: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
351: {
352:   PetscInt f;

356:   PetscSectionReset(s);

358:   s->numFields = numFields;
359:   PetscMalloc1(s->numFields, &s->numFieldComponents);
360:   PetscMalloc1(s->numFields, &s->fieldNames);
361:   PetscMalloc1(s->numFields, &s->compNames);
362:   PetscMalloc1(s->numFields, &s->field);
363:   for (f = 0; f < s->numFields; ++f) {
364:     char name[64];

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

368:     PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]);
369:     PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f);
370:     PetscStrallocpy(name, (char **)&s->fieldNames[f]);
371:     PetscSNPrintf(name, 64, "Component_0");
372:     PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
373:     PetscStrallocpy(name, (char **)&s->compNames[f][0]);
374:   }
375:   return 0;
376: }

378: /*@C
379:   PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`

381:   Not Collective

383:   Input Parameters:
384: + s     - the `PetscSection`
385: - field - the field number

387:   Output Parameter:
388: . fieldName - the field name

390:   Level: intermediate

392:   Note:
393:   Will error if the field number is out of range

395: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
396: @*/
397: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
398: {
401:   PetscSectionCheckValidField(field, s->numFields);
402:   *fieldName = s->fieldNames[field];
403:   return 0;
404: }

406: /*@C
407:   PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`

409:   Not Collective

411:   Input Parameters:
412: + s     - the `PetscSection`
413: . field - the field number
414: - fieldName - the field name

416:   Level: intermediate

418:   Note:
419:   Will error if the field number is out of range

421: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
422: @*/
423: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
424: {
427:   PetscSectionCheckValidField(field, s->numFields);
428:   PetscFree(s->fieldNames[field]);
429:   PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]);
430:   return 0;
431: }

433: /*@C
434:   PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`

436:   Not Collective

438:   Input Parameters:
439: + s     - the `PetscSection`
440: . field - the field number
441: - comp  - the component number

443:   Output Parameter:
444: . compName - the component name

446:   Level: intermediate

448:   Note:
449:   Will error if the field or component number do not exist

451: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
452:           `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
453: @*/
454: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
455: {
458:   PetscSectionCheckValidField(field, s->numFields);
459:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
460:   *compName = s->compNames[field][comp];
461:   return 0;
462: }

464: /*@C
465:   PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`

467:   Not Collective

469:   Input Parameters:
470: + s     - the PetscSection
471: . field - the field number
472: . comp  - the component number
473: - compName - the component name

475:   Level: intermediate

477:   Note:
478:   Will error if the field or component number do not exist

480: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
481:           `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
482: @*/
483: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
484: {
487:   PetscSectionCheckValidField(field, s->numFields);
488:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
489:   PetscFree(s->compNames[field][comp]);
490:   PetscStrallocpy(compName, (char **)&s->compNames[field][comp]);
491:   return 0;
492: }

494: /*@
495:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

497:   Not Collective

499:   Input Parameters:
500: + s - the `PetscSection`
501: - field - the field number

503:   Output Parameter:
504: . numComp - the number of field components

506:   Level: intermediate

508:   Developer Note:
509:   This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name

511: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`
512: @*/
513: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
514: {
517:   PetscSectionCheckValidField(field, 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: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionGetNumFields()`
535: @*/
536: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
537: {
538:   PetscInt c;

541:   PetscSectionCheckValidField(field, s->numFields);
542:   if (s->compNames) {
543:     for (c = 0; c < s->numFieldComponents[field]; ++c) PetscFree(s->compNames[field][c]);
544:     PetscFree(s->compNames[field]);
545:   }

547:   s->numFieldComponents[field] = numComp;
548:   if (numComp) {
549:     PetscMalloc1(numComp, (char ***)&s->compNames[field]);
550:     for (c = 0; c < numComp; ++c) {
551:       char name[64];

553:       PetscSNPrintf(name, 64, "%" PetscInt_FMT, c);
554:       PetscStrallocpy(name, (char **)&s->compNames[field][c]);
555:     }
556:   }
557:   return 0;
558: }

560: /*@
561:   PetscSectionGetChart - Returns the range [pStart, pEnd) in which points (indices) lie for this `PetscSection`

563:   Not Collective

565:   Input Parameter:
566: . s - the `PetscSection`

568:   Output Parameters:
569: + pStart - the first point
570: - pEnd - one past the last point

572:   Level: intermediate

574: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
575: @*/
576: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
577: {
579:   if (pStart) *pStart = s->pStart;
580:   if (pEnd) *pEnd = s->pEnd;
581:   return 0;
582: }

584: /*@
585:   PetscSectionSetChart - Sets the range [pStart, pEnd) in which points (indices) lie for this `PetscSection`

587:   Not Collective

589:   Input Parameters:
590: + s - the PetscSection
591: . pStart - the first point
592: - pEnd - one past the last point

594:   Level: intermediate

596: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`
597: @*/
598: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
599: {
600:   PetscInt f;

603:   if (pStart == s->pStart && pEnd == s->pEnd) return 0;
604:   /* Cannot Reset() because it destroys field information */
605:   s->setup = PETSC_FALSE;
606:   PetscSectionDestroy(&s->bc);
607:   PetscFree(s->bcIndices);
608:   PetscFree2(s->atlasDof, s->atlasOff);

610:   s->pStart = pStart;
611:   s->pEnd   = pEnd;
612:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
613:   PetscArrayzero(s->atlasDof, pEnd - pStart);
614:   for (f = 0; f < s->numFields; ++f) PetscSectionSetChart(s->field[f], pStart, pEnd);
615:   return 0;
616: }

618: /*@
619:   PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL that was set with `PetscSectionSetPermutation()`

621:   Not Collective

623:   Input Parameter:
624: . s - the `PetscSection`

626:   Output Parameters:
627: . perm - The permutation as an `IS`

629:   Level: intermediate

631: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
632: @*/
633: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
634: {
636:   if (perm) {
638:     *perm = s->perm;
639:   }
640:   return 0;
641: }

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

646:   Not Collective

648:   Input Parameters:
649: + s - the PetscSection
650: - perm - the permutation of points

652:   Level: intermediate

654:   Developer Note:
655:   What purpose does this permutation serve?

657: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetPermutation()`, `PetscSectionCreate()`
658: @*/
659: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
660: {
664:   if (s->perm != perm) {
665:     ISDestroy(&s->perm);
666:     if (perm) {
667:       s->perm = perm;
668:       PetscObjectReference((PetscObject)s->perm);
669:     }
670:   }
671:   return 0;
672: }

674: /*@
675:   PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major

677:   Not Collective

679:   Input Parameter:
680: . s - the `PetscSection`

682:   Output Parameter:
683: . pm - the flag for point major ordering

685:   Level: intermediate

687: .seealso: [PetscSection](sec_petscsection), `PetscSection`, PetscSectionSetPointMajor()`
688: @*/
689: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
690: {
693:   *pm = s->pointMajor;
694:   return 0;
695: }

697: /*@
698:   PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major

700:   Not Collective

702:   Input Parameters:
703: + s  - the `PetscSection`
704: - pm - the flag for point major ordering

706:   Level: intermediate

708: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`
709: @*/
710: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
711: {
714:   s->pointMajor = pm;
715:   return 0;
716: }

718: /*@
719:   PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
720:   The value is set with `PetscSectionSetIncludesConstraints()`

722:   Not Collective

724:   Input Parameter:
725: . s - the `PetscSection`

727:   Output Parameter:
728: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets

730:   Level: intermediate

732: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
733: @*/
734: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
735: {
738:   *includesConstraints = s->includesConstraints;
739:   return 0;
740: }

742: /*@
743:   PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets

745:   Not Collective

747:   Input Parameters:
748: + s  - the PetscSection
749: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets

751:   Level: intermediate

753: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
754: @*/
755: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
756: {
759:   s->includesConstraints = includesConstraints;
760:   return 0;
761: }

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

766:   Not Collective

768:   Input Parameters:
769: + s - the `PetscSection`
770: - point - the point

772:   Output Parameter:
773: . numDof - the number of dof

775:   Note:
776:   In a global section, this size will be negative for points not owned by this process.

778:   Level: intermediate

780: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
781: @*/
782: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
783: {
787:   PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
788:   *numDof = s->atlasDof[point - s->pStart];
789:   return 0;
790: }

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

795:   Not Collective

797:   Input Parameters:
798: + s - the `PetscSection`
799: . point - the point
800: - numDof - the number of dof

802:   Level: intermediate

804: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
805: @*/
806: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
807: {
809:   PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
810:   s->atlasDof[point - s->pStart] = numDof;
811:   PetscSectionInvalidateMaxDof_Internal(s);
812:   return 0;
813: }

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

818:   Not Collective

820:   Input Parameters:
821: + s - the `PetscSection`
822: . point - the point
823: - numDof - the number of additional dof

825:   Level: intermediate

827: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
828: @*/
829: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
830: {
833:   PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
834:   s->atlasDof[point - s->pStart] += numDof;
835:   PetscSectionInvalidateMaxDof_Internal(s);
836:   return 0;
837: }

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

842:   Not Collective

844:   Input Parameters:
845: + s - the `PetscSection`
846: . point - the point
847: - field - the field

849:   Output Parameter:
850: . numDof - the number of dof

852:   Level: intermediate

854: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
855: @*/
856: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
857: {
860:   PetscSectionCheckValidField(field, s->numFields);
861:   PetscSectionGetDof(s->field[field], point, numDof);
862:   return 0;
863: }

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

868:   Not Collective

870:   Input Parameters:
871: + s - the `PetscSection`
872: . point - the point
873: . field - the field
874: - numDof - the number of dof

876:   Level: intermediate

878: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
879: @*/
880: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
881: {
883:   PetscSectionCheckValidField(field, s->numFields);
884:   PetscSectionSetDof(s->field[field], point, numDof);
885:   return 0;
886: }

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

891:   Not Collective

893:   Input Parameters:
894: + s - the `PetscSection`
895: . point - the point
896: . field - the field
897: - numDof - the number of dof

899:   Level: intermediate

901: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
902: @*/
903: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
904: {
906:   PetscSectionCheckValidField(field, s->numFields);
907:   PetscSectionAddDof(s->field[field], point, numDof);
908:   return 0;
909: }

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

914:   Not Collective

916:   Input Parameters:
917: + s - the `PetscSection`
918: - point - the point

920:   Output Parameter:
921: . numDof - the number of dof which are fixed by constraints

923:   Level: intermediate

925: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
926: @*/
927: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
928: {
931:   if (s->bc) {
932:     PetscSectionGetDof(s->bc, point, numDof);
933:   } else *numDof = 0;
934:   return 0;
935: }

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

940:   Not Collective

942:   Input Parameters:
943: + s - the `PetscSection`
944: . point - the point
945: - numDof - the number of dof which are fixed by constraints

947:   Level: intermediate

949: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
950: @*/
951: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
952: {
954:   if (numDof) {
955:     PetscSectionCheckConstraints_Private(s);
956:     PetscSectionSetDof(s->bc, point, numDof);
957:   }
958:   return 0;
959: }

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

964:   Not Collective

966:   Input Parameters:
967: + s - the `PetscSection`
968: . point - the point
969: - numDof - the number of additional dof which are fixed by constraints

971:   Level: intermediate

973: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
974: @*/
975: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
976: {
978:   if (numDof) {
979:     PetscSectionCheckConstraints_Private(s);
980:     PetscSectionAddDof(s->bc, point, numDof);
981:   }
982:   return 0;
983: }

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

988:   Not Collective

990:   Input Parameters:
991: + s - the `PetscSection`
992: . point - the point
993: - field - the field

995:   Output Parameter:
996: . numDof - the number of dof which are fixed by constraints

998:   Level: intermediate

1000: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1001: @*/
1002: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1003: {
1006:   PetscSectionCheckValidField(field, s->numFields);
1007:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
1008:   return 0;
1009: }

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

1014:   Not Collective

1016:   Input Parameters:
1017: + s - the `PetscSection`
1018: . point - the point
1019: . field - the field
1020: - numDof - the number of dof which are fixed by constraints

1022:   Level: intermediate

1024: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1025: @*/
1026: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1027: {
1029:   PetscSectionCheckValidField(field, s->numFields);
1030:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
1031:   return 0;
1032: }

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

1037:   Not Collective

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

1045:   Level: intermediate

1047: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1048: @*/
1049: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1050: {
1052:   PetscSectionCheckValidField(field, s->numFields);
1053:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
1054:   return 0;
1055: }

1057: /*@
1058:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1060:   Not Collective

1062:   Input Parameter:
1063: . s - the `PetscSection`

1065:   Level: advanced

1067: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1068: @*/
1069: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1070: {
1072:   if (s->bc) {
1073:     const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;

1075:     PetscSectionSetUp(s->bc);
1076:     PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices);
1077:   }
1078:   return 0;
1079: }

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

1084:   Not Collective

1086:   Input Parameter:
1087: . s - the `PetscSection`

1089:   Level: intermediate

1091: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1092: @*/
1093: PetscErrorCode PetscSectionSetUp(PetscSection s)
1094: {
1095:   const PetscInt *pind   = NULL;
1096:   PetscInt        offset = 0, foff, p, f;

1099:   if (s->setup) return 0;
1100:   s->setup = PETSC_TRUE;
1101:   /* Set offsets and field offsets for all points */
1102:   /*   Assume that all fields have the same chart */
1104:   if (s->perm) ISGetIndices(s->perm, &pind);
1105:   if (s->pointMajor) {
1106:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1107:       const PetscInt q = pind ? pind[p] : p;

1109:       /* Set point offset */
1110:       s->atlasOff[q] = offset;
1111:       offset += s->atlasDof[q];
1112:       /* Set field offset */
1113:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1114:         PetscSection sf = s->field[f];

1116:         sf->atlasOff[q] = foff;
1117:         foff += sf->atlasDof[q];
1118:       }
1119:     }
1120:   } else {
1121:     /* Set field offsets for all points */
1122:     for (f = 0; f < s->numFields; ++f) {
1123:       PetscSection sf = s->field[f];

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

1128:         sf->atlasOff[q] = offset;
1129:         offset += sf->atlasDof[q];
1130:       }
1131:     }
1132:     /* Disable point offsets since these are unused */
1133:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1134:   }
1135:   if (s->perm) ISRestoreIndices(s->perm, &pind);
1136:   /* Setup BC sections */
1137:   PetscSectionSetUpBC(s);
1138:   for (f = 0; f < s->numFields; ++f) PetscSectionSetUpBC(s->field[f]);
1139:   return 0;
1140: }

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

1145:   Not Collective

1147:   Input Parameters:
1148: . s - the `PetscSection`

1150:   Output Parameter:
1151: . maxDof - the maximum dof

1153:   Level: intermediate

1155:   Note:
1156:   The returned number is up-to-date without need for `PetscSectionSetUp()`.

1158:   Developer Note:
1159:   The returned number is calculated lazily and stashed.

1161:   A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.

1163:   `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`

1165:   It should also be called every time atlasDof is modified directly.

1167: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1168: @*/
1169: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1170: {
1171:   PetscInt p;

1175:   if (s->maxDof == PETSC_MIN_INT) {
1176:     s->maxDof = 0;
1177:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1178:   }
1179:   *maxDof = s->maxDof;
1180:   return 0;
1181: }

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

1186:   Not Collective

1188:   Input Parameter:
1189: . s - the `PetscSection`

1191:   Output Parameter:
1192: . size - the size of an array which can hold all the dofs

1194:   Level: intermediate

1196: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1197: @*/
1198: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1199: {
1200:   PetscInt p, n = 0;

1204:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1205:   *size = n;
1206:   return 0;
1207: }

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

1212:   Not Collective

1214:   Input Parameter:
1215: . s - the `PetscSection`

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

1220:   Level: intermediate

1222: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1223: @*/
1224: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1225: {
1226:   PetscInt p, n = 0;

1230:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1231:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1232:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1233:   }
1234:   *size = n;
1235:   return 0;
1236: }

1238: /*@
1239:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1240:   the local section and a `PetscSF` describing the section point overlap.

1242:   Input Parameters:
1243: + s - The `PetscSection` for the local field layout
1244: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1245: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1246: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points

1248:   Output Parameter:
1249: . gsection - The `PetscSection` for the global field layout

1251:   Level: intermediate

1253:   Notes:
1254:   If we have a set of local sections defining the layout of a set of local vectors, and also a `PetscSF` to determine which section points are shared and the ownership, we can calculate a global section defining the parallel data layout, and the associated global vector.

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

1258: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1259: @*/
1260: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1261: {
1262:   PetscSection    gs;
1263:   const PetscInt *pind = NULL;
1264:   PetscInt       *recv = NULL, *neg = NULL;
1265:   PetscInt        pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1266:   PetscInt        numFields, f, numComponents;

1274:   PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs);
1275:   PetscSectionGetNumFields(s, &numFields);
1276:   if (numFields > 0) PetscSectionSetNumFields(gs, numFields);
1277:   PetscSectionGetChart(s, &pStart, &pEnd);
1278:   PetscSectionSetChart(gs, pStart, pEnd);
1279:   gs->includesConstraints = includeConstraints;
1280:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1281:   nlocal = nroots; /* The local/leaf space matches global/root space */
1282:   /* Must allocate for all points visible to SF, which may be more than this section */
1283:   if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1284:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
1287:     PetscMalloc2(nroots, &neg, nlocal, &recv);
1288:     PetscArrayzero(neg, nroots);
1289:   }
1290:   /* Mark all local points with negative dof */
1291:   for (p = pStart; p < pEnd; ++p) {
1292:     PetscSectionGetDof(s, p, &dof);
1293:     PetscSectionSetDof(gs, p, dof);
1294:     PetscSectionGetConstraintDof(s, p, &cdof);
1295:     if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(gs, p, cdof);
1296:     if (neg) neg[p] = -(dof + 1);
1297:   }
1298:   PetscSectionSetUpBC(gs);
1299:   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]);
1300:   if (nroots >= 0) {
1301:     PetscArrayzero(recv, nlocal);
1302:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1303:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1304:     for (p = pStart; p < pEnd; ++p) {
1305:       if (recv[p] < 0) {
1306:         gs->atlasDof[p - pStart] = recv[p];
1307:         PetscSectionGetDof(s, p, &dof);
1309:       }
1310:     }
1311:   }
1312:   /* Calculate new sizes, get process offset, and calculate point offsets */
1313:   if (s->perm) ISGetIndices(s->perm, &pind);
1314:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1315:     const PetscInt q = pind ? pind[p] : p;

1317:     cdof            = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1318:     gs->atlasOff[q] = off;
1319:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1320:   }
1321:   if (!localOffsets) {
1322:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1323:     globalOff -= off;
1324:   }
1325:   for (p = pStart, off = 0; p < pEnd; ++p) {
1326:     gs->atlasOff[p - pStart] += globalOff;
1327:     if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1328:   }
1329:   if (s->perm) ISRestoreIndices(s->perm, &pind);
1330:   /* Put in negative offsets for ghost points */
1331:   if (nroots >= 0) {
1332:     PetscArrayzero(recv, nlocal);
1333:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1334:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1335:     for (p = pStart; p < pEnd; ++p) {
1336:       if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1337:     }
1338:   }
1339:   PetscFree2(neg, recv);
1340:   /* Set field dofs/offsets/constraints */
1341:   for (f = 0; f < numFields; ++f) {
1342:     gs->field[f]->includesConstraints = includeConstraints;
1343:     PetscSectionGetFieldComponents(s, f, &numComponents);
1344:     PetscSectionSetFieldComponents(gs, f, numComponents);
1345:   }
1346:   for (p = pStart; p < pEnd; ++p) {
1347:     PetscSectionGetOffset(gs, p, &off);
1348:     for (f = 0, foff = off; f < numFields; ++f) {
1349:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1350:       if (!includeConstraints && cdof > 0) PetscSectionSetFieldConstraintDof(gs, p, f, cdof);
1351:       PetscSectionGetFieldDof(s, p, f, &dof);
1352:       PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof);
1353:       PetscSectionSetFieldOffset(gs, p, f, foff);
1354:       PetscSectionGetFieldConstraintDof(gs, p, f, &cdof);
1355:       foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1356:     }
1357:   }
1358:   for (f = 0; f < numFields; ++f) {
1359:     PetscSection gfs = gs->field[f];

1361:     PetscSectionSetUpBC(gfs);
1362:     if (gfs->bcIndices) PetscArraycpy(gfs->bcIndices, s->field[f]->bcIndices, gfs->bc->atlasOff[gfs->bc->pEnd - gfs->bc->pStart - 1] + gfs->bc->atlasDof[gfs->bc->pEnd - gfs->bc->pStart - 1]);
1363:   }
1364:   gs->setup = PETSC_TRUE;
1365:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1366:   *gsection = gs;
1367:   return 0;
1368: }

1370: /*@
1371:   PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the global field layout using
1372:   the local section and an `PetscSF` describing the section point overlap.

1374:   Input Parameters:
1375: + s - The `PetscSection` for the local field layout
1376: . sf - The `PetscSF` describing parallel layout of the section points
1377: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1378: . numExcludes - The number of exclusion ranges
1379: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs

1381:   Output Parameter:
1382: . gsection - The `PetscSection` for the global field layout

1384:   Level: advanced

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

1389: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1390: @*/
1391: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1392: {
1393:   const PetscInt *pind = NULL;
1394:   PetscInt       *neg = NULL, *tmpOff = NULL;
1395:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;

1400:   PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection);
1401:   PetscSectionGetChart(s, &pStart, &pEnd);
1402:   PetscSectionSetChart(*gsection, pStart, pEnd);
1403:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1404:   if (nroots >= 0) {
1406:     PetscCalloc1(nroots, &neg);
1407:     if (nroots > pEnd - pStart) {
1408:       PetscCalloc1(nroots, &tmpOff);
1409:     } else {
1410:       tmpOff = &(*gsection)->atlasDof[-pStart];
1411:     }
1412:   }
1413:   /* Mark ghost points with negative dof */
1414:   for (p = pStart; p < pEnd; ++p) {
1415:     for (e = 0; e < numExcludes; ++e) {
1416:       if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1417:         PetscSectionSetDof(*gsection, p, 0);
1418:         break;
1419:       }
1420:     }
1421:     if (e < numExcludes) continue;
1422:     PetscSectionGetDof(s, p, &dof);
1423:     PetscSectionSetDof(*gsection, p, dof);
1424:     PetscSectionGetConstraintDof(s, p, &cdof);
1425:     if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
1426:     if (neg) neg[p] = -(dof + 1);
1427:   }
1428:   PetscSectionSetUpBC(*gsection);
1429:   if (nroots >= 0) {
1430:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1431:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1432:     if (nroots > pEnd - pStart) {
1433:       for (p = pStart; p < pEnd; ++p) {
1434:         if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1435:       }
1436:     }
1437:   }
1438:   /* Calculate new sizes, get process offset, and calculate point offsets */
1439:   if (s->perm) ISGetIndices(s->perm, &pind);
1440:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1441:     const PetscInt q = pind ? pind[p] : p;

1443:     cdof                     = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1444:     (*gsection)->atlasOff[q] = off;
1445:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1446:   }
1447:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1448:   globalOff -= off;
1449:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1450:     (*gsection)->atlasOff[p] += globalOff;
1451:     if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1452:   }
1453:   if (s->perm) ISRestoreIndices(s->perm, &pind);
1454:   /* Put in negative offsets for ghost points */
1455:   if (nroots >= 0) {
1456:     if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1457:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1458:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1459:     if (nroots > pEnd - pStart) {
1460:       for (p = pStart; p < pEnd; ++p) {
1461:         if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1462:       }
1463:     }
1464:   }
1465:   if (nroots >= 0 && nroots > pEnd - pStart) PetscFree(tmpOff);
1466:   PetscFree(neg);
1467:   return 0;
1468: }

1470: /*@
1471:   PetscSectionGetPointLayout - Get the `PetscLayout` associated with the section points

1473:   Collective

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

1479:   Output Parameter:
1480: . layout - The point layout for the section

1482:   Level: advanced

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

1487: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1488: @*/
1489: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1490: {
1491:   PetscInt pStart, pEnd, p, localSize = 0;

1493:   PetscSectionGetChart(s, &pStart, &pEnd);
1494:   for (p = pStart; p < pEnd; ++p) {
1495:     PetscInt dof;

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

1507: /*@
1508:   PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs.

1510:   Collective

1512:   Input Parameters:
1513: + comm - The `MPI_Comm`
1514: - s    - The `PetscSection`

1516:   Output Parameter:
1517: . layout - The dof layout for the section

1519:   Level: advanced

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

1524: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1525: @*/
1526: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1527: {
1528:   PetscInt pStart, pEnd, p, localSize = 0;

1532:   PetscSectionGetChart(s, &pStart, &pEnd);
1533:   for (p = pStart; p < pEnd; ++p) {
1534:     PetscInt dof, cdof;

1536:     PetscSectionGetDof(s, p, &dof);
1537:     PetscSectionGetConstraintDof(s, p, &cdof);
1538:     if (dof - cdof > 0) localSize += dof - cdof;
1539:   }
1540:   PetscLayoutCreate(comm, layout);
1541:   PetscLayoutSetLocalSize(*layout, localSize);
1542:   PetscLayoutSetBlockSize(*layout, 1);
1543:   PetscLayoutSetUp(*layout);
1544:   return 0;
1545: }

1547: /*@
1548:   PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.

1550:   Not Collective

1552:   Input Parameters:
1553: + s - the `PetscSection`
1554: - point - the point

1556:   Output Parameter:
1557: . offset - the offset

1559:   Note:
1560:   In a global section, this offset will be negative for points not owned by this process.

1562:   Level: intermediate

1564: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`
1565: @*/
1566: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1567: {
1571:   *offset = s->atlasOff[point - s->pStart];
1572:   return 0;
1573: }

1575: /*@
1576:   PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.

1578:   Not Collective

1580:   Input Parameters:
1581: + s - the `PetscSection`
1582: . point - the point
1583: - offset - the offset

1585:   Level: intermediate

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

1590: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1591: @*/
1592: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1593: {
1596:   s->atlasOff[point - s->pStart] = offset;
1597:   return 0;
1598: }

1600: /*@
1601:   PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.

1603:   Not Collective

1605:   Input Parameters:
1606: + s - the `PetscSection`
1607: . point - the point
1608: - field - the field

1610:   Output Parameter:
1611: . offset - the offset

1613:   Note:
1614:   In a global section, this offset will be negative for points not owned by this process.

1616:   Level: intermediate

1618: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1619: @*/
1620: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1621: {
1624:   PetscSectionCheckValidField(field, s->numFields);
1625:   PetscSectionGetOffset(s->field[field], point, offset);
1626:   return 0;
1627: }

1629: /*@
1630:   PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.

1632:   Not Collective

1634:   Input Parameters:
1635: + s - the PetscSection
1636: . point - the point
1637: . field - the field
1638: - offset - the offset

1640:   Level: developer

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

1645: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1646: @*/
1647: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1648: {
1650:   PetscSectionCheckValidField(field, s->numFields);
1651:   PetscSectionSetOffset(s->field[field], point, offset);
1652:   return 0;
1653: }

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

1658:   Not Collective

1660:   Input Parameters:
1661: + s - the `PetscSection`
1662: . point - the point
1663: - field - the field

1665:   Output Parameter:
1666: . offset - the offset

1668:   Level: advanced

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

1674: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1675: @*/
1676: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1677: {
1678:   PetscInt off, foff;

1682:   PetscSectionCheckValidField(field, s->numFields);
1683:   PetscSectionGetOffset(s, point, &off);
1684:   PetscSectionGetOffset(s->field[field], point, &foff);
1685:   *offset = foff - off;
1686:   return 0;
1687: }

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

1692:   Not Collective

1694:   Input Parameter:
1695: . s - the `PetscSection`

1697:   Output Parameters:
1698: + start - the minimum offset
1699: - end   - one more than the maximum offset

1701:   Level: intermediate

1703: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1704: @*/
1705: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1706: {
1707:   PetscInt os = 0, oe = 0, pStart, pEnd, p;

1710:   if (s->atlasOff) {
1711:     os = s->atlasOff[0];
1712:     oe = s->atlasOff[0];
1713:   }
1714:   PetscSectionGetChart(s, &pStart, &pEnd);
1715:   for (p = 0; p < pEnd - pStart; ++p) {
1716:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1718:     if (off >= 0) {
1719:       os = PetscMin(os, off);
1720:       oe = PetscMax(oe, off + dof);
1721:     }
1722:   }
1723:   if (start) *start = os;
1724:   if (end) *end = oe;
1725:   return 0;
1726: }

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

1731:   Collective

1733:   Input Parameters:
1734: + s      - the `PetscSection`
1735: . len    - the number of subfields
1736: - fields - the subfield numbers

1738:   Output Parameter:
1739: . subs   - the subsection

1741:   Level: advanced

1743:   Notes:
1744:   The section offsets now refer to a new, smaller vector.

1746:   This will error if a fieldnumber is out of range

1748: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1749: @*/
1750: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1751: {
1752:   PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;

1754:   if (!len) return 0;
1758:   PetscSectionGetNumFields(s, &nF);
1760:   PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1761:   PetscSectionSetNumFields(*subs, len);
1762:   for (f = 0; f < len; ++f) {
1763:     const char *name    = NULL;
1764:     PetscInt    numComp = 0;

1766:     PetscSectionGetFieldName(s, fields[f], &name);
1767:     PetscSectionSetFieldName(*subs, f, name);
1768:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1769:     PetscSectionSetFieldComponents(*subs, f, numComp);
1770:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1771:       PetscSectionGetComponentName(s, fields[f], c, &name);
1772:       PetscSectionSetComponentName(*subs, f, c, name);
1773:     }
1774:   }
1775:   PetscSectionGetChart(s, &pStart, &pEnd);
1776:   PetscSectionSetChart(*subs, pStart, pEnd);
1777:   for (p = pStart; p < pEnd; ++p) {
1778:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1780:     for (f = 0; f < len; ++f) {
1781:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1782:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1783:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1784:       if (cfdof) PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);
1785:       dof += fdof;
1786:       cdof += cfdof;
1787:     }
1788:     PetscSectionSetDof(*subs, p, dof);
1789:     if (cdof) PetscSectionSetConstraintDof(*subs, p, cdof);
1790:     maxCdof = PetscMax(cdof, maxCdof);
1791:   }
1792:   PetscSectionSetUp(*subs);
1793:   if (maxCdof) {
1794:     PetscInt *indices;

1796:     PetscMalloc1(maxCdof, &indices);
1797:     for (p = pStart; p < pEnd; ++p) {
1798:       PetscInt cdof;

1800:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1801:       if (cdof) {
1802:         const PetscInt *oldIndices = NULL;
1803:         PetscInt        fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1805:         for (f = 0; f < len; ++f) {
1806:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1807:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1808:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1809:           PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices);
1810:           for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1811:           numConst += cfdof;
1812:           fOff += fdof;
1813:         }
1814:         PetscSectionSetConstraintIndices(*subs, p, indices);
1815:       }
1816:     }
1817:     PetscFree(indices);
1818:   }
1819:   return 0;
1820: }

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

1825:   Collective

1827:   Input Parameters:
1828: + s     - the input sections
1829: - len   - the number of input sections

1831:   Output Parameter:
1832: . supers - the supersection

1834:   Level: advanced

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

1839:   Developer Note:
1840:   Needs to explain how the sections are composed

1842: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
1843: @*/
1844: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1845: {
1846:   PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1848:   if (!len) return 0;
1849:   for (i = 0; i < len; ++i) {
1850:     PetscInt nf, pStarti, pEndi;

1852:     PetscSectionGetNumFields(s[i], &nf);
1853:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1854:     pStart = PetscMin(pStart, pStarti);
1855:     pEnd   = PetscMax(pEnd, pEndi);
1856:     Nf += nf;
1857:   }
1858:   PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers);
1859:   PetscSectionSetNumFields(*supers, Nf);
1860:   for (i = 0, f = 0; i < len; ++i) {
1861:     PetscInt nf, fi, ci;

1863:     PetscSectionGetNumFields(s[i], &nf);
1864:     for (fi = 0; fi < nf; ++fi, ++f) {
1865:       const char *name    = NULL;
1866:       PetscInt    numComp = 0;

1868:       PetscSectionGetFieldName(s[i], fi, &name);
1869:       PetscSectionSetFieldName(*supers, f, name);
1870:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1871:       PetscSectionSetFieldComponents(*supers, f, numComp);
1872:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1873:         PetscSectionGetComponentName(s[i], fi, ci, &name);
1874:         PetscSectionSetComponentName(*supers, f, ci, name);
1875:       }
1876:     }
1877:   }
1878:   PetscSectionSetChart(*supers, pStart, pEnd);
1879:   for (p = pStart; p < pEnd; ++p) {
1880:     PetscInt dof = 0, cdof = 0;

1882:     for (i = 0, f = 0; i < len; ++i) {
1883:       PetscInt nf, fi, pStarti, pEndi;
1884:       PetscInt fdof = 0, cfdof = 0;

1886:       PetscSectionGetNumFields(s[i], &nf);
1887:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1888:       if ((p < pStarti) || (p >= pEndi)) continue;
1889:       for (fi = 0; fi < nf; ++fi, ++f) {
1890:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1891:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1892:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1893:         if (cfdof) PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);
1894:         dof += fdof;
1895:         cdof += cfdof;
1896:       }
1897:     }
1898:     PetscSectionSetDof(*supers, p, dof);
1899:     if (cdof) PetscSectionSetConstraintDof(*supers, p, cdof);
1900:     maxCdof = PetscMax(cdof, maxCdof);
1901:   }
1902:   PetscSectionSetUp(*supers);
1903:   if (maxCdof) {
1904:     PetscInt *indices;

1906:     PetscMalloc1(maxCdof, &indices);
1907:     for (p = pStart; p < pEnd; ++p) {
1908:       PetscInt cdof;

1910:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1911:       if (cdof) {
1912:         PetscInt dof, numConst = 0, fOff = 0;

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

1918:           PetscSectionGetNumFields(s[i], &nf);
1919:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1920:           if ((p < pStarti) || (p >= pEndi)) continue;
1921:           for (fi = 0; fi < nf; ++fi, ++f) {
1922:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1923:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1924:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1925:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
1926:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1927:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
1928:             numConst += cfdof;
1929:           }
1930:           PetscSectionGetDof(s[i], p, &dof);
1931:           fOff += dof;
1932:         }
1933:         PetscSectionSetConstraintIndices(*supers, p, indices);
1934:       }
1935:     }
1936:     PetscFree(indices);
1937:   }
1938:   return 0;
1939: }

1941: PetscErrorCode PetscSectionCreateSubplexSection_Internal(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
1942: {
1943:   const PetscInt *points = NULL, *indices = NULL;
1944:   PetscInt        numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;

1949:   PetscSectionGetNumFields(s, &numFields);
1950:   PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1951:   if (numFields) PetscSectionSetNumFields(*subs, numFields);
1952:   for (f = 0; f < numFields; ++f) {
1953:     const char *name    = NULL;
1954:     PetscInt    numComp = 0;

1956:     PetscSectionGetFieldName(s, f, &name);
1957:     PetscSectionSetFieldName(*subs, f, name);
1958:     PetscSectionGetFieldComponents(s, f, &numComp);
1959:     PetscSectionSetFieldComponents(*subs, f, numComp);
1960:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
1961:       PetscSectionGetComponentName(s, f, c, &name);
1962:       PetscSectionSetComponentName(*subs, f, c, name);
1963:     }
1964:   }
1965:   /* For right now, we do not try to squeeze the subchart */
1966:   if (subpointMap) {
1967:     ISGetSize(subpointMap, &numSubpoints);
1968:     ISGetIndices(subpointMap, &points);
1969:   }
1970:   PetscSectionGetChart(s, &pStart, &pEnd);
1971:   if (renumberPoints) {
1972:     spStart = 0;
1973:     spEnd   = numSubpoints;
1974:   } else {
1975:     ISGetMinMax(subpointMap, &spStart, &spEnd);
1976:     ++spEnd;
1977:   }
1978:   PetscSectionSetChart(*subs, spStart, spEnd);
1979:   for (p = pStart; p < pEnd; ++p) {
1980:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1982:     PetscFindInt(p, numSubpoints, points, &subp);
1983:     if (subp < 0) continue;
1984:     if (!renumberPoints) subp = p;
1985:     for (f = 0; f < numFields; ++f) {
1986:       PetscSectionGetFieldDof(s, p, f, &fdof);
1987:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1988:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1989:       if (cfdof) PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);
1990:     }
1991:     PetscSectionGetDof(s, p, &dof);
1992:     PetscSectionSetDof(*subs, subp, dof);
1993:     PetscSectionGetConstraintDof(s, p, &cdof);
1994:     if (cdof) PetscSectionSetConstraintDof(*subs, subp, cdof);
1995:   }
1996:   PetscSectionSetUp(*subs);
1997:   /* Change offsets to original offsets */
1998:   for (p = pStart; p < pEnd; ++p) {
1999:     PetscInt off, foff = 0;

2001:     PetscFindInt(p, numSubpoints, points, &subp);
2002:     if (subp < 0) continue;
2003:     if (!renumberPoints) subp = p;
2004:     for (f = 0; f < numFields; ++f) {
2005:       PetscSectionGetFieldOffset(s, p, f, &foff);
2006:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
2007:     }
2008:     PetscSectionGetOffset(s, p, &off);
2009:     PetscSectionSetOffset(*subs, subp, off);
2010:   }
2011:   /* Copy constraint indices */
2012:   for (subp = spStart; subp < spEnd; ++subp) {
2013:     PetscInt cdof;

2015:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
2016:     if (cdof) {
2017:       for (f = 0; f < numFields; ++f) {
2018:         PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices);
2019:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2020:       }
2021:       PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices);
2022:       PetscSectionSetConstraintIndices(*subs, subp, indices);
2023:     }
2024:   }
2025:   if (subpointMap) ISRestoreIndices(subpointMap, &points);
2026:   return 0;
2027: }

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

2032:   Collective

2034:   Input Parameters:
2035: + s           - the `PetscSection`
2036: - subpointMap - a sorted list of points in the original mesh which are in the submesh

2038:   Output Parameter:
2039: . subs - the subsection

2041:   Level: advanced

2043:   Note:
2044:   The points are renumbered from 0, and the section offsets now refer to a new, smaller vector.

2046:   Developer Note:
2047:   The use of the term Submesh is confusing and unneeded

2049: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2050: @*/
2051: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2052: {
2053:   PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs);
2054:   return 0;
2055: }

2057: /*@
2058:   PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh

2060:   Collective

2062:   Input Parameters:
2063: + s           - the `PetscSection`
2064: - subpointMap - a sorted list of points in the original mesh which are in the subdomain

2066:   Output Parameter:
2067: . subs - the subsection

2069:   Level: advanced

2071:   Note:
2072:   The point numbers remain the same, but the section offsets now refer to a new, smaller vector.

2074:   Developer Notes:
2075:   It is unclear what the difference with `PetscSectionCreateSubmeshSection()` is.

2077:   The use of the term Subdomain is unneeded and confusing

2079: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2080: @*/
2081: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2082: {
2083:   PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs);
2084:   return 0;
2085: }

2087: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2088: {
2089:   PetscInt    p;
2090:   PetscMPIInt rank;

2092:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2093:   PetscViewerASCIIPushSynchronized(viewer);
2094:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2095:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2096:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2097:       PetscInt b;

2099:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2100:       if (s->bcIndices) {
2101:         for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
2102:       }
2103:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2104:     } else {
2105:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2106:     }
2107:   }
2108:   PetscViewerFlush(viewer);
2109:   PetscViewerASCIIPopSynchronized(viewer);
2110:   if (s->sym) {
2111:     PetscViewerASCIIPushTab(viewer);
2112:     PetscSectionSymView(s->sym, viewer);
2113:     PetscViewerASCIIPopTab(viewer);
2114:   }
2115:   return 0;
2116: }

2118: /*@C
2119:    PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database

2121:    Collective

2123:    Input Parameters:
2124: +  A - the PetscSection object to view
2125: .  obj - Optional object that provides the options prefix used for the options
2126: -  name - command line option

2128:    Level: intermediate

2130: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2131: @*/
2132: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2133: {
2135:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
2136:   return 0;
2137: }

2139: /*@C
2140:   PetscSectionView - Views a `PetscSection`

2142:   Collective

2144:   Input Parameters:
2145: + s - the `PetscSection` object to view
2146: - v - the viewer

2148:   Level: beginner

2150:   Note:
2151:   `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2152:   distribution independent data, such as dofs, offsets, constraint dofs,
2153:   and constraint indices. Points that have negative dofs, for instance,
2154:   are not saved as they represent points owned by other processes.
2155:   Point numbering and rank assignment is currently not stored.
2156:   The saved section can be loaded with `PetscSectionLoad()`.

2158: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2159: @*/
2160: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2161: {
2162:   PetscBool isascii, ishdf5;
2163:   PetscInt  f;

2166:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);
2168:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
2169:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2170:   if (isascii) {
2171:     PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer);
2172:     if (s->numFields) {
2173:       PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields);
2174:       for (f = 0; f < s->numFields; ++f) {
2175:         PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
2176:         PetscSectionView_ASCII(s->field[f], viewer);
2177:       }
2178:     } else {
2179:       PetscSectionView_ASCII(s, viewer);
2180:     }
2181:   } else if (ishdf5) {
2182: #if PetscDefined(HAVE_HDF5)
2183:     PetscSectionView_HDF5_Internal(s, viewer);
2184: #else
2185:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2186: #endif
2187:   }
2188:   return 0;
2189: }

2191: /*@C
2192:   PetscSectionLoad - Loads a `PetscSection`

2194:   Collective

2196:   Input Parameters:
2197: + s - the `PetscSection` object to load
2198: - v - the viewer

2200:   Level: beginner

2202:   Note:
2203:   `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2204:   a section saved with `PetscSectionView()`. The number of processes
2205:   used here (N) does not need to be the same as that used when saving.
2206:   After calling this function, the chart of s on rank i will be set
2207:   to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2208:   saved section points.

2210: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2211: @*/
2212: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2213: {
2214:   PetscBool ishdf5;

2218:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2219:   if (ishdf5) {
2220: #if PetscDefined(HAVE_HDF5)
2221:     PetscSectionLoad_HDF5_Internal(s, viewer);
2222:     return 0;
2223: #else
2224:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2225: #endif
2226:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2227: }

2229: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2230: {
2231:   PetscSectionClosurePermVal clVal;

2233:   if (!section->clHash) return 0;
2234:   kh_foreach_value(section->clHash, clVal, {
2235:     PetscFree(clVal.perm);
2236:     PetscFree(clVal.invPerm);
2237:   });
2238:   kh_destroy(ClPerm, section->clHash);
2239:   section->clHash = NULL;
2240:   return 0;
2241: }

2243: /*@
2244:   PetscSectionReset - Frees all section data.

2246:   Not Collective

2248:   Input Parameters:
2249: . s - the `PetscSection`

2251:   Level: beginner

2253: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2254: @*/
2255: PetscErrorCode PetscSectionReset(PetscSection s)
2256: {
2257:   PetscInt f, c;

2260:   for (f = 0; f < s->numFields; ++f) {
2261:     PetscSectionDestroy(&s->field[f]);
2262:     PetscFree(s->fieldNames[f]);
2263:     for (c = 0; c < s->numFieldComponents[f]; ++c) PetscFree(s->compNames[f][c]);
2264:     PetscFree(s->compNames[f]);
2265:   }
2266:   PetscFree(s->numFieldComponents);
2267:   PetscFree(s->fieldNames);
2268:   PetscFree(s->compNames);
2269:   PetscFree(s->field);
2270:   PetscSectionDestroy(&s->bc);
2271:   PetscFree(s->bcIndices);
2272:   PetscFree2(s->atlasDof, s->atlasOff);
2273:   PetscSectionDestroy(&s->clSection);
2274:   ISDestroy(&s->clPoints);
2275:   ISDestroy(&s->perm);
2276:   PetscSectionResetClosurePermutation(s);
2277:   PetscSectionSymDestroy(&s->sym);
2278:   PetscSectionDestroy(&s->clSection);
2279:   ISDestroy(&s->clPoints);
2280:   PetscSectionInvalidateMaxDof_Internal(s);
2281:   s->pStart    = -1;
2282:   s->pEnd      = -1;
2283:   s->maxDof    = 0;
2284:   s->setup     = PETSC_FALSE;
2285:   s->numFields = 0;
2286:   s->clObj     = NULL;
2287:   return 0;
2288: }

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

2293:   Not Collective

2295:   Input Parameters:
2296: . s - the `PetscSection`

2298:   Level: beginner

2300: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2301: @*/
2302: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2303: {
2304:   if (!*s) return 0;
2306:   if (--((PetscObject)(*s))->refct > 0) {
2307:     *s = NULL;
2308:     return 0;
2309:   }
2310:   PetscSectionReset(*s);
2311:   PetscHeaderDestroy(s);
2312:   return 0;
2313: }

2315: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2316: {
2317:   const PetscInt p = point - s->pStart;

2320:   *values = &baseArray[s->atlasOff[p]];
2321:   return 0;
2322: }

2324: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2325: {
2326:   PetscInt      *array;
2327:   const PetscInt p           = point - s->pStart;
2328:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2329:   PetscInt       cDim        = 0;

2332:   PetscSectionGetConstraintDof(s, p, &cDim);
2333:   array = &baseArray[s->atlasOff[p]];
2334:   if (!cDim) {
2335:     if (orientation >= 0) {
2336:       const PetscInt dim = s->atlasDof[p];
2337:       PetscInt       i;

2339:       if (mode == INSERT_VALUES) {
2340:         for (i = 0; i < dim; ++i) array[i] = values[i];
2341:       } else {
2342:         for (i = 0; i < dim; ++i) array[i] += values[i];
2343:       }
2344:     } else {
2345:       PetscInt offset = 0;
2346:       PetscInt j      = -1, field, i;

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

2351:         for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2352:         offset += dim;
2353:       }
2354:     }
2355:   } else {
2356:     if (orientation >= 0) {
2357:       const PetscInt  dim  = s->atlasDof[p];
2358:       PetscInt        cInd = 0, i;
2359:       const PetscInt *cDof;

2361:       PetscSectionGetConstraintIndices(s, point, &cDof);
2362:       if (mode == INSERT_VALUES) {
2363:         for (i = 0; i < dim; ++i) {
2364:           if ((cInd < cDim) && (i == cDof[cInd])) {
2365:             ++cInd;
2366:             continue;
2367:           }
2368:           array[i] = values[i];
2369:         }
2370:       } else {
2371:         for (i = 0; i < dim; ++i) {
2372:           if ((cInd < cDim) && (i == cDof[cInd])) {
2373:             ++cInd;
2374:             continue;
2375:           }
2376:           array[i] += values[i];
2377:         }
2378:       }
2379:     } else {
2380:       const PetscInt *cDof;
2381:       PetscInt        offset  = 0;
2382:       PetscInt        cOffset = 0;
2383:       PetscInt        j       = 0, field;

2385:       PetscSectionGetConstraintIndices(s, point, &cDof);
2386:       for (field = 0; field < s->numFields; ++field) {
2387:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2388:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2389:         const PetscInt sDim = dim - tDim;
2390:         PetscInt       cInd = 0, i, k;

2392:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2393:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2394:             ++cInd;
2395:             continue;
2396:           }
2397:           array[j] = values[k];
2398:         }
2399:         offset += dim;
2400:         cOffset += dim - tDim;
2401:       }
2402:     }
2403:   }
2404:   return 0;
2405: }

2407: /*@C
2408:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2410:   Not Collective

2412:   Input Parameter:
2413: . s - The `PetscSection`

2415:   Output Parameter:
2416: . hasConstraints - flag indicating that the section has constrained dofs

2418:   Level: intermediate

2420: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2421: @*/
2422: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2423: {
2426:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2427:   return 0;
2428: }

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

2433:   Not Collective

2435:   Input Parameters:
2436: + s     - The `PetscSection`
2437: - point - The point

2439:   Output Parameter:
2440: . indices - The constrained dofs

2442:   Level: intermediate

2444:   Fortran Note:
2445:   Call `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`

2447: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2448: @*/
2449: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2450: {
2452:   if (s->bc) {
2453:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2454:   } else *indices = NULL;
2455:   return 0;
2456: }

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

2461:   Not Collective

2463:   Input Parameters:
2464: + s     - The `PetscSection`
2465: . point - The point
2466: - indices - The constrained dofs

2468:   Level: intermediate

2470:   Fortran Note:
2471:   Use `PetscSectionSetConstraintIndicesF90()`

2473: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2474: @*/
2475: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2476: {
2478:   if (s->bc) {
2479:     const PetscInt dof  = s->atlasDof[point];
2480:     const PetscInt cdof = s->bc->atlasDof[point];
2481:     PetscInt       d;

2484:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2485:   }
2486:   return 0;
2487: }

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

2492:   Not Collective

2494:   Input Parameters:
2495: + s     - The `PetscSection`
2496: . field  - The field number
2497: - point - The point

2499:   Output Parameter:
2500: . indices - The constrained dofs sorted in ascending order

2502:   Level: intermediate

2504:   Note:
2505:   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()`.

2507:   Fortran Note:
2508:   Call `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`

2510: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2511: @*/
2512: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2513: {
2516:   PetscSectionCheckValidField(field, s->numFields);
2517:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2518:   return 0;
2519: }

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

2524:   Not Collective

2526:   Input Parameters:
2527: + s       - The `PetscSection`
2528: . point   - The point
2529: . field   - The field number
2530: - indices - The constrained dofs

2532:   Level: intermediate

2534:   Fortran Note:
2535:   Use `PetscSectionSetFieldConstraintIndicesF90()`

2537: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2538: @*/
2539: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2540: {
2542:   if (PetscDefined(USE_DEBUG)) {
2543:     PetscInt nfdof;

2545:     PetscSectionGetFieldConstraintDof(s, point, field, &nfdof);
2547:   }
2548:   PetscSectionCheckValidField(field, s->numFields);
2549:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2550:   return 0;
2551: }

2553: /*@
2554:   PetscSectionPermute - Reorder the section according to the input point permutation

2556:   Collective

2558:   Input Parameters:
2559: + section - The `PetscSection` object
2560: - perm - The point permutation, old point p becomes new point perm[p]

2562:   Output Parameter:
2563: . sectionNew - The permuted `PetscSection`

2565:   Level: intermediate

2567: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`
2568: @*/
2569: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2570: {
2571:   PetscSection    s = section, sNew;
2572:   const PetscInt *perm;
2573:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;

2578:   PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew);
2579:   PetscSectionGetNumFields(s, &numFields);
2580:   if (numFields) PetscSectionSetNumFields(sNew, numFields);
2581:   for (f = 0; f < numFields; ++f) {
2582:     const char *name;
2583:     PetscInt    numComp;

2585:     PetscSectionGetFieldName(s, f, &name);
2586:     PetscSectionSetFieldName(sNew, f, name);
2587:     PetscSectionGetFieldComponents(s, f, &numComp);
2588:     PetscSectionSetFieldComponents(sNew, f, numComp);
2589:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2590:       PetscSectionGetComponentName(s, f, c, &name);
2591:       PetscSectionSetComponentName(sNew, f, c, name);
2592:     }
2593:   }
2594:   ISGetLocalSize(permutation, &numPoints);
2595:   ISGetIndices(permutation, &perm);
2596:   PetscSectionGetChart(s, &pStart, &pEnd);
2597:   PetscSectionSetChart(sNew, pStart, pEnd);
2599:   for (p = pStart; p < pEnd; ++p) {
2600:     PetscInt dof, cdof;

2602:     PetscSectionGetDof(s, p, &dof);
2603:     PetscSectionSetDof(sNew, perm[p], dof);
2604:     PetscSectionGetConstraintDof(s, p, &cdof);
2605:     if (cdof) PetscSectionSetConstraintDof(sNew, perm[p], cdof);
2606:     for (f = 0; f < numFields; ++f) {
2607:       PetscSectionGetFieldDof(s, p, f, &dof);
2608:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2609:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2610:       if (cdof) PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);
2611:     }
2612:   }
2613:   PetscSectionSetUp(sNew);
2614:   for (p = pStart; p < pEnd; ++p) {
2615:     const PetscInt *cind;
2616:     PetscInt        cdof;

2618:     PetscSectionGetConstraintDof(s, p, &cdof);
2619:     if (cdof) {
2620:       PetscSectionGetConstraintIndices(s, p, &cind);
2621:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2622:     }
2623:     for (f = 0; f < numFields; ++f) {
2624:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2625:       if (cdof) {
2626:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2627:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2628:       }
2629:     }
2630:   }
2631:   ISRestoreIndices(permutation, &perm);
2632:   *sectionNew = sNew;
2633:   return 0;
2634: }

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

2639:   Collective

2641:   Input Parameters:
2642: + section   - The `PetscSection`
2643: . obj       - A `PetscObject` which serves as the key for this index
2644: . clSection - `PetscSection` giving the size of the closure of each point
2645: - clPoints  - `IS` giving the points in each closure

2647:   Level: advanced

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

2652: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2653: @*/
2654: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2655: {
2659:   if (section->clObj != obj) PetscSectionResetClosurePermutation(section);
2660:   section->clObj = obj;
2661:   PetscObjectReference((PetscObject)clSection);
2662:   PetscObjectReference((PetscObject)clPoints);
2663:   PetscSectionDestroy(&section->clSection);
2664:   ISDestroy(&section->clPoints);
2665:   section->clSection = clSection;
2666:   section->clPoints  = clPoints;
2667:   return 0;
2668: }

2670: /*@
2671:   PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`

2673:   Collective

2675:   Input Parameters:
2676: + section   - The `PetscSection`
2677: - obj       - A `PetscObject` which serves as the key for this index

2679:   Output Parameters:
2680: + clSection - `PetscSection` giving the size of the closure of each point
2681: - clPoints  - `IS` giving the points in each closure

2683:   Level: advanced

2685: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2686: @*/
2687: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2688: {
2689:   if (section->clObj == obj) {
2690:     if (clSection) *clSection = section->clSection;
2691:     if (clPoints) *clPoints = section->clPoints;
2692:   } else {
2693:     if (clSection) *clSection = NULL;
2694:     if (clPoints) *clPoints = NULL;
2695:   }
2696:   return 0;
2697: }

2699: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2700: {
2701:   PetscInt                    i;
2702:   khiter_t                    iter;
2703:   int                         new_entry;
2704:   PetscSectionClosurePermKey  key = {depth, clSize};
2705:   PetscSectionClosurePermVal *val;

2707:   if (section->clObj != obj) {
2708:     PetscSectionDestroy(&section->clSection);
2709:     ISDestroy(&section->clPoints);
2710:   }
2711:   section->clObj = obj;
2712:   if (!section->clHash) PetscClPermCreate(&section->clHash);
2713:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2714:   val  = &kh_val(section->clHash, iter);
2715:   if (!new_entry) {
2716:     PetscFree(val->perm);
2717:     PetscFree(val->invPerm);
2718:   }
2719:   if (mode == PETSC_COPY_VALUES) {
2720:     PetscMalloc1(clSize, &val->perm);
2721:     PetscArraycpy(val->perm, clPerm, clSize);
2722:   } else if (mode == PETSC_OWN_POINTER) {
2723:     val->perm = clPerm;
2724:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2725:   PetscMalloc1(clSize, &val->invPerm);
2726:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2727:   return 0;
2728: }

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

2733:   Not Collective

2735:   Input Parameters:
2736: + section - The `PetscSection`
2737: . obj     - A `PetscObject` which serves as the key for this index (usually a `DM`)
2738: . depth   - Depth of points on which to apply the given permutation
2739: - perm    - Permutation of the cell dof closure

2741:   Level: intermediate

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

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

2751: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2752: @*/
2753: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2754: {
2755:   const PetscInt *clPerm = NULL;
2756:   PetscInt        clSize = 0;

2758:   if (perm) {
2759:     ISGetLocalSize(perm, &clSize);
2760:     ISGetIndices(perm, &clPerm);
2761:   }
2762:   PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm);
2763:   if (perm) ISRestoreIndices(perm, &clPerm);
2764:   return 0;
2765: }

2767: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2768: {
2769:   if (section->clObj == obj) {
2770:     PetscSectionClosurePermKey k = {depth, size};
2771:     PetscSectionClosurePermVal v;
2772:     PetscClPermGet(section->clHash, k, &v);
2773:     if (perm) *perm = v.perm;
2774:   } else {
2775:     if (perm) *perm = NULL;
2776:   }
2777:   return 0;
2778: }

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

2783:   Not Collective

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

2791:   Output Parameter:
2792: . perm - The dof closure permutation

2794:   Level: intermediate

2796:   Note:
2797:   The user must destroy the `IS` that is returned.

2799: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2800: @*/
2801: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2802: {
2803:   const PetscInt *clPerm;

2805:   PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2806:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2807:   return 0;
2808: }

2810: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2811: {
2812:   if (section->clObj == obj && section->clHash) {
2813:     PetscSectionClosurePermKey k = {depth, size};
2814:     PetscSectionClosurePermVal v;
2815:     PetscClPermGet(section->clHash, k, &v);
2816:     if (perm) *perm = v.invPerm;
2817:   } else {
2818:     if (perm) *perm = NULL;
2819:   }
2820:   return 0;
2821: }

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

2826:   Not Collective

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

2834:   Output Parameters:
2835: . perm - The dof closure permutation

2837:   Level: intermediate

2839:   Note:
2840:   The user must destroy the `IS` that is returned.

2842: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2843: @*/
2844: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2845: {
2846:   const PetscInt *clPerm;

2848:   PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2849:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2850:   return 0;
2851: }

2853: /*@
2854:   PetscSectionGetField - Get the subsection associated with a single field

2856:   Input Parameters:
2857: + s     - The `PetscSection`
2858: - field - The field number

2860:   Output Parameter:
2861: . subs  - The subsection for the given field

2863:   Level: intermediate

2865:   Note:
2866:   Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`

2868: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
2869: @*/
2870: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2871: {
2874:   PetscSectionCheckValidField(field, s->numFields);
2875:   *subs = s->field[field];
2876:   return 0;
2877: }

2879: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2880: PetscFunctionList PetscSectionSymList = NULL;

2882: /*@
2883:   PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.

2885:   Collective

2887:   Input Parameter:
2888: . comm - the MPI communicator

2890:   Output Parameter:
2891: . sym - pointer to the new set of symmetries

2893:   Level: developer

2895: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
2896: @*/
2897: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2898: {
2900:   ISInitializePackage();
2901:   PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView);
2902:   return 0;
2903: }

2905: /*@C
2906:   PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.

2908:   Collective

2910:   Input Parameters:
2911: + sym    - The section symmetry object
2912: - method - The name of the section symmetry type

2914:   Level: developer

2916: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
2917: @*/
2918: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2919: {
2920:   PetscErrorCode (*r)(PetscSectionSym);
2921:   PetscBool match;

2924:   PetscObjectTypeCompare((PetscObject)sym, method, &match);
2925:   if (match) return 0;

2927:   PetscFunctionListFind(PetscSectionSymList, method, &r);
2929:   PetscTryTypeMethod(sym, destroy);
2930:   sym->ops->destroy = NULL;

2932:   (*r)(sym);
2933:   PetscObjectChangeTypeName((PetscObject)sym, method);
2934:   return 0;
2935: }

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

2940:   Not Collective

2942:   Input Parameter:
2943: . sym  - The section symmetry

2945:   Output Parameter:
2946: . type - The index set type name

2948:   Level: developer

2950: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
2951: @*/
2952: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2953: {
2956:   *type = ((PetscObject)sym)->type_name;
2957:   return 0;
2958: }

2960: /*@C
2961:   PetscSectionSymRegister - Registers a new section symmetry implementation

2963:   Not Collective

2965:   Input Parameters:
2966: + name        - The name of a new user-defined creation routine
2967: - create_func - The creation routine itself

2969:   Level: developer

2971:   Notes:
2972:   `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors

2974: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
2975: @*/
2976: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2977: {
2978:   ISInitializePackage();
2979:   PetscFunctionListAdd(&PetscSectionSymList, sname, function);
2980:   return 0;
2981: }

2983: /*@
2984:    PetscSectionSymDestroy - Destroys a section symmetry.

2986:    Collective

2988:    Input Parameters:
2989: .  sym - the section symmetry

2991:    Level: developer

2993: .seealso: [PetscSection](sec_petscsection),  `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSymDestroy()`
2994: @*/
2995: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2996: {
2997:   SymWorkLink link, next;

2999:   if (!*sym) return 0;
3001:   if (--((PetscObject)(*sym))->refct > 0) {
3002:     *sym = NULL;
3003:     return 0;
3004:   }
3005:   if ((*sym)->ops->destroy) (*(*sym)->ops->destroy)(*sym);
3007:   for (link = (*sym)->workin; link; link = next) {
3008:     PetscInt    **perms = (PetscInt **)link->perms;
3009:     PetscScalar **rots  = (PetscScalar **)link->rots;
3010:     PetscFree2(perms, rots);
3011:     next = link->next;
3012:     PetscFree(link);
3013:   }
3014:   (*sym)->workin = NULL;
3015:   PetscHeaderDestroy(sym);
3016:   return 0;
3017: }

3019: /*@C
3020:    PetscSectionSymView - Displays a section symmetry

3022:    Collective

3024:    Input Parameters:
3025: +  sym - the index set
3026: -  viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.

3028:    Level: developer

3030: .seealso:  `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3031: @*/
3032: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3033: {
3035:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer);
3038:   PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer);
3039:   PetscTryTypeMethod(sym, view, viewer);
3040:   return 0;
3041: }

3043: /*@
3044:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3046:   Collective

3048:   Input Parameters:
3049: + section - the section describing data layout
3050: - sym - the symmetry describing the affect of orientation on the access of the data

3052:   Level: developer

3054: .seealso: [PetscSection](sec_petscsection),  `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3055: @*/
3056: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3057: {
3059:   PetscSectionSymDestroy(&(section->sym));
3060:   if (sym) {
3063:     PetscObjectReference((PetscObject)sym);
3064:   }
3065:   section->sym = sym;
3066:   return 0;
3067: }

3069: /*@
3070:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3072:   Not Collective

3074:   Input Parameters:
3075: . section - the section describing data layout

3077:   Output Parameters:
3078: . sym - the symmetry describing the affect of orientation on the access of the data

3080:   Level: developer

3082: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3083: @*/
3084: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3085: {
3087:   *sym = section->sym;
3088:   return 0;
3089: }

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

3094:   Collective

3096:   Input Parameters:
3097: + section - the section describing data layout
3098: . field - the field number
3099: - sym - the symmetry describing the affect of orientation on the access of the data

3101:   Level: developer

3103: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3104: @*/
3105: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3106: {
3108:   PetscSectionCheckValidField(field, section->numFields);
3109:   PetscSectionSetSym(section->field[field], sym);
3110:   return 0;
3111: }

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

3116:   Collective

3118:   Input Parameters:
3119: + section - the section describing data layout
3120: - field - the field number

3122:   Output Parameters:
3123: . sym - the symmetry describing the affect of orientation on the access of the data

3125:   Level: developer

3127: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3128: @*/
3129: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3130: {
3132:   PetscSectionCheckValidField(field, section->numFields);
3133:   *sym = section->field[field]->sym;
3134:   return 0;
3135: }

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

3140:   Not Collective

3142:   Input Parameters:
3143: + section - the section
3144: . numPoints - the number of points
3145: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3146:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3147:     context, see DMPlexGetConeOrientation()).

3149:   Output Parameters:
3150: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3151: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3152:     identity).

3154:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3155: .vb
3156:      const PetscInt    **perms;
3157:      const PetscScalar **rots;
3158:      PetscInt            lOffset;

3160:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3161:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3162:        PetscInt           point = points[2*i], dof, sOffset;
3163:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3164:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3166:        PetscSectionGetDof(section,point,&dof);
3167:        PetscSectionGetOffset(section,point,&sOffset);

3169:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3170:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3171:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3172:        lOffset += dof;
3173:      }
3174:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3175: .ve

3177:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3178: .vb
3179:      const PetscInt    **perms;
3180:      const PetscScalar **rots;
3181:      PetscInt            lOffset;

3183:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3184:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3185:        PetscInt           point = points[2*i], dof, sOffset;
3186:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3187:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3189:        PetscSectionGetDof(section,point,&dof);
3190:        PetscSectionGetOffset(section,point,&sOff);

3192:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3193:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3194:        offset += dof;
3195:      }
3196:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3197: .ve

3199:   Level: developer

3201: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3202: @*/
3203: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3204: {
3205:   PetscSectionSym sym;

3209:   if (perms) *perms = NULL;
3210:   if (rots) *rots = NULL;
3211:   sym = section->sym;
3212:   if (sym && (perms || rots)) {
3213:     SymWorkLink link;

3215:     if (sym->workin) {
3216:       link        = sym->workin;
3217:       sym->workin = sym->workin->next;
3218:     } else {
3219:       PetscNew(&link);
3220:     }
3221:     if (numPoints > link->numPoints) {
3222:       PetscInt    **perms = (PetscInt **)link->perms;
3223:       PetscScalar **rots  = (PetscScalar **)link->rots;
3224:       PetscFree2(perms, rots);
3225:       PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots);
3226:       link->numPoints = numPoints;
3227:     }
3228:     link->next   = sym->workout;
3229:     sym->workout = link;
3230:     PetscArrayzero((PetscInt **)link->perms, numPoints);
3231:     PetscArrayzero((PetscInt **)link->rots, numPoints);
3232:     (*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots);
3233:     if (perms) *perms = link->perms;
3234:     if (rots) *rots = link->rots;
3235:   }
3236:   return 0;
3237: }

3239: /*@C
3240:   PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`

3242:   Not Collective

3244:   Input Parameters:
3245: + section - the section
3246: . numPoints - the number of points
3247: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3248:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3249:     context, see `DMPlexGetConeOrientation()`).

3251:   Output Parameters:
3252: + perms - The permutations for the given orientations: set to NULL at conclusion
3253: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3255:   Level: developer

3257: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3258: @*/
3259: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3260: {
3261:   PetscSectionSym sym;

3264:   sym = section->sym;
3265:   if (sym && (perms || rots)) {
3266:     SymWorkLink *p, link;

3268:     for (p = &sym->workout; (link = *p); p = &link->next) {
3269:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3270:         *p          = link->next;
3271:         link->next  = sym->workin;
3272:         sym->workin = link;
3273:         if (perms) *perms = NULL;
3274:         if (rots) *rots = NULL;
3275:         return 0;
3276:       }
3277:     }
3278:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3279:   }
3280:   return 0;
3281: }

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

3286:   Not Collective

3288:   Input Parameters:
3289: + section - the section
3290: . field - the field of the section
3291: . numPoints - the number of points
3292: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3293:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3294:     context, see `DMPlexGetConeOrientation()`).

3296:   Output Parameters:
3297: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3298: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3299:     identity).

3301:   Level: developer

3303: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3304: @*/
3305: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3306: {
3309:   PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots);
3310:   return 0;
3311: }

3313: /*@C
3314:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`

3316:   Not Collective

3318:   Input Parameters:
3319: + section - the section
3320: . field - the field number
3321: . numPoints - the number of points
3322: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3323:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3324:     context, see `DMPlexGetConeOrientation()`).

3326:   Output Parameters:
3327: + perms - The permutations for the given orientations: set to NULL at conclusion
3328: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3330:   Level: developer

3332: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3333: @*/
3334: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3335: {
3338:   PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots);
3339:   return 0;
3340: }

3342: /*@
3343:   PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible

3345:   Not Collective

3347:   Input Parameter:
3348: . sym - the `PetscSectionSym`

3350:   Output Parameter:
3351: . nsym - the equivalent symmetries

3353:   Level: developer

3355: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3356: @*/
3357: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3358: {
3361:   PetscTryTypeMethod(sym, copy, nsym);
3362:   return 0;
3363: }

3365: /*@
3366:   PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`

3368:   Collective

3370:   Input Parameters:
3371: + sym - the `PetscSectionSym`
3372: - migrationSF - the distribution map from roots to leaves

3374:   Output Parameters:
3375: . dsym - the redistributed symmetries

3377:   Level: developer

3379: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3380: @*/
3381: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3382: {
3386:   PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3387:   return 0;
3388: }

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

3393:   Not Collective

3395:   Input Parameter:
3396: . s - the global `PetscSection`

3398:   Output Parameters:
3399: . flg - the flag

3401:   Level: developer

3403: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3404: @*/
3405: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3406: {
3408:   *flg = s->useFieldOff;
3409:   return 0;
3410: }

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

3415:   Not Collective

3417:   Input Parameters:
3418: + s   - the global `PetscSection`
3419: - flg - the flag

3421:   Level: developer

3423: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3424: @*/
3425: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3426: {
3428:   s->useFieldOff = flg;
3429:   return 0;
3430: }

3432: #define PetscSectionExpandPoints_Loop(TYPE) \
3433:   { \
3434:     PetscInt i, n, o0, o1, size; \
3435:     TYPE    *a0 = (TYPE *)origArray, *a1; \
3436:     PetscSectionGetStorageSize(s, &size); \
3437:     PetscMalloc1(size, &a1); \
3438:     for (i = 0; i < npoints; i++) { \
3439:       PetscSectionGetOffset(origSection, points_[i], &o0); \
3440:       PetscSectionGetOffset(s, i, &o1); \
3441:       PetscSectionGetDof(s, i, &n); \
3442:       PetscMemcpy(&a1[o1], &a0[o0], n *unitsize); \
3443:     } \
3444:     *newArray = (void *)a1; \
3445:   }

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

3450:   Not Collective

3452:   Input Parameters:
3453: + origSection - the `PetscSection` describing the layout of the array
3454: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3455: . origArray - the array; its size must be equal to the storage size of origSection
3456: - points - `IS` with points to extract; its indices must lie in the chart of origSection

3458:   Output Parameters:
3459: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3460: - newArray - the array of the extracted DOFs; its size is the storage size of newSection

3462:   Level: developer

3464: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3465: @*/
3466: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3467: {
3468:   PetscSection    s;
3469:   const PetscInt *points_;
3470:   PetscInt        i, n, npoints, pStart, pEnd;
3471:   PetscMPIInt     unitsize;

3478:   MPI_Type_size(dataType, &unitsize);
3479:   ISGetLocalSize(points, &npoints);
3480:   ISGetIndices(points, &points_);
3481:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3482:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3483:   PetscSectionSetChart(s, 0, npoints);
3484:   for (i = 0; i < npoints; i++) {
3486:     PetscSectionGetDof(origSection, points_[i], &n);
3487:     PetscSectionSetDof(s, i, n);
3488:   }
3489:   PetscSectionSetUp(s);
3490:   if (newArray) {
3491:     if (dataType == MPIU_INT) {
3492:       PetscSectionExpandPoints_Loop(PetscInt);
3493:     } else if (dataType == MPIU_SCALAR) {
3494:       PetscSectionExpandPoints_Loop(PetscScalar);
3495:     } else if (dataType == MPIU_REAL) {
3496:       PetscSectionExpandPoints_Loop(PetscReal);
3497:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3498:   }
3499:   if (newSection) {
3500:     *newSection = s;
3501:   } else {
3502:     PetscSectionDestroy(&s);
3503:   }
3504:   ISRestoreIndices(points, &points_);
3505:   return 0;
3506: }