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 space 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: $       PetscSectionCreate(MPI_Comm,PetscSection *);
 24: $       PetscSectionSetNumFields(PetscSection, numFields);
 25: $       PetscSectionSetChart(PetscSection,low,high);
 26: $       PetscSectionSetDof(PetscSection,point,numdof);
 27: $       PetscSectionSetUp(PetscSection);
 28: $       PetscSectionGetOffset(PetscSection,point,PetscInt *);
 29: $       PetscSectionDestroy(PetscSection);

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

 34: .seealso: `PetscSection`, `PetscSectionDestroy()`
 35: @*/
 36: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
 37: {
 39:   ISInitializePackage();

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

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

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

 69:   Collective

 71:   Input Parameter:
 72: . section - the PetscSection

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

 77:   Level: intermediate

 79: .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
 80: @*/
 81: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 82: {
 83:   PetscSectionSym sym;
 84:   IS              perm;
 85:   PetscInt        numFields, f, c, pStart, pEnd, p;

 89:   PetscSectionReset(newSection);
 90:   PetscSectionGetNumFields(section, &numFields);
 91:   if (numFields) PetscSectionSetNumFields(newSection, numFields);
 92:   for (f = 0; f < numFields; ++f) {
 93:     const char *fieldName = NULL, *compName = NULL;
 94:     PetscInt    numComp = 0;

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

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

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

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

156:   Collective

158:   Input Parameter:
159: . section - the PetscSection

161:   Output Parameter:
162: . newSection - the copy

164:   Level: beginner

166: .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
167: @*/
168: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
169: {
172:   PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection);
173:   PetscSectionCopy(section, *newSection);
174:   return 0;
175: }

177: /*@
178:   PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database

180:   Collective

182:   Input Parameter:
183: . section - the PetscSection

185:   Options Database:
186: . -petscsection_point_major - PETSC_TRUE for point-major order

188:   Level: intermediate

190: .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
191: @*/
192: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
193: {
195:   PetscObjectOptionsBegin((PetscObject)s);
196:   PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
197:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
198:   PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject);
199:   PetscOptionsEnd();
200:   PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view");
201:   return 0;
202: }

204: /*@
205:   PetscSectionCompare - Compares two sections

207:   Collective

209:   Input Parameters:
210: + s1 - the first PetscSection
211: - s2 - the second PetscSection

213:   Output Parameter:
214: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise

216:   Level: intermediate

218:   Notes:
219:   Field names are disregarded.

221: .seealso: `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
222: @*/
223: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
224: {
225:   PetscInt        pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
226:   const PetscInt *idx1, *idx2;
227:   IS              perm1, perm2;
228:   PetscBool       flg;
229:   PetscMPIInt     mflg;

234:   flg = PETSC_FALSE;

236:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg);
237:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
238:     *congruent = PETSC_FALSE;
239:     return 0;
240:   }

242:   PetscSectionGetChart(s1, &pStart, &pEnd);
243:   PetscSectionGetChart(s2, &n1, &n2);
244:   if (pStart != n1 || pEnd != n2) goto not_congruent;

246:   PetscSectionGetPermutation(s1, &perm1);
247:   PetscSectionGetPermutation(s2, &perm2);
248:   if (perm1 && perm2) {
249:     ISEqual(perm1, perm2, congruent);
250:     if (!(*congruent)) goto not_congruent;
251:   } else if (perm1 != perm2) goto not_congruent;

253:   for (p = pStart; p < pEnd; ++p) {
254:     PetscSectionGetOffset(s1, p, &n1);
255:     PetscSectionGetOffset(s2, p, &n2);
256:     if (n1 != n2) goto not_congruent;

258:     PetscSectionGetDof(s1, p, &n1);
259:     PetscSectionGetDof(s2, p, &n2);
260:     if (n1 != n2) goto not_congruent;

262:     PetscSectionGetConstraintDof(s1, p, &ncdof);
263:     PetscSectionGetConstraintDof(s2, p, &n2);
264:     if (ncdof != n2) goto not_congruent;

266:     PetscSectionGetConstraintIndices(s1, p, &idx1);
267:     PetscSectionGetConstraintIndices(s2, p, &idx2);
268:     PetscArraycmp(idx1, idx2, ncdof, congruent);
269:     if (!(*congruent)) goto not_congruent;
270:   }

272:   PetscSectionGetNumFields(s1, &nfields);
273:   PetscSectionGetNumFields(s2, &n2);
274:   if (nfields != n2) goto not_congruent;

276:   for (f = 0; f < nfields; ++f) {
277:     PetscSectionGetFieldComponents(s1, f, &n1);
278:     PetscSectionGetFieldComponents(s2, f, &n2);
279:     if (n1 != n2) goto not_congruent;

281:     for (p = pStart; p < pEnd; ++p) {
282:       PetscSectionGetFieldOffset(s1, p, f, &n1);
283:       PetscSectionGetFieldOffset(s2, p, f, &n2);
284:       if (n1 != n2) goto not_congruent;

286:       PetscSectionGetFieldDof(s1, p, f, &n1);
287:       PetscSectionGetFieldDof(s2, p, f, &n2);
288:       if (n1 != n2) goto not_congruent;

290:       PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
291:       PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
292:       if (nfcdof != n2) goto not_congruent;

294:       PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
295:       PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
296:       PetscArraycmp(idx1, idx2, nfcdof, congruent);
297:       if (!(*congruent)) goto not_congruent;
298:     }
299:   }

301:   flg = PETSC_TRUE;
302: not_congruent:
303:   MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1));
304:   return 0;
305: }

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

310:   Not Collective

312:   Input Parameter:
313: . s - the PetscSection

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

318:   Level: intermediate

320: .seealso: `PetscSectionSetNumFields()`
321: @*/
322: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
323: {
326:   *numFields = s->numFields;
327:   return 0;
328: }

330: /*@
331:   PetscSectionSetNumFields - Sets the number of fields.

333:   Not Collective

335:   Input Parameters:
336: + s - the PetscSection
337: - numFields - the number of fields

339:   Level: intermediate

341: .seealso: `PetscSectionGetNumFields()`
342: @*/
343: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
344: {
345:   PetscInt f;

349:   PetscSectionReset(s);

351:   s->numFields = numFields;
352:   PetscMalloc1(s->numFields, &s->numFieldComponents);
353:   PetscMalloc1(s->numFields, &s->fieldNames);
354:   PetscMalloc1(s->numFields, &s->compNames);
355:   PetscMalloc1(s->numFields, &s->field);
356:   for (f = 0; f < s->numFields; ++f) {
357:     char name[64];

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

361:     PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]);
362:     PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f);
363:     PetscStrallocpy(name, (char **)&s->fieldNames[f]);
364:     PetscSNPrintf(name, 64, "Component_0");
365:     PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
366:     PetscStrallocpy(name, (char **)&s->compNames[f][0]);
367:   }
368:   return 0;
369: }

371: /*@C
372:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

374:   Not Collective

376:   Input Parameters:
377: + s     - the PetscSection
378: - field - the field number

380:   Output Parameter:
381: . fieldName - the field name

383:   Level: intermediate

385: .seealso: `PetscSectionSetFieldName()`
386: @*/
387: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
388: {
391:   PetscSectionCheckValidField(field, s->numFields);
392:   *fieldName = s->fieldNames[field];
393:   return 0;
394: }

396: /*@C
397:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

399:   Not Collective

401:   Input Parameters:
402: + s     - the PetscSection
403: . field - the field number
404: - fieldName - the field name

406:   Level: intermediate

408: .seealso: `PetscSectionGetFieldName()`
409: @*/
410: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
411: {
414:   PetscSectionCheckValidField(field, s->numFields);
415:   PetscFree(s->fieldNames[field]);
416:   PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]);
417:   return 0;
418: }

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

423:   Not Collective

425:   Input Parameters:
426: + s     - the PetscSection
427: . field - the field number
428: . comp  - the component number
429: - compName - the component name

431:   Level: intermediate

433: .seealso: `PetscSectionSetComponentName()`
434: @*/
435: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
436: {
439:   PetscSectionCheckValidField(field, s->numFields);
440:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
441:   *compName = s->compNames[field][comp];
442:   return 0;
443: }

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

448:   Not Collective

450:   Input Parameters:
451: + s     - the PetscSection
452: . field - the field number
453: . comp  - the component number
454: - compName - the component name

456:   Level: intermediate

458: .seealso: `PetscSectionGetComponentName()`
459: @*/
460: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
461: {
464:   PetscSectionCheckValidField(field, s->numFields);
465:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
466:   PetscFree(s->compNames[field][comp]);
467:   PetscStrallocpy(compName, (char **)&s->compNames[field][comp]);
468:   return 0;
469: }

471: /*@
472:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

474:   Not Collective

476:   Input Parameters:
477: + s - the PetscSection
478: - field - the field number

480:   Output Parameter:
481: . numComp - the number of field components

483:   Level: intermediate

485: .seealso: `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`
486: @*/
487: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
488: {
491:   PetscSectionCheckValidField(field, s->numFields);
492:   *numComp = s->numFieldComponents[field];
493:   return 0;
494: }

496: /*@
497:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

499:   Not Collective

501:   Input Parameters:
502: + s - the PetscSection
503: . field - the field number
504: - numComp - the number of field components

506:   Level: intermediate

508: .seealso: `PetscSectionGetFieldComponents()`, `PetscSectionGetNumFields()`
509: @*/
510: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
511: {
512:   PetscInt c;

515:   PetscSectionCheckValidField(field, s->numFields);
516:   if (s->compNames) {
517:     for (c = 0; c < s->numFieldComponents[field]; ++c) PetscFree(s->compNames[field][c]);
518:     PetscFree(s->compNames[field]);
519:   }

521:   s->numFieldComponents[field] = numComp;
522:   if (numComp) {
523:     PetscMalloc1(numComp, (char ***)&s->compNames[field]);
524:     for (c = 0; c < numComp; ++c) {
525:       char name[64];

527:       PetscSNPrintf(name, 64, "%" PetscInt_FMT, c);
528:       PetscStrallocpy(name, (char **)&s->compNames[field][c]);
529:     }
530:   }
531:   return 0;
532: }

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

537:   Not Collective

539:   Input Parameter:
540: . s - the PetscSection

542:   Output Parameters:
543: + pStart - the first point
544: - pEnd - one past the last point

546:   Level: intermediate

548: .seealso: `PetscSectionSetChart()`, `PetscSectionCreate()`
549: @*/
550: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
551: {
553:   if (pStart) *pStart = s->pStart;
554:   if (pEnd) *pEnd = s->pEnd;
555:   return 0;
556: }

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

561:   Not Collective

563:   Input Parameters:
564: + s - the PetscSection
565: . pStart - the first point
566: - pEnd - one past the last point

568:   Level: intermediate

570: .seealso: `PetscSectionGetChart()`, `PetscSectionCreate()`
571: @*/
572: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
573: {
574:   PetscInt f;

577:   if (pStart == s->pStart && pEnd == s->pEnd) return 0;
578:   /* Cannot Reset() because it destroys field information */
579:   s->setup = PETSC_FALSE;
580:   PetscSectionDestroy(&s->bc);
581:   PetscFree(s->bcIndices);
582:   PetscFree2(s->atlasDof, s->atlasOff);

584:   s->pStart = pStart;
585:   s->pEnd   = pEnd;
586:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
587:   PetscArrayzero(s->atlasDof, pEnd - pStart);
588:   for (f = 0; f < s->numFields; ++f) PetscSectionSetChart(s->field[f], pStart, pEnd);
589:   return 0;
590: }

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

595:   Not Collective

597:   Input Parameter:
598: . s - the PetscSection

600:   Output Parameters:
601: . perm - The permutation as an IS

603:   Level: intermediate

605: .seealso: `PetscSectionSetPermutation()`, `PetscSectionCreate()`
606: @*/
607: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
608: {
610:   if (perm) {
612:     *perm = s->perm;
613:   }
614:   return 0;
615: }

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

620:   Not Collective

622:   Input Parameters:
623: + s - the PetscSection
624: - perm - the permutation of points

626:   Level: intermediate

628: .seealso: `PetscSectionGetPermutation()`, `PetscSectionCreate()`
629: @*/
630: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
631: {
635:   if (s->perm != perm) {
636:     ISDestroy(&s->perm);
637:     if (perm) {
638:       s->perm = perm;
639:       PetscObjectReference((PetscObject)s->perm);
640:     }
641:   }
642:   return 0;
643: }

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

648:   Not Collective

650:   Input Parameter:
651: . s - the PetscSection

653:   Output Parameter:
654: . pm - the flag for point major ordering

656:   Level: intermediate

658: .seealso: `PetscSectionSetPointMajor()`
659: @*/
660: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
661: {
664:   *pm = s->pointMajor;
665:   return 0;
666: }

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

671:   Not Collective

673:   Input Parameters:
674: + s  - the PetscSection
675: - pm - the flag for point major ordering

677:   Level: intermediate

679: .seealso: `PetscSectionGetPointMajor()`
680: @*/
681: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
682: {
685:   s->pointMajor = pm;
686:   return 0;
687: }

689: /*@
690:   PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets

692:   Not Collective

694:   Input Parameter:
695: . s - the PetscSection

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

700:   Level: intermediate

702: .seealso: `PetscSectionSetIncludesConstraints()`
703: @*/
704: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
705: {
708:   *includesConstraints = s->includesConstraints;
709:   return 0;
710: }

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

715:   Not Collective

717:   Input Parameters:
718: + s  - the PetscSection
719: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets

721:   Level: intermediate

723: .seealso: `PetscSectionGetIncludesConstraints()`
724: @*/
725: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
726: {
729:   s->includesConstraints = includesConstraints;
730:   return 0;
731: }

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

736:   Not Collective

738:   Input Parameters:
739: + s - the PetscSection
740: - point - the point

742:   Output Parameter:
743: . numDof - the number of dof

745:   Level: intermediate

747: .seealso: `PetscSectionSetDof()`, `PetscSectionCreate()`
748: @*/
749: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
750: {
754:   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);
755:   *numDof = s->atlasDof[point - s->pStart];
756:   return 0;
757: }

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

762:   Not Collective

764:   Input Parameters:
765: + s - the PetscSection
766: . point - the point
767: - numDof - the number of dof

769:   Level: intermediate

771: .seealso: `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
772: @*/
773: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
774: {
776:   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);
777:   s->atlasDof[point - s->pStart] = numDof;
778:   PetscSectionInvalidateMaxDof_Internal(s);
779:   return 0;
780: }

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

785:   Not Collective

787:   Input Parameters:
788: + s - the PetscSection
789: . point - the point
790: - numDof - the number of additional dof

792:   Level: intermediate

794: .seealso: `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
795: @*/
796: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
797: {
800:   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);
801:   s->atlasDof[point - s->pStart] += numDof;
802:   PetscSectionInvalidateMaxDof_Internal(s);
803:   return 0;
804: }

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

809:   Not Collective

811:   Input Parameters:
812: + s - the PetscSection
813: . point - the point
814: - field - the field

816:   Output Parameter:
817: . numDof - the number of dof

819:   Level: intermediate

821: .seealso: `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
822: @*/
823: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
824: {
827:   PetscSectionCheckValidField(field, s->numFields);
828:   PetscSectionGetDof(s->field[field], point, numDof);
829:   return 0;
830: }

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

835:   Not Collective

837:   Input Parameters:
838: + s - the PetscSection
839: . point - the point
840: . field - the field
841: - numDof - the number of dof

843:   Level: intermediate

845: .seealso: `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
846: @*/
847: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
848: {
850:   PetscSectionCheckValidField(field, s->numFields);
851:   PetscSectionSetDof(s->field[field], point, numDof);
852:   return 0;
853: }

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

858:   Not Collective

860:   Input Parameters:
861: + s - the PetscSection
862: . point - the point
863: . field - the field
864: - numDof - the number of dof

866:   Level: intermediate

868: .seealso: `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
869: @*/
870: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
871: {
873:   PetscSectionCheckValidField(field, s->numFields);
874:   PetscSectionAddDof(s->field[field], point, numDof);
875:   return 0;
876: }

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

881:   Not Collective

883:   Input Parameters:
884: + s - the PetscSection
885: - point - the point

887:   Output Parameter:
888: . numDof - the number of dof which are fixed by constraints

890:   Level: intermediate

892: .seealso: `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
893: @*/
894: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
895: {
898:   if (s->bc) {
899:     PetscSectionGetDof(s->bc, point, numDof);
900:   } else *numDof = 0;
901:   return 0;
902: }

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

907:   Not Collective

909:   Input Parameters:
910: + s - the PetscSection
911: . point - the point
912: - numDof - the number of dof which are fixed by constraints

914:   Level: intermediate

916: .seealso: `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
917: @*/
918: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
919: {
921:   if (numDof) {
922:     PetscSectionCheckConstraints_Private(s);
923:     PetscSectionSetDof(s->bc, point, numDof);
924:   }
925:   return 0;
926: }

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

931:   Not Collective

933:   Input Parameters:
934: + s - the PetscSection
935: . point - the point
936: - numDof - the number of additional dof which are fixed by constraints

938:   Level: intermediate

940: .seealso: `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
941: @*/
942: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
943: {
945:   if (numDof) {
946:     PetscSectionCheckConstraints_Private(s);
947:     PetscSectionAddDof(s->bc, point, numDof);
948:   }
949:   return 0;
950: }

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

955:   Not Collective

957:   Input Parameters:
958: + s - the PetscSection
959: . point - the point
960: - field - the field

962:   Output Parameter:
963: . numDof - the number of dof which are fixed by constraints

965:   Level: intermediate

967: .seealso: `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
968: @*/
969: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
970: {
973:   PetscSectionCheckValidField(field, s->numFields);
974:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
975:   return 0;
976: }

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

981:   Not Collective

983:   Input Parameters:
984: + s - the PetscSection
985: . point - the point
986: . field - the field
987: - numDof - the number of dof which are fixed by constraints

989:   Level: intermediate

991: .seealso: `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
992: @*/
993: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
994: {
996:   PetscSectionCheckValidField(field, s->numFields);
997:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
998:   return 0;
999: }

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

1004:   Not Collective

1006:   Input Parameters:
1007: + s - the PetscSection
1008: . point - the point
1009: . field - the field
1010: - numDof - the number of additional dof which are fixed by constraints

1012:   Level: intermediate

1014: .seealso: `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1015: @*/
1016: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1017: {
1019:   PetscSectionCheckValidField(field, s->numFields);
1020:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
1021:   return 0;
1022: }

1024: /*@
1025:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1027:   Not Collective

1029:   Input Parameter:
1030: . s - the PetscSection

1032:   Level: advanced

1034: .seealso: `PetscSectionSetUp()`, `PetscSectionCreate()`
1035: @*/
1036: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1037: {
1039:   if (s->bc) {
1040:     const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;

1042:     PetscSectionSetUp(s->bc);
1043:     PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices);
1044:   }
1045:   return 0;
1046: }

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

1051:   Not Collective

1053:   Input Parameter:
1054: . s - the PetscSection

1056:   Level: intermediate

1058: .seealso: `PetscSectionCreate()`
1059: @*/
1060: PetscErrorCode PetscSectionSetUp(PetscSection s)
1061: {
1062:   const PetscInt *pind   = NULL;
1063:   PetscInt        offset = 0, foff, p, f;

1066:   if (s->setup) return 0;
1067:   s->setup = PETSC_TRUE;
1068:   /* Set offsets and field offsets for all points */
1069:   /*   Assume that all fields have the same chart */
1071:   if (s->perm) ISGetIndices(s->perm, &pind);
1072:   if (s->pointMajor) {
1073:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1074:       const PetscInt q = pind ? pind[p] : p;

1076:       /* Set point offset */
1077:       s->atlasOff[q] = offset;
1078:       offset += s->atlasDof[q];
1079:       /* Set field offset */
1080:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1081:         PetscSection sf = s->field[f];

1083:         sf->atlasOff[q] = foff;
1084:         foff += sf->atlasDof[q];
1085:       }
1086:     }
1087:   } else {
1088:     /* Set field offsets for all points */
1089:     for (f = 0; f < s->numFields; ++f) {
1090:       PetscSection sf = s->field[f];

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

1095:         sf->atlasOff[q] = offset;
1096:         offset += sf->atlasDof[q];
1097:       }
1098:     }
1099:     /* Disable point offsets since these are unused */
1100:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1101:   }
1102:   if (s->perm) ISRestoreIndices(s->perm, &pind);
1103:   /* Setup BC sections */
1104:   PetscSectionSetUpBC(s);
1105:   for (f = 0; f < s->numFields; ++f) PetscSectionSetUpBC(s->field[f]);
1106:   return 0;
1107: }

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

1112:   Not Collective

1114:   Input Parameters:
1115: . s - the PetscSection

1117:   Output Parameter:
1118: . maxDof - the maximum dof

1120:   Level: intermediate

1122:   Notes:
1123:   The returned number is up-to-date without need for PetscSectionSetUp().

1125:   Developer Notes:
1126:   The returned number is calculated lazily and stashed.
1127:   A call to PetscSectionInvalidateMaxDof_Internal() invalidates the stashed value.
1128:   PetscSectionInvalidateMaxDof_Internal() is called in PetscSectionSetDof(), PetscSectionAddDof() and PetscSectionReset().
1129:   It should also be called every time atlasDof is modified directly.

1131: .seealso: `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1132: @*/
1133: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1134: {
1135:   PetscInt p;

1139:   if (s->maxDof == PETSC_MIN_INT) {
1140:     s->maxDof = 0;
1141:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1142:   }
1143:   *maxDof = s->maxDof;
1144:   return 0;
1145: }

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

1150:   Not Collective

1152:   Input Parameter:
1153: . s - the PetscSection

1155:   Output Parameter:
1156: . size - the size of an array which can hold all the dofs

1158:   Level: intermediate

1160: .seealso: `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1161: @*/
1162: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1163: {
1164:   PetscInt p, n = 0;

1168:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1169:   *size = n;
1170:   return 0;
1171: }

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

1176:   Not Collective

1178:   Input Parameter:
1179: . s - the PetscSection

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

1184:   Level: intermediate

1186: .seealso: `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1187: @*/
1188: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1189: {
1190:   PetscInt p, n = 0;

1194:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1195:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1196:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1197:   }
1198:   *size = n;
1199:   return 0;
1200: }

1202: /*@
1203:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1204:   the local section and an SF describing the section point overlap.

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

1212:   Output Parameter:
1213: . gsection - The PetscSection for the global field layout

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

1217:   Level: intermediate

1219: .seealso: `PetscSectionCreate()`
1220: @*/
1221: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1222: {
1223:   PetscSection    gs;
1224:   const PetscInt *pind = NULL;
1225:   PetscInt       *recv = NULL, *neg = NULL;
1226:   PetscInt        pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1227:   PetscInt        numFields, f, numComponents;

1235:   PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs);
1236:   PetscSectionGetNumFields(s, &numFields);
1237:   if (numFields > 0) PetscSectionSetNumFields(gs, numFields);
1238:   PetscSectionGetChart(s, &pStart, &pEnd);
1239:   PetscSectionSetChart(gs, pStart, pEnd);
1240:   gs->includesConstraints = includeConstraints;
1241:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1242:   nlocal = nroots; /* The local/leaf space matches global/root space */
1243:   /* Must allocate for all points visible to SF, which may be more than this section */
1244:   if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1245:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
1248:     PetscMalloc2(nroots, &neg, nlocal, &recv);
1249:     PetscArrayzero(neg, nroots);
1250:   }
1251:   /* Mark all local points with negative dof */
1252:   for (p = pStart; p < pEnd; ++p) {
1253:     PetscSectionGetDof(s, p, &dof);
1254:     PetscSectionSetDof(gs, p, dof);
1255:     PetscSectionGetConstraintDof(s, p, &cdof);
1256:     if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(gs, p, cdof);
1257:     if (neg) neg[p] = -(dof + 1);
1258:   }
1259:   PetscSectionSetUpBC(gs);
1260:   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]);
1261:   if (nroots >= 0) {
1262:     PetscArrayzero(recv, nlocal);
1263:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1264:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1265:     for (p = pStart; p < pEnd; ++p) {
1266:       if (recv[p] < 0) {
1267:         gs->atlasDof[p - pStart] = recv[p];
1268:         PetscSectionGetDof(s, p, &dof);
1270:       }
1271:     }
1272:   }
1273:   /* Calculate new sizes, get process offset, and calculate point offsets */
1274:   if (s->perm) ISGetIndices(s->perm, &pind);
1275:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1276:     const PetscInt q = pind ? pind[p] : p;

1278:     cdof            = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1279:     gs->atlasOff[q] = off;
1280:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1281:   }
1282:   if (!localOffsets) {
1283:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1284:     globalOff -= off;
1285:   }
1286:   for (p = pStart, off = 0; p < pEnd; ++p) {
1287:     gs->atlasOff[p - pStart] += globalOff;
1288:     if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1289:   }
1290:   if (s->perm) ISRestoreIndices(s->perm, &pind);
1291:   /* Put in negative offsets for ghost points */
1292:   if (nroots >= 0) {
1293:     PetscArrayzero(recv, nlocal);
1294:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1295:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1296:     for (p = pStart; p < pEnd; ++p) {
1297:       if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1298:     }
1299:   }
1300:   PetscFree2(neg, recv);
1301:   /* Set field dofs/offsets/constraints */
1302:   for (f = 0; f < numFields; ++f) {
1303:     gs->field[f]->includesConstraints = includeConstraints;
1304:     PetscSectionGetFieldComponents(s, f, &numComponents);
1305:     PetscSectionSetFieldComponents(gs, f, numComponents);
1306:   }
1307:   for (p = pStart; p < pEnd; ++p) {
1308:     PetscSectionGetOffset(gs, p, &off);
1309:     for (f = 0, foff = off; f < numFields; ++f) {
1310:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1311:       if (!includeConstraints && cdof > 0) PetscSectionSetFieldConstraintDof(gs, p, f, cdof);
1312:       PetscSectionGetFieldDof(s, p, f, &dof);
1313:       PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof);
1314:       PetscSectionSetFieldOffset(gs, p, f, foff);
1315:       PetscSectionGetFieldConstraintDof(gs, p, f, &cdof);
1316:       foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1317:     }
1318:   }
1319:   for (f = 0; f < numFields; ++f) {
1320:     PetscSection gfs = gs->field[f];

1322:     PetscSectionSetUpBC(gfs);
1323:     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]);
1324:   }
1325:   gs->setup = PETSC_TRUE;
1326:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1327:   *gsection = gs;
1328:   return 0;
1329: }

1331: /*@
1332:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1333:   the local section and an SF describing the section point overlap.

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

1342:   Output Parameter:
1343: . gsection - The PetscSection for the global field layout

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

1347:   Level: advanced

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

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

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

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

1433:   Collective on comm

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

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

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

1444:   Level: advanced

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

1452:   PetscSectionGetChart(s, &pStart, &pEnd);
1453:   for (p = pStart; p < pEnd; ++p) {
1454:     PetscInt dof;

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

1466: /*@
1467:   PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.

1469:   Collective on comm

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

1475:   Output Parameter:
1476: . layout - The dof layout for the section

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

1480:   Level: advanced

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

1490:   PetscSectionGetChart(s, &pStart, &pEnd);
1491:   for (p = pStart; p < pEnd; ++p) {
1492:     PetscInt dof, cdof;

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

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

1508:   Not Collective

1510:   Input Parameters:
1511: + s - the PetscSection
1512: - point - the point

1514:   Output Parameter:
1515: . offset - the offset

1517:   Level: intermediate

1519: .seealso: `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`
1520: @*/
1521: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1522: {
1526:   *offset = s->atlasOff[point - s->pStart];
1527:   return 0;
1528: }

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

1533:   Not Collective

1535:   Input Parameters:
1536: + s - the PetscSection
1537: . point - the point
1538: - offset - the offset

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

1542:   Level: intermediate

1544: .seealso: `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1545: @*/
1546: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1547: {
1550:   s->atlasOff[point - s->pStart] = offset;
1551:   return 0;
1552: }

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

1557:   Not Collective

1559:   Input Parameters:
1560: + s - the PetscSection
1561: . point - the point
1562: - field - the field

1564:   Output Parameter:
1565: . offset - the offset

1567:   Level: intermediate

1569: .seealso: `PetscSectionGetOffset()`, `PetscSectionCreate()`
1570: @*/
1571: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1572: {
1575:   PetscSectionCheckValidField(field, s->numFields);
1576:   PetscSectionGetOffset(s->field[field], point, offset);
1577:   return 0;
1578: }

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

1583:   Not Collective

1585:   Input Parameters:
1586: + s - the PetscSection
1587: . point - the point
1588: . field - the field
1589: - offset - the offset

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

1593:   Level: intermediate

1595: .seealso: `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1596: @*/
1597: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1598: {
1600:   PetscSectionCheckValidField(field, s->numFields);
1601:   PetscSectionSetOffset(s->field[field], point, offset);
1602:   return 0;
1603: }

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

1608:   Not Collective

1610:   Input Parameters:
1611: + s - the PetscSection
1612: . point - the point
1613: - field - the field

1615:   Output Parameter:
1616: . offset - the offset

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

1621:   Level: advanced

1623: .seealso: `PetscSectionGetOffset()`, `PetscSectionCreate()`
1624: @*/
1625: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1626: {
1627:   PetscInt off, foff;

1631:   PetscSectionCheckValidField(field, s->numFields);
1632:   PetscSectionGetOffset(s, point, &off);
1633:   PetscSectionGetOffset(s->field[field], point, &foff);
1634:   *offset = foff - off;
1635:   return 0;
1636: }

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

1641:   Not Collective

1643:   Input Parameter:
1644: . s - the PetscSection

1646:   Output Parameters:
1647: + start - the minimum offset
1648: - end   - one more than the maximum offset

1650:   Level: intermediate

1652: .seealso: `PetscSectionGetOffset()`, `PetscSectionCreate()`
1653: @*/
1654: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1655: {
1656:   PetscInt os = 0, oe = 0, pStart, pEnd, p;

1659:   if (s->atlasOff) {
1660:     os = s->atlasOff[0];
1661:     oe = s->atlasOff[0];
1662:   }
1663:   PetscSectionGetChart(s, &pStart, &pEnd);
1664:   for (p = 0; p < pEnd - pStart; ++p) {
1665:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1667:     if (off >= 0) {
1668:       os = PetscMin(os, off);
1669:       oe = PetscMax(oe, off + dof);
1670:     }
1671:   }
1672:   if (start) *start = os;
1673:   if (end) *end = oe;
1674:   return 0;
1675: }

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

1680:   Collective

1682:   Input Parameters:
1683: + s      - the PetscSection
1684: . len    - the number of subfields
1685: - fields - the subfield numbers

1687:   Output Parameter:
1688: . subs   - the subsection

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

1692:   Level: advanced

1694: .seealso: `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1695: @*/
1696: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1697: {
1698:   PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;

1700:   if (!len) return 0;
1704:   PetscSectionGetNumFields(s, &nF);
1706:   PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1707:   PetscSectionSetNumFields(*subs, len);
1708:   for (f = 0; f < len; ++f) {
1709:     const char *name    = NULL;
1710:     PetscInt    numComp = 0;

1712:     PetscSectionGetFieldName(s, fields[f], &name);
1713:     PetscSectionSetFieldName(*subs, f, name);
1714:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
1715:     PetscSectionSetFieldComponents(*subs, f, numComp);
1716:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1717:       PetscSectionGetComponentName(s, fields[f], c, &name);
1718:       PetscSectionSetComponentName(*subs, f, c, name);
1719:     }
1720:   }
1721:   PetscSectionGetChart(s, &pStart, &pEnd);
1722:   PetscSectionSetChart(*subs, pStart, pEnd);
1723:   for (p = pStart; p < pEnd; ++p) {
1724:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1726:     for (f = 0; f < len; ++f) {
1727:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1728:       PetscSectionSetFieldDof(*subs, p, f, fdof);
1729:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1730:       if (cfdof) PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);
1731:       dof += fdof;
1732:       cdof += cfdof;
1733:     }
1734:     PetscSectionSetDof(*subs, p, dof);
1735:     if (cdof) PetscSectionSetConstraintDof(*subs, p, cdof);
1736:     maxCdof = PetscMax(cdof, maxCdof);
1737:   }
1738:   PetscSectionSetUp(*subs);
1739:   if (maxCdof) {
1740:     PetscInt *indices;

1742:     PetscMalloc1(maxCdof, &indices);
1743:     for (p = pStart; p < pEnd; ++p) {
1744:       PetscInt cdof;

1746:       PetscSectionGetConstraintDof(*subs, p, &cdof);
1747:       if (cdof) {
1748:         const PetscInt *oldIndices = NULL;
1749:         PetscInt        fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1751:         for (f = 0; f < len; ++f) {
1752:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1753:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1754:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1755:           PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices);
1756:           for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1757:           numConst += cfdof;
1758:           fOff += fdof;
1759:         }
1760:         PetscSectionSetConstraintIndices(*subs, p, indices);
1761:       }
1762:     }
1763:     PetscFree(indices);
1764:   }
1765:   return 0;
1766: }

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

1771:   Collective

1773:   Input Parameters:
1774: + s     - the input sections
1775: - len   - the number of input sections

1777:   Output Parameter:
1778: . supers - the supersection

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

1782:   Level: advanced

1784: .seealso: `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
1785: @*/
1786: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1787: {
1788:   PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

1790:   if (!len) return 0;
1791:   for (i = 0; i < len; ++i) {
1792:     PetscInt nf, pStarti, pEndi;

1794:     PetscSectionGetNumFields(s[i], &nf);
1795:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
1796:     pStart = PetscMin(pStart, pStarti);
1797:     pEnd   = PetscMax(pEnd, pEndi);
1798:     Nf += nf;
1799:   }
1800:   PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers);
1801:   PetscSectionSetNumFields(*supers, Nf);
1802:   for (i = 0, f = 0; i < len; ++i) {
1803:     PetscInt nf, fi, ci;

1805:     PetscSectionGetNumFields(s[i], &nf);
1806:     for (fi = 0; fi < nf; ++fi, ++f) {
1807:       const char *name    = NULL;
1808:       PetscInt    numComp = 0;

1810:       PetscSectionGetFieldName(s[i], fi, &name);
1811:       PetscSectionSetFieldName(*supers, f, name);
1812:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
1813:       PetscSectionSetFieldComponents(*supers, f, numComp);
1814:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1815:         PetscSectionGetComponentName(s[i], fi, ci, &name);
1816:         PetscSectionSetComponentName(*supers, f, ci, name);
1817:       }
1818:     }
1819:   }
1820:   PetscSectionSetChart(*supers, pStart, pEnd);
1821:   for (p = pStart; p < pEnd; ++p) {
1822:     PetscInt dof = 0, cdof = 0;

1824:     for (i = 0, f = 0; i < len; ++i) {
1825:       PetscInt nf, fi, pStarti, pEndi;
1826:       PetscInt fdof = 0, cfdof = 0;

1828:       PetscSectionGetNumFields(s[i], &nf);
1829:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
1830:       if ((p < pStarti) || (p >= pEndi)) continue;
1831:       for (fi = 0; fi < nf; ++fi, ++f) {
1832:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1833:         PetscSectionAddFieldDof(*supers, p, f, fdof);
1834:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1835:         if (cfdof) PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);
1836:         dof += fdof;
1837:         cdof += cfdof;
1838:       }
1839:     }
1840:     PetscSectionSetDof(*supers, p, dof);
1841:     if (cdof) PetscSectionSetConstraintDof(*supers, p, cdof);
1842:     maxCdof = PetscMax(cdof, maxCdof);
1843:   }
1844:   PetscSectionSetUp(*supers);
1845:   if (maxCdof) {
1846:     PetscInt *indices;

1848:     PetscMalloc1(maxCdof, &indices);
1849:     for (p = pStart; p < pEnd; ++p) {
1850:       PetscInt cdof;

1852:       PetscSectionGetConstraintDof(*supers, p, &cdof);
1853:       if (cdof) {
1854:         PetscInt dof, numConst = 0, fOff = 0;

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

1860:           PetscSectionGetNumFields(s[i], &nf);
1861:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
1862:           if ((p < pStarti) || (p >= pEndi)) continue;
1863:           for (fi = 0; fi < nf; ++fi, ++f) {
1864:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1865:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1866:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1867:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
1868:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1869:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
1870:             numConst += cfdof;
1871:           }
1872:           PetscSectionGetDof(s[i], p, &dof);
1873:           fOff += dof;
1874:         }
1875:         PetscSectionSetConstraintIndices(*supers, p, indices);
1876:       }
1877:     }
1878:     PetscFree(indices);
1879:   }
1880:   return 0;
1881: }

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

1891:   PetscSectionGetNumFields(s, &numFields);
1892:   PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1893:   if (numFields) PetscSectionSetNumFields(*subs, numFields);
1894:   for (f = 0; f < numFields; ++f) {
1895:     const char *name    = NULL;
1896:     PetscInt    numComp = 0;

1898:     PetscSectionGetFieldName(s, f, &name);
1899:     PetscSectionSetFieldName(*subs, f, name);
1900:     PetscSectionGetFieldComponents(s, f, &numComp);
1901:     PetscSectionSetFieldComponents(*subs, f, numComp);
1902:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
1903:       PetscSectionGetComponentName(s, f, c, &name);
1904:       PetscSectionSetComponentName(*subs, f, c, name);
1905:     }
1906:   }
1907:   /* For right now, we do not try to squeeze the subchart */
1908:   if (subpointMap) {
1909:     ISGetSize(subpointMap, &numSubpoints);
1910:     ISGetIndices(subpointMap, &points);
1911:   }
1912:   PetscSectionGetChart(s, &pStart, &pEnd);
1913:   if (renumberPoints) {
1914:     spStart = 0;
1915:     spEnd   = numSubpoints;
1916:   } else {
1917:     ISGetMinMax(subpointMap, &spStart, &spEnd);
1918:     ++spEnd;
1919:   }
1920:   PetscSectionSetChart(*subs, spStart, spEnd);
1921:   for (p = pStart; p < pEnd; ++p) {
1922:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

1924:     PetscFindInt(p, numSubpoints, points, &subp);
1925:     if (subp < 0) continue;
1926:     if (!renumberPoints) subp = p;
1927:     for (f = 0; f < numFields; ++f) {
1928:       PetscSectionGetFieldDof(s, p, f, &fdof);
1929:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
1930:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1931:       if (cfdof) PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);
1932:     }
1933:     PetscSectionGetDof(s, p, &dof);
1934:     PetscSectionSetDof(*subs, subp, dof);
1935:     PetscSectionGetConstraintDof(s, p, &cdof);
1936:     if (cdof) PetscSectionSetConstraintDof(*subs, subp, cdof);
1937:   }
1938:   PetscSectionSetUp(*subs);
1939:   /* Change offsets to original offsets */
1940:   for (p = pStart; p < pEnd; ++p) {
1941:     PetscInt off, foff = 0;

1943:     PetscFindInt(p, numSubpoints, points, &subp);
1944:     if (subp < 0) continue;
1945:     if (!renumberPoints) subp = p;
1946:     for (f = 0; f < numFields; ++f) {
1947:       PetscSectionGetFieldOffset(s, p, f, &foff);
1948:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
1949:     }
1950:     PetscSectionGetOffset(s, p, &off);
1951:     PetscSectionSetOffset(*subs, subp, off);
1952:   }
1953:   /* Copy constraint indices */
1954:   for (subp = spStart; subp < spEnd; ++subp) {
1955:     PetscInt cdof;

1957:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
1958:     if (cdof) {
1959:       for (f = 0; f < numFields; ++f) {
1960:         PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices);
1961:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1962:       }
1963:       PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices);
1964:       PetscSectionSetConstraintIndices(*subs, subp, indices);
1965:     }
1966:   }
1967:   if (subpointMap) ISRestoreIndices(subpointMap, &points);
1968:   return 0;
1969: }

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

1974:   Collective on s

1976:   Input Parameters:
1977: + s           - the PetscSection
1978: - subpointMap - a sorted list of points in the original mesh which are in the submesh

1980:   Output Parameter:
1981: . subs - the subsection

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

1986:   Level: advanced

1988: .seealso: `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
1989: @*/
1990: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1991: {
1992:   PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs);
1993:   return 0;
1994: }

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

1999:   Collective on s

2001:   Input Parameters:
2002: + s           - the PetscSection
2003: - subpointMap - a sorted list of points in the original mesh which are in the subdomain

2005:   Output Parameter:
2006: . subs - the subsection

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

2011:   Level: advanced

2013: .seealso: `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2014: @*/
2015: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2016: {
2017:   PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs);
2018:   return 0;
2019: }

2021: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2022: {
2023:   PetscInt    p;
2024:   PetscMPIInt rank;

2026:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2027:   PetscViewerASCIIPushSynchronized(viewer);
2028:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2029:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2030:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2031:       PetscInt b;

2033:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2034:       if (s->bcIndices) {
2035:         for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
2036:       }
2037:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2038:     } else {
2039:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2040:     }
2041:   }
2042:   PetscViewerFlush(viewer);
2043:   PetscViewerASCIIPopSynchronized(viewer);
2044:   if (s->sym) {
2045:     PetscViewerASCIIPushTab(viewer);
2046:     PetscSectionSymView(s->sym, viewer);
2047:     PetscViewerASCIIPopTab(viewer);
2048:   }
2049:   return 0;
2050: }

2052: /*@C
2053:    PetscSectionViewFromOptions - View from Options

2055:    Collective

2057:    Input Parameters:
2058: +  A - the PetscSection object to view
2059: .  obj - Optional object
2060: -  name - command line option

2062:    Level: intermediate
2063: .seealso: `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`
2064: @*/
2065: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2066: {
2068:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
2069:   return 0;
2070: }

2072: /*@C
2073:   PetscSectionView - Views a PetscSection

2075:   Collective

2077:   Input Parameters:
2078: + s - the PetscSection object to view
2079: - v - the viewer

2081:   Note:
2082:   PetscSectionView(), when viewer is of type PETSCVIEWERHDF5, only saves
2083:   distribution independent data, such as dofs, offsets, constraint dofs,
2084:   and constraint indices. Points that have negative dofs, for instance,
2085:   are not saved as they represent points owned by other processes.
2086:   Point numbering and rank assignment is currently not stored.
2087:   The saved section can be loaded with PetscSectionLoad().

2089:   Level: beginner

2091: .seealso `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`
2092: @*/
2093: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2094: {
2095:   PetscBool isascii, ishdf5;
2096:   PetscInt  f;

2099:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);
2101:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
2102:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2103:   if (isascii) {
2104:     PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer);
2105:     if (s->numFields) {
2106:       PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields);
2107:       for (f = 0; f < s->numFields; ++f) {
2108:         PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
2109:         PetscSectionView_ASCII(s->field[f], viewer);
2110:       }
2111:     } else {
2112:       PetscSectionView_ASCII(s, viewer);
2113:     }
2114:   } else if (ishdf5) {
2115: #if PetscDefined(HAVE_HDF5)
2116:     PetscSectionView_HDF5_Internal(s, viewer);
2117: #else
2118:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2119: #endif
2120:   }
2121:   return 0;
2122: }

2124: /*@C
2125:   PetscSectionLoad - Loads a PetscSection

2127:   Collective

2129:   Input Parameters:
2130: + s - the PetscSection object to load
2131: - v - the viewer

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

2141:   Level: beginner

2143: .seealso `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2144: @*/
2145: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2146: {
2147:   PetscBool ishdf5;

2151:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2152:   if (ishdf5) {
2153: #if PetscDefined(HAVE_HDF5)
2154:     PetscSectionLoad_HDF5_Internal(s, viewer);
2155:     return 0;
2156: #else
2157:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2158: #endif
2159:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2160: }

2162: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2163: {
2164:   PetscSectionClosurePermVal clVal;

2166:   if (!section->clHash) return 0;
2167:   kh_foreach_value(section->clHash, clVal, {
2168:     PetscFree(clVal.perm);
2169:     PetscFree(clVal.invPerm);
2170:   });
2171:   kh_destroy(ClPerm, section->clHash);
2172:   section->clHash = NULL;
2173:   return 0;
2174: }

2176: /*@
2177:   PetscSectionReset - Frees all section data.

2179:   Not Collective

2181:   Input Parameters:
2182: . s - the PetscSection

2184:   Level: beginner

2186: .seealso: `PetscSection`, `PetscSectionCreate()`
2187: @*/
2188: PetscErrorCode PetscSectionReset(PetscSection s)
2189: {
2190:   PetscInt f, c;

2193:   for (f = 0; f < s->numFields; ++f) {
2194:     PetscSectionDestroy(&s->field[f]);
2195:     PetscFree(s->fieldNames[f]);
2196:     for (c = 0; c < s->numFieldComponents[f]; ++c) PetscFree(s->compNames[f][c]);
2197:     PetscFree(s->compNames[f]);
2198:   }
2199:   PetscFree(s->numFieldComponents);
2200:   PetscFree(s->fieldNames);
2201:   PetscFree(s->compNames);
2202:   PetscFree(s->field);
2203:   PetscSectionDestroy(&s->bc);
2204:   PetscFree(s->bcIndices);
2205:   PetscFree2(s->atlasDof, s->atlasOff);
2206:   PetscSectionDestroy(&s->clSection);
2207:   ISDestroy(&s->clPoints);
2208:   ISDestroy(&s->perm);
2209:   PetscSectionResetClosurePermutation(s);
2210:   PetscSectionSymDestroy(&s->sym);
2211:   PetscSectionDestroy(&s->clSection);
2212:   ISDestroy(&s->clPoints);
2213:   PetscSectionInvalidateMaxDof_Internal(s);
2214:   s->pStart    = -1;
2215:   s->pEnd      = -1;
2216:   s->maxDof    = 0;
2217:   s->setup     = PETSC_FALSE;
2218:   s->numFields = 0;
2219:   s->clObj     = NULL;
2220:   return 0;
2221: }

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

2226:   Not Collective

2228:   Input Parameters:
2229: . s - the PetscSection

2231:   Level: beginner

2233: .seealso: `PetscSection`, `PetscSectionCreate()`
2234: @*/
2235: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2236: {
2237:   if (!*s) return 0;
2239:   if (--((PetscObject)(*s))->refct > 0) {
2240:     *s = NULL;
2241:     return 0;
2242:   }
2243:   PetscSectionReset(*s);
2244:   PetscHeaderDestroy(s);
2245:   return 0;
2246: }

2248: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2249: {
2250:   const PetscInt p = point - s->pStart;

2253:   *values = &baseArray[s->atlasOff[p]];
2254:   return 0;
2255: }

2257: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2258: {
2259:   PetscInt      *array;
2260:   const PetscInt p           = point - s->pStart;
2261:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2262:   PetscInt       cDim        = 0;

2265:   PetscSectionGetConstraintDof(s, p, &cDim);
2266:   array = &baseArray[s->atlasOff[p]];
2267:   if (!cDim) {
2268:     if (orientation >= 0) {
2269:       const PetscInt dim = s->atlasDof[p];
2270:       PetscInt       i;

2272:       if (mode == INSERT_VALUES) {
2273:         for (i = 0; i < dim; ++i) array[i] = values[i];
2274:       } else {
2275:         for (i = 0; i < dim; ++i) array[i] += values[i];
2276:       }
2277:     } else {
2278:       PetscInt offset = 0;
2279:       PetscInt j      = -1, field, i;

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

2284:         for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2285:         offset += dim;
2286:       }
2287:     }
2288:   } else {
2289:     if (orientation >= 0) {
2290:       const PetscInt  dim  = s->atlasDof[p];
2291:       PetscInt        cInd = 0, i;
2292:       const PetscInt *cDof;

2294:       PetscSectionGetConstraintIndices(s, point, &cDof);
2295:       if (mode == INSERT_VALUES) {
2296:         for (i = 0; i < dim; ++i) {
2297:           if ((cInd < cDim) && (i == cDof[cInd])) {
2298:             ++cInd;
2299:             continue;
2300:           }
2301:           array[i] = values[i];
2302:         }
2303:       } else {
2304:         for (i = 0; i < dim; ++i) {
2305:           if ((cInd < cDim) && (i == cDof[cInd])) {
2306:             ++cInd;
2307:             continue;
2308:           }
2309:           array[i] += values[i];
2310:         }
2311:       }
2312:     } else {
2313:       const PetscInt *cDof;
2314:       PetscInt        offset  = 0;
2315:       PetscInt        cOffset = 0;
2316:       PetscInt        j       = 0, field;

2318:       PetscSectionGetConstraintIndices(s, point, &cDof);
2319:       for (field = 0; field < s->numFields; ++field) {
2320:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2321:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2322:         const PetscInt sDim = dim - tDim;
2323:         PetscInt       cInd = 0, i, k;

2325:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2326:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2327:             ++cInd;
2328:             continue;
2329:           }
2330:           array[j] = values[k];
2331:         }
2332:         offset += dim;
2333:         cOffset += dim - tDim;
2334:       }
2335:     }
2336:   }
2337:   return 0;
2338: }

2340: /*@C
2341:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

2343:   Not Collective

2345:   Input Parameter:
2346: . s - The PetscSection

2348:   Output Parameter:
2349: . hasConstraints - flag indicating that the section has constrained dofs

2351:   Level: intermediate

2353: .seealso: `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2354: @*/
2355: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2356: {
2359:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2360:   return 0;
2361: }

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

2366:   Not Collective

2368:   Input Parameters:
2369: + s     - The PetscSection
2370: - point - The point

2372:   Output Parameter:
2373: . indices - The constrained dofs

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

2377:   Level: intermediate

2379: .seealso: `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2380: @*/
2381: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2382: {
2384:   if (s->bc) {
2385:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2386:   } else *indices = NULL;
2387:   return 0;
2388: }

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

2393:   Not Collective

2395:   Input Parameters:
2396: + s     - The PetscSection
2397: . point - The point
2398: - indices - The constrained dofs

2400:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

2402:   Level: intermediate

2404: .seealso: `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2405: @*/
2406: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2407: {
2409:   if (s->bc) {
2410:     const PetscInt dof  = s->atlasDof[point];
2411:     const PetscInt cdof = s->bc->atlasDof[point];
2412:     PetscInt       d;

2415:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2416:   }
2417:   return 0;
2418: }

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

2423:   Not Collective

2425:   Input Parameters:
2426: + s     - The PetscSection
2427: . field  - The field number
2428: - point - The point

2430:   Output Parameter:
2431: . indices - The constrained dofs sorted in ascending order

2433:   Notes:
2434:   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().

2436:   Fortran Note:
2437:   In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()

2439:   Level: intermediate

2441: .seealso: `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2442: @*/
2443: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2444: {
2447:   PetscSectionCheckValidField(field, s->numFields);
2448:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
2449:   return 0;
2450: }

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

2455:   Not Collective

2457:   Input Parameters:
2458: + s       - The PetscSection
2459: . point   - The point
2460: . field   - The field number
2461: - indices - The constrained dofs

2463:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

2465:   Level: intermediate

2467: .seealso: `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2468: @*/
2469: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2470: {
2472:   if (PetscDefined(USE_DEBUG)) {
2473:     PetscInt nfdof;

2475:     PetscSectionGetFieldConstraintDof(s, point, field, &nfdof);
2477:   }
2478:   PetscSectionCheckValidField(field, s->numFields);
2479:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
2480:   return 0;
2481: }

2483: /*@
2484:   PetscSectionPermute - Reorder the section according to the input point permutation

2486:   Collective

2488:   Input Parameters:
2489: + section - The PetscSection object
2490: - perm - The point permutation, old point p becomes new point perm[p]

2492:   Output Parameter:
2493: . sectionNew - The permuted PetscSection

2495:   Level: intermediate

2497: .seealso: `MatPermute()`
2498: @*/
2499: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2500: {
2501:   PetscSection    s = section, sNew;
2502:   const PetscInt *perm;
2503:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;

2508:   PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew);
2509:   PetscSectionGetNumFields(s, &numFields);
2510:   if (numFields) PetscSectionSetNumFields(sNew, numFields);
2511:   for (f = 0; f < numFields; ++f) {
2512:     const char *name;
2513:     PetscInt    numComp;

2515:     PetscSectionGetFieldName(s, f, &name);
2516:     PetscSectionSetFieldName(sNew, f, name);
2517:     PetscSectionGetFieldComponents(s, f, &numComp);
2518:     PetscSectionSetFieldComponents(sNew, f, numComp);
2519:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2520:       PetscSectionGetComponentName(s, f, c, &name);
2521:       PetscSectionSetComponentName(sNew, f, c, name);
2522:     }
2523:   }
2524:   ISGetLocalSize(permutation, &numPoints);
2525:   ISGetIndices(permutation, &perm);
2526:   PetscSectionGetChart(s, &pStart, &pEnd);
2527:   PetscSectionSetChart(sNew, pStart, pEnd);
2529:   for (p = pStart; p < pEnd; ++p) {
2530:     PetscInt dof, cdof;

2532:     PetscSectionGetDof(s, p, &dof);
2533:     PetscSectionSetDof(sNew, perm[p], dof);
2534:     PetscSectionGetConstraintDof(s, p, &cdof);
2535:     if (cdof) PetscSectionSetConstraintDof(sNew, perm[p], cdof);
2536:     for (f = 0; f < numFields; ++f) {
2537:       PetscSectionGetFieldDof(s, p, f, &dof);
2538:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2539:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2540:       if (cdof) PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);
2541:     }
2542:   }
2543:   PetscSectionSetUp(sNew);
2544:   for (p = pStart; p < pEnd; ++p) {
2545:     const PetscInt *cind;
2546:     PetscInt        cdof;

2548:     PetscSectionGetConstraintDof(s, p, &cdof);
2549:     if (cdof) {
2550:       PetscSectionGetConstraintIndices(s, p, &cind);
2551:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2552:     }
2553:     for (f = 0; f < numFields; ++f) {
2554:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2555:       if (cdof) {
2556:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2557:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2558:       }
2559:     }
2560:   }
2561:   ISRestoreIndices(permutation, &perm);
2562:   *sectionNew = sNew;
2563:   return 0;
2564: }

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

2569:   Collective

2571:   Input Parameters:
2572: + section   - The PetscSection
2573: . obj       - A PetscObject which serves as the key for this index
2574: . clSection - Section giving the size of the closure of each point
2575: - clPoints  - IS giving the points in each closure

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

2579:   Level: advanced

2581: .seealso: `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2582: @*/
2583: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2584: {
2588:   if (section->clObj != obj) PetscSectionResetClosurePermutation(section);
2589:   section->clObj = obj;
2590:   PetscObjectReference((PetscObject)clSection);
2591:   PetscObjectReference((PetscObject)clPoints);
2592:   PetscSectionDestroy(&section->clSection);
2593:   ISDestroy(&section->clPoints);
2594:   section->clSection = clSection;
2595:   section->clPoints  = clPoints;
2596:   return 0;
2597: }

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

2602:   Collective

2604:   Input Parameters:
2605: + section   - The PetscSection
2606: - obj       - A PetscObject which serves as the key for this index

2608:   Output Parameters:
2609: + clSection - Section giving the size of the closure of each point
2610: - clPoints  - IS giving the points in each closure

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

2614:   Level: advanced

2616: .seealso: `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2617: @*/
2618: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2619: {
2620:   if (section->clObj == obj) {
2621:     if (clSection) *clSection = section->clSection;
2622:     if (clPoints) *clPoints = section->clPoints;
2623:   } else {
2624:     if (clSection) *clSection = NULL;
2625:     if (clPoints) *clPoints = NULL;
2626:   }
2627:   return 0;
2628: }

2630: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2631: {
2632:   PetscInt                    i;
2633:   khiter_t                    iter;
2634:   int                         new_entry;
2635:   PetscSectionClosurePermKey  key = {depth, clSize};
2636:   PetscSectionClosurePermVal *val;

2638:   if (section->clObj != obj) {
2639:     PetscSectionDestroy(&section->clSection);
2640:     ISDestroy(&section->clPoints);
2641:   }
2642:   section->clObj = obj;
2643:   if (!section->clHash) PetscClPermCreate(&section->clHash);
2644:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2645:   val  = &kh_val(section->clHash, iter);
2646:   if (!new_entry) {
2647:     PetscFree(val->perm);
2648:     PetscFree(val->invPerm);
2649:   }
2650:   if (mode == PETSC_COPY_VALUES) {
2651:     PetscMalloc1(clSize, &val->perm);
2652:     PetscArraycpy(val->perm, clPerm, clSize);
2653:   } else if (mode == PETSC_OWN_POINTER) {
2654:     val->perm = clPerm;
2655:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2656:   PetscMalloc1(clSize, &val->invPerm);
2657:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2658:   return 0;
2659: }

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

2664:   Not Collective

2666:   Input Parameters:
2667: + section - The PetscSection
2668: . obj     - A PetscObject which serves as the key for this index (usually a DM)
2669: . depth   - Depth of points on which to apply the given permutation
2670: - perm    - Permutation of the cell dof closure

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

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

2680:   Level: intermediate

2682: .seealso: `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2683: @*/
2684: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2685: {
2686:   const PetscInt *clPerm = NULL;
2687:   PetscInt        clSize = 0;

2689:   if (perm) {
2690:     ISGetLocalSize(perm, &clSize);
2691:     ISGetIndices(perm, &clPerm);
2692:   }
2693:   PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm);
2694:   if (perm) ISRestoreIndices(perm, &clPerm);
2695:   return 0;
2696: }

2698: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2699: {
2700:   if (section->clObj == obj) {
2701:     PetscSectionClosurePermKey k = {depth, size};
2702:     PetscSectionClosurePermVal v;
2703:     PetscClPermGet(section->clHash, k, &v);
2704:     if (perm) *perm = v.perm;
2705:   } else {
2706:     if (perm) *perm = NULL;
2707:   }
2708:   return 0;
2709: }

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

2714:   Not Collective

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

2722:   Output Parameter:
2723: . perm - The dof closure permutation

2725:   Note:
2726:   The user must destroy the IS that is returned.

2728:   Level: intermediate

2730: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2731: @*/
2732: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2733: {
2734:   const PetscInt *clPerm;

2736:   PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2737:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2738:   return 0;
2739: }

2741: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2742: {
2743:   if (section->clObj == obj && section->clHash) {
2744:     PetscSectionClosurePermKey k = {depth, size};
2745:     PetscSectionClosurePermVal v;
2746:     PetscClPermGet(section->clHash, k, &v);
2747:     if (perm) *perm = v.invPerm;
2748:   } else {
2749:     if (perm) *perm = NULL;
2750:   }
2751:   return 0;
2752: }

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

2757:   Not Collective

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

2765:   Output Parameters:
2766: . perm - The dof closure permutation

2768:   Note:
2769:   The user must destroy the IS that is returned.

2771:   Level: intermediate

2773: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2774: @*/
2775: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2776: {
2777:   const PetscInt *clPerm;

2779:   PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2780:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2781:   return 0;
2782: }

2784: /*@
2785:   PetscSectionGetField - Get the subsection associated with a single field

2787:   Input Parameters:
2788: + s     - The PetscSection
2789: - field - The field number

2791:   Output Parameter:
2792: . subs  - The subsection for the given field

2794:   Level: intermediate

2796: .seealso: `PetscSectionSetNumFields()`
2797: @*/
2798: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2799: {
2802:   PetscSectionCheckValidField(field, s->numFields);
2803:   *subs = s->field[field];
2804:   return 0;
2805: }

2807: PetscClassId      PETSC_SECTION_SYM_CLASSID;
2808: PetscFunctionList PetscSectionSymList = NULL;

2810: /*@
2811:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

2813:   Collective

2815:   Input Parameter:
2816: . comm - the MPI communicator

2818:   Output Parameter:
2819: . sym - pointer to the new set of symmetries

2821:   Level: developer

2823: .seealso: `PetscSectionSym`, `PetscSectionSymDestroy()`
2824: @*/
2825: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2826: {
2828:   ISInitializePackage();
2829:   PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView);
2830:   return 0;
2831: }

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

2836:   Collective

2838:   Input Parameters:
2839: + sym    - The section symmetry object
2840: - method - The name of the section symmetry type

2842:   Level: developer

2844: .seealso: `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
2845: @*/
2846: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2847: {
2848:   PetscErrorCode (*r)(PetscSectionSym);
2849:   PetscBool match;

2852:   PetscObjectTypeCompare((PetscObject)sym, method, &match);
2853:   if (match) return 0;

2855:   PetscFunctionListFind(PetscSectionSymList, method, &r);
2857:   PetscTryTypeMethod(sym, destroy);
2858:   sym->ops->destroy = NULL;

2860:   (*r)(sym);
2861:   PetscObjectChangeTypeName((PetscObject)sym, method);
2862:   return 0;
2863: }

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

2868:   Not Collective

2870:   Input Parameter:
2871: . sym  - The section symmetry

2873:   Output Parameter:
2874: . type - The index set type name

2876:   Level: developer

2878: .seealso: `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
2879: @*/
2880: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2881: {
2884:   *type = ((PetscObject)sym)->type_name;
2885:   return 0;
2886: }

2888: /*@C
2889:   PetscSectionSymRegister - Adds a new section symmetry implementation

2891:   Not Collective

2893:   Input Parameters:
2894: + name        - The name of a new user-defined creation routine
2895: - create_func - The creation routine itself

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

2900:   Level: developer

2902: .seealso: `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
2903: @*/
2904: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2905: {
2906:   ISInitializePackage();
2907:   PetscFunctionListAdd(&PetscSectionSymList, sname, function);
2908:   return 0;
2909: }

2911: /*@
2912:    PetscSectionSymDestroy - Destroys a section symmetry.

2914:    Collective

2916:    Input Parameters:
2917: .  sym - the section symmetry

2919:    Level: developer

2921: .seealso: `PetscSectionSymCreate()`, `PetscSectionSymDestroy()`
2922: @*/
2923: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2924: {
2925:   SymWorkLink link, next;

2927:   if (!*sym) return 0;
2929:   if (--((PetscObject)(*sym))->refct > 0) {
2930:     *sym = NULL;
2931:     return 0;
2932:   }
2933:   if ((*sym)->ops->destroy) (*(*sym)->ops->destroy)(*sym);
2935:   for (link = (*sym)->workin; link; link = next) {
2936:     PetscInt    **perms = (PetscInt **)link->perms;
2937:     PetscScalar **rots  = (PetscScalar **)link->rots;
2938:     PetscFree2(perms, rots);
2939:     next = link->next;
2940:     PetscFree(link);
2941:   }
2942:   (*sym)->workin = NULL;
2943:   PetscHeaderDestroy(sym);
2944:   return 0;
2945: }

2947: /*@C
2948:    PetscSectionSymView - Displays a section symmetry

2950:    Collective

2952:    Input Parameters:
2953: +  sym - the index set
2954: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

2956:    Level: developer

2958: .seealso: `PetscViewerASCIIOpen()`
2959: @*/
2960: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
2961: {
2963:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer);
2966:   PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer);
2967:   PetscTryTypeMethod(sym, view, viewer);
2968:   return 0;
2969: }

2971: /*@
2972:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

2974:   Collective

2976:   Input Parameters:
2977: + section - the section describing data layout
2978: - sym - the symmetry describing the affect of orientation on the access of the data

2980:   Level: developer

2982: .seealso: `PetscSectionGetSym()`, `PetscSectionSymCreate()`
2983: @*/
2984: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2985: {
2987:   PetscSectionSymDestroy(&(section->sym));
2988:   if (sym) {
2991:     PetscObjectReference((PetscObject)sym);
2992:   }
2993:   section->sym = sym;
2994:   return 0;
2995: }

2997: /*@
2998:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3000:   Not Collective

3002:   Input Parameters:
3003: . section - the section describing data layout

3005:   Output Parameters:
3006: . sym - the symmetry describing the affect of orientation on the access of the data

3008:   Level: developer

3010: .seealso: `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3011: @*/
3012: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3013: {
3015:   *sym = section->sym;
3016:   return 0;
3017: }

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

3022:   Collective

3024:   Input Parameters:
3025: + section - the section describing data layout
3026: . field - the field number
3027: - sym - the symmetry describing the affect of orientation on the access of the data

3029:   Level: developer

3031: .seealso: `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3032: @*/
3033: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3034: {
3036:   PetscSectionCheckValidField(field, section->numFields);
3037:   PetscSectionSetSym(section->field[field], sym);
3038:   return 0;
3039: }

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

3044:   Collective

3046:   Input Parameters:
3047: + section - the section describing data layout
3048: - field - the field number

3050:   Output Parameters:
3051: . sym - the symmetry describing the affect of orientation on the access of the data

3053:   Level: developer

3055: .seealso: `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3056: @*/
3057: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3058: {
3060:   PetscSectionCheckValidField(field, section->numFields);
3061:   *sym = section->field[field]->sym;
3062:   return 0;
3063: }

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

3068:   Not Collective

3070:   Input Parameters:
3071: + section - the section
3072: . numPoints - the number of points
3073: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3074:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3075:     context, see DMPlexGetConeOrientation()).

3077:   Output Parameters:
3078: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3079: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3080:     identity).

3082:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3083: .vb
3084:      const PetscInt    **perms;
3085:      const PetscScalar **rots;
3086:      PetscInt            lOffset;

3088:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3089:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3090:        PetscInt           point = points[2*i], dof, sOffset;
3091:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3092:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3094:        PetscSectionGetDof(section,point,&dof);
3095:        PetscSectionGetOffset(section,point,&sOffset);

3097:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3098:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3099:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3100:        lOffset += dof;
3101:      }
3102:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3103: .ve

3105:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3106: .vb
3107:      const PetscInt    **perms;
3108:      const PetscScalar **rots;
3109:      PetscInt            lOffset;

3111:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3112:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3113:        PetscInt           point = points[2*i], dof, sOffset;
3114:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3115:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3117:        PetscSectionGetDof(section,point,&dof);
3118:        PetscSectionGetOffset(section,point,&sOff);

3120:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3121:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3122:        offset += dof;
3123:      }
3124:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3125: .ve

3127:   Level: developer

3129: .seealso: `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3130: @*/
3131: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3132: {
3133:   PetscSectionSym sym;

3137:   if (perms) *perms = NULL;
3138:   if (rots) *rots = NULL;
3139:   sym = section->sym;
3140:   if (sym && (perms || rots)) {
3141:     SymWorkLink link;

3143:     if (sym->workin) {
3144:       link        = sym->workin;
3145:       sym->workin = sym->workin->next;
3146:     } else {
3147:       PetscNew(&link);
3148:     }
3149:     if (numPoints > link->numPoints) {
3150:       PetscInt    **perms = (PetscInt **)link->perms;
3151:       PetscScalar **rots  = (PetscScalar **)link->rots;
3152:       PetscFree2(perms, rots);
3153:       PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots);
3154:       link->numPoints = numPoints;
3155:     }
3156:     link->next   = sym->workout;
3157:     sym->workout = link;
3158:     PetscArrayzero((PetscInt **)link->perms, numPoints);
3159:     PetscArrayzero((PetscInt **)link->rots, numPoints);
3160:     (*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots);
3161:     if (perms) *perms = link->perms;
3162:     if (rots) *rots = link->rots;
3163:   }
3164:   return 0;
3165: }

3167: /*@C
3168:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

3170:   Not Collective

3172:   Input Parameters:
3173: + section - the section
3174: . numPoints - the number of points
3175: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3176:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3177:     context, see DMPlexGetConeOrientation()).

3179:   Output Parameters:
3180: + perms - The permutations for the given orientations: set to NULL at conclusion
3181: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3183:   Level: developer

3185: .seealso: `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3186: @*/
3187: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3188: {
3189:   PetscSectionSym sym;

3192:   sym = section->sym;
3193:   if (sym && (perms || rots)) {
3194:     SymWorkLink *p, link;

3196:     for (p = &sym->workout; (link = *p); p = &link->next) {
3197:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3198:         *p          = link->next;
3199:         link->next  = sym->workin;
3200:         sym->workin = link;
3201:         if (perms) *perms = NULL;
3202:         if (rots) *rots = NULL;
3203:         return 0;
3204:       }
3205:     }
3206:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3207:   }
3208:   return 0;
3209: }

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

3214:   Not Collective

3216:   Input Parameters:
3217: + section - the section
3218: . field - the field of the section
3219: . numPoints - the number of points
3220: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3221:     arbitrary integer: its interpretation is up to sym.  Orientations are used by DM: for their interpretation in that
3222:     context, see DMPlexGetConeOrientation()).

3224:   Output Parameters:
3225: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3226: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3227:     identity).

3229:   Level: developer

3231: .seealso: `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3232: @*/
3233: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3234: {
3237:   PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots);
3238:   return 0;
3239: }

3241: /*@C
3242:   PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()

3244:   Not Collective

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

3254:   Output Parameters:
3255: + perms - The permutations for the given orientations: set to NULL at conclusion
3256: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

3258:   Level: developer

3260: .seealso: `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3261: @*/
3262: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3263: {
3266:   PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots);
3267:   return 0;
3268: }

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

3273:   Not Collective

3275:   Input Parameter:
3276: . sym - the PetscSectionSym

3278:   Output Parameter:
3279: . nsym - the equivalent symmetries

3281:   Level: developer

3283: .seealso: `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3284: @*/
3285: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3286: {
3289:   PetscTryTypeMethod(sym, copy, nsym);
3290:   return 0;
3291: }

3293: /*@
3294:   PetscSectionSymDistribute - Distribute the symmetries in accordance with the input SF

3296:   Collective

3298:   Input Parameters:
3299: + sym - the PetscSectionSym
3300: - migrationSF - the distribution map from roots to leaves

3302:   Output Parameters:
3303: . dsym - the redistributed symmetries

3305:   Level: developer

3307: .seealso: `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3308: @*/
3309: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3310: {
3314:   PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3315:   return 0;
3316: }

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

3321:   Not Collective

3323:   Input Parameter:
3324: . s - the global PetscSection

3326:   Output Parameters:
3327: . flg - the flag

3329:   Level: developer

3331: .seealso: `PetscSectionSetChart()`, `PetscSectionCreate()`
3332: @*/
3333: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3334: {
3336:   *flg = s->useFieldOff;
3337:   return 0;
3338: }

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

3343:   Not Collective

3345:   Input Parameters:
3346: + s   - the global PetscSection
3347: - flg - the flag

3349:   Level: developer

3351: .seealso: `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3352: @*/
3353: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3354: {
3356:   s->useFieldOff = flg;
3357:   return 0;
3358: }

3360: #define PetscSectionExpandPoints_Loop(TYPE) \
3361:   { \
3362:     PetscInt i, n, o0, o1, size; \
3363:     TYPE    *a0 = (TYPE *)origArray, *a1; \
3364:     PetscSectionGetStorageSize(s, &size); \
3365:     PetscMalloc1(size, &a1); \
3366:     for (i = 0; i < npoints; i++) { \
3367:       PetscSectionGetOffset(origSection, points_[i], &o0); \
3368:       PetscSectionGetOffset(s, i, &o1); \
3369:       PetscSectionGetDof(s, i, &n); \
3370:       PetscMemcpy(&a1[o1], &a0[o0], n *unitsize); \
3371:     } \
3372:     *newArray = (void *)a1; \
3373:   }

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

3378:   Not Collective

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

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

3390:   Level: developer

3392: .seealso: `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3393: @*/
3394: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3395: {
3396:   PetscSection    s;
3397:   const PetscInt *points_;
3398:   PetscInt        i, n, npoints, pStart, pEnd;
3399:   PetscMPIInt     unitsize;

3406:   MPI_Type_size(dataType, &unitsize);
3407:   ISGetLocalSize(points, &npoints);
3408:   ISGetIndices(points, &points_);
3409:   PetscSectionGetChart(origSection, &pStart, &pEnd);
3410:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3411:   PetscSectionSetChart(s, 0, npoints);
3412:   for (i = 0; i < npoints; i++) {
3414:     PetscSectionGetDof(origSection, points_[i], &n);
3415:     PetscSectionSetDof(s, i, n);
3416:   }
3417:   PetscSectionSetUp(s);
3418:   if (newArray) {
3419:     if (dataType == MPIU_INT) {
3420:       PetscSectionExpandPoints_Loop(PetscInt);
3421:     } else if (dataType == MPIU_SCALAR) {
3422:       PetscSectionExpandPoints_Loop(PetscScalar);
3423:     } else if (dataType == MPIU_REAL) {
3424:       PetscSectionExpandPoints_Loop(PetscReal);
3425:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3426:   }
3427:   if (newSection) {
3428:     *newSection = s;
3429:   } else {
3430:     PetscSectionDestroy(&s);
3431:   }
3432:   ISRestoreIndices(points, &points_);
3433:   return 0;
3434: }