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.

 35: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionDestroy()`
 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:   Ouput 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:   Level: intermediate

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

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

792:   Not Collective

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

799:   Level: intermediate

801: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
802: @*/
803: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
804: {
806:   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);
807:   s->atlasDof[point - s->pStart] = numDof;
808:   PetscSectionInvalidateMaxDof_Internal(s);
809:   return 0;
810: }

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

815:   Not Collective

817:   Input Parameters:
818: + s - the `PetscSection`
819: . point - the point
820: - numDof - the number of additional dof

822:   Level: intermediate

824: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
825: @*/
826: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
827: {
830:   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);
831:   s->atlasDof[point - s->pStart] += numDof;
832:   PetscSectionInvalidateMaxDof_Internal(s);
833:   return 0;
834: }

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

839:   Not Collective

841:   Input Parameters:
842: + s - the `PetscSection`
843: . point - the point
844: - field - the field

846:   Output Parameter:
847: . numDof - the number of dof

849:   Level: intermediate

851: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
852: @*/
853: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
854: {
857:   PetscSectionCheckValidField(field, s->numFields);
858:   PetscSectionGetDof(s->field[field], point, numDof);
859:   return 0;
860: }

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

865:   Not Collective

867:   Input Parameters:
868: + s - the `PetscSection`
869: . point - the point
870: . field - the field
871: - numDof - the number of dof

873:   Level: intermediate

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

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

888:   Not Collective

890:   Input Parameters:
891: + s - the `PetscSection`
892: . point - the point
893: . field - the field
894: - numDof - the number of dof

896:   Level: intermediate

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

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

911:   Not Collective

913:   Input Parameters:
914: + s - the `PetscSection`
915: - point - the point

917:   Output Parameter:
918: . numDof - the number of dof which are fixed by constraints

920:   Level: intermediate

922: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
923: @*/
924: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
925: {
928:   if (s->bc) {
929:     PetscSectionGetDof(s->bc, point, numDof);
930:   } else *numDof = 0;
931:   return 0;
932: }

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

937:   Not Collective

939:   Input Parameters:
940: + s - the `PetscSection`
941: . point - the point
942: - numDof - the number of dof which are fixed by constraints

944:   Level: intermediate

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

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

961:   Not Collective

963:   Input Parameters:
964: + s - the `PetscSection`
965: . point - the point
966: - numDof - the number of additional dof which are fixed by constraints

968:   Level: intermediate

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

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

985:   Not Collective

987:   Input Parameters:
988: + s - the `PetscSection`
989: . point - the point
990: - field - the field

992:   Output Parameter:
993: . numDof - the number of dof which are fixed by constraints

995:   Level: intermediate

997: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
998: @*/
999: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1000: {
1003:   PetscSectionCheckValidField(field, s->numFields);
1004:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
1005:   return 0;
1006: }

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

1011:   Not Collective

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

1019:   Level: intermediate

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

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

1034:   Not Collective

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

1042:   Level: intermediate

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

1054: /*@
1055:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1057:   Not Collective

1059:   Input Parameter:
1060: . s - the `PetscSection`

1062:   Level: advanced

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

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

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

1081:   Not Collective

1083:   Input Parameter:
1084: . s - the `PetscSection`

1086:   Level: intermediate

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

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

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

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

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

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

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

1142:   Not Collective

1144:   Input Parameters:
1145: . s - the `PetscSection`

1147:   Output Parameter:
1148: . maxDof - the maximum dof

1150:   Level: intermediate

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

1155:   Developer Note:
1156:   The returned number is calculated lazily and stashed.

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

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

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

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

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

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

1183:   Not Collective

1185:   Input Parameter:
1186: . s - the `PetscSection`

1188:   Output Parameter:
1189: . size - the size of an array which can hold all the dofs

1191:   Level: intermediate

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

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

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

1209:   Not Collective

1211:   Input Parameter:
1212: . s - the `PetscSection`

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

1217:   Level: intermediate

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

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

1235: /*@
1236:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1237:   the local section and a `PetscSF` describing the section point overlap.

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

1245:   Output Parameter:
1246: . gsection - The `PetscSection` for the global field layout

1248:   Level: intermediate

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

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

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

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

1356:     PetscSectionSetUpBC(gfs);
1357:     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]);
1358:   }
1359:   gs->setup = PETSC_TRUE;
1360:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1361:   *gsection = gs;
1362:   return 0;
1363: }

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

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

1376:   Output Parameter:
1377: . gsection - The `PetscSection` for the global field layout

1379:   Level: advanced

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

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

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

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

1465: /*@
1466:   PetscSectionGetPointLayout - Get the `PetscLayout` associated with the section points

1468:   Collective

1470:   Input Parameters:
1471: + comm - The `MPI_Comm`
1472: - s    - The `PetscSection`

1474:   Output Parameter:
1475: . layout - The point layout for the section

1477:   Level: advanced

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

1482: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1483: @*/
1484: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1485: {
1486:   PetscInt pStart, pEnd, p, localSize = 0;

1488:   PetscSectionGetChart(s, &pStart, &pEnd);
1489:   for (p = pStart; p < pEnd; ++p) {
1490:     PetscInt dof;

1492:     PetscSectionGetDof(s, p, &dof);
1493:     if (dof >= 0) ++localSize;
1494:   }
1495:   PetscLayoutCreate(comm, layout);
1496:   PetscLayoutSetLocalSize(*layout, localSize);
1497:   PetscLayoutSetBlockSize(*layout, 1);
1498:   PetscLayoutSetUp(*layout);
1499:   return 0;
1500: }

1502: /*@
1503:   PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs.

1505:   Collective

1507:   Input Parameters:
1508: + comm - The `MPI_Comm`
1509: - s    - The `PetscSection`

1511:   Output Parameter:
1512: . layout - The dof layout for the section

1514:   Level: advanced

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

1519: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1520: @*/
1521: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1522: {
1523:   PetscInt pStart, pEnd, p, localSize = 0;

1527:   PetscSectionGetChart(s, &pStart, &pEnd);
1528:   for (p = pStart; p < pEnd; ++p) {
1529:     PetscInt dof, cdof;

1531:     PetscSectionGetDof(s, p, &dof);
1532:     PetscSectionGetConstraintDof(s, p, &cdof);
1533:     if (dof - cdof > 0) localSize += dof - cdof;
1534:   }
1535:   PetscLayoutCreate(comm, layout);
1536:   PetscLayoutSetLocalSize(*layout, localSize);
1537:   PetscLayoutSetBlockSize(*layout, 1);
1538:   PetscLayoutSetUp(*layout);
1539:   return 0;
1540: }

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

1545:   Not Collective

1547:   Input Parameters:
1548: + s - the `PetscSection`
1549: - point - the point

1551:   Output Parameter:
1552: . offset - the offset

1554:   Level: intermediate

1556: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`
1557: @*/
1558: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1559: {
1563:   *offset = s->atlasOff[point - s->pStart];
1564:   return 0;
1565: }

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

1570:   Not Collective

1572:   Input Parameters:
1573: + s - the `PetscSection`
1574: . point - the point
1575: - offset - the offset

1577:   Level: intermediate

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

1582: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1583: @*/
1584: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1585: {
1588:   s->atlasOff[point - s->pStart] = offset;
1589:   return 0;
1590: }

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

1595:   Not Collective

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

1602:   Output Parameter:
1603: . offset - the offset

1605:   Level: intermediate

1607: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1608: @*/
1609: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1610: {
1613:   PetscSectionCheckValidField(field, s->numFields);
1614:   PetscSectionGetOffset(s->field[field], point, offset);
1615:   return 0;
1616: }

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

1621:   Not Collective

1623:   Input Parameters:
1624: + s - the PetscSection
1625: . point - the point
1626: . field - the field
1627: - offset - the offset

1629:   Level: developer

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

1634: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1635: @*/
1636: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1637: {
1639:   PetscSectionCheckValidField(field, s->numFields);
1640:   PetscSectionSetOffset(s->field[field], point, offset);
1641:   return 0;
1642: }

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

1647:   Not Collective

1649:   Input Parameters:
1650: + s - the `PetscSection`
1651: . point - the point
1652: - field - the field

1654:   Output Parameter:
1655: . offset - the offset

1657:   Level: advanced

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

1663: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1664: @*/
1665: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1666: {
1667:   PetscInt off, foff;

1671:   PetscSectionCheckValidField(field, s->numFields);
1672:   PetscSectionGetOffset(s, point, &off);
1673:   PetscSectionGetOffset(s->field[field], point, &foff);
1674:   *offset = foff - off;
1675:   return 0;
1676: }

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

1681:   Not Collective

1683:   Input Parameter:
1684: . s - the `PetscSection`

1686:   Output Parameters:
1687: + start - the minimum offset
1688: - end   - one more than the maximum offset

1690:   Level: intermediate

1692: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1693: @*/
1694: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1695: {
1696:   PetscInt os = 0, oe = 0, pStart, pEnd, p;

1699:   if (s->atlasOff) {
1700:     os = s->atlasOff[0];
1701:     oe = s->atlasOff[0];
1702:   }
1703:   PetscSectionGetChart(s, &pStart, &pEnd);
1704:   for (p = 0; p < pEnd - pStart; ++p) {
1705:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1707:     if (off >= 0) {
1708:       os = PetscMin(os, off);
1709:       oe = PetscMax(oe, off + dof);
1710:     }
1711:   }
1712:   if (start) *start = os;
1713:   if (end) *end = oe;
1714:   return 0;
1715: }

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

1720:   Collective

1722:   Input Parameters:
1723: + s      - the `PetscSection`
1724: . len    - the number of subfields
1725: - fields - the subfield numbers

1727:   Output Parameter:
1728: . subs   - the subsection

1730:   Level: advanced

1732:   Notes:
1733:   The section offsets now refer to a new, smaller vector.

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

1737: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1738: @*/
1739: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1740: {
1741:   PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;

1743:   if (!len) return 0;
1747:   PetscSectionGetNumFields(s, &nF);
1749:   PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1750:   PetscSectionSetNumFields(*subs, len);
1751:   for (f = 0; f < len; ++f) {
1752:     const char *name    = NULL;
1753:     PetscInt    numComp = 0;

1755:     PetscSectionGetFieldName(s, fields[f], &name);
1756:     PetscSectionSetFieldName(*subs, f, name);
1757:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1758:     PetscSectionSetFieldComponents(*subs, f, numComp);
1759:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1760:       PetscSectionGetComponentName(s, fields[f], c, &name);
1761:       PetscSectionSetComponentName(*subs, f, c, name);
1762:     }
1763:   }
1764:   PetscSectionGetChart(s, &pStart, &pEnd);
1765:   PetscSectionSetChart(*subs, pStart, pEnd);
1766:   for (p = pStart; p < pEnd; ++p) {
1767:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1769:     for (f = 0; f < len; ++f) {
1770:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1771:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1772:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1773:       if (cfdof) PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);
1774:       dof += fdof;
1775:       cdof += cfdof;
1776:     }
1777:     PetscSectionSetDof(*subs, p, dof);
1778:     if (cdof) PetscSectionSetConstraintDof(*subs, p, cdof);
1779:     maxCdof = PetscMax(cdof, maxCdof);
1780:   }
1781:   PetscSectionSetUp(*subs);
1782:   if (maxCdof) {
1783:     PetscInt *indices;

1785:     PetscMalloc1(maxCdof, &indices);
1786:     for (p = pStart; p < pEnd; ++p) {
1787:       PetscInt cdof;

1789:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1790:       if (cdof) {
1791:         const PetscInt *oldIndices = NULL;
1792:         PetscInt        fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1794:         for (f = 0; f < len; ++f) {
1795:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1796:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1797:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1798:           PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices);
1799:           for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1800:           numConst += cfdof;
1801:           fOff += fdof;
1802:         }
1803:         PetscSectionSetConstraintIndices(*subs, p, indices);
1804:       }
1805:     }
1806:     PetscFree(indices);
1807:   }
1808:   return 0;
1809: }

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

1814:   Collective

1816:   Input Parameters:
1817: + s     - the input sections
1818: - len   - the number of input sections

1820:   Output Parameter:
1821: . supers - the supersection

1823:   Level: advanced

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

1828:   Developer Note:
1829:   Needs to explain how the sections are composed

1831: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
1832: @*/
1833: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1834: {
1835:   PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1837:   if (!len) return 0;
1838:   for (i = 0; i < len; ++i) {
1839:     PetscInt nf, pStarti, pEndi;

1841:     PetscSectionGetNumFields(s[i], &nf);
1842:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1843:     pStart = PetscMin(pStart, pStarti);
1844:     pEnd   = PetscMax(pEnd, pEndi);
1845:     Nf += nf;
1846:   }
1847:   PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers);
1848:   PetscSectionSetNumFields(*supers, Nf);
1849:   for (i = 0, f = 0; i < len; ++i) {
1850:     PetscInt nf, fi, ci;

1852:     PetscSectionGetNumFields(s[i], &nf);
1853:     for (fi = 0; fi < nf; ++fi, ++f) {
1854:       const char *name    = NULL;
1855:       PetscInt    numComp = 0;

1857:       PetscSectionGetFieldName(s[i], fi, &name);
1858:       PetscSectionSetFieldName(*supers, f, name);
1859:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1860:       PetscSectionSetFieldComponents(*supers, f, numComp);
1861:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1862:         PetscSectionGetComponentName(s[i], fi, ci, &name);
1863:         PetscSectionSetComponentName(*supers, f, ci, name);
1864:       }
1865:     }
1866:   }
1867:   PetscSectionSetChart(*supers, pStart, pEnd);
1868:   for (p = pStart; p < pEnd; ++p) {
1869:     PetscInt dof = 0, cdof = 0;

1871:     for (i = 0, f = 0; i < len; ++i) {
1872:       PetscInt nf, fi, pStarti, pEndi;
1873:       PetscInt fdof = 0, cfdof = 0;

1875:       PetscSectionGetNumFields(s[i], &nf);
1876:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1877:       if ((p < pStarti) || (p >= pEndi)) continue;
1878:       for (fi = 0; fi < nf; ++fi, ++f) {
1879:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1880:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1881:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1882:         if (cfdof) PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);
1883:         dof += fdof;
1884:         cdof += cfdof;
1885:       }
1886:     }
1887:     PetscSectionSetDof(*supers, p, dof);
1888:     if (cdof) PetscSectionSetConstraintDof(*supers, p, cdof);
1889:     maxCdof = PetscMax(cdof, maxCdof);
1890:   }
1891:   PetscSectionSetUp(*supers);
1892:   if (maxCdof) {
1893:     PetscInt *indices;

1895:     PetscMalloc1(maxCdof, &indices);
1896:     for (p = pStart; p < pEnd; ++p) {
1897:       PetscInt cdof;

1899:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1900:       if (cdof) {
1901:         PetscInt dof, numConst = 0, fOff = 0;

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

1907:           PetscSectionGetNumFields(s[i], &nf);
1908:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1909:           if ((p < pStarti) || (p >= pEndi)) continue;
1910:           for (fi = 0; fi < nf; ++fi, ++f) {
1911:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1912:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1913:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1914:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
1915:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1916:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
1917:             numConst += cfdof;
1918:           }
1919:           PetscSectionGetDof(s[i], p, &dof);
1920:           fOff += dof;
1921:         }
1922:         PetscSectionSetConstraintIndices(*supers, p, indices);
1923:       }
1924:     }
1925:     PetscFree(indices);
1926:   }
1927:   return 0;
1928: }

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

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

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

1971:     PetscFindInt(p, numSubpoints, points, &subp);
1972:     if (subp < 0) continue;
1973:     if (!renumberPoints) subp = p;
1974:     for (f = 0; f < numFields; ++f) {
1975:       PetscSectionGetFieldDof(s, p, f, &fdof);
1976:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1977:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1978:       if (cfdof) PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);
1979:     }
1980:     PetscSectionGetDof(s, p, &dof);
1981:     PetscSectionSetDof(*subs, subp, dof);
1982:     PetscSectionGetConstraintDof(s, p, &cdof);
1983:     if (cdof) PetscSectionSetConstraintDof(*subs, subp, cdof);
1984:   }
1985:   PetscSectionSetUp(*subs);
1986:   /* Change offsets to original offsets */
1987:   for (p = pStart; p < pEnd; ++p) {
1988:     PetscInt off, foff = 0;

1990:     PetscFindInt(p, numSubpoints, points, &subp);
1991:     if (subp < 0) continue;
1992:     if (!renumberPoints) subp = p;
1993:     for (f = 0; f < numFields; ++f) {
1994:       PetscSectionGetFieldOffset(s, p, f, &foff);
1995:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1996:     }
1997:     PetscSectionGetOffset(s, p, &off);
1998:     PetscSectionSetOffset(*subs, subp, off);
1999:   }
2000:   /* Copy constraint indices */
2001:   for (subp = spStart; subp < spEnd; ++subp) {
2002:     PetscInt cdof;

2004:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
2005:     if (cdof) {
2006:       for (f = 0; f < numFields; ++f) {
2007:         PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices);
2008:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2009:       }
2010:       PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices);
2011:       PetscSectionSetConstraintIndices(*subs, subp, indices);
2012:     }
2013:   }
2014:   if (subpointMap) ISRestoreIndices(subpointMap, &points);
2015:   return 0;
2016: }

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

2021:   Collective on s

2023:   Input Parameters:
2024: + s           - the `PetscSection`
2025: - subpointMap - a sorted list of points in the original mesh which are in the submesh

2027:   Output Parameter:
2028: . subs - the subsection

2030:   Level: advanced

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

2035:   Developer Note:
2036:   The use of the term Submesh is confusing and unneeded

2038: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2039: @*/
2040: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2041: {
2042:   PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs);
2043:   return 0;
2044: }

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

2049:   Collective on s

2051:   Input Parameters:
2052: + s           - the `PetscSection`
2053: - subpointMap - a sorted list of points in the original mesh which are in the subdomain

2055:   Output Parameter:
2056: . subs - the subsection

2058:   Level: advanced

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

2063:   Developer Notes:
2064:   It is unclear what the difference with `PetscSectionCreateSubmeshSection()` is.

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

2068: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2069: @*/
2070: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2071: {
2072:   PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs);
2073:   return 0;
2074: }

2076: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2077: {
2078:   PetscInt    p;
2079:   PetscMPIInt rank;

2081:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2082:   PetscViewerASCIIPushSynchronized(viewer);
2083:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2084:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2085:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2086:       PetscInt b;

2088:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2089:       if (s->bcIndices) {
2090:         for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
2091:       }
2092:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2093:     } else {
2094:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2095:     }
2096:   }
2097:   PetscViewerFlush(viewer);
2098:   PetscViewerASCIIPopSynchronized(viewer);
2099:   if (s->sym) {
2100:     PetscViewerASCIIPushTab(viewer);
2101:     PetscSectionSymView(s->sym, viewer);
2102:     PetscViewerASCIIPopTab(viewer);
2103:   }
2104:   return 0;
2105: }

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

2110:    Collective

2112:    Input Parameters:
2113: +  A - the PetscSection object to view
2114: .  obj - Optional object that provides the options prefix used for the options
2115: -  name - command line option

2117:    Level: intermediate

2119: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2120: @*/
2121: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2122: {
2124:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
2125:   return 0;
2126: }

2128: /*@C
2129:   PetscSectionView - Views a `PetscSection`

2131:   Collective

2133:   Input Parameters:
2134: + s - the `PetscSection` object to view
2135: - v - the viewer

2137:   Level: beginner

2139:   Note:
2140:   `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2141:   distribution independent data, such as dofs, offsets, constraint dofs,
2142:   and constraint indices. Points that have negative dofs, for instance,
2143:   are not saved as they represent points owned by other processes.
2144:   Point numbering and rank assignment is currently not stored.
2145:   The saved section can be loaded with `PetscSectionLoad()`.

2147: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2148: @*/
2149: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2150: {
2151:   PetscBool isascii, ishdf5;
2152:   PetscInt  f;

2155:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);
2157:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
2158:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2159:   if (isascii) {
2160:     PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer);
2161:     if (s->numFields) {
2162:       PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields);
2163:       for (f = 0; f < s->numFields; ++f) {
2164:         PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
2165:         PetscSectionView_ASCII(s->field[f], viewer);
2166:       }
2167:     } else {
2168:       PetscSectionView_ASCII(s, viewer);
2169:     }
2170:   } else if (ishdf5) {
2171: #if PetscDefined(HAVE_HDF5)
2172:     PetscSectionView_HDF5_Internal(s, viewer);
2173: #else
2174:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2175: #endif
2176:   }
2177:   return 0;
2178: }

2180: /*@C
2181:   PetscSectionLoad - Loads a `PetscSection`

2183:   Collective

2185:   Input Parameters:
2186: + s - the `PetscSection` object to load
2187: - v - the viewer

2189:   Level: beginner

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

2199: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2200: @*/
2201: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2202: {
2203:   PetscBool ishdf5;

2207:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2208:   if (ishdf5) {
2209: #if PetscDefined(HAVE_HDF5)
2210:     PetscSectionLoad_HDF5_Internal(s, viewer);
2211:     return 0;
2212: #else
2213:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2214: #endif
2215:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2216: }

2218: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2219: {
2220:   PetscSectionClosurePermVal clVal;

2222:   if (!section->clHash) return 0;
2223:   kh_foreach_value(section->clHash, clVal, {
2224:     PetscFree(clVal.perm);
2225:     PetscFree(clVal.invPerm);
2226:   });
2227:   kh_destroy(ClPerm, section->clHash);
2228:   section->clHash = NULL;
2229:   return 0;
2230: }

2232: /*@
2233:   PetscSectionReset - Frees all section data.

2235:   Not Collective

2237:   Input Parameters:
2238: . s - the `PetscSection`

2240:   Level: beginner

2242: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2243: @*/
2244: PetscErrorCode PetscSectionReset(PetscSection s)
2245: {
2246:   PetscInt f, c;

2249:   for (f = 0; f < s->numFields; ++f) {
2250:     PetscSectionDestroy(&s->field[f]);
2251:     PetscFree(s->fieldNames[f]);
2252:     for (c = 0; c < s->numFieldComponents[f]; ++c) PetscFree(s->compNames[f][c]);
2253:     PetscFree(s->compNames[f]);
2254:   }
2255:   PetscFree(s->numFieldComponents);
2256:   PetscFree(s->fieldNames);
2257:   PetscFree(s->compNames);
2258:   PetscFree(s->field);
2259:   PetscSectionDestroy(&s->bc);
2260:   PetscFree(s->bcIndices);
2261:   PetscFree2(s->atlasDof, s->atlasOff);
2262:   PetscSectionDestroy(&s->clSection);
2263:   ISDestroy(&s->clPoints);
2264:   ISDestroy(&s->perm);
2265:   PetscSectionResetClosurePermutation(s);
2266:   PetscSectionSymDestroy(&s->sym);
2267:   PetscSectionDestroy(&s->clSection);
2268:   ISDestroy(&s->clPoints);
2269:   PetscSectionInvalidateMaxDof_Internal(s);
2270:   s->pStart    = -1;
2271:   s->pEnd      = -1;
2272:   s->maxDof    = 0;
2273:   s->setup     = PETSC_FALSE;
2274:   s->numFields = 0;
2275:   s->clObj     = NULL;
2276:   return 0;
2277: }

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

2282:   Not Collective

2284:   Input Parameters:
2285: . s - the `PetscSection`

2287:   Level: beginner

2289: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2290: @*/
2291: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2292: {
2293:   if (!*s) return 0;
2295:   if (--((PetscObject)(*s))->refct > 0) {
2296:     *s = NULL;
2297:     return 0;
2298:   }
2299:   PetscSectionReset(*s);
2300:   PetscHeaderDestroy(s);
2301:   return 0;
2302: }

2304: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2305: {
2306:   const PetscInt p = point - s->pStart;

2309:   *values = &baseArray[s->atlasOff[p]];
2310:   return 0;
2311: }

2313: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2314: {
2315:   PetscInt      *array;
2316:   const PetscInt p           = point - s->pStart;
2317:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2318:   PetscInt       cDim        = 0;

2321:   PetscSectionGetConstraintDof(s, p, &cDim);
2322:   array = &baseArray[s->atlasOff[p]];
2323:   if (!cDim) {
2324:     if (orientation >= 0) {
2325:       const PetscInt dim = s->atlasDof[p];
2326:       PetscInt       i;

2328:       if (mode == INSERT_VALUES) {
2329:         for (i = 0; i < dim; ++i) array[i] = values[i];
2330:       } else {
2331:         for (i = 0; i < dim; ++i) array[i] += values[i];
2332:       }
2333:     } else {
2334:       PetscInt offset = 0;
2335:       PetscInt j      = -1, field, i;

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

2340:         for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2341:         offset += dim;
2342:       }
2343:     }
2344:   } else {
2345:     if (orientation >= 0) {
2346:       const PetscInt  dim  = s->atlasDof[p];
2347:       PetscInt        cInd = 0, i;
2348:       const PetscInt *cDof;

2350:       PetscSectionGetConstraintIndices(s, point, &cDof);
2351:       if (mode == INSERT_VALUES) {
2352:         for (i = 0; i < dim; ++i) {
2353:           if ((cInd < cDim) && (i == cDof[cInd])) {
2354:             ++cInd;
2355:             continue;
2356:           }
2357:           array[i] = values[i];
2358:         }
2359:       } else {
2360:         for (i = 0; i < dim; ++i) {
2361:           if ((cInd < cDim) && (i == cDof[cInd])) {
2362:             ++cInd;
2363:             continue;
2364:           }
2365:           array[i] += values[i];
2366:         }
2367:       }
2368:     } else {
2369:       const PetscInt *cDof;
2370:       PetscInt        offset  = 0;
2371:       PetscInt        cOffset = 0;
2372:       PetscInt        j       = 0, field;

2374:       PetscSectionGetConstraintIndices(s, point, &cDof);
2375:       for (field = 0; field < s->numFields; ++field) {
2376:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2377:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2378:         const PetscInt sDim = dim - tDim;
2379:         PetscInt       cInd = 0, i, k;

2381:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2382:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2383:             ++cInd;
2384:             continue;
2385:           }
2386:           array[j] = values[k];
2387:         }
2388:         offset += dim;
2389:         cOffset += dim - tDim;
2390:       }
2391:     }
2392:   }
2393:   return 0;
2394: }

2396: /*@C
2397:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2399:   Not Collective

2401:   Input Parameter:
2402: . s - The `PetscSection`

2404:   Output Parameter:
2405: . hasConstraints - flag indicating that the section has constrained dofs

2407:   Level: intermediate

2409: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2410: @*/
2411: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2412: {
2415:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2416:   return 0;
2417: }

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

2422:   Not Collective

2424:   Input Parameters:
2425: + s     - The `PetscSection`
2426: - point - The point

2428:   Output Parameter:
2429: . indices - The constrained dofs

2431:   Level: intermediate

2433:   Fortran Note:
2434:   Call `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`

2436: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2437: @*/
2438: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2439: {
2441:   if (s->bc) {
2442:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2443:   } else *indices = NULL;
2444:   return 0;
2445: }

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

2450:   Not Collective

2452:   Input Parameters:
2453: + s     - The `PetscSection`
2454: . point - The point
2455: - indices - The constrained dofs

2457:   Level: intermediate

2459:   Fortran Note:
2460:   Use `PetscSectionSetConstraintIndicesF90()`

2462: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2463: @*/
2464: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2465: {
2467:   if (s->bc) {
2468:     const PetscInt dof  = s->atlasDof[point];
2469:     const PetscInt cdof = s->bc->atlasDof[point];
2470:     PetscInt       d;

2473:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2474:   }
2475:   return 0;
2476: }

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

2481:   Not Collective

2483:   Input Parameters:
2484: + s     - The `PetscSection`
2485: . field  - The field number
2486: - point - The point

2488:   Output Parameter:
2489: . indices - The constrained dofs sorted in ascending order

2491:   Level: intermediate

2493:   Note:
2494:   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()`.

2496:   Fortran Note:
2497:   Call `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`

2499: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2500: @*/
2501: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2502: {
2505:   PetscSectionCheckValidField(field, s->numFields);
2506:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2507:   return 0;
2508: }

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

2513:   Not Collective

2515:   Input Parameters:
2516: + s       - The `PetscSection`
2517: . point   - The point
2518: . field   - The field number
2519: - indices - The constrained dofs

2521:   Level: intermediate

2523:   Fortran Note:
2524:   Use `PetscSectionSetFieldConstraintIndicesF90()`

2526: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2527: @*/
2528: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2529: {
2531:   if (PetscDefined(USE_DEBUG)) {
2532:     PetscInt nfdof;

2534:     PetscSectionGetFieldConstraintDof(s, point, field, &nfdof);
2536:   }
2537:   PetscSectionCheckValidField(field, s->numFields);
2538:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2539:   return 0;
2540: }

2542: /*@
2543:   PetscSectionPermute - Reorder the section according to the input point permutation

2545:   Collective

2547:   Input Parameters:
2548: + section - The `PetscSection` object
2549: - perm - The point permutation, old point p becomes new point perm[p]

2551:   Output Parameter:
2552: . sectionNew - The permuted `PetscSection`

2554:   Level: intermediate

2556: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`
2557: @*/
2558: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2559: {
2560:   PetscSection    s = section, sNew;
2561:   const PetscInt *perm;
2562:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;

2567:   PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew);
2568:   PetscSectionGetNumFields(s, &numFields);
2569:   if (numFields) PetscSectionSetNumFields(sNew, numFields);
2570:   for (f = 0; f < numFields; ++f) {
2571:     const char *name;
2572:     PetscInt    numComp;

2574:     PetscSectionGetFieldName(s, f, &name);
2575:     PetscSectionSetFieldName(sNew, f, name);
2576:     PetscSectionGetFieldComponents(s, f, &numComp);
2577:     PetscSectionSetFieldComponents(sNew, f, numComp);
2578:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2579:       PetscSectionGetComponentName(s, f, c, &name);
2580:       PetscSectionSetComponentName(sNew, f, c, name);
2581:     }
2582:   }
2583:   ISGetLocalSize(permutation, &numPoints);
2584:   ISGetIndices(permutation, &perm);
2585:   PetscSectionGetChart(s, &pStart, &pEnd);
2586:   PetscSectionSetChart(sNew, pStart, pEnd);
2588:   for (p = pStart; p < pEnd; ++p) {
2589:     PetscInt dof, cdof;

2591:     PetscSectionGetDof(s, p, &dof);
2592:     PetscSectionSetDof(sNew, perm[p], dof);
2593:     PetscSectionGetConstraintDof(s, p, &cdof);
2594:     if (cdof) PetscSectionSetConstraintDof(sNew, perm[p], cdof);
2595:     for (f = 0; f < numFields; ++f) {
2596:       PetscSectionGetFieldDof(s, p, f, &dof);
2597:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2598:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2599:       if (cdof) PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);
2600:     }
2601:   }
2602:   PetscSectionSetUp(sNew);
2603:   for (p = pStart; p < pEnd; ++p) {
2604:     const PetscInt *cind;
2605:     PetscInt        cdof;

2607:     PetscSectionGetConstraintDof(s, p, &cdof);
2608:     if (cdof) {
2609:       PetscSectionGetConstraintIndices(s, p, &cind);
2610:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2611:     }
2612:     for (f = 0; f < numFields; ++f) {
2613:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2614:       if (cdof) {
2615:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2616:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2617:       }
2618:     }
2619:   }
2620:   ISRestoreIndices(permutation, &perm);
2621:   *sectionNew = sNew;
2622:   return 0;
2623: }

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

2628:   Collective

2630:   Input Parameters:
2631: + section   - The `PetscSection`
2632: . obj       - A `PetscObject` which serves as the key for this index
2633: . clSection - `PetscSection` giving the size of the closure of each point
2634: - clPoints  - `IS` giving the points in each closure

2636:   Level: advanced

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

2641: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2642: @*/
2643: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2644: {
2648:   if (section->clObj != obj) PetscSectionResetClosurePermutation(section);
2649:   section->clObj = obj;
2650:   PetscObjectReference((PetscObject)clSection);
2651:   PetscObjectReference((PetscObject)clPoints);
2652:   PetscSectionDestroy(&section->clSection);
2653:   ISDestroy(&section->clPoints);
2654:   section->clSection = clSection;
2655:   section->clPoints  = clPoints;
2656:   return 0;
2657: }

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

2662:   Collective

2664:   Input Parameters:
2665: + section   - The `PetscSection`
2666: - obj       - A `PetscObject` which serves as the key for this index

2668:   Output Parameters:
2669: + clSection - `PetscSection` giving the size of the closure of each point
2670: - clPoints  - `IS` giving the points in each closure

2672:   Level: advanced

2674: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2675: @*/
2676: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2677: {
2678:   if (section->clObj == obj) {
2679:     if (clSection) *clSection = section->clSection;
2680:     if (clPoints) *clPoints = section->clPoints;
2681:   } else {
2682:     if (clSection) *clSection = NULL;
2683:     if (clPoints) *clPoints = NULL;
2684:   }
2685:   return 0;
2686: }

2688: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2689: {
2690:   PetscInt                    i;
2691:   khiter_t                    iter;
2692:   int                         new_entry;
2693:   PetscSectionClosurePermKey  key = {depth, clSize};
2694:   PetscSectionClosurePermVal *val;

2696:   if (section->clObj != obj) {
2697:     PetscSectionDestroy(&section->clSection);
2698:     ISDestroy(&section->clPoints);
2699:   }
2700:   section->clObj = obj;
2701:   if (!section->clHash) PetscClPermCreate(&section->clHash);
2702:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2703:   val  = &kh_val(section->clHash, iter);
2704:   if (!new_entry) {
2705:     PetscFree(val->perm);
2706:     PetscFree(val->invPerm);
2707:   }
2708:   if (mode == PETSC_COPY_VALUES) {
2709:     PetscMalloc1(clSize, &val->perm);
2710:     PetscArraycpy(val->perm, clPerm, clSize);
2711:   } else if (mode == PETSC_OWN_POINTER) {
2712:     val->perm = clPerm;
2713:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2714:   PetscMalloc1(clSize, &val->invPerm);
2715:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2716:   return 0;
2717: }

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

2722:   Not Collective

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

2730:   Level: intermediate

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

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

2740: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2741: @*/
2742: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2743: {
2744:   const PetscInt *clPerm = NULL;
2745:   PetscInt        clSize = 0;

2747:   if (perm) {
2748:     ISGetLocalSize(perm, &clSize);
2749:     ISGetIndices(perm, &clPerm);
2750:   }
2751:   PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm);
2752:   if (perm) ISRestoreIndices(perm, &clPerm);
2753:   return 0;
2754: }

2756: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2757: {
2758:   if (section->clObj == obj) {
2759:     PetscSectionClosurePermKey k = {depth, size};
2760:     PetscSectionClosurePermVal v;
2761:     PetscClPermGet(section->clHash, k, &v);
2762:     if (perm) *perm = v.perm;
2763:   } else {
2764:     if (perm) *perm = NULL;
2765:   }
2766:   return 0;
2767: }

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

2772:   Not Collective

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

2780:   Output Parameter:
2781: . perm - The dof closure permutation

2783:   Level: intermediate

2785:   Note:
2786:   The user must destroy the `IS` that is returned.

2788: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2789: @*/
2790: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2791: {
2792:   const PetscInt *clPerm;

2794:   PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2795:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2796:   return 0;
2797: }

2799: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2800: {
2801:   if (section->clObj == obj && section->clHash) {
2802:     PetscSectionClosurePermKey k = {depth, size};
2803:     PetscSectionClosurePermVal v;
2804:     PetscClPermGet(section->clHash, k, &v);
2805:     if (perm) *perm = v.invPerm;
2806:   } else {
2807:     if (perm) *perm = NULL;
2808:   }
2809:   return 0;
2810: }

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

2815:   Not Collective

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

2823:   Output Parameters:
2824: . perm - The dof closure permutation

2826:   Level: intermediate

2828:   Note:
2829:   The user must destroy the `IS` that is returned.

2831: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2832: @*/
2833: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2834: {
2835:   const PetscInt *clPerm;

2837:   PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2838:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2839:   return 0;
2840: }

2842: /*@
2843:   PetscSectionGetField - Get the subsection associated with a single field

2845:   Input Parameters:
2846: + s     - The `PetscSection`
2847: - field - The field number

2849:   Output Parameter:
2850: . subs  - The subsection for the given field

2852:   Level: intermediate

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

2857: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
2858: @*/
2859: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2860: {
2863:   PetscSectionCheckValidField(field, s->numFields);
2864:   *subs = s->field[field];
2865:   return 0;
2866: }

2868: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2869: PetscFunctionList PetscSectionSymList = NULL;

2871: /*@
2872:   PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.

2874:   Collective

2876:   Input Parameter:
2877: . comm - the MPI communicator

2879:   Output Parameter:
2880: . sym - pointer to the new set of symmetries

2882:   Level: developer

2884: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
2885: @*/
2886: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2887: {
2889:   ISInitializePackage();
2890:   PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView);
2891:   return 0;
2892: }

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

2897:   Collective

2899:   Input Parameters:
2900: + sym    - The section symmetry object
2901: - method - The name of the section symmetry type

2903:   Level: developer

2905: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
2906: @*/
2907: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2908: {
2909:   PetscErrorCode (*r)(PetscSectionSym);
2910:   PetscBool match;

2913:   PetscObjectTypeCompare((PetscObject)sym, method, &match);
2914:   if (match) return 0;

2916:   PetscFunctionListFind(PetscSectionSymList, method, &r);
2918:   PetscTryTypeMethod(sym, destroy);
2919:   sym->ops->destroy = NULL;

2921:   (*r)(sym);
2922:   PetscObjectChangeTypeName((PetscObject)sym, method);
2923:   return 0;
2924: }

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

2929:   Not Collective

2931:   Input Parameter:
2932: . sym  - The section symmetry

2934:   Output Parameter:
2935: . type - The index set type name

2937:   Level: developer

2939: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
2940: @*/
2941: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2942: {
2945:   *type = ((PetscObject)sym)->type_name;
2946:   return 0;
2947: }

2949: /*@C
2950:   PetscSectionSymRegister - Registers a new section symmetry implementation

2952:   Not Collective

2954:   Input Parameters:
2955: + name        - The name of a new user-defined creation routine
2956: - create_func - The creation routine itself

2958:   Level: developer

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

2963: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
2964: @*/
2965: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2966: {
2967:   ISInitializePackage();
2968:   PetscFunctionListAdd(&PetscSectionSymList, sname, function);
2969:   return 0;
2970: }

2972: /*@
2973:    PetscSectionSymDestroy - Destroys a section symmetry.

2975:    Collective

2977:    Input Parameters:
2978: .  sym - the section symmetry

2980:    Level: developer

2982: .seealso: [PetscSection](sec_petscsection),  `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSymDestroy()`
2983: @*/
2984: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2985: {
2986:   SymWorkLink link, next;

2988:   if (!*sym) return 0;
2990:   if (--((PetscObject)(*sym))->refct > 0) {
2991:     *sym = NULL;
2992:     return 0;
2993:   }
2994:   if ((*sym)->ops->destroy) (*(*sym)->ops->destroy)(*sym);
2996:   for (link = (*sym)->workin; link; link = next) {
2997:     PetscInt    **perms = (PetscInt **)link->perms;
2998:     PetscScalar **rots  = (PetscScalar **)link->rots;
2999:     PetscFree2(perms, rots);
3000:     next = link->next;
3001:     PetscFree(link);
3002:   }
3003:   (*sym)->workin = NULL;
3004:   PetscHeaderDestroy(sym);
3005:   return 0;
3006: }

3008: /*@C
3009:    PetscSectionSymView - Displays a section symmetry

3011:    Collective

3013:    Input Parameters:
3014: +  sym - the index set
3015: -  viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.

3017:    Level: developer

3019: .seealso:  `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3020: @*/
3021: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3022: {
3024:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer);
3027:   PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer);
3028:   PetscTryTypeMethod(sym, view, viewer);
3029:   return 0;
3030: }

3032: /*@
3033:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3035:   Collective

3037:   Input Parameters:
3038: + section - the section describing data layout
3039: - sym - the symmetry describing the affect of orientation on the access of the data

3041:   Level: developer

3043: .seealso: [PetscSection](sec_petscsection),  `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3044: @*/
3045: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3046: {
3048:   PetscSectionSymDestroy(&(section->sym));
3049:   if (sym) {
3052:     PetscObjectReference((PetscObject)sym);
3053:   }
3054:   section->sym = sym;
3055:   return 0;
3056: }

3058: /*@
3059:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3061:   Not Collective

3063:   Input Parameters:
3064: . section - the section describing data layout

3066:   Output Parameters:
3067: . sym - the symmetry describing the affect of orientation on the access of the data

3069:   Level: developer

3071: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3072: @*/
3073: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3074: {
3076:   *sym = section->sym;
3077:   return 0;
3078: }

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

3083:   Collective

3085:   Input Parameters:
3086: + section - the section describing data layout
3087: . field - the field number
3088: - sym - the symmetry describing the affect of orientation on the access of the data

3090:   Level: developer

3092: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3093: @*/
3094: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3095: {
3097:   PetscSectionCheckValidField(field, section->numFields);
3098:   PetscSectionSetSym(section->field[field], sym);
3099:   return 0;
3100: }

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

3105:   Collective

3107:   Input Parameters:
3108: + section - the section describing data layout
3109: - field - the field number

3111:   Output Parameters:
3112: . sym - the symmetry describing the affect of orientation on the access of the data

3114:   Level: developer

3116: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3117: @*/
3118: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3119: {
3121:   PetscSectionCheckValidField(field, section->numFields);
3122:   *sym = section->field[field]->sym;
3123:   return 0;
3124: }

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

3129:   Not Collective

3131:   Input Parameters:
3132: + section - the section
3133: . numPoints - the number of points
3134: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3135:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3136:     context, see DMPlexGetConeOrientation()).

3138:   Output Parameters:
3139: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3140: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3141:     identity).

3143:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3144: .vb
3145:      const PetscInt    **perms;
3146:      const PetscScalar **rots;
3147:      PetscInt            lOffset;

3149:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3150:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3151:        PetscInt           point = points[2*i], dof, sOffset;
3152:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3153:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3155:        PetscSectionGetDof(section,point,&dof);
3156:        PetscSectionGetOffset(section,point,&sOffset);

3158:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3159:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3160:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3161:        lOffset += dof;
3162:      }
3163:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3164: .ve

3166:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3167: .vb
3168:      const PetscInt    **perms;
3169:      const PetscScalar **rots;
3170:      PetscInt            lOffset;

3172:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3173:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3174:        PetscInt           point = points[2*i], dof, sOffset;
3175:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3176:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3178:        PetscSectionGetDof(section,point,&dof);
3179:        PetscSectionGetOffset(section,point,&sOff);

3181:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3182:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3183:        offset += dof;
3184:      }
3185:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3186: .ve

3188:   Level: developer

3190: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3191: @*/
3192: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3193: {
3194:   PetscSectionSym sym;

3198:   if (perms) *perms = NULL;
3199:   if (rots) *rots = NULL;
3200:   sym = section->sym;
3201:   if (sym && (perms || rots)) {
3202:     SymWorkLink link;

3204:     if (sym->workin) {
3205:       link        = sym->workin;
3206:       sym->workin = sym->workin->next;
3207:     } else {
3208:       PetscNew(&link);
3209:     }
3210:     if (numPoints > link->numPoints) {
3211:       PetscInt    **perms = (PetscInt **)link->perms;
3212:       PetscScalar **rots  = (PetscScalar **)link->rots;
3213:       PetscFree2(perms, rots);
3214:       PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots);
3215:       link->numPoints = numPoints;
3216:     }
3217:     link->next   = sym->workout;
3218:     sym->workout = link;
3219:     PetscArrayzero((PetscInt **)link->perms, numPoints);
3220:     PetscArrayzero((PetscInt **)link->rots, numPoints);
3221:     (*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots);
3222:     if (perms) *perms = link->perms;
3223:     if (rots) *rots = link->rots;
3224:   }
3225:   return 0;
3226: }

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

3231:   Not Collective

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

3240:   Output Parameters:
3241: + perms - The permutations for the given orientations: set to NULL at conclusion
3242: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3244:   Level: developer

3246: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3247: @*/
3248: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3249: {
3250:   PetscSectionSym sym;

3253:   sym = section->sym;
3254:   if (sym && (perms || rots)) {
3255:     SymWorkLink *p, link;

3257:     for (p = &sym->workout; (link = *p); p = &link->next) {
3258:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3259:         *p          = link->next;
3260:         link->next  = sym->workin;
3261:         sym->workin = link;
3262:         if (perms) *perms = NULL;
3263:         if (rots) *rots = NULL;
3264:         return 0;
3265:       }
3266:     }
3267:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3268:   }
3269:   return 0;
3270: }

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

3275:   Not Collective

3277:   Input Parameters:
3278: + section - the section
3279: . field - the field of the section
3280: . numPoints - the number of points
3281: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3282:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3283:     context, see `DMPlexGetConeOrientation()`).

3285:   Output Parameters:
3286: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3287: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3288:     identity).

3290:   Level: developer

3292: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3293: @*/
3294: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3295: {
3298:   PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots);
3299:   return 0;
3300: }

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

3305:   Not Collective

3307:   Input Parameters:
3308: + section - the section
3309: . field - the field number
3310: . numPoints - the number of points
3311: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3312:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3313:     context, see `DMPlexGetConeOrientation()`).

3315:   Output Parameters:
3316: + perms - The permutations for the given orientations: set to NULL at conclusion
3317: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3319:   Level: developer

3321: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3322: @*/
3323: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3324: {
3327:   PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots);
3328:   return 0;
3329: }

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

3334:   Not Collective

3336:   Input Parameter:
3337: . sym - the `PetscSectionSym`

3339:   Output Parameter:
3340: . nsym - the equivalent symmetries

3342:   Level: developer

3344: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3345: @*/
3346: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3347: {
3350:   PetscTryTypeMethod(sym, copy, nsym);
3351:   return 0;
3352: }

3354: /*@
3355:   PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`

3357:   Collective

3359:   Input Parameters:
3360: + sym - the `PetscSectionSym`
3361: - migrationSF - the distribution map from roots to leaves

3363:   Output Parameters:
3364: . dsym - the redistributed symmetries

3366:   Level: developer

3368: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3369: @*/
3370: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3371: {
3375:   PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3376:   return 0;
3377: }

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

3382:   Not Collective

3384:   Input Parameter:
3385: . s - the global `PetscSection`

3387:   Output Parameters:
3388: . flg - the flag

3390:   Level: developer

3392: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3393: @*/
3394: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3395: {
3397:   *flg = s->useFieldOff;
3398:   return 0;
3399: }

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

3404:   Not Collective

3406:   Input Parameters:
3407: + s   - the global `PetscSection`
3408: - flg - the flag

3410:   Level: developer

3412: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3413: @*/
3414: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3415: {
3417:   s->useFieldOff = flg;
3418:   return 0;
3419: }

3421: #define PetscSectionExpandPoints_Loop(TYPE) \
3422:   { \
3423:     PetscInt i, n, o0, o1, size; \
3424:     TYPE    *a0 = (TYPE *)origArray, *a1; \
3425:     PetscSectionGetStorageSize(s, &size); \
3426:     PetscMalloc1(size, &a1); \
3427:     for (i = 0; i < npoints; i++) { \
3428:       PetscSectionGetOffset(origSection, points_[i], &o0); \
3429:       PetscSectionGetOffset(s, i, &o1); \
3430:       PetscSectionGetDof(s, i, &n); \
3431:       PetscMemcpy(&a1[o1], &a0[o0], n *unitsize); \
3432:     } \
3433:     *newArray = (void *)a1; \
3434:   }

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

3439:   Not Collective

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

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

3451:   Level: developer

3453: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3454: @*/
3455: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3456: {
3457:   PetscSection    s;
3458:   const PetscInt *points_;
3459:   PetscInt        i, n, npoints, pStart, pEnd;
3460:   PetscMPIInt     unitsize;

3467:   MPI_Type_size(dataType, &unitsize);
3468:   ISGetLocalSize(points, &npoints);
3469:   ISGetIndices(points, &points_);
3470:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3471:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3472:   PetscSectionSetChart(s, 0, npoints);
3473:   for (i = 0; i < npoints; i++) {
3475:     PetscSectionGetDof(origSection, points_[i], &n);
3476:     PetscSectionSetDof(s, i, n);
3477:   }
3478:   PetscSectionSetUp(s);
3479:   if (newArray) {
3480:     if (dataType == MPIU_INT) {
3481:       PetscSectionExpandPoints_Loop(PetscInt);
3482:     } else if (dataType == MPIU_SCALAR) {
3483:       PetscSectionExpandPoints_Loop(PetscScalar);
3484:     } else if (dataType == MPIU_REAL) {
3485:       PetscSectionExpandPoints_Loop(PetscReal);
3486:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3487:   }
3488:   if (newSection) {
3489:     *newSection = s;
3490:   } else {
3491:     PetscSectionDestroy(&s);
3492:   }
3493:   ISRestoreIndices(points, &points_);
3494:   return 0;
3495: }