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 a `PetscSection` and sets the map contents to the default.

 13:   Collective

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

 19:   Level: beginner

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

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

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

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

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

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

 71:   Collective

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

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

 79:   Level: intermediate

 81:   Developer Notes:
 82:   What exactly does shallow mean in this context?

 84: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
 85: @*/
 86: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
 87: {
 88:   PetscFunctionBegin;
 91:   PetscCall(PetscSectionCopy_Internal(section, newSection, NULL));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: PetscErrorCode PetscSectionCopy_Internal(PetscSection section, PetscSection newSection, PetscBT constrained_dofs)
 96: {
 97:   PetscSectionSym sym;
 98:   IS              perm;
 99:   PetscInt        numFields, f, c, pStart, pEnd, p;

101:   PetscFunctionBegin;
104:   PetscCall(PetscSectionReset(newSection));
105:   PetscCall(PetscSectionGetNumFields(section, &numFields));
106:   if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields));
107:   for (f = 0; f < numFields; ++f) {
108:     const char *fieldName = NULL, *compName = NULL;
109:     PetscInt    numComp = 0;

111:     PetscCall(PetscSectionGetFieldName(section, f, &fieldName));
112:     PetscCall(PetscSectionSetFieldName(newSection, f, fieldName));
113:     PetscCall(PetscSectionGetFieldComponents(section, f, &numComp));
114:     PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp));
115:     for (c = 0; c < numComp; ++c) {
116:       PetscCall(PetscSectionGetComponentName(section, f, c, &compName));
117:       PetscCall(PetscSectionSetComponentName(newSection, f, c, compName));
118:     }
119:     PetscCall(PetscSectionGetFieldSym(section, f, &sym));
120:     PetscCall(PetscSectionSetFieldSym(newSection, f, sym));
121:   }
122:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
123:   PetscCall(PetscSectionSetChart(newSection, pStart, pEnd));
124:   PetscCall(PetscSectionGetPermutation(section, &perm));
125:   PetscCall(PetscSectionSetPermutation(newSection, perm));
126:   PetscCall(PetscSectionGetSym(section, &sym));
127:   PetscCall(PetscSectionSetSym(newSection, sym));
128:   for (p = pStart; p < pEnd; ++p) {
129:     PetscInt  dof, cdof, fcdof = 0;
130:     PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));

132:     PetscCall(PetscSectionGetDof(section, p, &dof));
133:     PetscCall(PetscSectionSetDof(newSection, p, dof));
134:     if (force_constrained) cdof = dof;
135:     else PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
136:     if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof));
137:     for (f = 0; f < numFields; ++f) {
138:       PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
139:       PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof));
140:       if (cdof) {
141:         if (force_constrained) fcdof = dof;
142:         else PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
143:         if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof));
144:       }
145:     }
146:   }
147:   PetscCall(PetscSectionSetUp(newSection));
148:   for (p = pStart; p < pEnd; ++p) {
149:     PetscInt        off, cdof, fcdof = 0;
150:     const PetscInt *cInd;
151:     PetscBool       force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));

153:     /* Must set offsets in case they do not agree with the prefix sums */
154:     PetscCall(PetscSectionGetOffset(section, p, &off));
155:     PetscCall(PetscSectionSetOffset(newSection, p, off));
156:     PetscCall(PetscSectionGetConstraintDof(newSection, p, &cdof));
157:     if (cdof) {
158:       if (force_constrained) cInd = NULL;
159:       else PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd));
160:       PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd));
161:       for (f = 0; f < numFields; ++f) {
162:         PetscCall(PetscSectionGetFieldOffset(section, p, f, &off));
163:         PetscCall(PetscSectionSetFieldOffset(newSection, p, f, off));
164:         PetscCall(PetscSectionGetFieldConstraintDof(newSection, p, f, &fcdof));
165:         if (fcdof) {
166:           if (force_constrained) cInd = NULL;
167:           else PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd));
168:           PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd));
169:         }
170:       }
171:     }
172:   }
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

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

179:   Collective

181:   Input Parameter:
182: . section - the `PetscSection`

184:   Output Parameter:
185: . newSection - the copy

187:   Level: beginner

189:   Developer Notes:
190:   With standard PETSc terminology this should be called `PetscSectionDuplicate()`

192: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
193: @*/
194: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
195: {
196:   PetscFunctionBegin;
198:   PetscAssertPointer(newSection, 2);
199:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
200:   PetscCall(PetscSectionCopy(section, *newSection));
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

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

207:   Collective

209:   Input Parameter:
210: . s - the `PetscSection`

212:   Options Database Key:
213: . -petscsection_point_major - `PETSC_TRUE` for point-major order

215:   Level: intermediate

217: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
218: @*/
219: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
220: {
221:   PetscFunctionBegin;
223:   PetscObjectOptionsBegin((PetscObject)s);
224:   PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL));
225:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
226:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject));
227:   PetscOptionsEnd();
228:   PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view"));
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: /*@
233:   PetscSectionCompare - Compares two sections

235:   Collective

237:   Input Parameters:
238: + s1 - the first `PetscSection`
239: - s2 - the second `PetscSection`

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

244:   Level: intermediate

246:   Note:
247:   Field names are disregarded.

249: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
250: @*/
251: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
252: {
253:   PetscInt        pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
254:   const PetscInt *idx1, *idx2;
255:   IS              perm1, perm2;
256:   PetscBool       flg;
257:   PetscMPIInt     mflg;

259:   PetscFunctionBegin;
262:   PetscAssertPointer(congruent, 3);
263:   flg = PETSC_FALSE;

265:   PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
266:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
267:     *congruent = PETSC_FALSE;
268:     PetscFunctionReturn(PETSC_SUCCESS);
269:   }

271:   PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd));
272:   PetscCall(PetscSectionGetChart(s2, &n1, &n2));
273:   if (pStart != n1 || pEnd != n2) goto not_congruent;

275:   PetscCall(PetscSectionGetPermutation(s1, &perm1));
276:   PetscCall(PetscSectionGetPermutation(s2, &perm2));
277:   if (perm1 && perm2) {
278:     PetscCall(ISEqual(perm1, perm2, congruent));
279:     if (!*congruent) goto not_congruent;
280:   } else if (perm1 != perm2) goto not_congruent;

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

287:     PetscCall(PetscSectionGetDof(s1, p, &n1));
288:     PetscCall(PetscSectionGetDof(s2, p, &n2));
289:     if (n1 != n2) goto not_congruent;

291:     PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
292:     PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
293:     if (ncdof != n2) goto not_congruent;

295:     PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1));
296:     PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2));
297:     PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent));
298:     if (!*congruent) goto not_congruent;
299:   }

301:   PetscCall(PetscSectionGetNumFields(s1, &nfields));
302:   PetscCall(PetscSectionGetNumFields(s2, &n2));
303:   if (nfields != n2) goto not_congruent;

305:   for (f = 0; f < nfields; ++f) {
306:     PetscCall(PetscSectionGetFieldComponents(s1, f, &n1));
307:     PetscCall(PetscSectionGetFieldComponents(s2, f, &n2));
308:     if (n1 != n2) goto not_congruent;

310:     for (p = pStart; p < pEnd; ++p) {
311:       PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1));
312:       PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2));
313:       if (n1 != n2) goto not_congruent;

315:       PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
316:       PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
317:       if (n1 != n2) goto not_congruent;

319:       PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
320:       PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
321:       if (nfcdof != n2) goto not_congruent;

323:       PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1));
324:       PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2));
325:       PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent));
326:       if (!*congruent) goto not_congruent;
327:     }
328:   }

330:   flg = PETSC_TRUE;
331: not_congruent:
332:   PetscCall(MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

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

339:   Not Collective

341:   Input Parameter:
342: . s - the `PetscSection`

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

347:   Level: intermediate

349: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
350: @*/
351: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
352: {
353:   PetscFunctionBegin;
355:   PetscAssertPointer(numFields, 2);
356:   *numFields = s->numFields;
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: /*@
361:   PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`

363:   Not Collective

365:   Input Parameters:
366: + s         - the `PetscSection`
367: - numFields - the number of fields

369:   Level: intermediate

371:   Notes:
372:   Calling this destroys all the information in the `PetscSection` including the chart.

374:   You must call `PetscSectionSetChart()` after calling this.

376: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
377: @*/
378: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
379: {
380:   PetscInt f;

382:   PetscFunctionBegin;
384:   PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields);
385:   PetscCall(PetscSectionReset(s));

387:   s->numFields = numFields;
388:   PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents));
389:   PetscCall(PetscMalloc1(s->numFields, &s->fieldNames));
390:   PetscCall(PetscMalloc1(s->numFields, &s->compNames));
391:   PetscCall(PetscMalloc1(s->numFields, &s->field));
392:   for (f = 0; f < s->numFields; ++f) {
393:     char name[64];

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

397:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
398:     PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
399:     PetscCall(PetscStrallocpy(name, (char **)&s->fieldNames[f]));
400:     PetscCall(PetscSNPrintf(name, 64, "Component_0"));
401:     PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
402:     PetscCall(PetscStrallocpy(name, (char **)&s->compNames[f][0]));
403:   }
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }

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

410:   Not Collective

412:   Input Parameters:
413: + s     - the `PetscSection`
414: - field - the field number

416:   Output Parameter:
417: . fieldName - the field name

419:   Level: intermediate

421:   Note:
422:   Will error if the field number is out of range

424: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
425: @*/
426: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
427: {
428:   PetscFunctionBegin;
430:   PetscAssertPointer(fieldName, 3);
431:   PetscSectionCheckValidField(field, s->numFields);
432:   *fieldName = s->fieldNames[field];
433:   PetscFunctionReturn(PETSC_SUCCESS);
434: }

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

439:   Not Collective

441:   Input Parameters:
442: + s         - the `PetscSection`
443: . field     - the field number
444: - fieldName - the field name

446:   Level: intermediate

448:   Note:
449:   Will error if the field number is out of range

451: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
452: @*/
453: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
454: {
455:   PetscFunctionBegin;
457:   if (fieldName) PetscAssertPointer(fieldName, 3);
458:   PetscSectionCheckValidField(field, s->numFields);
459:   PetscCall(PetscFree(s->fieldNames[field]));
460:   PetscCall(PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]));
461:   PetscFunctionReturn(PETSC_SUCCESS);
462: }

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

467:   Not Collective

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

474:   Output Parameter:
475: . compName - the component name

477:   Level: intermediate

479:   Note:
480:   Will error if the field or component number do not exist

482:   Developer Notes:
483:   The function name should have Field in it since they are field components.

485: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
486:           `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
487: @*/
488: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
489: {
490:   PetscFunctionBegin;
492:   PetscAssertPointer(compName, 4);
493:   PetscSectionCheckValidField(field, s->numFields);
494:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
495:   *compName = s->compNames[field][comp];
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }

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

502:   Not Collective

504:   Input Parameters:
505: + s        - the `PetscSection`
506: . field    - the field number
507: . comp     - the component number
508: - compName - the component name

510:   Level: advanced

512:   Note:
513:   Will error if the field or component number do not exist

515:   Developer Notes:
516:   The function name should have Field in it since they are field components.

518: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
519:           `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
520: @*/
521: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
522: {
523:   PetscFunctionBegin;
525:   if (compName) PetscAssertPointer(compName, 4);
526:   PetscSectionCheckValidField(field, s->numFields);
527:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
528:   PetscCall(PetscFree(s->compNames[field][comp]));
529:   PetscCall(PetscStrallocpy(compName, (char **)&s->compNames[field][comp]));
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: /*@
534:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

536:   Not Collective

538:   Input Parameters:
539: + s     - the `PetscSection`
540: - field - the field number

542:   Output Parameter:
543: . numComp - the number of field components

545:   Level: advanced

547:   Developer Notes:
548:   This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name

550: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
551:           `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
552: @*/
553: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
554: {
555:   PetscFunctionBegin;
557:   PetscAssertPointer(numComp, 3);
558:   PetscSectionCheckValidField(field, s->numFields);
559:   *numComp = s->numFieldComponents[field];
560:   PetscFunctionReturn(PETSC_SUCCESS);
561: }

563: /*@
564:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

566:   Not Collective

568:   Input Parameters:
569: + s       - the `PetscSection`
570: . field   - the field number
571: - numComp - the number of field components

573:   Level: advanced

575:   Note:
576:   This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
577:   components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
578:   at a point in a discretization. For example, if in three dimensions the field is velocity, it will have 3 components, u, v, and w but
579:   an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.

581:   The value set with this function are not needed or used in `PetscSectionSetUp()`.

583:   Developer Notes:
584:   This function is misnamed. There is a Num in `PetscSectionSetNumFields()` but not in this name

586: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
587:           `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
588: @*/
589: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
590: {
591:   PetscInt c;

593:   PetscFunctionBegin;
595:   PetscSectionCheckValidField(field, s->numFields);
596:   if (s->compNames) {
597:     for (c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
598:     PetscCall(PetscFree(s->compNames[field]));
599:   }

601:   s->numFieldComponents[field] = numComp;
602:   if (numComp) {
603:     PetscCall(PetscMalloc1(numComp, (char ***)&s->compNames[field]));
604:     for (c = 0; c < numComp; ++c) {
605:       char name[64];

607:       PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
608:       PetscCall(PetscStrallocpy(name, (char **)&s->compNames[field][c]));
609:     }
610:   }
611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

614: /*@
615:   PetscSectionGetChart - Returns the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process

617:   Not Collective

619:   Input Parameter:
620: . s - the `PetscSection`

622:   Output Parameters:
623: + pStart - the first point
624: - pEnd   - one past the last point

626:   Level: intermediate

628: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
629: @*/
630: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
631: {
632:   PetscFunctionBegin;
634:   if (pStart) *pStart = s->pStart;
635:   if (pEnd) *pEnd = s->pEnd;
636:   PetscFunctionReturn(PETSC_SUCCESS);
637: }

639: /*@
640:   PetscSectionSetChart - Sets the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process

642:   Not Collective

644:   Input Parameters:
645: + s      - the `PetscSection`
646: . pStart - the first point
647: - pEnd   - one past the last point, `pStart` $ \le $ `pEnd`

649:   Level: intermediate

651:   Notes:
652:   The charts on different MPI processes may (and often do) overlap

654:   If you intend to use `PetscSectionSetNumFields()` it must be called before this call.

656:   The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.

658: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
659: @*/
660: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
661: {
662:   PetscInt f;

664:   PetscFunctionBegin;
666:   PetscCheck(pEnd >= pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Chart pEnd %" PetscInt_FMT " cannot be smaller than chart pStart %" PetscInt_FMT, pEnd, pStart);
667:   if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
668:   /* Cannot Reset() because it destroys field information */
669:   s->setup = PETSC_FALSE;
670:   PetscCall(PetscSectionDestroy(&s->bc));
671:   PetscCall(PetscFree(s->bcIndices));
672:   PetscCall(PetscFree2(s->atlasDof, s->atlasOff));

674:   s->pStart = pStart;
675:   s->pEnd   = pEnd;
676:   PetscCall(PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff));
677:   PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
678:   for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
679:   PetscFunctionReturn(PETSC_SUCCESS);
680: }

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

685:   Not Collective

687:   Input Parameter:
688: . s - the `PetscSection`

690:   Output Parameter:
691: . perm - The permutation as an `IS`

693:   Level: intermediate

695: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
696: @*/
697: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
698: {
699:   PetscFunctionBegin;
701:   if (perm) {
702:     PetscAssertPointer(perm, 2);
703:     *perm = s->perm;
704:   }
705:   PetscFunctionReturn(PETSC_SUCCESS);
706: }

708: /*@
709:   PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information

711:   Not Collective

713:   Input Parameters:
714: + s    - the `PetscSection`
715: - perm - the permutation of points

717:   Level: intermediate

719:   Notes:
720:   The permutation must be provided before `PetscSectionSetUp()`.

722:   The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed

724:   Compare to `PetscSectionPermute()`

726: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
727: @*/
728: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
729: {
730:   PetscFunctionBegin;
733:   PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
734:   if (s->perm != perm) {
735:     PetscCall(ISDestroy(&s->perm));
736:     if (perm) {
737:       s->perm = perm;
738:       PetscCall(PetscObjectReference((PetscObject)s->perm));
739:     }
740:   }
741:   PetscFunctionReturn(PETSC_SUCCESS);
742: }

744: /*@C
745:   PetscSectionGetBlockStarts - Returns a table indicating which points start new blocks

747:   Not Collective

749:   Input Parameter:
750: . s - the `PetscSection`

752:   Output Parameter:
753: . blockStarts - The `PetscBT` with a 1 for each point that begins a block

755:   Notes:
756:   The table is on [0, `pEnd` - `pStart`).

758:   This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.

760:   Level: intermediate

762: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
763: @*/
764: PetscErrorCode PetscSectionGetBlockStarts(PetscSection s, PetscBT *blockStarts)
765: {
766:   PetscFunctionBegin;
768:   if (blockStarts) {
769:     PetscAssertPointer(blockStarts, 2);
770:     *blockStarts = s->blockStarts;
771:   }
772:   PetscFunctionReturn(PETSC_SUCCESS);
773: }

775: /*@C
776:   PetscSectionSetBlockStarts - Sets a table indicating which points start new blocks

778:   Not Collective

780:   Input Parameters:
781: + s           - the `PetscSection`
782: - blockStarts - The `PetscBT` with a 1 for each point that begins a block

784:   Level: intermediate

786:   Notes:
787:   The table is on [0, `pEnd` - `pStart`). PETSc takes ownership of the `PetscBT` when it is passed in and will destroy it. The user should not destroy it.

789:   This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.

791: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
792: @*/
793: PetscErrorCode PetscSectionSetBlockStarts(PetscSection s, PetscBT blockStarts)
794: {
795:   PetscFunctionBegin;
797:   if (s->blockStarts != blockStarts) {
798:     PetscCall(PetscBTDestroy(&s->blockStarts));
799:     s->blockStarts = blockStarts;
800:   }
801:   PetscFunctionReturn(PETSC_SUCCESS);
802: }

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

807:   Not Collective

809:   Input Parameter:
810: . s - the `PetscSection`

812:   Output Parameter:
813: . pm - the flag for point major ordering

815:   Level: intermediate

817: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
818: @*/
819: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
820: {
821:   PetscFunctionBegin;
823:   PetscAssertPointer(pm, 2);
824:   *pm = s->pointMajor;
825:   PetscFunctionReturn(PETSC_SUCCESS);
826: }

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

831:   Not Collective

833:   Input Parameters:
834: + s  - the `PetscSection`
835: - pm - the flag for point major ordering

837:   Level: intermediate

839:   Note:
840:   Field-major order is not recommended unless you are managing the entire problem yourself, since many higher-level functions in PETSc depend on point-major order.

842:   Point major order means the degrees of freedom are stored as follows
843: .vb
844:     all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
845:     for each point
846:        the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
847: .ve

849:   Field major order means the degrees of freedom are stored as follows
850: .vb
851:     all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
852:     for each field (started with unnamed default field)
853:       the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
854: .ve

856: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
857: @*/
858: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
859: {
860:   PetscFunctionBegin;
862:   PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
863:   s->pointMajor = pm;
864:   PetscFunctionReturn(PETSC_SUCCESS);
865: }

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

871:   Not Collective

873:   Input Parameter:
874: . s - the `PetscSection`

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

879:   Level: intermediate

881: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
882: @*/
883: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
884: {
885:   PetscFunctionBegin;
887:   PetscAssertPointer(includesConstraints, 2);
888:   *includesConstraints = s->includesConstraints;
889:   PetscFunctionReturn(PETSC_SUCCESS);
890: }

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

895:   Not Collective

897:   Input Parameters:
898: + s                   - the `PetscSection`
899: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets

901:   Level: intermediate

903: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
904: @*/
905: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
906: {
907:   PetscFunctionBegin;
909:   PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
910:   s->includesConstraints = includesConstraints;
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }

914: /*@
915:   PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.

917:   Not Collective

919:   Input Parameters:
920: + s     - the `PetscSection`
921: - point - the point

923:   Output Parameter:
924: . numDof - the number of dof

926:   Level: intermediate

928:   Notes:
929:   In a global section, this size will be negative for points not owned by this process.

931:   This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point

933: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
934: @*/
935: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
936: {
937:   PetscFunctionBeginHot;
939:   PetscAssertPointer(numDof, 3);
940:   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);
941:   *numDof = s->atlasDof[point - s->pStart];
942:   PetscFunctionReturn(PETSC_SUCCESS);
943: }

945: /*@
946:   PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.

948:   Not Collective

950:   Input Parameters:
951: + s      - the `PetscSection`
952: . point  - the point
953: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process

955:   Level: intermediate

957:   Note:
958:   This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point

960: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
961: @*/
962: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
963: {
964:   PetscFunctionBegin;
966:   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);
967:   s->atlasDof[point - s->pStart] = numDof;
968:   PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

972: /*@
973:   PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.

975:   Not Collective

977:   Input Parameters:
978: + s      - the `PetscSection`
979: . point  - the point
980: - numDof - the number of additional dof

982:   Level: intermediate

984:   Note:
985:   This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point

987: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
988: @*/
989: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
990: {
991:   PetscFunctionBeginHot;
993:   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);
994:   PetscCheck(numDof >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numDof %" PetscInt_FMT " should not be negative", numDof);
995:   s->atlasDof[point - s->pStart] += numDof;
996:   PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
997:   PetscFunctionReturn(PETSC_SUCCESS);
998: }

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

1003:   Not Collective

1005:   Input Parameters:
1006: + s     - the `PetscSection`
1007: . point - the point
1008: - field - the field

1010:   Output Parameter:
1011: . numDof - the number of dof

1013:   Level: intermediate

1015: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
1016: @*/
1017: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1018: {
1019:   PetscFunctionBegin;
1021:   PetscAssertPointer(numDof, 4);
1022:   PetscSectionCheckValidField(field, s->numFields);
1023:   PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
1024:   PetscFunctionReturn(PETSC_SUCCESS);
1025: }

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

1030:   Not Collective

1032:   Input Parameters:
1033: + s      - the `PetscSection`
1034: . point  - the point
1035: . field  - the field
1036: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process

1038:   Level: intermediate

1040:   Note:
1041:   When setting the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1042:   the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`

1044:   This is equivalent to
1045: .vb
1046:      PetscSection fs;
1047:      PetscSectionGetField(s,field,&fs)
1048:      PetscSectionSetDof(fs,numDof)
1049: .ve

1051: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
1052: @*/
1053: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1054: {
1055:   PetscFunctionBegin;
1057:   PetscSectionCheckValidField(field, s->numFields);
1058:   PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
1059:   PetscFunctionReturn(PETSC_SUCCESS);
1060: }

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

1065:   Not Collective

1067:   Input Parameters:
1068: + s      - the `PetscSection`
1069: . point  - the point
1070: . field  - the field
1071: - numDof - the number of dof

1073:   Level: intermediate

1075:   Notes:
1076:   When adding to the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1077:   the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`

1079:   This is equivalent to
1080: .vb
1081:      PetscSection fs;
1082:      PetscSectionGetField(s,field,&fs)
1083:      PetscSectionAddDof(fs,numDof)
1084: .ve

1086: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1087: @*/
1088: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1089: {
1090:   PetscFunctionBegin;
1092:   PetscSectionCheckValidField(field, s->numFields);
1093:   PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1094:   PetscFunctionReturn(PETSC_SUCCESS);
1095: }

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

1100:   Not Collective

1102:   Input Parameters:
1103: + s     - the `PetscSection`
1104: - point - the point

1106:   Output Parameter:
1107: . numDof - the number of dof which are fixed by constraints

1109:   Level: intermediate

1111: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1112: @*/
1113: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1114: {
1115:   PetscFunctionBegin;
1117:   PetscAssertPointer(numDof, 3);
1118:   if (s->bc) {
1119:     PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1120:   } else *numDof = 0;
1121:   PetscFunctionReturn(PETSC_SUCCESS);
1122: }

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

1127:   Not Collective

1129:   Input Parameters:
1130: + s      - the `PetscSection`
1131: . point  - the point
1132: - numDof - the number of dof which are fixed by constraints

1134:   Level: intermediate

1136: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1137: @*/
1138: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1139: {
1140:   PetscFunctionBegin;
1142:   if (numDof) {
1143:     PetscCall(PetscSectionCheckConstraints_Private(s));
1144:     PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1145:   }
1146:   PetscFunctionReturn(PETSC_SUCCESS);
1147: }

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

1152:   Not Collective

1154:   Input Parameters:
1155: + s      - the `PetscSection`
1156: . point  - the point
1157: - numDof - the number of additional dof which are fixed by constraints

1159:   Level: intermediate

1161: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1162: @*/
1163: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1164: {
1165:   PetscFunctionBegin;
1167:   if (numDof) {
1168:     PetscCall(PetscSectionCheckConstraints_Private(s));
1169:     PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1170:   }
1171:   PetscFunctionReturn(PETSC_SUCCESS);
1172: }

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

1177:   Not Collective

1179:   Input Parameters:
1180: + s     - the `PetscSection`
1181: . point - the point
1182: - field - the field

1184:   Output Parameter:
1185: . numDof - the number of dof which are fixed by constraints

1187:   Level: intermediate

1189: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1190: @*/
1191: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1192: {
1193:   PetscFunctionBegin;
1195:   PetscAssertPointer(numDof, 4);
1196:   PetscSectionCheckValidField(field, s->numFields);
1197:   PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1198:   PetscFunctionReturn(PETSC_SUCCESS);
1199: }

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

1204:   Not Collective

1206:   Input Parameters:
1207: + s      - the `PetscSection`
1208: . point  - the point
1209: . field  - the field
1210: - numDof - the number of dof which are fixed by constraints

1212:   Level: intermediate

1214: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1215: @*/
1216: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1217: {
1218:   PetscFunctionBegin;
1220:   PetscSectionCheckValidField(field, s->numFields);
1221:   PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1222:   PetscFunctionReturn(PETSC_SUCCESS);
1223: }

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

1228:   Not Collective

1230:   Input Parameters:
1231: + s      - the `PetscSection`
1232: . point  - the point
1233: . field  - the field
1234: - numDof - the number of additional dof which are fixed by constraints

1236:   Level: intermediate

1238: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1239: @*/
1240: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1241: {
1242:   PetscFunctionBegin;
1244:   PetscSectionCheckValidField(field, s->numFields);
1245:   PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1246:   PetscFunctionReturn(PETSC_SUCCESS);
1247: }

1249: /*@
1250:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1252:   Not Collective

1254:   Input Parameter:
1255: . s - the `PetscSection`

1257:   Level: advanced

1259: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1260: @*/
1261: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1262: {
1263:   PetscFunctionBegin;
1265:   if (s->bc) {
1266:     const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;

1268:     PetscCall(PetscSectionSetUp(s->bc));
1269:     PetscCall(PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices));
1270:   }
1271:   PetscFunctionReturn(PETSC_SUCCESS);
1272: }

1274: /*@
1275:   PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point in preparation for use of the `PetscSection`

1277:   Not Collective

1279:   Input Parameter:
1280: . s - the `PetscSection`

1282:   Level: intermediate

1284:   Notes:
1285:   If used, `PetscSectionSetPermutation()` must be called before this routine.

1287:   `PetscSectionSetPointMajor()`, cannot be called after this routine.

1289: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1290: @*/
1291: PetscErrorCode PetscSectionSetUp(PetscSection s)
1292: {
1293:   const PetscInt *pind   = NULL;
1294:   PetscInt        offset = 0, foff, p, f;

1296:   PetscFunctionBegin;
1298:   if (s->setup) PetscFunctionReturn(PETSC_SUCCESS);
1299:   s->setup = PETSC_TRUE;
1300:   /* Set offsets and field offsets for all points */
1301:   /*   Assume that all fields have the same chart */
1302:   PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1303:   if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1304:   if (s->pointMajor) {
1305:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1306:       const PetscInt q = pind ? pind[p] : p;

1308:       /* Set point offset */
1309:       s->atlasOff[q] = offset;
1310:       offset += s->atlasDof[q];
1311:       /* Set field offset */
1312:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1313:         PetscSection sf = s->field[f];

1315:         sf->atlasOff[q] = foff;
1316:         foff += sf->atlasDof[q];
1317:       }
1318:     }
1319:   } else {
1320:     /* Set field offsets for all points */
1321:     for (f = 0; f < s->numFields; ++f) {
1322:       PetscSection sf = s->field[f];

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

1327:         sf->atlasOff[q] = offset;
1328:         offset += sf->atlasDof[q];
1329:       }
1330:     }
1331:     /* Disable point offsets since these are unused */
1332:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1333:   }
1334:   if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1335:   /* Setup BC sections */
1336:   PetscCall(PetscSectionSetUpBC(s));
1337:   for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1338:   PetscFunctionReturn(PETSC_SUCCESS);
1339: }

1341: /*@
1342:   PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the `PetscSection`

1344:   Not Collective

1346:   Input Parameter:
1347: . s - the `PetscSection`

1349:   Output Parameter:
1350: . maxDof - the maximum dof

1352:   Level: intermediate

1354:   Notes:
1355:   The returned number is up-to-date without need for `PetscSectionSetUp()`.

1357:   This is the maximum over all points of the sum of the number of dof in the unnamed default field plus all named fields. This is equivalent to
1358:   the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process

1360:   Developer Notes:
1361:   The returned number is calculated lazily and stashed.

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

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

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

1369: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1370: @*/
1371: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1372: {
1373:   PetscInt p;

1375:   PetscFunctionBegin;
1377:   PetscAssertPointer(maxDof, 2);
1378:   if (s->maxDof == PETSC_MIN_INT) {
1379:     s->maxDof = 0;
1380:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1381:   }
1382:   *maxDof = s->maxDof;
1383:   PetscFunctionReturn(PETSC_SUCCESS);
1384: }

1386: /*@
1387:   PetscSectionGetStorageSize - Return the size of an array or local `Vec` capable of holding all the degrees of freedom defined in a `PetscSection`

1389:   Not Collective

1391:   Input Parameter:
1392: . s - the `PetscSection`

1394:   Output Parameter:
1395: . size - the size of an array which can hold all the dofs

1397:   Level: intermediate

1399: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1400: @*/
1401: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1402: {
1403:   PetscInt p, n = 0;

1405:   PetscFunctionBegin;
1407:   PetscAssertPointer(size, 2);
1408:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1409:   *size = n;
1410:   PetscFunctionReturn(PETSC_SUCCESS);
1411: }

1413: /*@
1414:   PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom in a `PetscSection`

1416:   Not Collective

1418:   Input Parameter:
1419: . s - the `PetscSection`

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

1424:   Level: intermediate

1426: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1427: @*/
1428: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1429: {
1430:   PetscInt p, n = 0;

1432:   PetscFunctionBegin;
1434:   PetscAssertPointer(size, 2);
1435:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1436:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1437:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1438:   }
1439:   *size = n;
1440:   PetscFunctionReturn(PETSC_SUCCESS);
1441: }

1443: /*@
1444:   PetscSectionCreateGlobalSection - Create a parallel section describing the global layout using
1445:   a local (sequential) `PetscSection` on each MPI process and a `PetscSF` describing the section point overlap.

1447:   Input Parameters:
1448: + s                  - The `PetscSection` for the local field layout
1449: . sf                 - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1450: . usePermutation     - By default this is `PETSC_TRUE`, meaning any permutation of the local section is transferred to the global section
1451: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1452: - localOffsets       - If `PETSC_TRUE`, use local rather than global offsets for the points

1454:   Output Parameter:
1455: . gsection - The `PetscSection` for the global field layout

1457:   Level: intermediate

1459:   Notes:
1460:   On each MPI process `gsection` inherits the chart of the `s` on that process.

1462:   This sets negative sizes and offsets to points not owned by this process as defined by `sf` but that are within the local value of the chart of `gsection`.
1463:   In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).

1465: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1466: @*/
1467: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool usePermutation, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1468: {
1469:   PetscSection    gs;
1470:   const PetscInt *pind = NULL;
1471:   PetscInt       *recv = NULL, *neg = NULL;
1472:   PetscInt        pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1473:   PetscInt        numFields, f, numComponents;

1475:   PetscFunctionBegin;
1481:   PetscAssertPointer(gsection, 6);
1482:   PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1483:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1484:   PetscCall(PetscSectionGetNumFields(s, &numFields));
1485:   if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1486:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1487:   PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1488:   gs->includesConstraints = includeConstraints;
1489:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1490:   nlocal = nroots; /* The local/leaf space matches global/root space */
1491:   /* Must allocate for all points visible to SF, which may be more than this section */
1492:   if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1493:     PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1494:     PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1495:     PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1496:     PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1497:     PetscCall(PetscArrayzero(neg, nroots));
1498:   }
1499:   /* Mark all local points with negative dof */
1500:   for (p = pStart; p < pEnd; ++p) {
1501:     PetscCall(PetscSectionGetDof(s, p, &dof));
1502:     PetscCall(PetscSectionSetDof(gs, p, dof));
1503:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1504:     if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1505:     if (neg) neg[p] = -(dof + 1);
1506:   }
1507:   PetscCall(PetscSectionSetUpBC(gs));
1508:   if (gs->bcIndices) PetscCall(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]));
1509:   if (nroots >= 0) {
1510:     PetscCall(PetscArrayzero(recv, nlocal));
1511:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1512:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1513:     for (p = pStart; p < pEnd; ++p) {
1514:       if (recv[p] < 0) {
1515:         gs->atlasDof[p - pStart] = recv[p];
1516:         PetscCall(PetscSectionGetDof(s, p, &dof));
1517:         PetscCheck(-(recv[p] + 1) == dof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %" PetscInt_FMT " for point %" PetscInt_FMT " is not the unconstrained %" PetscInt_FMT, -(recv[p] + 1), p, dof);
1518:       }
1519:     }
1520:   }
1521:   /* Calculate new sizes, get process offset, and calculate point offsets */
1522:   if (usePermutation && s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1523:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1524:     const PetscInt q = pind ? pind[p] : p;

1526:     cdof            = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1527:     gs->atlasOff[q] = off;
1528:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1529:   }
1530:   if (!localOffsets) {
1531:     PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1532:     globalOff -= off;
1533:   }
1534:   for (p = pStart, off = 0; p < pEnd; ++p) {
1535:     gs->atlasOff[p - pStart] += globalOff;
1536:     if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1537:   }
1538:   if (usePermutation && s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1539:   /* Put in negative offsets for ghost points */
1540:   if (nroots >= 0) {
1541:     PetscCall(PetscArrayzero(recv, nlocal));
1542:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1543:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1544:     for (p = pStart; p < pEnd; ++p) {
1545:       if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1546:     }
1547:   }
1548:   PetscCall(PetscFree2(neg, recv));
1549:   /* Set field dofs/offsets/constraints */
1550:   for (f = 0; f < numFields; ++f) {
1551:     const char *name;

1553:     gs->field[f]->includesConstraints = includeConstraints;
1554:     PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1555:     PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1556:     PetscCall(PetscSectionGetFieldName(s, f, &name));
1557:     PetscCall(PetscSectionSetFieldName(gs, f, name));
1558:   }
1559:   for (p = pStart; p < pEnd; ++p) {
1560:     PetscCall(PetscSectionGetOffset(gs, p, &off));
1561:     for (f = 0, foff = off; f < numFields; ++f) {
1562:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1563:       if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1564:       PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1565:       PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1566:       PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff));
1567:       PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1568:       foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1569:     }
1570:   }
1571:   for (f = 0; f < numFields; ++f) {
1572:     PetscSection gfs = gs->field[f];

1574:     PetscCall(PetscSectionSetUpBC(gfs));
1575:     if (gfs->bcIndices) PetscCall(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]));
1576:   }
1577:   gs->setup = PETSC_TRUE;
1578:   PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1579:   *gsection = gs;
1580:   PetscFunctionReturn(PETSC_SUCCESS);
1581: }

1583: /*@
1584:   PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1585:   a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.

1587:   Input Parameters:
1588: + s                  - The `PetscSection` for the local field layout
1589: . sf                 - The `PetscSF` describing parallel layout of the section points
1590: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1591: . numExcludes        - The number of exclusion ranges, this must have the same value on all MPI processes
1592: - excludes           - An array [start_0, end_0, start_1, end_1, ...] where there are `numExcludes` pairs and must have the same values on all MPI processes

1594:   Output Parameter:
1595: . gsection - The `PetscSection` for the global field layout

1597:   Level: advanced

1599:   Notes:
1600:   On each MPI process `gsection` inherits the chart of the `s` on that process.

1602:   This sets negative sizes and offsets to points not owned by this process as defined by `sf` but that are within the local value of the chart of `gsection`.
1603:   In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).

1605:   This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`

1607:   Developer Notes:
1608:   This is a terrible function name

1610: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1611: @*/
1612: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1613: {
1614:   const PetscInt *pind = NULL;
1615:   PetscInt       *neg = NULL, *tmpOff = NULL;
1616:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;

1618:   PetscFunctionBegin;
1621:   PetscAssertPointer(gsection, 6);
1622:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1623:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1624:   PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1625:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1626:   if (nroots >= 0) {
1627:     PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1628:     PetscCall(PetscCalloc1(nroots, &neg));
1629:     if (nroots > pEnd - pStart) {
1630:       PetscCall(PetscCalloc1(nroots, &tmpOff));
1631:     } else {
1632:       tmpOff = &(*gsection)->atlasDof[-pStart];
1633:     }
1634:   }
1635:   /* Mark ghost points with negative dof */
1636:   for (p = pStart; p < pEnd; ++p) {
1637:     for (e = 0; e < numExcludes; ++e) {
1638:       if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1639:         PetscCall(PetscSectionSetDof(*gsection, p, 0));
1640:         break;
1641:       }
1642:     }
1643:     if (e < numExcludes) continue;
1644:     PetscCall(PetscSectionGetDof(s, p, &dof));
1645:     PetscCall(PetscSectionSetDof(*gsection, p, dof));
1646:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1647:     if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1648:     if (neg) neg[p] = -(dof + 1);
1649:   }
1650:   PetscCall(PetscSectionSetUpBC(*gsection));
1651:   if (nroots >= 0) {
1652:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1653:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1654:     if (nroots > pEnd - pStart) {
1655:       for (p = pStart; p < pEnd; ++p) {
1656:         if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1657:       }
1658:     }
1659:   }
1660:   /* Calculate new sizes, get process offset, and calculate point offsets */
1661:   if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1662:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1663:     const PetscInt q = pind ? pind[p] : p;

1665:     cdof                     = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1666:     (*gsection)->atlasOff[q] = off;
1667:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1668:   }
1669:   PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1670:   globalOff -= off;
1671:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1672:     (*gsection)->atlasOff[p] += globalOff;
1673:     if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1674:   }
1675:   if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1676:   /* Put in negative offsets for ghost points */
1677:   if (nroots >= 0) {
1678:     if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1679:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1680:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1681:     if (nroots > pEnd - pStart) {
1682:       for (p = pStart; p < pEnd; ++p) {
1683:         if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1684:       }
1685:     }
1686:   }
1687:   if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1688:   PetscCall(PetscFree(neg));
1689:   PetscFunctionReturn(PETSC_SUCCESS);
1690: }

1692: /*@
1693:   PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart

1695:   Collective

1697:   Input Parameters:
1698: + comm - The `MPI_Comm`
1699: - s    - The `PetscSection`

1701:   Output Parameter:
1702: . layout - The point layout for the data that defines the section

1704:   Level: advanced

1706:   Notes:
1707:   `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1708:   degrees of freedom).

1710:   This count includes constrained degrees of freedom

1712:   This is usually called on the default global section.

1714:   Example:
1715: .vb
1716:      The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1717:      The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1718: .ve

1720:   Developer Notes:
1721:   I find the names of these two functions extremely non-informative

1723: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1724: @*/
1725: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1726: {
1727:   PetscInt pStart, pEnd, p, localSize = 0;

1729:   PetscFunctionBegin;
1730:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1731:   for (p = pStart; p < pEnd; ++p) {
1732:     PetscInt dof;

1734:     PetscCall(PetscSectionGetDof(s, p, &dof));
1735:     if (dof >= 0) ++localSize;
1736:   }
1737:   PetscCall(PetscLayoutCreate(comm, layout));
1738:   PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1739:   PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1740:   PetscCall(PetscLayoutSetUp(*layout));
1741:   PetscFunctionReturn(PETSC_SUCCESS);
1742: }

1744: /*@
1745:   PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`

1747:   Collective

1749:   Input Parameters:
1750: + comm - The `MPI_Comm`
1751: - s    - The `PetscSection`

1753:   Output Parameter:
1754: . layout - The dof layout for the section

1756:   Level: advanced

1758:   Notes:
1759:   `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1760:   including the constrained degrees of freedom

1762:   This is usually called for the default global section.

1764:   Example:
1765: .vb
1766:      The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1767:      The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1768: .ve

1770: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1771: @*/
1772: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1773: {
1774:   PetscInt pStart, pEnd, p, localSize = 0;

1776:   PetscFunctionBegin;
1778:   PetscAssertPointer(layout, 3);
1779:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1780:   for (p = pStart; p < pEnd; ++p) {
1781:     PetscInt dof, cdof;

1783:     PetscCall(PetscSectionGetDof(s, p, &dof));
1784:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1785:     if (dof - cdof > 0) localSize += dof - cdof;
1786:   }
1787:   PetscCall(PetscLayoutCreate(comm, layout));
1788:   PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1789:   PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1790:   PetscCall(PetscLayoutSetUp(*layout));
1791:   PetscFunctionReturn(PETSC_SUCCESS);
1792: }

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

1797:   Not Collective

1799:   Input Parameters:
1800: + s     - the `PetscSection`
1801: - point - the point

1803:   Output Parameter:
1804: . offset - the offset

1806:   Level: intermediate

1808:   Notes:
1809:   In a global section, `offset` will be negative for points not owned by this process.

1811:   This is for the unnamed default field in the `PetscSection` not the named fields

1813:   The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`

1815: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1816: @*/
1817: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1818: {
1819:   PetscFunctionBegin;
1821:   PetscAssertPointer(offset, 3);
1822:   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);
1823:   *offset = s->atlasOff[point - s->pStart];
1824:   PetscFunctionReturn(PETSC_SUCCESS);
1825: }

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

1830:   Not Collective

1832:   Input Parameters:
1833: + s      - the `PetscSection`
1834: . point  - the point
1835: - offset - the offset, these values may be negative indicating the values are off process

1837:   Level: developer

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

1842: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1843: @*/
1844: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1845: {
1846:   PetscFunctionBegin;
1848:   PetscCheck(!(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);
1849:   s->atlasOff[point - s->pStart] = offset;
1850:   PetscFunctionReturn(PETSC_SUCCESS);
1851: }

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

1856:   Not Collective

1858:   Input Parameters:
1859: + s     - the `PetscSection`
1860: . point - the point
1861: - field - the field

1863:   Output Parameter:
1864: . offset - the offset

1866:   Level: intermediate

1868:   Notes:
1869:   In a global section, `offset` will be negative for points not owned by this process.

1871:   The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`

1873: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1874: @*/
1875: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1876: {
1877:   PetscFunctionBegin;
1879:   PetscAssertPointer(offset, 4);
1880:   PetscSectionCheckValidField(field, s->numFields);
1881:   PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1882:   PetscFunctionReturn(PETSC_SUCCESS);
1883: }

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

1888:   Not Collective

1890:   Input Parameters:
1891: + s      - the `PetscSection`
1892: . point  - the point
1893: . field  - the field
1894: - offset - the offset, these values may be negative indicating the values are off process

1896:   Level: developer

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

1901: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1902: @*/
1903: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1904: {
1905:   PetscFunctionBegin;
1907:   PetscSectionCheckValidField(field, s->numFields);
1908:   PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1909:   PetscFunctionReturn(PETSC_SUCCESS);
1910: }

1912: /*@
1913:   PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1914:   unnamed default field's first dof

1916:   Not Collective

1918:   Input Parameters:
1919: + s     - the `PetscSection`
1920: . point - the point
1921: - field - the field

1923:   Output Parameter:
1924: . offset - the offset

1926:   Level: advanced

1928:   Note:
1929:   This ignores constraints

1931:   Example:
1932: .vb
1933:   if PetscSectionSetPointMajor(s,PETSC_TRUE)
1934:   The unnamed default field has 3 dof at `point`
1935:   Field 0 has 2 dof at `point`
1936:   Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1937: .ve

1939: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1940: @*/
1941: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1942: {
1943:   PetscInt off, foff;

1945:   PetscFunctionBegin;
1947:   PetscAssertPointer(offset, 4);
1948:   PetscSectionCheckValidField(field, s->numFields);
1949:   PetscCall(PetscSectionGetOffset(s, point, &off));
1950:   PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1951:   *offset = foff - off;
1952:   PetscFunctionReturn(PETSC_SUCCESS);
1953: }

1955: /*@
1956:   PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`

1958:   Not Collective

1960:   Input Parameter:
1961: . s - the `PetscSection`

1963:   Output Parameters:
1964: + start - the minimum offset
1965: - end   - one more than the maximum offset

1967:   Level: intermediate

1969: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1970: @*/
1971: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1972: {
1973:   PetscInt os = 0, oe = 0, pStart, pEnd, p;

1975:   PetscFunctionBegin;
1977:   if (s->atlasOff) {
1978:     os = s->atlasOff[0];
1979:     oe = s->atlasOff[0];
1980:   }
1981:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1982:   for (p = 0; p < pEnd - pStart; ++p) {
1983:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1985:     if (off >= 0) {
1986:       os = PetscMin(os, off);
1987:       oe = PetscMax(oe, off + dof);
1988:     }
1989:   }
1990:   if (start) *start = os;
1991:   if (end) *end = oe;
1992:   PetscFunctionReturn(PETSC_SUCCESS);
1993: }

1995: /*@
1996:   PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields

1998:   Collective

2000:   Input Parameters:
2001: + s      - the `PetscSection`
2002: . len    - the number of subfields
2003: - fields - the subfield numbers

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

2008:   Level: advanced

2010:   Notes:
2011:   The chart of `subs` is the same as the chart of `s`

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

2015: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2016: @*/
2017: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
2018: {
2019:   PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;

2021:   PetscFunctionBegin;
2022:   if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2024:   PetscAssertPointer(fields, 3);
2025:   PetscAssertPointer(subs, 4);
2026:   PetscCall(PetscSectionGetNumFields(s, &nF));
2027:   PetscCheck(len <= nF, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of fields %" PetscInt_FMT, len, nF);
2028:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2029:   PetscCall(PetscSectionSetNumFields(*subs, len));
2030:   for (f = 0; f < len; ++f) {
2031:     const char     *name    = NULL;
2032:     PetscInt        numComp = 0;
2033:     PetscSectionSym sym;

2035:     PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
2036:     PetscCall(PetscSectionSetFieldName(*subs, f, name));
2037:     PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
2038:     PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2039:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
2040:       PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
2041:       PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2042:     }
2043:     PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
2044:     PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
2045:   }
2046:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2047:   PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2048:   for (p = pStart; p < pEnd; ++p) {
2049:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

2051:     for (f = 0; f < len; ++f) {
2052:       PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2053:       PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
2054:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2055:       if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
2056:       dof += fdof;
2057:       cdof += cfdof;
2058:     }
2059:     PetscCall(PetscSectionSetDof(*subs, p, dof));
2060:     if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
2061:     maxCdof = PetscMax(cdof, maxCdof);
2062:   }
2063:   PetscCall(PetscSectionSetUp(*subs));
2064:   if (maxCdof) {
2065:     PetscInt *indices;

2067:     PetscCall(PetscMalloc1(maxCdof, &indices));
2068:     for (p = pStart; p < pEnd; ++p) {
2069:       PetscInt cdof;

2071:       PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
2072:       if (cdof) {
2073:         const PetscInt *oldIndices = NULL;
2074:         PetscInt        fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

2076:         for (f = 0; f < len; ++f) {
2077:           PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2078:           PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2079:           PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
2080:           PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
2081:           for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
2082:           numConst += cfdof;
2083:           fOff += fdof;
2084:         }
2085:         PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2086:       }
2087:     }
2088:     PetscCall(PetscFree(indices));
2089:   }
2090:   PetscFunctionReturn(PETSC_SUCCESS);
2091: }

2093: /*@
2094:   PetscSectionCreateComponentSubsection - Create a new, smaller `PetscSection` composed of only selected components

2096:   Collective

2098:   Input Parameters:
2099: + s     - the `PetscSection`
2100: . len   - the number of components
2101: - comps - the component numbers

2103:   Output Parameter:
2104: . subs - the subsection

2106:   Level: advanced

2108:   Notes:
2109:   The chart of `subs` is the same as the chart of `s`

2111:   This will error if the section has more than one field, or if a component number is out of range

2113: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2114: @*/
2115: PetscErrorCode PetscSectionCreateComponentSubsection(PetscSection s, PetscInt len, const PetscInt comps[], PetscSection *subs)
2116: {
2117:   PetscSectionSym sym;
2118:   const char     *name = NULL;
2119:   PetscInt        Nf, pStart, pEnd;

2121:   PetscFunctionBegin;
2122:   if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2124:   PetscAssertPointer(comps, 3);
2125:   PetscAssertPointer(subs, 4);
2126:   PetscCall(PetscSectionGetNumFields(s, &Nf));
2127:   PetscCheck(Nf == 1, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "This method can only handle one field, not %" PetscInt_FMT, Nf);
2128:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2129:   PetscCall(PetscSectionSetNumFields(*subs, 1));
2130:   PetscCall(PetscSectionGetFieldName(s, 0, &name));
2131:   PetscCall(PetscSectionSetFieldName(*subs, 0, name));
2132:   PetscCall(PetscSectionSetFieldComponents(*subs, 0, len));
2133:   PetscCall(PetscSectionGetFieldSym(s, 0, &sym));
2134:   PetscCall(PetscSectionSetFieldSym(*subs, 0, sym));
2135:   for (PetscInt c = 0; c < len; ++c) {
2136:     PetscCall(PetscSectionGetComponentName(s, 0, comps[c], &name));
2137:     PetscCall(PetscSectionSetComponentName(*subs, 0, c, name));
2138:   }
2139:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2140:   PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2141:   for (PetscInt p = pStart; p < pEnd; ++p) {
2142:     PetscInt dof, cdof, cfdof;

2144:     PetscCall(PetscSectionGetDof(s, p, &dof));
2145:     if (!dof) continue;
2146:     PetscCall(PetscSectionGetFieldConstraintDof(s, p, 0, &cfdof));
2147:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2148:     PetscCheck(!cdof && !cfdof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Component selection does not work with constraints");
2149:     PetscCall(PetscSectionSetFieldDof(*subs, p, 0, len));
2150:     PetscCall(PetscSectionSetDof(*subs, p, len));
2151:   }
2152:   PetscCall(PetscSectionSetUp(*subs));
2153:   PetscFunctionReturn(PETSC_SUCCESS);
2154: }

2156: /*@
2157:   PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s

2159:   Collective

2161:   Input Parameters:
2162: + s   - the input sections
2163: - len - the number of input sections

2165:   Output Parameter:
2166: . supers - the supersection

2168:   Level: advanced

2170:   Notes:
2171:   The section offsets now refer to a new, larger vector.

2173:   Developer Notes:
2174:   Needs to explain how the sections are composed

2176: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2177: @*/
2178: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2179: {
2180:   PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

2182:   PetscFunctionBegin;
2183:   if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2184:   for (i = 0; i < len; ++i) {
2185:     PetscInt nf, pStarti, pEndi;

2187:     PetscCall(PetscSectionGetNumFields(s[i], &nf));
2188:     PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2189:     pStart = PetscMin(pStart, pStarti);
2190:     pEnd   = PetscMax(pEnd, pEndi);
2191:     Nf += nf;
2192:   }
2193:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2194:   PetscCall(PetscSectionSetNumFields(*supers, Nf));
2195:   for (i = 0, f = 0; i < len; ++i) {
2196:     PetscInt nf, fi, ci;

2198:     PetscCall(PetscSectionGetNumFields(s[i], &nf));
2199:     for (fi = 0; fi < nf; ++fi, ++f) {
2200:       const char *name    = NULL;
2201:       PetscInt    numComp = 0;

2203:       PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2204:       PetscCall(PetscSectionSetFieldName(*supers, f, name));
2205:       PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2206:       PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2207:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2208:         PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2209:         PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2210:       }
2211:     }
2212:   }
2213:   PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2214:   for (p = pStart; p < pEnd; ++p) {
2215:     PetscInt dof = 0, cdof = 0;

2217:     for (i = 0, f = 0; i < len; ++i) {
2218:       PetscInt nf, fi, pStarti, pEndi;
2219:       PetscInt fdof = 0, cfdof = 0;

2221:       PetscCall(PetscSectionGetNumFields(s[i], &nf));
2222:       PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2223:       if ((p < pStarti) || (p >= pEndi)) continue;
2224:       for (fi = 0; fi < nf; ++fi, ++f) {
2225:         PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2226:         PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2227:         PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2228:         if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2229:         dof += fdof;
2230:         cdof += cfdof;
2231:       }
2232:     }
2233:     PetscCall(PetscSectionSetDof(*supers, p, dof));
2234:     if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2235:     maxCdof = PetscMax(cdof, maxCdof);
2236:   }
2237:   PetscCall(PetscSectionSetUp(*supers));
2238:   if (maxCdof) {
2239:     PetscInt *indices;

2241:     PetscCall(PetscMalloc1(maxCdof, &indices));
2242:     for (p = pStart; p < pEnd; ++p) {
2243:       PetscInt cdof;

2245:       PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2246:       if (cdof) {
2247:         PetscInt dof, numConst = 0, fOff = 0;

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

2253:           PetscCall(PetscSectionGetNumFields(s[i], &nf));
2254:           PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2255:           if ((p < pStarti) || (p >= pEndi)) continue;
2256:           for (fi = 0; fi < nf; ++fi, ++f) {
2257:             PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2258:             PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2259:             PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2260:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2261:             PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2262:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2263:             numConst += cfdof;
2264:           }
2265:           PetscCall(PetscSectionGetDof(s[i], p, &dof));
2266:           fOff += dof;
2267:         }
2268:         PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2269:       }
2270:     }
2271:     PetscCall(PetscFree(indices));
2272:   }
2273:   PetscFunctionReturn(PETSC_SUCCESS);
2274: }

2276: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
2277: {
2278:   const PetscInt *points = NULL, *indices = NULL;
2279:   PetscInt        numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;

2281:   PetscFunctionBegin;
2284:   PetscAssertPointer(subs, 4);
2285:   PetscCall(PetscSectionGetNumFields(s, &numFields));
2286:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2287:   if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2288:   for (f = 0; f < numFields; ++f) {
2289:     const char *name    = NULL;
2290:     PetscInt    numComp = 0;

2292:     PetscCall(PetscSectionGetFieldName(s, f, &name));
2293:     PetscCall(PetscSectionSetFieldName(*subs, f, name));
2294:     PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2295:     PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2296:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2297:       PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2298:       PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2299:     }
2300:   }
2301:   /* For right now, we do not try to squeeze the subchart */
2302:   if (subpointMap) {
2303:     PetscCall(ISGetSize(subpointMap, &numSubpoints));
2304:     PetscCall(ISGetIndices(subpointMap, &points));
2305:   }
2306:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2307:   if (renumberPoints) {
2308:     spStart = 0;
2309:     spEnd   = numSubpoints;
2310:   } else {
2311:     PetscCall(ISGetMinMax(subpointMap, &spStart, &spEnd));
2312:     ++spEnd;
2313:   }
2314:   PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2315:   for (p = pStart; p < pEnd; ++p) {
2316:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

2318:     PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2319:     if (subp < 0) continue;
2320:     if (!renumberPoints) subp = p;
2321:     for (f = 0; f < numFields; ++f) {
2322:       PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2323:       PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2324:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2325:       if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2326:     }
2327:     PetscCall(PetscSectionGetDof(s, p, &dof));
2328:     PetscCall(PetscSectionSetDof(*subs, subp, dof));
2329:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2330:     if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2331:   }
2332:   PetscCall(PetscSectionSetUp(*subs));
2333:   /* Change offsets to original offsets */
2334:   for (p = pStart; p < pEnd; ++p) {
2335:     PetscInt off, foff = 0;

2337:     PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2338:     if (subp < 0) continue;
2339:     if (!renumberPoints) subp = p;
2340:     for (f = 0; f < numFields; ++f) {
2341:       PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2342:       PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2343:     }
2344:     PetscCall(PetscSectionGetOffset(s, p, &off));
2345:     PetscCall(PetscSectionSetOffset(*subs, subp, off));
2346:   }
2347:   /* Copy constraint indices */
2348:   for (subp = spStart; subp < spEnd; ++subp) {
2349:     PetscInt cdof;

2351:     PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2352:     if (cdof) {
2353:       for (f = 0; f < numFields; ++f) {
2354:         PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2355:         PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2356:       }
2357:       PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2358:       PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2359:     }
2360:   }
2361:   if (subpointMap) PetscCall(ISRestoreIndices(subpointMap, &points));
2362:   PetscFunctionReturn(PETSC_SUCCESS);
2363: }

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

2368:   Collective

2370:   Input Parameters:
2371: + s           - the `PetscSection`
2372: - subpointMap - a sorted list of points in the original mesh which are in the submesh

2374:   Output Parameter:
2375: . subs - the subsection

2377:   Level: advanced

2379:   Notes:
2380:   The points are renumbered from 0, and the section offsets now refer to a new, smaller vector. That is the chart of `subs` is `[0,sizeof(subpointmap))`

2382:   Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before

2384:   Developer Notes:
2385:   The use of the term Submesh is confusing and needs clarification, it is not specific to meshes. It appears to be just a subset of the chart of the original `PetscSection`

2387: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2388: @*/
2389: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2390: {
2391:   PetscFunctionBegin;
2392:   PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_TRUE, subs));
2393:   PetscFunctionReturn(PETSC_SUCCESS);
2394: }

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

2399:   Collective

2401:   Input Parameters:
2402: + s           - the `PetscSection`
2403: - subpointMap - a sorted list of points in the original mesh which are in the subdomain

2405:   Output Parameter:
2406: . subs - the subsection

2408:   Level: advanced

2410:   Notes:
2411:   The point numbers remain the same as in the larger `PetscSection`, but the section offsets now refer to a new, smaller vector. The chart of `subs`
2412:   is `[min(subpointMap),max(subpointMap)+1)`

2414:   Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero

2416:   Developer Notes:
2417:   The use of the term Subdomain is unneeded and needs clarification, it is not specific to meshes. It appears to be just a subset of the chart of the original `PetscSection`

2419: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2420: @*/
2421: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2422: {
2423:   PetscFunctionBegin;
2424:   PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2425:   PetscFunctionReturn(PETSC_SUCCESS);
2426: }

2428: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2429: {
2430:   PetscInt    p;
2431:   PetscMPIInt rank;

2433:   PetscFunctionBegin;
2434:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2435:   PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2436:   PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2437:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2438:     if (s->bc && s->bc->atlasDof[p] > 0) {
2439:       PetscInt b;
2440:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2441:       if (s->bcIndices) {
2442:         for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2443:       }
2444:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2445:     } else {
2446:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2447:     }
2448:   }
2449:   PetscCall(PetscViewerFlush(viewer));
2450:   PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2451:   if (s->sym) {
2452:     PetscCall(PetscViewerASCIIPushTab(viewer));
2453:     PetscCall(PetscSectionSymView(s->sym, viewer));
2454:     PetscCall(PetscViewerASCIIPopTab(viewer));
2455:   }
2456:   PetscFunctionReturn(PETSC_SUCCESS);
2457: }

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

2462:   Collective

2464:   Input Parameters:
2465: + A    - the `PetscSection` object to view
2466: . obj  - Optional object that provides the options prefix used for the options
2467: - name - command line option

2469:   Level: intermediate

2471:   Note:
2472:   See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`

2474: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2475: @*/
2476: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2477: {
2478:   PetscFunctionBegin;
2480:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2481:   PetscFunctionReturn(PETSC_SUCCESS);
2482: }

2484: /*@C
2485:   PetscSectionView - Views a `PetscSection`

2487:   Collective

2489:   Input Parameters:
2490: + s      - the `PetscSection` object to view
2491: - viewer - the viewer

2493:   Level: beginner

2495:   Note:
2496:   `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2497:   distribution independent data, such as dofs, offsets, constraint dofs,
2498:   and constraint indices. Points that have negative dofs, for instance,
2499:   are not saved as they represent points owned by other processes.
2500:   Point numbering and rank assignment is currently not stored.
2501:   The saved section can be loaded with `PetscSectionLoad()`.

2503: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2504: @*/
2505: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2506: {
2507:   PetscBool isascii, ishdf5;
2508:   PetscInt  f;

2510:   PetscFunctionBegin;
2512:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2514:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2515:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2516:   if (isascii) {
2517:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2518:     if (s->numFields) {
2519:       PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2520:       for (f = 0; f < s->numFields; ++f) {
2521:         PetscCall(PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " \"%s\" with %" PetscInt_FMT " components\n", f, s->fieldNames[f], s->numFieldComponents[f]));
2522:         PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2523:       }
2524:     } else {
2525:       PetscCall(PetscSectionView_ASCII(s, viewer));
2526:     }
2527:   } else if (ishdf5) {
2528: #if PetscDefined(HAVE_HDF5)
2529:     PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2530: #else
2531:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2532: #endif
2533:   }
2534:   PetscFunctionReturn(PETSC_SUCCESS);
2535: }

2537: /*@C
2538:   PetscSectionLoad - Loads a `PetscSection`

2540:   Collective

2542:   Input Parameters:
2543: + s      - the `PetscSection` object to load
2544: - viewer - the viewer

2546:   Level: beginner

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

2556: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2557: @*/
2558: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2559: {
2560:   PetscBool ishdf5;

2562:   PetscFunctionBegin;
2565:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2566:   if (ishdf5) {
2567: #if PetscDefined(HAVE_HDF5)
2568:     PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2569:     PetscFunctionReturn(PETSC_SUCCESS);
2570: #else
2571:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2572: #endif
2573:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2574: }

2576: /*@
2577:   PetscSectionResetClosurePermutation - Remove any existing closure permutation

2579:   Input Parameter:
2580: . section - The `PetscSection`

2582:   Level: intermediate

2584: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2585: @*/
2586: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2587: {
2588:   PetscSectionClosurePermVal clVal;

2590:   PetscFunctionBegin;
2591:   if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2592:   kh_foreach_value(section->clHash, clVal, {
2593:     PetscCall(PetscFree(clVal.perm));
2594:     PetscCall(PetscFree(clVal.invPerm));
2595:   });
2596:   kh_destroy(ClPerm, section->clHash);
2597:   section->clHash = NULL;
2598:   PetscFunctionReturn(PETSC_SUCCESS);
2599: }

2601: /*@
2602:   PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.

2604:   Not Collective

2606:   Input Parameter:
2607: . s - the `PetscSection`

2609:   Level: beginner

2611: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2612: @*/
2613: PetscErrorCode PetscSectionReset(PetscSection s)
2614: {
2615:   PetscInt f, c;

2617:   PetscFunctionBegin;
2619:   for (f = 0; f < s->numFields; ++f) {
2620:     PetscCall(PetscSectionDestroy(&s->field[f]));
2621:     PetscCall(PetscFree(s->fieldNames[f]));
2622:     for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2623:     PetscCall(PetscFree(s->compNames[f]));
2624:   }
2625:   PetscCall(PetscFree(s->numFieldComponents));
2626:   PetscCall(PetscFree(s->fieldNames));
2627:   PetscCall(PetscFree(s->compNames));
2628:   PetscCall(PetscFree(s->field));
2629:   PetscCall(PetscSectionDestroy(&s->bc));
2630:   PetscCall(PetscFree(s->bcIndices));
2631:   PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2632:   PetscCall(PetscSectionDestroy(&s->clSection));
2633:   PetscCall(ISDestroy(&s->clPoints));
2634:   PetscCall(ISDestroy(&s->perm));
2635:   PetscCall(PetscBTDestroy(&s->blockStarts));
2636:   PetscCall(PetscSectionResetClosurePermutation(s));
2637:   PetscCall(PetscSectionSymDestroy(&s->sym));
2638:   PetscCall(PetscSectionDestroy(&s->clSection));
2639:   PetscCall(ISDestroy(&s->clPoints));
2640:   PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2641:   s->pStart    = -1;
2642:   s->pEnd      = -1;
2643:   s->maxDof    = 0;
2644:   s->setup     = PETSC_FALSE;
2645:   s->numFields = 0;
2646:   s->clObj     = NULL;
2647:   PetscFunctionReturn(PETSC_SUCCESS);
2648: }

2650: /*@
2651:   PetscSectionDestroy - Frees a `PetscSection`

2653:   Not Collective

2655:   Input Parameter:
2656: . s - the `PetscSection`

2658:   Level: beginner

2660: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2661: @*/
2662: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2663: {
2664:   PetscFunctionBegin;
2665:   if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2667:   if (--((PetscObject)*s)->refct > 0) {
2668:     *s = NULL;
2669:     PetscFunctionReturn(PETSC_SUCCESS);
2670:   }
2671:   PetscCall(PetscSectionReset(*s));
2672:   PetscCall(PetscHeaderDestroy(s));
2673:   PetscFunctionReturn(PETSC_SUCCESS);
2674: }

2676: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2677: {
2678:   const PetscInt p = point - s->pStart;

2680:   PetscFunctionBegin;
2682:   *values = &baseArray[s->atlasOff[p]];
2683:   PetscFunctionReturn(PETSC_SUCCESS);
2684: }

2686: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2687: {
2688:   PetscInt      *array;
2689:   const PetscInt p           = point - s->pStart;
2690:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2691:   PetscInt       cDim        = 0;

2693:   PetscFunctionBegin;
2695:   PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2696:   array = &baseArray[s->atlasOff[p]];
2697:   if (!cDim) {
2698:     if (orientation >= 0) {
2699:       const PetscInt dim = s->atlasDof[p];
2700:       PetscInt       i;

2702:       if (mode == INSERT_VALUES) {
2703:         for (i = 0; i < dim; ++i) array[i] = values ? values[i] : i;
2704:       } else {
2705:         for (i = 0; i < dim; ++i) array[i] += values[i];
2706:       }
2707:     } else {
2708:       PetscInt offset = 0;
2709:       PetscInt j      = -1, field, i;

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

2714:         for (i = dim - 1; i >= 0; --i) array[++j] = values ? values[i + offset] : i + offset;
2715:         offset += dim;
2716:       }
2717:     }
2718:   } else {
2719:     if (orientation >= 0) {
2720:       const PetscInt  dim  = s->atlasDof[p];
2721:       PetscInt        cInd = 0, i;
2722:       const PetscInt *cDof;

2724:       PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2725:       if (mode == INSERT_VALUES) {
2726:         for (i = 0; i < dim; ++i) {
2727:           if ((cInd < cDim) && (i == cDof[cInd])) {
2728:             ++cInd;
2729:             continue;
2730:           }
2731:           array[i] = values ? values[i] : i;
2732:         }
2733:       } else {
2734:         for (i = 0; i < dim; ++i) {
2735:           if ((cInd < cDim) && (i == cDof[cInd])) {
2736:             ++cInd;
2737:             continue;
2738:           }
2739:           array[i] += values[i];
2740:         }
2741:       }
2742:     } else {
2743:       const PetscInt *cDof;
2744:       PetscInt        offset  = 0;
2745:       PetscInt        cOffset = 0;
2746:       PetscInt        j       = 0, field;

2748:       PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2749:       for (field = 0; field < s->numFields; ++field) {
2750:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2751:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2752:         const PetscInt sDim = dim - tDim;
2753:         PetscInt       cInd = 0, i, k;

2755:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2756:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2757:             ++cInd;
2758:             continue;
2759:           }
2760:           array[j] = values ? values[k] : k;
2761:         }
2762:         offset += dim;
2763:         cOffset += dim - tDim;
2764:       }
2765:     }
2766:   }
2767:   PetscFunctionReturn(PETSC_SUCCESS);
2768: }

2770: /*@C
2771:   PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs

2773:   Not Collective

2775:   Input Parameter:
2776: . s - The `PetscSection`

2778:   Output Parameter:
2779: . hasConstraints - flag indicating that the section has constrained dofs

2781:   Level: intermediate

2783: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2784: @*/
2785: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2786: {
2787:   PetscFunctionBegin;
2789:   PetscAssertPointer(hasConstraints, 2);
2790:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2791:   PetscFunctionReturn(PETSC_SUCCESS);
2792: }

2794: /*@C
2795:   PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point

2797:   Not Collective

2799:   Input Parameters:
2800: + s     - The `PetscSection`
2801: - point - The point

2803:   Output Parameter:
2804: . indices - The constrained dofs

2806:   Level: intermediate

2808:   Fortran Notes:
2809:   Use `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`

2811: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2812: @*/
2813: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2814: {
2815:   PetscFunctionBegin;
2817:   if (s->bc) {
2818:     PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2819:   } else *indices = NULL;
2820:   PetscFunctionReturn(PETSC_SUCCESS);
2821: }

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

2826:   Not Collective

2828:   Input Parameters:
2829: + s       - The `PetscSection`
2830: . point   - The point
2831: - indices - The constrained dofs

2833:   Level: intermediate

2835:   Fortran Notes:
2836:   Use `PetscSectionSetConstraintIndicesF90()`

2838: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2839: @*/
2840: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2841: {
2842:   PetscFunctionBegin;
2844:   if (s->bc) {
2845:     const PetscInt dof  = s->atlasDof[point];
2846:     const PetscInt cdof = s->bc->atlasDof[point];
2847:     PetscInt       d;

2849:     if (indices)
2850:       for (d = 0; d < cdof; ++d) PetscCheck(indices[d] < dof, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " dof %" PetscInt_FMT ", invalid constraint index[%" PetscInt_FMT "]: %" PetscInt_FMT, point, dof, d, indices[d]);
2851:     PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
2852:   }
2853:   PetscFunctionReturn(PETSC_SUCCESS);
2854: }

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

2859:   Not Collective

2861:   Input Parameters:
2862: + s     - The `PetscSection`
2863: . field - The field number
2864: - point - The point

2866:   Output Parameter:
2867: . indices - The constrained dofs sorted in ascending order

2869:   Level: intermediate

2871:   Note:
2872:   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()`.

2874:   Fortran Notes:
2875:   Use `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`

2877: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2878: @*/
2879: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2880: {
2881:   PetscFunctionBegin;
2883:   PetscAssertPointer(indices, 4);
2884:   PetscSectionCheckValidField(field, s->numFields);
2885:   PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
2886:   PetscFunctionReturn(PETSC_SUCCESS);
2887: }

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

2892:   Not Collective

2894:   Input Parameters:
2895: + s       - The `PetscSection`
2896: . point   - The point
2897: . field   - The field number
2898: - indices - The constrained dofs

2900:   Level: intermediate

2902:   Fortran Notes:
2903:   Use `PetscSectionSetFieldConstraintIndicesF90()`

2905: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2906: @*/
2907: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2908: {
2909:   PetscFunctionBegin;
2911:   PetscSectionCheckValidField(field, s->numFields);
2912:   PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
2913:   PetscFunctionReturn(PETSC_SUCCESS);
2914: }

2916: /*@
2917:   PetscSectionPermute - Reorder the section according to the input point permutation

2919:   Collective

2921:   Input Parameters:
2922: + section     - The `PetscSection` object
2923: - permutation - The point permutation, old point p becomes new point perm[p]

2925:   Output Parameter:
2926: . sectionNew - The permuted `PetscSection`

2928:   Level: intermediate

2930:   Note:
2931:   The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`

2933:   Compare to `PetscSectionSetPermutation()`

2935: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
2936: @*/
2937: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2938: {
2939:   PetscSection    s = section, sNew;
2940:   const PetscInt *perm;
2941:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;

2943:   PetscFunctionBegin;
2946:   PetscAssertPointer(sectionNew, 3);
2947:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
2948:   PetscCall(PetscSectionGetNumFields(s, &numFields));
2949:   if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
2950:   for (f = 0; f < numFields; ++f) {
2951:     const char *name;
2952:     PetscInt    numComp;

2954:     PetscCall(PetscSectionGetFieldName(s, f, &name));
2955:     PetscCall(PetscSectionSetFieldName(sNew, f, name));
2956:     PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2957:     PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
2958:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2959:       PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2960:       PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
2961:     }
2962:   }
2963:   PetscCall(ISGetLocalSize(permutation, &numPoints));
2964:   PetscCall(ISGetIndices(permutation, &perm));
2965:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2966:   PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
2967:   PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
2968:   for (p = pStart; p < pEnd; ++p) {
2969:     PetscInt dof, cdof;

2971:     PetscCall(PetscSectionGetDof(s, p, &dof));
2972:     PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
2973:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2974:     if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
2975:     for (f = 0; f < numFields; ++f) {
2976:       PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
2977:       PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
2978:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2979:       if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
2980:     }
2981:   }
2982:   PetscCall(PetscSectionSetUp(sNew));
2983:   for (p = pStart; p < pEnd; ++p) {
2984:     const PetscInt *cind;
2985:     PetscInt        cdof;

2987:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2988:     if (cdof) {
2989:       PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
2990:       PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
2991:     }
2992:     for (f = 0; f < numFields; ++f) {
2993:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2994:       if (cdof) {
2995:         PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
2996:         PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
2997:       }
2998:     }
2999:   }
3000:   PetscCall(ISRestoreIndices(permutation, &perm));
3001:   *sectionNew = sNew;
3002:   PetscFunctionReturn(PETSC_SUCCESS);
3003: }

3005: /*@
3006:   PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.

3008:   Collective

3010:   Input Parameters:
3011: + section   - The `PetscSection`
3012: . obj       - A `PetscObject` which serves as the key for this index
3013: . clSection - `PetscSection` giving the size of the closure of each point
3014: - clPoints  - `IS` giving the points in each closure

3016:   Level: advanced

3018:   Note:
3019:   This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.

3021:   Developer Notes:
3022:   The information provided here is completely opaque

3024: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
3025: @*/
3026: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
3027: {
3028:   PetscFunctionBegin;
3032:   if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
3033:   section->clObj = obj;
3034:   PetscCall(PetscObjectReference((PetscObject)clSection));
3035:   PetscCall(PetscObjectReference((PetscObject)clPoints));
3036:   PetscCall(PetscSectionDestroy(&section->clSection));
3037:   PetscCall(ISDestroy(&section->clPoints));
3038:   section->clSection = clSection;
3039:   section->clPoints  = clPoints;
3040:   PetscFunctionReturn(PETSC_SUCCESS);
3041: }

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

3046:   Collective

3048:   Input Parameters:
3049: + section - The `PetscSection`
3050: - obj     - A `PetscObject` which serves as the key for this index

3052:   Output Parameters:
3053: + clSection - `PetscSection` giving the size of the closure of each point
3054: - clPoints  - `IS` giving the points in each closure

3056:   Level: advanced

3058: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3059: @*/
3060: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
3061: {
3062:   PetscFunctionBegin;
3063:   if (section->clObj == obj) {
3064:     if (clSection) *clSection = section->clSection;
3065:     if (clPoints) *clPoints = section->clPoints;
3066:   } else {
3067:     if (clSection) *clSection = NULL;
3068:     if (clPoints) *clPoints = NULL;
3069:   }
3070:   PetscFunctionReturn(PETSC_SUCCESS);
3071: }

3073: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
3074: {
3075:   PetscInt                    i;
3076:   khiter_t                    iter;
3077:   int                         new_entry;
3078:   PetscSectionClosurePermKey  key = {depth, clSize};
3079:   PetscSectionClosurePermVal *val;

3081:   PetscFunctionBegin;
3082:   if (section->clObj != obj) {
3083:     PetscCall(PetscSectionDestroy(&section->clSection));
3084:     PetscCall(ISDestroy(&section->clPoints));
3085:   }
3086:   section->clObj = obj;
3087:   if (!section->clHash) PetscCall(PetscClPermCreate(&section->clHash));
3088:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
3089:   val  = &kh_val(section->clHash, iter);
3090:   if (!new_entry) {
3091:     PetscCall(PetscFree(val->perm));
3092:     PetscCall(PetscFree(val->invPerm));
3093:   }
3094:   if (mode == PETSC_COPY_VALUES) {
3095:     PetscCall(PetscMalloc1(clSize, &val->perm));
3096:     PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
3097:   } else if (mode == PETSC_OWN_POINTER) {
3098:     val->perm = clPerm;
3099:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
3100:   PetscCall(PetscMalloc1(clSize, &val->invPerm));
3101:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
3102:   PetscFunctionReturn(PETSC_SUCCESS);
3103: }

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

3108:   Not Collective

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

3116:   Level: intermediate

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

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

3126: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
3127: @*/
3128: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
3129: {
3130:   const PetscInt *clPerm = NULL;
3131:   PetscInt        clSize = 0;

3133:   PetscFunctionBegin;
3134:   if (perm) {
3135:     PetscCall(ISGetLocalSize(perm, &clSize));
3136:     PetscCall(ISGetIndices(perm, &clPerm));
3137:   }
3138:   PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
3139:   if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
3140:   PetscFunctionReturn(PETSC_SUCCESS);
3141: }

3143: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3144: {
3145:   PetscFunctionBegin;
3146:   if (section->clObj == obj) {
3147:     PetscSectionClosurePermKey k = {depth, size};
3148:     PetscSectionClosurePermVal v;

3150:     PetscCall(PetscClPermGet(section->clHash, k, &v));
3151:     if (perm) *perm = v.perm;
3152:   } else {
3153:     if (perm) *perm = NULL;
3154:   }
3155:   PetscFunctionReturn(PETSC_SUCCESS);
3156: }

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

3161:   Not Collective

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

3169:   Output Parameter:
3170: . perm - The dof closure permutation

3172:   Level: intermediate

3174:   Note:
3175:   The user must destroy the `IS` that is returned.

3177: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3178: @*/
3179: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3180: {
3181:   const PetscInt *clPerm = NULL;

3183:   PetscFunctionBegin;
3184:   PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3185:   PetscCheck(clPerm, PetscObjectComm((PetscObject)obj), PETSC_ERR_ARG_WRONG, "There is no closure permutation associated with this object for depth %" PetscInt_FMT " of size %" PetscInt_FMT, depth, clSize);
3186:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3187:   PetscFunctionReturn(PETSC_SUCCESS);
3188: }

3190: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3191: {
3192:   PetscFunctionBegin;
3193:   if (section->clObj == obj && section->clHash) {
3194:     PetscSectionClosurePermKey k = {depth, size};
3195:     PetscSectionClosurePermVal v;
3196:     PetscCall(PetscClPermGet(section->clHash, k, &v));
3197:     if (perm) *perm = v.invPerm;
3198:   } else {
3199:     if (perm) *perm = NULL;
3200:   }
3201:   PetscFunctionReturn(PETSC_SUCCESS);
3202: }

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

3207:   Not Collective

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

3215:   Output Parameter:
3216: . perm - The dof closure permutation

3218:   Level: intermediate

3220:   Note:
3221:   The user must destroy the `IS` that is returned.

3223: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3224: @*/
3225: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3226: {
3227:   const PetscInt *clPerm = NULL;

3229:   PetscFunctionBegin;
3230:   PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3231:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3232:   PetscFunctionReturn(PETSC_SUCCESS);
3233: }

3235: /*@
3236:   PetscSectionGetField - Get the `PetscSection` associated with a single field

3238:   Input Parameters:
3239: + s     - The `PetscSection`
3240: - field - The field number

3242:   Output Parameter:
3243: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set

3245:   Level: intermediate

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

3250: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3251: @*/
3252: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3253: {
3254:   PetscFunctionBegin;
3256:   PetscAssertPointer(subs, 3);
3257:   PetscSectionCheckValidField(field, s->numFields);
3258:   *subs = s->field[field];
3259:   PetscFunctionReturn(PETSC_SUCCESS);
3260: }

3262: PetscClassId      PETSC_SECTION_SYM_CLASSID;
3263: PetscFunctionList PetscSectionSymList = NULL;

3265: /*@
3266:   PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.

3268:   Collective

3270:   Input Parameter:
3271: . comm - the MPI communicator

3273:   Output Parameter:
3274: . sym - pointer to the new set of symmetries

3276:   Level: developer

3278: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3279: @*/
3280: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3281: {
3282:   PetscFunctionBegin;
3283:   PetscAssertPointer(sym, 2);
3284:   PetscCall(ISInitializePackage());
3285:   PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3286:   PetscFunctionReturn(PETSC_SUCCESS);
3287: }

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

3292:   Collective

3294:   Input Parameters:
3295: + sym    - The section symmetry object
3296: - method - The name of the section symmetry type

3298:   Level: developer

3300: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3301: @*/
3302: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3303: {
3304:   PetscErrorCode (*r)(PetscSectionSym);
3305:   PetscBool match;

3307:   PetscFunctionBegin;
3309:   PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3310:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

3312:   PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3313:   PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3314:   PetscTryTypeMethod(sym, destroy);
3315:   sym->ops->destroy = NULL;

3317:   PetscCall((*r)(sym));
3318:   PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3319:   PetscFunctionReturn(PETSC_SUCCESS);
3320: }

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

3325:   Not Collective

3327:   Input Parameter:
3328: . sym - The section symmetry

3330:   Output Parameter:
3331: . type - The index set type name

3333:   Level: developer

3335: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3336: @*/
3337: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3338: {
3339:   PetscFunctionBegin;
3341:   PetscAssertPointer(type, 2);
3342:   *type = ((PetscObject)sym)->type_name;
3343:   PetscFunctionReturn(PETSC_SUCCESS);
3344: }

3346: /*@C
3347:   PetscSectionSymRegister - Registers a new section symmetry implementation

3349:   Not Collective

3351:   Input Parameters:
3352: + sname    - The name of a new user-defined creation routine
3353: - function - The creation routine itself

3355:   Level: developer

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

3360: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3361: @*/
3362: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3363: {
3364:   PetscFunctionBegin;
3365:   PetscCall(ISInitializePackage());
3366:   PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3367:   PetscFunctionReturn(PETSC_SUCCESS);
3368: }

3370: /*@
3371:   PetscSectionSymDestroy - Destroys a section symmetry.

3373:   Collective

3375:   Input Parameter:
3376: . sym - the section symmetry

3378:   Level: developer

3380: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3381: @*/
3382: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3383: {
3384:   SymWorkLink link, next;

3386:   PetscFunctionBegin;
3387:   if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3389:   if (--((PetscObject)*sym)->refct > 0) {
3390:     *sym = NULL;
3391:     PetscFunctionReturn(PETSC_SUCCESS);
3392:   }
3393:   PetscTryTypeMethod(*sym, destroy);
3394:   PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3395:   for (link = (*sym)->workin; link; link = next) {
3396:     PetscInt    **perms = (PetscInt **)link->perms;
3397:     PetscScalar **rots  = (PetscScalar **)link->rots;
3398:     PetscCall(PetscFree2(perms, rots));
3399:     next = link->next;
3400:     PetscCall(PetscFree(link));
3401:   }
3402:   (*sym)->workin = NULL;
3403:   PetscCall(PetscHeaderDestroy(sym));
3404:   PetscFunctionReturn(PETSC_SUCCESS);
3405: }

3407: /*@C
3408:   PetscSectionSymView - Displays a section symmetry

3410:   Collective

3412:   Input Parameters:
3413: + sym    - the index set
3414: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.

3416:   Level: developer

3418: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3419: @*/
3420: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3421: {
3422:   PetscFunctionBegin;
3424:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3426:   PetscCheckSameComm(sym, 1, viewer, 2);
3427:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3428:   PetscTryTypeMethod(sym, view, viewer);
3429:   PetscFunctionReturn(PETSC_SUCCESS);
3430: }

3432: /*@
3433:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3435:   Collective

3437:   Input Parameters:
3438: + section - the section describing data layout
3439: - sym     - the symmetry describing the affect of orientation on the access of the data

3441:   Level: developer

3443: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3444: @*/
3445: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3446: {
3447:   PetscFunctionBegin;
3449:   PetscCall(PetscSectionSymDestroy(&section->sym));
3450:   if (sym) {
3452:     PetscCheckSameComm(section, 1, sym, 2);
3453:     PetscCall(PetscObjectReference((PetscObject)sym));
3454:   }
3455:   section->sym = sym;
3456:   PetscFunctionReturn(PETSC_SUCCESS);
3457: }

3459: /*@
3460:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3462:   Not Collective

3464:   Input Parameter:
3465: . section - the section describing data layout

3467:   Output Parameter:
3468: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`

3470:   Level: developer

3472: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3473: @*/
3474: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3475: {
3476:   PetscFunctionBegin;
3478:   *sym = section->sym;
3479:   PetscFunctionReturn(PETSC_SUCCESS);
3480: }

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

3485:   Collective

3487:   Input Parameters:
3488: + section - the section describing data layout
3489: . field   - the field number
3490: - sym     - the symmetry describing the affect of orientation on the access of the data

3492:   Level: developer

3494: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3495: @*/
3496: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3497: {
3498:   PetscFunctionBegin;
3500:   PetscSectionCheckValidField(field, section->numFields);
3501:   PetscCall(PetscSectionSetSym(section->field[field], sym));
3502:   PetscFunctionReturn(PETSC_SUCCESS);
3503: }

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

3508:   Collective

3510:   Input Parameters:
3511: + section - the section describing data layout
3512: - field   - the field number

3514:   Output Parameter:
3515: . sym - the symmetry describing the affect of orientation on the access of the data

3517:   Level: developer

3519: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3520: @*/
3521: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3522: {
3523:   PetscFunctionBegin;
3525:   PetscSectionCheckValidField(field, section->numFields);
3526:   *sym = section->field[field]->sym;
3527:   PetscFunctionReturn(PETSC_SUCCESS);
3528: }

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

3533:   Not Collective

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

3542:   Output Parameters:
3543: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3544: - rots  - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3545:     identity).

3547:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3548: .vb
3549:      const PetscInt    **perms;
3550:      const PetscScalar **rots;
3551:      PetscInt            lOffset;

3553:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3554:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3555:        PetscInt           point = points[2*i], dof, sOffset;
3556:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3557:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3559:        PetscSectionGetDof(section,point,&dof);
3560:        PetscSectionGetOffset(section,point,&sOffset);

3562:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3563:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3564:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3565:        lOffset += dof;
3566:      }
3567:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3568: .ve

3570:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3571: .vb
3572:      const PetscInt    **perms;
3573:      const PetscScalar **rots;
3574:      PetscInt            lOffset;

3576:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3577:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3578:        PetscInt           point = points[2*i], dof, sOffset;
3579:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3580:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3582:        PetscSectionGetDof(section,point,&dof);
3583:        PetscSectionGetOffset(section,point,&sOff);

3585:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3586:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3587:        offset += dof;
3588:      }
3589:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3590: .ve

3592:   Level: developer

3594:   Notes:
3595:   `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`

3597:   Use `PetscSectionRestorePointSyms()` when finished with the data

3599: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3600: @*/
3601: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3602: {
3603:   PetscSectionSym sym;

3605:   PetscFunctionBegin;
3607:   if (numPoints) PetscAssertPointer(points, 3);
3608:   if (perms) *perms = NULL;
3609:   if (rots) *rots = NULL;
3610:   sym = section->sym;
3611:   if (sym && (perms || rots)) {
3612:     SymWorkLink link;

3614:     if (sym->workin) {
3615:       link        = sym->workin;
3616:       sym->workin = sym->workin->next;
3617:     } else {
3618:       PetscCall(PetscNew(&link));
3619:     }
3620:     if (numPoints > link->numPoints) {
3621:       PetscInt    **perms = (PetscInt **)link->perms;
3622:       PetscScalar **rots  = (PetscScalar **)link->rots;
3623:       PetscCall(PetscFree2(perms, rots));
3624:       PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3625:       link->numPoints = numPoints;
3626:     }
3627:     link->next   = sym->workout;
3628:     sym->workout = link;
3629:     PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3630:     PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3631:     PetscUseTypeMethod(sym, getpoints, section, numPoints, points, link->perms, link->rots);
3632:     if (perms) *perms = link->perms;
3633:     if (rots) *rots = link->rots;
3634:   }
3635:   PetscFunctionReturn(PETSC_SUCCESS);
3636: }

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

3641:   Not Collective

3643:   Input Parameters:
3644: + section   - the section
3645: . numPoints - the number of points
3646: . points    - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3647:               arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3648:               context, see `DMPlexGetConeOrientation()`).
3649: . perms     - The permutations for the given orientations: set to `NULL` at conclusion
3650: - rots      - The field rotations symmetries for the given orientations: set to `NULL` at conclusion

3652:   Level: developer

3654: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3655: @*/
3656: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3657: {
3658:   PetscSectionSym sym;

3660:   PetscFunctionBegin;
3662:   sym = section->sym;
3663:   if (sym && (perms || rots)) {
3664:     SymWorkLink *p, link;

3666:     for (p = &sym->workout; (link = *p); p = &link->next) {
3667:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3668:         *p          = link->next;
3669:         link->next  = sym->workin;
3670:         sym->workin = link;
3671:         if (perms) *perms = NULL;
3672:         if (rots) *rots = NULL;
3673:         PetscFunctionReturn(PETSC_SUCCESS);
3674:       }
3675:     }
3676:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3677:   }
3678:   PetscFunctionReturn(PETSC_SUCCESS);
3679: }

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

3684:   Not Collective

3686:   Input Parameters:
3687: + section   - the section
3688: . field     - the field of the section
3689: . numPoints - the number of points
3690: - points    - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3691:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3692:     context, see `DMPlexGetConeOrientation()`).

3694:   Output Parameters:
3695: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3696: - rots  - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3697:     identity).

3699:   Level: developer

3701:   Notes:
3702:   `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`

3704:   Use `PetscSectionRestoreFieldPointSyms()` when finished with the data

3706: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3707: @*/
3708: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3709: {
3710:   PetscFunctionBegin;
3712:   PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields);
3713:   PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3714:   PetscFunctionReturn(PETSC_SUCCESS);
3715: }

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

3720:   Not Collective

3722:   Input Parameters:
3723: + section   - the section
3724: . field     - the field number
3725: . numPoints - the number of points
3726: . points    - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3727:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3728:     context, see `DMPlexGetConeOrientation()`).
3729: . perms     - The permutations for the given orientations: set to NULL at conclusion
3730: - rots      - The field rotations symmetries for the given orientations: set to NULL at conclusion

3732:   Level: developer

3734: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3735: @*/
3736: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3737: {
3738:   PetscFunctionBegin;
3740:   PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields);
3741:   PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3742:   PetscFunctionReturn(PETSC_SUCCESS);
3743: }

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

3748:   Not Collective

3750:   Input Parameter:
3751: . sym - the `PetscSectionSym`

3753:   Output Parameter:
3754: . nsym - the equivalent symmetries

3756:   Level: developer

3758: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3759: @*/
3760: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3761: {
3762:   PetscFunctionBegin;
3765:   PetscTryTypeMethod(sym, copy, nsym);
3766:   PetscFunctionReturn(PETSC_SUCCESS);
3767: }

3769: /*@
3770:   PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`

3772:   Collective

3774:   Input Parameters:
3775: + sym         - the `PetscSectionSym`
3776: - migrationSF - the distribution map from roots to leaves

3778:   Output Parameter:
3779: . dsym - the redistributed symmetries

3781:   Level: developer

3783: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3784: @*/
3785: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3786: {
3787:   PetscFunctionBegin;
3790:   PetscAssertPointer(dsym, 3);
3791:   PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3792:   PetscFunctionReturn(PETSC_SUCCESS);
3793: }

3795: /*@
3796:   PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset

3798:   Not Collective

3800:   Input Parameter:
3801: . s - the global `PetscSection`

3803:   Output Parameter:
3804: . flg - the flag

3806:   Level: developer

3808: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3809: @*/
3810: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3811: {
3812:   PetscFunctionBegin;
3814:   *flg = s->useFieldOff;
3815:   PetscFunctionReturn(PETSC_SUCCESS);
3816: }

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

3821:   Not Collective

3823:   Input Parameters:
3824: + s   - the global `PetscSection`
3825: - flg - the flag

3827:   Level: developer

3829: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3830: @*/
3831: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3832: {
3833:   PetscFunctionBegin;
3835:   s->useFieldOff = flg;
3836:   PetscFunctionReturn(PETSC_SUCCESS);
3837: }

3839: #define PetscSectionExpandPoints_Loop(TYPE) \
3840:   do { \
3841:     PetscInt i, n, o0, o1, size; \
3842:     TYPE    *a0 = (TYPE *)origArray, *a1; \
3843:     PetscCall(PetscSectionGetStorageSize(s, &size)); \
3844:     PetscCall(PetscMalloc1(size, &a1)); \
3845:     for (i = 0; i < npoints; i++) { \
3846:       PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3847:       PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3848:       PetscCall(PetscSectionGetDof(s, i, &n)); \
3849:       PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \
3850:     } \
3851:     *newArray = (void *)a1; \
3852:   } while (0)

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

3857:   Not Collective

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

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

3869:   Level: developer

3871: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3872: @*/
3873: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3874: {
3875:   PetscSection    s;
3876:   const PetscInt *points_;
3877:   PetscInt        i, n, npoints, pStart, pEnd;
3878:   PetscMPIInt     unitsize;

3880:   PetscFunctionBegin;
3882:   PetscAssertPointer(origArray, 3);
3884:   if (newSection) PetscAssertPointer(newSection, 5);
3885:   if (newArray) PetscAssertPointer(newArray, 6);
3886:   PetscCallMPI(MPI_Type_size(dataType, &unitsize));
3887:   PetscCall(ISGetLocalSize(points, &npoints));
3888:   PetscCall(ISGetIndices(points, &points_));
3889:   PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
3890:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
3891:   PetscCall(PetscSectionSetChart(s, 0, npoints));
3892:   for (i = 0; i < npoints; i++) {
3893:     PetscCheck(points_[i] >= pStart && points_[i] < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " (index %" PetscInt_FMT ") in input IS out of input section's chart", points_[i], i);
3894:     PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
3895:     PetscCall(PetscSectionSetDof(s, i, n));
3896:   }
3897:   PetscCall(PetscSectionSetUp(s));
3898:   if (newArray) {
3899:     if (dataType == MPIU_INT) {
3900:       PetscSectionExpandPoints_Loop(PetscInt);
3901:     } else if (dataType == MPIU_SCALAR) {
3902:       PetscSectionExpandPoints_Loop(PetscScalar);
3903:     } else if (dataType == MPIU_REAL) {
3904:       PetscSectionExpandPoints_Loop(PetscReal);
3905:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3906:   }
3907:   if (newSection) {
3908:     *newSection = s;
3909:   } else {
3910:     PetscCall(PetscSectionDestroy(&s));
3911:   }
3912:   PetscCall(ISRestoreIndices(points, &points_));
3913:   PetscFunctionReturn(PETSC_SUCCESS);
3914: }