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:   PetscSectionSym sym;
 89:   IS              perm;
 90:   PetscInt        numFields, f, c, pStart, pEnd, p;

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

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

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

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

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

162:   Collective

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

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

170:   Level: beginner

172:   Developer Notes:
173:   With standard PETSc terminology this should be called `PetscSectionDuplicate()`

175: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
176: @*/
177: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
178: {
179:   PetscFunctionBegin;
181:   PetscAssertPointer(newSection, 2);
182:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
183:   PetscCall(PetscSectionCopy(section, *newSection));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

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

190:   Collective

192:   Input Parameter:
193: . s - the `PetscSection`

195:   Options Database Key:
196: . -petscsection_point_major - `PETSC_TRUE` for point-major order

198:   Level: intermediate

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

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

218:   Collective

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

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

227:   Level: intermediate

229:   Note:
230:   Field names are disregarded.

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

242:   PetscFunctionBegin;
245:   PetscAssertPointer(congruent, 3);
246:   flg = PETSC_FALSE;

248:   PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
249:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
250:     *congruent = PETSC_FALSE;
251:     PetscFunctionReturn(PETSC_SUCCESS);
252:   }

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

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

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

270:     PetscCall(PetscSectionGetDof(s1, p, &n1));
271:     PetscCall(PetscSectionGetDof(s2, p, &n2));
272:     if (n1 != n2) goto not_congruent;

274:     PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
275:     PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
276:     if (ncdof != n2) goto not_congruent;

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

284:   PetscCall(PetscSectionGetNumFields(s1, &nfields));
285:   PetscCall(PetscSectionGetNumFields(s2, &n2));
286:   if (nfields != n2) goto not_congruent;

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

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

298:       PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
299:       PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
300:       if (n1 != n2) goto not_congruent;

302:       PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
303:       PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
304:       if (nfcdof != n2) goto not_congruent;

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

313:   flg = PETSC_TRUE;
314: not_congruent:
315:   PetscCall(MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }

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

322:   Not Collective

324:   Input Parameter:
325: . s - the `PetscSection`

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

330:   Level: intermediate

332: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
333: @*/
334: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
335: {
336:   PetscFunctionBegin;
338:   PetscAssertPointer(numFields, 2);
339:   *numFields = s->numFields;
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

343: /*@
344:   PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`

346:   Not Collective

348:   Input Parameters:
349: + s         - the `PetscSection`
350: - numFields - the number of fields

352:   Level: intermediate

354:   Notes:
355:   Calling this destroys all the information in the `PetscSection` including the chart.

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

359: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
360: @*/
361: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
362: {
363:   PetscInt f;

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

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

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

380:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
381:     PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
382:     PetscCall(PetscStrallocpy(name, (char **)&s->fieldNames[f]));
383:     PetscCall(PetscSNPrintf(name, 64, "Component_0"));
384:     PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
385:     PetscCall(PetscStrallocpy(name, (char **)&s->compNames[f][0]));
386:   }
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

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

393:   Not Collective

395:   Input Parameters:
396: + s     - the `PetscSection`
397: - field - the field number

399:   Output Parameter:
400: . fieldName - the field name

402:   Level: intermediate

404:   Note:
405:   Will error if the field number is out of range

407: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
408: @*/
409: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
410: {
411:   PetscFunctionBegin;
413:   PetscAssertPointer(fieldName, 3);
414:   PetscSectionCheckValidField(field, s->numFields);
415:   *fieldName = s->fieldNames[field];
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

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

422:   Not Collective

424:   Input Parameters:
425: + s         - the `PetscSection`
426: . field     - the field number
427: - fieldName - the field name

429:   Level: intermediate

431:   Note:
432:   Will error if the field number is out of range

434: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
435: @*/
436: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
437: {
438:   PetscFunctionBegin;
440:   if (fieldName) PetscAssertPointer(fieldName, 3);
441:   PetscSectionCheckValidField(field, s->numFields);
442:   PetscCall(PetscFree(s->fieldNames[field]));
443:   PetscCall(PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]));
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

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

450:   Not Collective

452:   Input Parameters:
453: + s     - the `PetscSection`
454: . field - the field number
455: - comp  - the component number

457:   Output Parameter:
458: . compName - the component name

460:   Level: intermediate

462:   Note:
463:   Will error if the field or component number do not exist

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

468: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
469:           `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
470: @*/
471: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
472: {
473:   PetscFunctionBegin;
475:   PetscAssertPointer(compName, 4);
476:   PetscSectionCheckValidField(field, s->numFields);
477:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
478:   *compName = s->compNames[field][comp];
479:   PetscFunctionReturn(PETSC_SUCCESS);
480: }

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

485:   Not Collective

487:   Input Parameters:
488: + s        - the `PetscSection`
489: . field    - the field number
490: . comp     - the component number
491: - compName - the component name

493:   Level: advanced

495:   Note:
496:   Will error if the field or component number do not exist

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

501: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
502:           `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
503: @*/
504: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
505: {
506:   PetscFunctionBegin;
508:   if (compName) PetscAssertPointer(compName, 4);
509:   PetscSectionCheckValidField(field, s->numFields);
510:   PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
511:   PetscCall(PetscFree(s->compNames[field][comp]));
512:   PetscCall(PetscStrallocpy(compName, (char **)&s->compNames[field][comp]));
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: /*@
517:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

519:   Not Collective

521:   Input Parameters:
522: + s     - the `PetscSection`
523: - field - the field number

525:   Output Parameter:
526: . numComp - the number of field components

528:   Level: advanced

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

533: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
534:           `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
535: @*/
536: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
537: {
538:   PetscFunctionBegin;
540:   PetscAssertPointer(numComp, 3);
541:   PetscSectionCheckValidField(field, s->numFields);
542:   *numComp = s->numFieldComponents[field];
543:   PetscFunctionReturn(PETSC_SUCCESS);
544: }

546: /*@
547:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

549:   Not Collective

551:   Input Parameters:
552: + s       - the `PetscSection`
553: . field   - the field number
554: - numComp - the number of field components

556:   Level: advanced

558:   Note:
559:   This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
560:   components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
561:   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
562:   an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.

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

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

569: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
570:           `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
571: @*/
572: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
573: {
574:   PetscInt c;

576:   PetscFunctionBegin;
578:   PetscSectionCheckValidField(field, s->numFields);
579:   if (s->compNames) {
580:     for (c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
581:     PetscCall(PetscFree(s->compNames[field]));
582:   }

584:   s->numFieldComponents[field] = numComp;
585:   if (numComp) {
586:     PetscCall(PetscMalloc1(numComp, (char ***)&s->compNames[field]));
587:     for (c = 0; c < numComp; ++c) {
588:       char name[64];

590:       PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
591:       PetscCall(PetscStrallocpy(name, (char **)&s->compNames[field][c]));
592:     }
593:   }
594:   PetscFunctionReturn(PETSC_SUCCESS);
595: }

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

600:   Not Collective

602:   Input Parameter:
603: . s - the `PetscSection`

605:   Output Parameters:
606: + pStart - the first point
607: - pEnd   - one past the last point

609:   Level: intermediate

611: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
612: @*/
613: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
614: {
615:   PetscFunctionBegin;
617:   if (pStart) *pStart = s->pStart;
618:   if (pEnd) *pEnd = s->pEnd;
619:   PetscFunctionReturn(PETSC_SUCCESS);
620: }

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

625:   Not Collective

627:   Input Parameters:
628: + s      - the `PetscSection`
629: . pStart - the first point
630: - pEnd   - one past the last point, `pStart` $ \le $ `pEnd`

632:   Level: intermediate

634:   Notes:
635:   The charts on different MPI processes may (and often do) overlap

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

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

641: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
642: @*/
643: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
644: {
645:   PetscInt f;

647:   PetscFunctionBegin;
649:   if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
650:   /* Cannot Reset() because it destroys field information */
651:   s->setup = PETSC_FALSE;
652:   PetscCall(PetscSectionDestroy(&s->bc));
653:   PetscCall(PetscFree(s->bcIndices));
654:   PetscCall(PetscFree2(s->atlasDof, s->atlasOff));

656:   s->pStart = pStart;
657:   s->pEnd   = pEnd;
658:   PetscCall(PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff));
659:   PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
660:   for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
661:   PetscFunctionReturn(PETSC_SUCCESS);
662: }

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

667:   Not Collective

669:   Input Parameter:
670: . s - the `PetscSection`

672:   Output Parameter:
673: . perm - The permutation as an `IS`

675:   Level: intermediate

677: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
678: @*/
679: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
680: {
681:   PetscFunctionBegin;
683:   if (perm) {
684:     PetscAssertPointer(perm, 2);
685:     *perm = s->perm;
686:   }
687:   PetscFunctionReturn(PETSC_SUCCESS);
688: }

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

693:   Not Collective

695:   Input Parameters:
696: + s    - the `PetscSection`
697: - perm - the permutation of points

699:   Level: intermediate

701:   Notes:
702:   The permutation must be provided before `PetscSectionSetUp()`.

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

706:   Compart to `PetscSectionPermute()`

708: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
709: @*/
710: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
711: {
712:   PetscFunctionBegin;
715:   PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
716:   if (s->perm != perm) {
717:     PetscCall(ISDestroy(&s->perm));
718:     if (perm) {
719:       s->perm = perm;
720:       PetscCall(PetscObjectReference((PetscObject)s->perm));
721:     }
722:   }
723:   PetscFunctionReturn(PETSC_SUCCESS);
724: }

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

729:   Not Collective

731:   Input Parameter:
732: . s - the `PetscSection`

734:   Output Parameter:
735: . pm - the flag for point major ordering

737:   Level: intermediate

739: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
740: @*/
741: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
742: {
743:   PetscFunctionBegin;
745:   PetscAssertPointer(pm, 2);
746:   *pm = s->pointMajor;
747:   PetscFunctionReturn(PETSC_SUCCESS);
748: }

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

753:   Not Collective

755:   Input Parameters:
756: + s  - the `PetscSection`
757: - pm - the flag for point major ordering

759:   Level: intermediate

761:   Note:
762:   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.

764:   Point major order means the degrees of freedom are stored as follows
765: .vb
766:     all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
767:     for each point
768:        the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
769: .ve

771:   Field major order means the degrees of freedom are stored as follows
772: .vb
773:     all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
774:     for each field (started with unnamed default field)
775:       the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
776: .ve

778: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
779: @*/
780: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
781: {
782:   PetscFunctionBegin;
784:   PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
785:   s->pointMajor = pm;
786:   PetscFunctionReturn(PETSC_SUCCESS);
787: }

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

793:   Not Collective

795:   Input Parameter:
796: . s - the `PetscSection`

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

801:   Level: intermediate

803: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
804: @*/
805: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
806: {
807:   PetscFunctionBegin;
809:   PetscAssertPointer(includesConstraints, 2);
810:   *includesConstraints = s->includesConstraints;
811:   PetscFunctionReturn(PETSC_SUCCESS);
812: }

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

817:   Not Collective

819:   Input Parameters:
820: + s                   - the `PetscSection`
821: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets

823:   Level: intermediate

825: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
826: @*/
827: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
828: {
829:   PetscFunctionBegin;
831:   PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
832:   s->includesConstraints = includesConstraints;
833:   PetscFunctionReturn(PETSC_SUCCESS);
834: }

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

839:   Not Collective

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

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

848:   Level: intermediate

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

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

855: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
856: @*/
857: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
858: {
859:   PetscFunctionBeginHot;
861:   PetscAssertPointer(numDof, 3);
862:   if (PetscDefined(USE_DEBUG)) 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);
863:   *numDof = s->atlasDof[point - s->pStart];
864:   PetscFunctionReturn(PETSC_SUCCESS);
865: }

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

870:   Not Collective

872:   Input Parameters:
873: + s      - the `PetscSection`
874: . point  - the point
875: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process

877:   Level: intermediate

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

882: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
883: @*/
884: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
885: {
886:   PetscFunctionBegin;
888:   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);
889:   s->atlasDof[point - s->pStart] = numDof;
890:   PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }

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

897:   Not Collective

899:   Input Parameters:
900: + s      - the `PetscSection`
901: . point  - the point
902: - numDof - the number of additional dof

904:   Level: intermediate

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

909: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
910: @*/
911: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
912: {
913:   PetscFunctionBeginHot;
915:   if (PetscDefined(USE_DEBUG)) 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);
916:   s->atlasDof[point - s->pStart] += numDof;
917:   PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
918:   PetscFunctionReturn(PETSC_SUCCESS);
919: }

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

924:   Not Collective

926:   Input Parameters:
927: + s     - the `PetscSection`
928: . point - the point
929: - field - the field

931:   Output Parameter:
932: . numDof - the number of dof

934:   Level: intermediate

936: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
937: @*/
938: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
939: {
940:   PetscFunctionBegin;
942:   PetscAssertPointer(numDof, 4);
943:   PetscSectionCheckValidField(field, s->numFields);
944:   PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
945:   PetscFunctionReturn(PETSC_SUCCESS);
946: }

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

951:   Not Collective

953:   Input Parameters:
954: + s      - the `PetscSection`
955: . point  - the point
956: . field  - the field
957: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process

959:   Level: intermediate

961:   Note:
962:   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
963:   the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`

965:   This is equivalent to
966: .vb
967:      PetscSection fs;
968:      PetscSectionGetField(s,field,&fs)
969:      PetscSectionSetDof(fs,numDof)
970: .ve

972: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
973: @*/
974: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
975: {
976:   PetscFunctionBegin;
978:   PetscSectionCheckValidField(field, s->numFields);
979:   PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
980:   PetscFunctionReturn(PETSC_SUCCESS);
981: }

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

986:   Not Collective

988:   Input Parameters:
989: + s      - the `PetscSection`
990: . point  - the point
991: . field  - the field
992: - numDof - the number of dof

994:   Level: intermediate

996:   Notes:
997:   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
998:   the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`

1000:   This is equivalent to
1001: .vb
1002:      PetscSection fs;
1003:      PetscSectionGetField(s,field,&fs)
1004:      PetscSectionAddDof(fs,numDof)
1005: .ve

1007: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1008: @*/
1009: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1010: {
1011:   PetscFunctionBegin;
1013:   PetscSectionCheckValidField(field, s->numFields);
1014:   PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1015:   PetscFunctionReturn(PETSC_SUCCESS);
1016: }

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

1021:   Not Collective

1023:   Input Parameters:
1024: + s     - the `PetscSection`
1025: - point - the point

1027:   Output Parameter:
1028: . numDof - the number of dof which are fixed by constraints

1030:   Level: intermediate

1032: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1033: @*/
1034: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1035: {
1036:   PetscFunctionBegin;
1038:   PetscAssertPointer(numDof, 3);
1039:   if (s->bc) {
1040:     PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1041:   } else *numDof = 0;
1042:   PetscFunctionReturn(PETSC_SUCCESS);
1043: }

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

1048:   Not Collective

1050:   Input Parameters:
1051: + s      - the `PetscSection`
1052: . point  - the point
1053: - numDof - the number of dof which are fixed by constraints

1055:   Level: intermediate

1057: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1058: @*/
1059: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1060: {
1061:   PetscFunctionBegin;
1063:   if (numDof) {
1064:     PetscCall(PetscSectionCheckConstraints_Private(s));
1065:     PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1066:   }
1067:   PetscFunctionReturn(PETSC_SUCCESS);
1068: }

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

1073:   Not Collective

1075:   Input Parameters:
1076: + s      - the `PetscSection`
1077: . point  - the point
1078: - numDof - the number of additional dof which are fixed by constraints

1080:   Level: intermediate

1082: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1083: @*/
1084: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1085: {
1086:   PetscFunctionBegin;
1088:   if (numDof) {
1089:     PetscCall(PetscSectionCheckConstraints_Private(s));
1090:     PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1091:   }
1092:   PetscFunctionReturn(PETSC_SUCCESS);
1093: }

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

1098:   Not Collective

1100:   Input Parameters:
1101: + s     - the `PetscSection`
1102: . point - the point
1103: - field - the field

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

1108:   Level: intermediate

1110: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1111: @*/
1112: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1113: {
1114:   PetscFunctionBegin;
1116:   PetscAssertPointer(numDof, 4);
1117:   PetscSectionCheckValidField(field, s->numFields);
1118:   PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1119:   PetscFunctionReturn(PETSC_SUCCESS);
1120: }

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

1125:   Not Collective

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

1133:   Level: intermediate

1135: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1136: @*/
1137: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1138: {
1139:   PetscFunctionBegin;
1141:   PetscSectionCheckValidField(field, s->numFields);
1142:   PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1143:   PetscFunctionReturn(PETSC_SUCCESS);
1144: }

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

1149:   Not Collective

1151:   Input Parameters:
1152: + s      - the `PetscSection`
1153: . point  - the point
1154: . field  - the field
1155: - numDof - the number of additional dof which are fixed by constraints

1157:   Level: intermediate

1159: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1160: @*/
1161: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1162: {
1163:   PetscFunctionBegin;
1165:   PetscSectionCheckValidField(field, s->numFields);
1166:   PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1167:   PetscFunctionReturn(PETSC_SUCCESS);
1168: }

1170: /*@
1171:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

1173:   Not Collective

1175:   Input Parameter:
1176: . s - the `PetscSection`

1178:   Level: advanced

1180: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1181: @*/
1182: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1183: {
1184:   PetscFunctionBegin;
1186:   if (s->bc) {
1187:     const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;

1189:     PetscCall(PetscSectionSetUp(s->bc));
1190:     PetscCall(PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices));
1191:   }
1192:   PetscFunctionReturn(PETSC_SUCCESS);
1193: }

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

1198:   Not Collective

1200:   Input Parameter:
1201: . s - the `PetscSection`

1203:   Level: intermediate

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

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

1210: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1211: @*/
1212: PetscErrorCode PetscSectionSetUp(PetscSection s)
1213: {
1214:   const PetscInt *pind   = NULL;
1215:   PetscInt        offset = 0, foff, p, f;

1217:   PetscFunctionBegin;
1219:   if (s->setup) PetscFunctionReturn(PETSC_SUCCESS);
1220:   s->setup = PETSC_TRUE;
1221:   /* Set offsets and field offsets for all points */
1222:   /*   Assume that all fields have the same chart */
1223:   PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1224:   if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1225:   if (s->pointMajor) {
1226:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
1227:       const PetscInt q = pind ? pind[p] : p;

1229:       /* Set point offset */
1230:       s->atlasOff[q] = offset;
1231:       offset += s->atlasDof[q];
1232:       /* Set field offset */
1233:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1234:         PetscSection sf = s->field[f];

1236:         sf->atlasOff[q] = foff;
1237:         foff += sf->atlasDof[q];
1238:       }
1239:     }
1240:   } else {
1241:     /* Set field offsets for all points */
1242:     for (f = 0; f < s->numFields; ++f) {
1243:       PetscSection sf = s->field[f];

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

1248:         sf->atlasOff[q] = offset;
1249:         offset += sf->atlasDof[q];
1250:       }
1251:     }
1252:     /* Disable point offsets since these are unused */
1253:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1254:   }
1255:   if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1256:   /* Setup BC sections */
1257:   PetscCall(PetscSectionSetUpBC(s));
1258:   for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1259:   PetscFunctionReturn(PETSC_SUCCESS);
1260: }

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

1265:   Not Collective

1267:   Input Parameter:
1268: . s - the `PetscSection`

1270:   Output Parameter:
1271: . maxDof - the maximum dof

1273:   Level: intermediate

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

1278:   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
1279:   the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process

1281:   Developer Notes:
1282:   The returned number is calculated lazily and stashed.

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

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

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

1290: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1291: @*/
1292: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1293: {
1294:   PetscInt p;

1296:   PetscFunctionBegin;
1298:   PetscAssertPointer(maxDof, 2);
1299:   if (s->maxDof == PETSC_MIN_INT) {
1300:     s->maxDof = 0;
1301:     for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1302:   }
1303:   *maxDof = s->maxDof;
1304:   PetscFunctionReturn(PETSC_SUCCESS);
1305: }

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

1310:   Not Collective

1312:   Input Parameter:
1313: . s - the `PetscSection`

1315:   Output Parameter:
1316: . size - the size of an array which can hold all the dofs

1318:   Level: intermediate

1320: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1321: @*/
1322: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1323: {
1324:   PetscInt p, n = 0;

1326:   PetscFunctionBegin;
1328:   PetscAssertPointer(size, 2);
1329:   for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1330:   *size = n;
1331:   PetscFunctionReturn(PETSC_SUCCESS);
1332: }

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

1337:   Not Collective

1339:   Input Parameter:
1340: . s - the `PetscSection`

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

1345:   Level: intermediate

1347: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1348: @*/
1349: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1350: {
1351:   PetscInt p, n = 0;

1353:   PetscFunctionBegin;
1355:   PetscAssertPointer(size, 2);
1356:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
1357:     const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1358:     n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1359:   }
1360:   *size = n;
1361:   PetscFunctionReturn(PETSC_SUCCESS);
1362: }

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

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

1374:   Output Parameter:
1375: . gsection - The `PetscSection` for the global field layout

1377:   Level: intermediate

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

1382:   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`.
1383:   In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).

1385: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1386: @*/
1387: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1388: {
1389:   PetscSection    gs;
1390:   const PetscInt *pind = NULL;
1391:   PetscInt       *recv = NULL, *neg = NULL;
1392:   PetscInt        pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1393:   PetscInt        numFields, f, numComponents;

1395:   PetscFunctionBegin;
1400:   PetscAssertPointer(gsection, 5);
1401:   PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1402:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1403:   PetscCall(PetscSectionGetNumFields(s, &numFields));
1404:   if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1405:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1406:   PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1407:   gs->includesConstraints = includeConstraints;
1408:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1409:   nlocal = nroots; /* The local/leaf space matches global/root space */
1410:   /* Must allocate for all points visible to SF, which may be more than this section */
1411:   if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1412:     PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1413:     PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1414:     PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1415:     PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1416:     PetscCall(PetscArrayzero(neg, nroots));
1417:   }
1418:   /* Mark all local points with negative dof */
1419:   for (p = pStart; p < pEnd; ++p) {
1420:     PetscCall(PetscSectionGetDof(s, p, &dof));
1421:     PetscCall(PetscSectionSetDof(gs, p, dof));
1422:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1423:     if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1424:     if (neg) neg[p] = -(dof + 1);
1425:   }
1426:   PetscCall(PetscSectionSetUpBC(gs));
1427:   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]));
1428:   if (nroots >= 0) {
1429:     PetscCall(PetscArrayzero(recv, nlocal));
1430:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1431:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1432:     for (p = pStart; p < pEnd; ++p) {
1433:       if (recv[p] < 0) {
1434:         gs->atlasDof[p - pStart] = recv[p];
1435:         PetscCall(PetscSectionGetDof(s, p, &dof));
1436:         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);
1437:       }
1438:     }
1439:   }
1440:   /* Calculate new sizes, get process offset, and calculate point offsets */
1441:   if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1442:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1443:     const PetscInt q = pind ? pind[p] : p;

1445:     cdof            = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1446:     gs->atlasOff[q] = off;
1447:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1448:   }
1449:   if (!localOffsets) {
1450:     PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1451:     globalOff -= off;
1452:   }
1453:   for (p = pStart, off = 0; p < pEnd; ++p) {
1454:     gs->atlasOff[p - pStart] += globalOff;
1455:     if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1456:   }
1457:   if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1458:   /* Put in negative offsets for ghost points */
1459:   if (nroots >= 0) {
1460:     PetscCall(PetscArrayzero(recv, nlocal));
1461:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1462:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1463:     for (p = pStart; p < pEnd; ++p) {
1464:       if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1465:     }
1466:   }
1467:   PetscCall(PetscFree2(neg, recv));
1468:   /* Set field dofs/offsets/constraints */
1469:   for (f = 0; f < numFields; ++f) {
1470:     gs->field[f]->includesConstraints = includeConstraints;
1471:     PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1472:     PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1473:   }
1474:   for (p = pStart; p < pEnd; ++p) {
1475:     PetscCall(PetscSectionGetOffset(gs, p, &off));
1476:     for (f = 0, foff = off; f < numFields; ++f) {
1477:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1478:       if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1479:       PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1480:       PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1481:       PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff));
1482:       PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1483:       foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1484:     }
1485:   }
1486:   for (f = 0; f < numFields; ++f) {
1487:     PetscSection gfs = gs->field[f];

1489:     PetscCall(PetscSectionSetUpBC(gfs));
1490:     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]));
1491:   }
1492:   gs->setup = PETSC_TRUE;
1493:   PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1494:   *gsection = gs;
1495:   PetscFunctionReturn(PETSC_SUCCESS);
1496: }

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

1502:   Input Parameters:
1503: + s                  - The `PetscSection` for the local field layout
1504: . sf                 - The `PetscSF` describing parallel layout of the section points
1505: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1506: . numExcludes        - The number of exclusion ranges, this must have the same value on all MPI processes
1507: - 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

1509:   Output Parameter:
1510: . gsection - The `PetscSection` for the global field layout

1512:   Level: advanced

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

1517:   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`.
1518:   In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).

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

1522:   Developer Notes:
1523:   This is a terrible function name

1525: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1526: @*/
1527: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1528: {
1529:   const PetscInt *pind = NULL;
1530:   PetscInt       *neg = NULL, *tmpOff = NULL;
1531:   PetscInt        pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;

1533:   PetscFunctionBegin;
1536:   PetscAssertPointer(gsection, 6);
1537:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1538:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1539:   PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1540:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1541:   if (nroots >= 0) {
1542:     PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1543:     PetscCall(PetscCalloc1(nroots, &neg));
1544:     if (nroots > pEnd - pStart) {
1545:       PetscCall(PetscCalloc1(nroots, &tmpOff));
1546:     } else {
1547:       tmpOff = &(*gsection)->atlasDof[-pStart];
1548:     }
1549:   }
1550:   /* Mark ghost points with negative dof */
1551:   for (p = pStart; p < pEnd; ++p) {
1552:     for (e = 0; e < numExcludes; ++e) {
1553:       if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1554:         PetscCall(PetscSectionSetDof(*gsection, p, 0));
1555:         break;
1556:       }
1557:     }
1558:     if (e < numExcludes) continue;
1559:     PetscCall(PetscSectionGetDof(s, p, &dof));
1560:     PetscCall(PetscSectionSetDof(*gsection, p, dof));
1561:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1562:     if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1563:     if (neg) neg[p] = -(dof + 1);
1564:   }
1565:   PetscCall(PetscSectionSetUpBC(*gsection));
1566:   if (nroots >= 0) {
1567:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1568:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1569:     if (nroots > pEnd - pStart) {
1570:       for (p = pStart; p < pEnd; ++p) {
1571:         if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1572:       }
1573:     }
1574:   }
1575:   /* Calculate new sizes, get process offset, and calculate point offsets */
1576:   if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1577:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1578:     const PetscInt q = pind ? pind[p] : p;

1580:     cdof                     = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1581:     (*gsection)->atlasOff[q] = off;
1582:     off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1583:   }
1584:   PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1585:   globalOff -= off;
1586:   for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1587:     (*gsection)->atlasOff[p] += globalOff;
1588:     if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1589:   }
1590:   if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1591:   /* Put in negative offsets for ghost points */
1592:   if (nroots >= 0) {
1593:     if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1594:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1595:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1596:     if (nroots > pEnd - pStart) {
1597:       for (p = pStart; p < pEnd; ++p) {
1598:         if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1599:       }
1600:     }
1601:   }
1602:   if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1603:   PetscCall(PetscFree(neg));
1604:   PetscFunctionReturn(PETSC_SUCCESS);
1605: }

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

1610:   Collective

1612:   Input Parameters:
1613: + comm - The `MPI_Comm`
1614: - s    - The `PetscSection`

1616:   Output Parameter:
1617: . layout - The point layout for the data that defines the section

1619:   Level: advanced

1621:   Notes:
1622:   `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1623:   degrees of freedom).

1625:   This count includes constrained degrees of freedom

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

1629:   Example:
1630: .vb
1631:      The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1632:      The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1633: .ve

1635:   Developer Notes:
1636:   I find the names of these two functions extremely non-informative

1638: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1639: @*/
1640: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1641: {
1642:   PetscInt pStart, pEnd, p, localSize = 0;

1644:   PetscFunctionBegin;
1645:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1646:   for (p = pStart; p < pEnd; ++p) {
1647:     PetscInt dof;

1649:     PetscCall(PetscSectionGetDof(s, p, &dof));
1650:     if (dof >= 0) ++localSize;
1651:   }
1652:   PetscCall(PetscLayoutCreate(comm, layout));
1653:   PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1654:   PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1655:   PetscCall(PetscLayoutSetUp(*layout));
1656:   PetscFunctionReturn(PETSC_SUCCESS);
1657: }

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

1662:   Collective

1664:   Input Parameters:
1665: + comm - The `MPI_Comm`
1666: - s    - The `PetscSection`

1668:   Output Parameter:
1669: . layout - The dof layout for the section

1671:   Level: advanced

1673:   Notes:
1674:   `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1675:   including the constrained degrees of freedom

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

1679:   Example:
1680: .vb
1681:      The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1682:      The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1683: .ve

1685: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1686: @*/
1687: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1688: {
1689:   PetscInt pStart, pEnd, p, localSize = 0;

1691:   PetscFunctionBegin;
1693:   PetscAssertPointer(layout, 3);
1694:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1695:   for (p = pStart; p < pEnd; ++p) {
1696:     PetscInt dof, cdof;

1698:     PetscCall(PetscSectionGetDof(s, p, &dof));
1699:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1700:     if (dof - cdof > 0) localSize += dof - cdof;
1701:   }
1702:   PetscCall(PetscLayoutCreate(comm, layout));
1703:   PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1704:   PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1705:   PetscCall(PetscLayoutSetUp(*layout));
1706:   PetscFunctionReturn(PETSC_SUCCESS);
1707: }

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

1712:   Not Collective

1714:   Input Parameters:
1715: + s     - the `PetscSection`
1716: - point - the point

1718:   Output Parameter:
1719: . offset - the offset

1721:   Level: intermediate

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

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

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

1730: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1731: @*/
1732: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1733: {
1734:   PetscFunctionBegin;
1736:   PetscAssertPointer(offset, 3);
1737:   if (PetscDefined(USE_DEBUG)) 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);
1738:   *offset = s->atlasOff[point - s->pStart];
1739:   PetscFunctionReturn(PETSC_SUCCESS);
1740: }

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

1745:   Not Collective

1747:   Input Parameters:
1748: + s      - the `PetscSection`
1749: . point  - the point
1750: - offset - the offset, these values may be negative indicating the values are off process

1752:   Level: developer

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

1757: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1758: @*/
1759: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1760: {
1761:   PetscFunctionBegin;
1763:   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);
1764:   s->atlasOff[point - s->pStart] = offset;
1765:   PetscFunctionReturn(PETSC_SUCCESS);
1766: }

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

1771:   Not Collective

1773:   Input Parameters:
1774: + s     - the `PetscSection`
1775: . point - the point
1776: - field - the field

1778:   Output Parameter:
1779: . offset - the offset

1781:   Level: intermediate

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

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

1788: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1789: @*/
1790: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1791: {
1792:   PetscFunctionBegin;
1794:   PetscAssertPointer(offset, 4);
1795:   PetscSectionCheckValidField(field, s->numFields);
1796:   PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1797:   PetscFunctionReturn(PETSC_SUCCESS);
1798: }

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

1803:   Not Collective

1805:   Input Parameters:
1806: + s      - the `PetscSection`
1807: . point  - the point
1808: . field  - the field
1809: - offset - the offset, these values may be negative indicating the values are off process

1811:   Level: developer

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

1816: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1817: @*/
1818: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1819: {
1820:   PetscFunctionBegin;
1822:   PetscSectionCheckValidField(field, s->numFields);
1823:   PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1824:   PetscFunctionReturn(PETSC_SUCCESS);
1825: }

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

1831:   Not Collective

1833:   Input Parameters:
1834: + s     - the `PetscSection`
1835: . point - the point
1836: - field - the field

1838:   Output Parameter:
1839: . offset - the offset

1841:   Level: advanced

1843:   Note:
1844:   This ignores constraints

1846:   Example:
1847: .vb
1848:   if PetscSectionSetPointMajor(s,PETSC_TRUE)
1849:   The unnamed default field has 3 dof at `point`
1850:   Field 0 has 2 dof at `point`
1851:   Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1852: .ve

1854: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1855: @*/
1856: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1857: {
1858:   PetscInt off, foff;

1860:   PetscFunctionBegin;
1862:   PetscAssertPointer(offset, 4);
1863:   PetscSectionCheckValidField(field, s->numFields);
1864:   PetscCall(PetscSectionGetOffset(s, point, &off));
1865:   PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1866:   *offset = foff - off;
1867:   PetscFunctionReturn(PETSC_SUCCESS);
1868: }

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

1873:   Not Collective

1875:   Input Parameter:
1876: . s - the `PetscSection`

1878:   Output Parameters:
1879: + start - the minimum offset
1880: - end   - one more than the maximum offset

1882:   Level: intermediate

1884: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1885: @*/
1886: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1887: {
1888:   PetscInt os = 0, oe = 0, pStart, pEnd, p;

1890:   PetscFunctionBegin;
1892:   if (s->atlasOff) {
1893:     os = s->atlasOff[0];
1894:     oe = s->atlasOff[0];
1895:   }
1896:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1897:   for (p = 0; p < pEnd - pStart; ++p) {
1898:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

1900:     if (off >= 0) {
1901:       os = PetscMin(os, off);
1902:       oe = PetscMax(oe, off + dof);
1903:     }
1904:   }
1905:   if (start) *start = os;
1906:   if (end) *end = oe;
1907:   PetscFunctionReturn(PETSC_SUCCESS);
1908: }

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

1913:   Collective

1915:   Input Parameters:
1916: + s      - the `PetscSection`
1917: . len    - the number of subfields
1918: - fields - the subfield numbers

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

1923:   Level: advanced

1925:   Notes:
1926:   The chart of `subs` is the same as the chart of `s`

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

1930: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1931: @*/
1932: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1933: {
1934:   PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;

1936:   PetscFunctionBegin;
1937:   if (!len) PetscFunctionReturn(PETSC_SUCCESS);
1939:   PetscAssertPointer(fields, 3);
1940:   PetscAssertPointer(subs, 4);
1941:   PetscCall(PetscSectionGetNumFields(s, &nF));
1942:   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);
1943:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
1944:   PetscCall(PetscSectionSetNumFields(*subs, len));
1945:   for (f = 0; f < len; ++f) {
1946:     const char     *name    = NULL;
1947:     PetscInt        numComp = 0;
1948:     PetscSectionSym sym;

1950:     PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
1951:     PetscCall(PetscSectionSetFieldName(*subs, f, name));
1952:     PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
1953:     PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
1954:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1955:       PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
1956:       PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
1957:     }
1958:     PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
1959:     PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
1960:   }
1961:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1962:   PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
1963:   for (p = pStart; p < pEnd; ++p) {
1964:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

1966:     for (f = 0; f < len; ++f) {
1967:       PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
1968:       PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
1969:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
1970:       if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
1971:       dof += fdof;
1972:       cdof += cfdof;
1973:     }
1974:     PetscCall(PetscSectionSetDof(*subs, p, dof));
1975:     if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
1976:     maxCdof = PetscMax(cdof, maxCdof);
1977:   }
1978:   PetscCall(PetscSectionSetUp(*subs));
1979:   if (maxCdof) {
1980:     PetscInt *indices;

1982:     PetscCall(PetscMalloc1(maxCdof, &indices));
1983:     for (p = pStart; p < pEnd; ++p) {
1984:       PetscInt cdof;

1986:       PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
1987:       if (cdof) {
1988:         const PetscInt *oldIndices = NULL;
1989:         PetscInt        fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

1991:         for (f = 0; f < len; ++f) {
1992:           PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
1993:           PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
1994:           PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
1995:           PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
1996:           for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1997:           numConst += cfdof;
1998:           fOff += fdof;
1999:         }
2000:         PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2001:       }
2002:     }
2003:     PetscCall(PetscFree(indices));
2004:   }
2005:   PetscFunctionReturn(PETSC_SUCCESS);
2006: }

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

2011:   Collective

2013:   Input Parameters:
2014: + s   - the input sections
2015: - len - the number of input sections

2017:   Output Parameter:
2018: . supers - the supersection

2020:   Level: advanced

2022:   Notes:
2023:   The section offsets now refer to a new, larger vector.

2025:   Developer Notes:
2026:   Needs to explain how the sections are composed

2028: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2029: @*/
2030: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2031: {
2032:   PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

2034:   PetscFunctionBegin;
2035:   if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2036:   for (i = 0; i < len; ++i) {
2037:     PetscInt nf, pStarti, pEndi;

2039:     PetscCall(PetscSectionGetNumFields(s[i], &nf));
2040:     PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2041:     pStart = PetscMin(pStart, pStarti);
2042:     pEnd   = PetscMax(pEnd, pEndi);
2043:     Nf += nf;
2044:   }
2045:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2046:   PetscCall(PetscSectionSetNumFields(*supers, Nf));
2047:   for (i = 0, f = 0; i < len; ++i) {
2048:     PetscInt nf, fi, ci;

2050:     PetscCall(PetscSectionGetNumFields(s[i], &nf));
2051:     for (fi = 0; fi < nf; ++fi, ++f) {
2052:       const char *name    = NULL;
2053:       PetscInt    numComp = 0;

2055:       PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2056:       PetscCall(PetscSectionSetFieldName(*supers, f, name));
2057:       PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2058:       PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2059:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2060:         PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2061:         PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2062:       }
2063:     }
2064:   }
2065:   PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2066:   for (p = pStart; p < pEnd; ++p) {
2067:     PetscInt dof = 0, cdof = 0;

2069:     for (i = 0, f = 0; i < len; ++i) {
2070:       PetscInt nf, fi, pStarti, pEndi;
2071:       PetscInt fdof = 0, cfdof = 0;

2073:       PetscCall(PetscSectionGetNumFields(s[i], &nf));
2074:       PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2075:       if ((p < pStarti) || (p >= pEndi)) continue;
2076:       for (fi = 0; fi < nf; ++fi, ++f) {
2077:         PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2078:         PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2079:         PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2080:         if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2081:         dof += fdof;
2082:         cdof += cfdof;
2083:       }
2084:     }
2085:     PetscCall(PetscSectionSetDof(*supers, p, dof));
2086:     if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2087:     maxCdof = PetscMax(cdof, maxCdof);
2088:   }
2089:   PetscCall(PetscSectionSetUp(*supers));
2090:   if (maxCdof) {
2091:     PetscInt *indices;

2093:     PetscCall(PetscMalloc1(maxCdof, &indices));
2094:     for (p = pStart; p < pEnd; ++p) {
2095:       PetscInt cdof;

2097:       PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2098:       if (cdof) {
2099:         PetscInt dof, numConst = 0, fOff = 0;

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

2105:           PetscCall(PetscSectionGetNumFields(s[i], &nf));
2106:           PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2107:           if ((p < pStarti) || (p >= pEndi)) continue;
2108:           for (fi = 0; fi < nf; ++fi, ++f) {
2109:             PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2110:             PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2111:             PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2112:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2113:             PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2114:             for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2115:             numConst += cfdof;
2116:           }
2117:           PetscCall(PetscSectionGetDof(s[i], p, &dof));
2118:           fOff += dof;
2119:         }
2120:         PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2121:       }
2122:     }
2123:     PetscCall(PetscFree(indices));
2124:   }
2125:   PetscFunctionReturn(PETSC_SUCCESS);
2126: }

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

2133:   PetscFunctionBegin;
2136:   PetscAssertPointer(subs, 4);
2137:   PetscCall(PetscSectionGetNumFields(s, &numFields));
2138:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2139:   if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2140:   for (f = 0; f < numFields; ++f) {
2141:     const char *name    = NULL;
2142:     PetscInt    numComp = 0;

2144:     PetscCall(PetscSectionGetFieldName(s, f, &name));
2145:     PetscCall(PetscSectionSetFieldName(*subs, f, name));
2146:     PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2147:     PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2148:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2149:       PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2150:       PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2151:     }
2152:   }
2153:   /* For right now, we do not try to squeeze the subchart */
2154:   if (subpointMap) {
2155:     PetscCall(ISGetSize(subpointMap, &numSubpoints));
2156:     PetscCall(ISGetIndices(subpointMap, &points));
2157:   }
2158:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2159:   if (renumberPoints) {
2160:     spStart = 0;
2161:     spEnd   = numSubpoints;
2162:   } else {
2163:     PetscCall(ISGetMinMax(subpointMap, &spStart, &spEnd));
2164:     ++spEnd;
2165:   }
2166:   PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2167:   for (p = pStart; p < pEnd; ++p) {
2168:     PetscInt dof, cdof, fdof = 0, cfdof = 0;

2170:     PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2171:     if (subp < 0) continue;
2172:     if (!renumberPoints) subp = p;
2173:     for (f = 0; f < numFields; ++f) {
2174:       PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2175:       PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2176:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2177:       if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2178:     }
2179:     PetscCall(PetscSectionGetDof(s, p, &dof));
2180:     PetscCall(PetscSectionSetDof(*subs, subp, dof));
2181:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2182:     if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2183:   }
2184:   PetscCall(PetscSectionSetUp(*subs));
2185:   /* Change offsets to original offsets */
2186:   for (p = pStart; p < pEnd; ++p) {
2187:     PetscInt off, foff = 0;

2189:     PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2190:     if (subp < 0) continue;
2191:     if (!renumberPoints) subp = p;
2192:     for (f = 0; f < numFields; ++f) {
2193:       PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2194:       PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2195:     }
2196:     PetscCall(PetscSectionGetOffset(s, p, &off));
2197:     PetscCall(PetscSectionSetOffset(*subs, subp, off));
2198:   }
2199:   /* Copy constraint indices */
2200:   for (subp = spStart; subp < spEnd; ++subp) {
2201:     PetscInt cdof;

2203:     PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2204:     if (cdof) {
2205:       for (f = 0; f < numFields; ++f) {
2206:         PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2207:         PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2208:       }
2209:       PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2210:       PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2211:     }
2212:   }
2213:   if (subpointMap) PetscCall(ISRestoreIndices(subpointMap, &points));
2214:   PetscFunctionReturn(PETSC_SUCCESS);
2215: }

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

2220:   Collective

2222:   Input Parameters:
2223: + s           - the `PetscSection`
2224: - subpointMap - a sorted list of points in the original mesh which are in the submesh

2226:   Output Parameter:
2227: . subs - the subsection

2229:   Level: advanced

2231:   Notes:
2232:   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))`

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

2236:   Developer Notes:
2237:   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`

2239: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2240: @*/
2241: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2242: {
2243:   PetscFunctionBegin;
2244:   PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_TRUE, subs));
2245:   PetscFunctionReturn(PETSC_SUCCESS);
2246: }

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

2251:   Collective

2253:   Input Parameters:
2254: + s           - the `PetscSection`
2255: - subpointMap - a sorted list of points in the original mesh which are in the subdomain

2257:   Output Parameter:
2258: . subs - the subsection

2260:   Level: advanced

2262:   Notes:
2263:   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`
2264:   is `[min(subpointMap),max(subpointMap)+1)`

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

2268:   Developer Notes:
2269:   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`

2271: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2272: @*/
2273: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2274: {
2275:   PetscFunctionBegin;
2276:   PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2277:   PetscFunctionReturn(PETSC_SUCCESS);
2278: }

2280: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2281: {
2282:   PetscInt    p;
2283:   PetscMPIInt rank;

2285:   PetscFunctionBegin;
2286:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2287:   PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2288:   PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2289:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
2290:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2291:       PetscInt b;

2293:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2294:       if (s->bcIndices) {
2295:         for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2296:       }
2297:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2298:     } else {
2299:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2300:     }
2301:   }
2302:   PetscCall(PetscViewerFlush(viewer));
2303:   PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2304:   if (s->sym) {
2305:     PetscCall(PetscViewerASCIIPushTab(viewer));
2306:     PetscCall(PetscSectionSymView(s->sym, viewer));
2307:     PetscCall(PetscViewerASCIIPopTab(viewer));
2308:   }
2309:   PetscFunctionReturn(PETSC_SUCCESS);
2310: }

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

2315:   Collective

2317:   Input Parameters:
2318: + A    - the `PetscSection` object to view
2319: . obj  - Optional object that provides the options prefix used for the options
2320: - name - command line option

2322:   Level: intermediate

2324:   Note:
2325:   See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`

2327: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2328: @*/
2329: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2330: {
2331:   PetscFunctionBegin;
2333:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2334:   PetscFunctionReturn(PETSC_SUCCESS);
2335: }

2337: /*@C
2338:   PetscSectionView - Views a `PetscSection`

2340:   Collective

2342:   Input Parameters:
2343: + s      - the `PetscSection` object to view
2344: - viewer - the viewer

2346:   Level: beginner

2348:   Note:
2349:   `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2350:   distribution independent data, such as dofs, offsets, constraint dofs,
2351:   and constraint indices. Points that have negative dofs, for instance,
2352:   are not saved as they represent points owned by other processes.
2353:   Point numbering and rank assignment is currently not stored.
2354:   The saved section can be loaded with `PetscSectionLoad()`.

2356: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2357: @*/
2358: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2359: {
2360:   PetscBool isascii, ishdf5;
2361:   PetscInt  f;

2363:   PetscFunctionBegin;
2365:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2367:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2368:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2369:   if (isascii) {
2370:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2371:     if (s->numFields) {
2372:       PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2373:       for (f = 0; f < s->numFields; ++f) {
2374:         PetscCall(PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
2375:         PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2376:       }
2377:     } else {
2378:       PetscCall(PetscSectionView_ASCII(s, viewer));
2379:     }
2380:   } else if (ishdf5) {
2381: #if PetscDefined(HAVE_HDF5)
2382:     PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2383: #else
2384:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2385: #endif
2386:   }
2387:   PetscFunctionReturn(PETSC_SUCCESS);
2388: }

2390: /*@C
2391:   PetscSectionLoad - Loads a `PetscSection`

2393:   Collective

2395:   Input Parameters:
2396: + s      - the `PetscSection` object to load
2397: - viewer - the viewer

2399:   Level: beginner

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

2409: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2410: @*/
2411: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2412: {
2413:   PetscBool ishdf5;

2415:   PetscFunctionBegin;
2418:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2419:   if (ishdf5) {
2420: #if PetscDefined(HAVE_HDF5)
2421:     PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2422:     PetscFunctionReturn(PETSC_SUCCESS);
2423: #else
2424:     SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2425: #endif
2426:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2427: }

2429: /*@
2430:   PetscSectionResetClosurePermutation - Remove any existing closure permutation

2432:   Input Parameter:
2433: . section - The `PetscSection`

2435:   Level: intermediate

2437: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2438: @*/
2439: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2440: {
2441:   PetscSectionClosurePermVal clVal;

2443:   PetscFunctionBegin;
2444:   if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2445:   kh_foreach_value(section->clHash, clVal, {
2446:     PetscCall(PetscFree(clVal.perm));
2447:     PetscCall(PetscFree(clVal.invPerm));
2448:   });
2449:   kh_destroy(ClPerm, section->clHash);
2450:   section->clHash = NULL;
2451:   PetscFunctionReturn(PETSC_SUCCESS);
2452: }

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

2457:   Not Collective

2459:   Input Parameter:
2460: . s - the `PetscSection`

2462:   Level: beginner

2464: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2465: @*/
2466: PetscErrorCode PetscSectionReset(PetscSection s)
2467: {
2468:   PetscInt f, c;

2470:   PetscFunctionBegin;
2472:   for (f = 0; f < s->numFields; ++f) {
2473:     PetscCall(PetscSectionDestroy(&s->field[f]));
2474:     PetscCall(PetscFree(s->fieldNames[f]));
2475:     for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2476:     PetscCall(PetscFree(s->compNames[f]));
2477:   }
2478:   PetscCall(PetscFree(s->numFieldComponents));
2479:   PetscCall(PetscFree(s->fieldNames));
2480:   PetscCall(PetscFree(s->compNames));
2481:   PetscCall(PetscFree(s->field));
2482:   PetscCall(PetscSectionDestroy(&s->bc));
2483:   PetscCall(PetscFree(s->bcIndices));
2484:   PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2485:   PetscCall(PetscSectionDestroy(&s->clSection));
2486:   PetscCall(ISDestroy(&s->clPoints));
2487:   PetscCall(ISDestroy(&s->perm));
2488:   PetscCall(PetscSectionResetClosurePermutation(s));
2489:   PetscCall(PetscSectionSymDestroy(&s->sym));
2490:   PetscCall(PetscSectionDestroy(&s->clSection));
2491:   PetscCall(ISDestroy(&s->clPoints));
2492:   PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2493:   s->pStart    = -1;
2494:   s->pEnd      = -1;
2495:   s->maxDof    = 0;
2496:   s->setup     = PETSC_FALSE;
2497:   s->numFields = 0;
2498:   s->clObj     = NULL;
2499:   PetscFunctionReturn(PETSC_SUCCESS);
2500: }

2502: /*@
2503:   PetscSectionDestroy - Frees a `PetscSection`

2505:   Not Collective

2507:   Input Parameter:
2508: . s - the `PetscSection`

2510:   Level: beginner

2512: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2513: @*/
2514: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2515: {
2516:   PetscFunctionBegin;
2517:   if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2519:   if (--((PetscObject)(*s))->refct > 0) {
2520:     *s = NULL;
2521:     PetscFunctionReturn(PETSC_SUCCESS);
2522:   }
2523:   PetscCall(PetscSectionReset(*s));
2524:   PetscCall(PetscHeaderDestroy(s));
2525:   PetscFunctionReturn(PETSC_SUCCESS);
2526: }

2528: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2529: {
2530:   const PetscInt p = point - s->pStart;

2532:   PetscFunctionBegin;
2534:   *values = &baseArray[s->atlasOff[p]];
2535:   PetscFunctionReturn(PETSC_SUCCESS);
2536: }

2538: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2539: {
2540:   PetscInt      *array;
2541:   const PetscInt p           = point - s->pStart;
2542:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2543:   PetscInt       cDim        = 0;

2545:   PetscFunctionBegin;
2547:   PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2548:   array = &baseArray[s->atlasOff[p]];
2549:   if (!cDim) {
2550:     if (orientation >= 0) {
2551:       const PetscInt dim = s->atlasDof[p];
2552:       PetscInt       i;

2554:       if (mode == INSERT_VALUES) {
2555:         for (i = 0; i < dim; ++i) array[i] = values[i];
2556:       } else {
2557:         for (i = 0; i < dim; ++i) array[i] += values[i];
2558:       }
2559:     } else {
2560:       PetscInt offset = 0;
2561:       PetscInt j      = -1, field, i;

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

2566:         for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2567:         offset += dim;
2568:       }
2569:     }
2570:   } else {
2571:     if (orientation >= 0) {
2572:       const PetscInt  dim  = s->atlasDof[p];
2573:       PetscInt        cInd = 0, i;
2574:       const PetscInt *cDof;

2576:       PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2577:       if (mode == INSERT_VALUES) {
2578:         for (i = 0; i < dim; ++i) {
2579:           if ((cInd < cDim) && (i == cDof[cInd])) {
2580:             ++cInd;
2581:             continue;
2582:           }
2583:           array[i] = values[i];
2584:         }
2585:       } else {
2586:         for (i = 0; i < dim; ++i) {
2587:           if ((cInd < cDim) && (i == cDof[cInd])) {
2588:             ++cInd;
2589:             continue;
2590:           }
2591:           array[i] += values[i];
2592:         }
2593:       }
2594:     } else {
2595:       const PetscInt *cDof;
2596:       PetscInt        offset  = 0;
2597:       PetscInt        cOffset = 0;
2598:       PetscInt        j       = 0, field;

2600:       PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2601:       for (field = 0; field < s->numFields; ++field) {
2602:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
2603:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2604:         const PetscInt sDim = dim - tDim;
2605:         PetscInt       cInd = 0, i, k;

2607:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2608:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2609:             ++cInd;
2610:             continue;
2611:           }
2612:           array[j] = values[k];
2613:         }
2614:         offset += dim;
2615:         cOffset += dim - tDim;
2616:       }
2617:     }
2618:   }
2619:   PetscFunctionReturn(PETSC_SUCCESS);
2620: }

2622: /*@C
2623:   PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs

2625:   Not Collective

2627:   Input Parameter:
2628: . s - The `PetscSection`

2630:   Output Parameter:
2631: . hasConstraints - flag indicating that the section has constrained dofs

2633:   Level: intermediate

2635: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2636: @*/
2637: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2638: {
2639:   PetscFunctionBegin;
2641:   PetscAssertPointer(hasConstraints, 2);
2642:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2643:   PetscFunctionReturn(PETSC_SUCCESS);
2644: }

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

2649:   Not Collective

2651:   Input Parameters:
2652: + s     - The `PetscSection`
2653: - point - The point

2655:   Output Parameter:
2656: . indices - The constrained dofs

2658:   Level: intermediate

2660:   Fortran Notes:
2661:   Use `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`

2663: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2664: @*/
2665: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2666: {
2667:   PetscFunctionBegin;
2669:   if (s->bc) {
2670:     PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2671:   } else *indices = NULL;
2672:   PetscFunctionReturn(PETSC_SUCCESS);
2673: }

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

2678:   Not Collective

2680:   Input Parameters:
2681: + s       - The `PetscSection`
2682: . point   - The point
2683: - indices - The constrained dofs

2685:   Level: intermediate

2687:   Fortran Notes:
2688:   Use `PetscSectionSetConstraintIndicesF90()`

2690: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2691: @*/
2692: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2693: {
2694:   PetscFunctionBegin;
2696:   if (s->bc) {
2697:     const PetscInt dof  = s->atlasDof[point];
2698:     const PetscInt cdof = s->bc->atlasDof[point];
2699:     PetscInt       d;

2701:     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]);
2702:     PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
2703:   }
2704:   PetscFunctionReturn(PETSC_SUCCESS);
2705: }

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

2710:   Not Collective

2712:   Input Parameters:
2713: + s     - The `PetscSection`
2714: . field - The field number
2715: - point - The point

2717:   Output Parameter:
2718: . indices - The constrained dofs sorted in ascending order

2720:   Level: intermediate

2722:   Note:
2723:   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()`.

2725:   Fortran Notes:
2726:   Use `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`

2728: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2729: @*/
2730: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2731: {
2732:   PetscFunctionBegin;
2734:   PetscAssertPointer(indices, 4);
2735:   PetscSectionCheckValidField(field, s->numFields);
2736:   PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
2737:   PetscFunctionReturn(PETSC_SUCCESS);
2738: }

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

2743:   Not Collective

2745:   Input Parameters:
2746: + s       - The `PetscSection`
2747: . point   - The point
2748: . field   - The field number
2749: - indices - The constrained dofs

2751:   Level: intermediate

2753:   Fortran Notes:
2754:   Use `PetscSectionSetFieldConstraintIndicesF90()`

2756: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2757: @*/
2758: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2759: {
2760:   PetscFunctionBegin;
2762:   if (PetscDefined(USE_DEBUG)) {
2763:     PetscInt nfdof;

2765:     PetscCall(PetscSectionGetFieldConstraintDof(s, point, field, &nfdof));
2766:     if (nfdof) PetscAssertPointer(indices, 4);
2767:   }
2768:   PetscSectionCheckValidField(field, s->numFields);
2769:   PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
2770:   PetscFunctionReturn(PETSC_SUCCESS);
2771: }

2773: /*@
2774:   PetscSectionPermute - Reorder the section according to the input point permutation

2776:   Collective

2778:   Input Parameters:
2779: + section     - The `PetscSection` object
2780: - permutation - The point permutation, old point p becomes new point perm[p]

2782:   Output Parameter:
2783: . sectionNew - The permuted `PetscSection`

2785:   Level: intermediate

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

2790:   Compare to `PetscSectionSetPermutation()`

2792: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
2793: @*/
2794: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2795: {
2796:   PetscSection    s = section, sNew;
2797:   const PetscInt *perm;
2798:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;

2800:   PetscFunctionBegin;
2803:   PetscAssertPointer(sectionNew, 3);
2804:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
2805:   PetscCall(PetscSectionGetNumFields(s, &numFields));
2806:   if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
2807:   for (f = 0; f < numFields; ++f) {
2808:     const char *name;
2809:     PetscInt    numComp;

2811:     PetscCall(PetscSectionGetFieldName(s, f, &name));
2812:     PetscCall(PetscSectionSetFieldName(sNew, f, name));
2813:     PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2814:     PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
2815:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
2816:       PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2817:       PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
2818:     }
2819:   }
2820:   PetscCall(ISGetLocalSize(permutation, &numPoints));
2821:   PetscCall(ISGetIndices(permutation, &perm));
2822:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2823:   PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
2824:   PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
2825:   for (p = pStart; p < pEnd; ++p) {
2826:     PetscInt dof, cdof;

2828:     PetscCall(PetscSectionGetDof(s, p, &dof));
2829:     PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
2830:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2831:     if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
2832:     for (f = 0; f < numFields; ++f) {
2833:       PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
2834:       PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
2835:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2836:       if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
2837:     }
2838:   }
2839:   PetscCall(PetscSectionSetUp(sNew));
2840:   for (p = pStart; p < pEnd; ++p) {
2841:     const PetscInt *cind;
2842:     PetscInt        cdof;

2844:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2845:     if (cdof) {
2846:       PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
2847:       PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
2848:     }
2849:     for (f = 0; f < numFields; ++f) {
2850:       PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2851:       if (cdof) {
2852:         PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
2853:         PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
2854:       }
2855:     }
2856:   }
2857:   PetscCall(ISRestoreIndices(permutation, &perm));
2858:   *sectionNew = sNew;
2859:   PetscFunctionReturn(PETSC_SUCCESS);
2860: }

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

2865:   Collective

2867:   Input Parameters:
2868: + section   - The `PetscSection`
2869: . obj       - A `PetscObject` which serves as the key for this index
2870: . clSection - `PetscSection` giving the size of the closure of each point
2871: - clPoints  - `IS` giving the points in each closure

2873:   Level: advanced

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

2878:   Developer Notes:
2879:   The information provided here is completely opaque

2881: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2882: @*/
2883: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2884: {
2885:   PetscFunctionBegin;
2889:   if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
2890:   section->clObj = obj;
2891:   PetscCall(PetscObjectReference((PetscObject)clSection));
2892:   PetscCall(PetscObjectReference((PetscObject)clPoints));
2893:   PetscCall(PetscSectionDestroy(&section->clSection));
2894:   PetscCall(ISDestroy(&section->clPoints));
2895:   section->clSection = clSection;
2896:   section->clPoints  = clPoints;
2897:   PetscFunctionReturn(PETSC_SUCCESS);
2898: }

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

2903:   Collective

2905:   Input Parameters:
2906: + section - The `PetscSection`
2907: - obj     - A `PetscObject` which serves as the key for this index

2909:   Output Parameters:
2910: + clSection - `PetscSection` giving the size of the closure of each point
2911: - clPoints  - `IS` giving the points in each closure

2913:   Level: advanced

2915: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2916: @*/
2917: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2918: {
2919:   PetscFunctionBegin;
2920:   if (section->clObj == obj) {
2921:     if (clSection) *clSection = section->clSection;
2922:     if (clPoints) *clPoints = section->clPoints;
2923:   } else {
2924:     if (clSection) *clSection = NULL;
2925:     if (clPoints) *clPoints = NULL;
2926:   }
2927:   PetscFunctionReturn(PETSC_SUCCESS);
2928: }

2930: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2931: {
2932:   PetscInt                    i;
2933:   khiter_t                    iter;
2934:   int                         new_entry;
2935:   PetscSectionClosurePermKey  key = {depth, clSize};
2936:   PetscSectionClosurePermVal *val;

2938:   PetscFunctionBegin;
2939:   if (section->clObj != obj) {
2940:     PetscCall(PetscSectionDestroy(&section->clSection));
2941:     PetscCall(ISDestroy(&section->clPoints));
2942:   }
2943:   section->clObj = obj;
2944:   if (!section->clHash) PetscCall(PetscClPermCreate(&section->clHash));
2945:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2946:   val  = &kh_val(section->clHash, iter);
2947:   if (!new_entry) {
2948:     PetscCall(PetscFree(val->perm));
2949:     PetscCall(PetscFree(val->invPerm));
2950:   }
2951:   if (mode == PETSC_COPY_VALUES) {
2952:     PetscCall(PetscMalloc1(clSize, &val->perm));
2953:     PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
2954:   } else if (mode == PETSC_OWN_POINTER) {
2955:     val->perm = clPerm;
2956:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2957:   PetscCall(PetscMalloc1(clSize, &val->invPerm));
2958:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2959:   PetscFunctionReturn(PETSC_SUCCESS);
2960: }

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

2965:   Not Collective

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

2973:   Level: intermediate

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

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

2983: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2984: @*/
2985: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2986: {
2987:   const PetscInt *clPerm = NULL;
2988:   PetscInt        clSize = 0;

2990:   PetscFunctionBegin;
2991:   if (perm) {
2992:     PetscCall(ISGetLocalSize(perm, &clSize));
2993:     PetscCall(ISGetIndices(perm, &clPerm));
2994:   }
2995:   PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
2996:   if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
2997:   PetscFunctionReturn(PETSC_SUCCESS);
2998: }

3000: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3001: {
3002:   PetscFunctionBegin;
3003:   if (section->clObj == obj) {
3004:     PetscSectionClosurePermKey k = {depth, size};
3005:     PetscSectionClosurePermVal v;

3007:     PetscCall(PetscClPermGet(section->clHash, k, &v));
3008:     if (perm) *perm = v.perm;
3009:   } else {
3010:     if (perm) *perm = NULL;
3011:   }
3012:   PetscFunctionReturn(PETSC_SUCCESS);
3013: }

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

3018:   Not Collective

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

3026:   Output Parameter:
3027: . perm - The dof closure permutation

3029:   Level: intermediate

3031:   Note:
3032:   The user must destroy the `IS` that is returned.

3034: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3035: @*/
3036: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3037: {
3038:   const PetscInt *clPerm = NULL;

3040:   PetscFunctionBegin;
3041:   PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3042:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3043:   PetscFunctionReturn(PETSC_SUCCESS);
3044: }

3046: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3047: {
3048:   PetscFunctionBegin;
3049:   if (section->clObj == obj && section->clHash) {
3050:     PetscSectionClosurePermKey k = {depth, size};
3051:     PetscSectionClosurePermVal v;
3052:     PetscCall(PetscClPermGet(section->clHash, k, &v));
3053:     if (perm) *perm = v.invPerm;
3054:   } else {
3055:     if (perm) *perm = NULL;
3056:   }
3057:   PetscFunctionReturn(PETSC_SUCCESS);
3058: }

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

3063:   Not Collective

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

3071:   Output Parameter:
3072: . perm - The dof closure permutation

3074:   Level: intermediate

3076:   Note:
3077:   The user must destroy the `IS` that is returned.

3079: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3080: @*/
3081: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3082: {
3083:   const PetscInt *clPerm = NULL;

3085:   PetscFunctionBegin;
3086:   PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3087:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3088:   PetscFunctionReturn(PETSC_SUCCESS);
3089: }

3091: /*@
3092:   PetscSectionGetField - Get the `PetscSection` associated with a single field

3094:   Input Parameters:
3095: + s     - The `PetscSection`
3096: - field - The field number

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

3101:   Level: intermediate

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

3106: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3107: @*/
3108: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3109: {
3110:   PetscFunctionBegin;
3112:   PetscAssertPointer(subs, 3);
3113:   PetscSectionCheckValidField(field, s->numFields);
3114:   *subs = s->field[field];
3115:   PetscFunctionReturn(PETSC_SUCCESS);
3116: }

3118: PetscClassId      PETSC_SECTION_SYM_CLASSID;
3119: PetscFunctionList PetscSectionSymList = NULL;

3121: /*@
3122:   PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.

3124:   Collective

3126:   Input Parameter:
3127: . comm - the MPI communicator

3129:   Output Parameter:
3130: . sym - pointer to the new set of symmetries

3132:   Level: developer

3134: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3135: @*/
3136: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3137: {
3138:   PetscFunctionBegin;
3139:   PetscAssertPointer(sym, 2);
3140:   PetscCall(ISInitializePackage());
3141:   PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3142:   PetscFunctionReturn(PETSC_SUCCESS);
3143: }

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

3148:   Collective

3150:   Input Parameters:
3151: + sym    - The section symmetry object
3152: - method - The name of the section symmetry type

3154:   Level: developer

3156: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3157: @*/
3158: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3159: {
3160:   PetscErrorCode (*r)(PetscSectionSym);
3161:   PetscBool match;

3163:   PetscFunctionBegin;
3165:   PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3166:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

3168:   PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3169:   PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3170:   PetscTryTypeMethod(sym, destroy);
3171:   sym->ops->destroy = NULL;

3173:   PetscCall((*r)(sym));
3174:   PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3175:   PetscFunctionReturn(PETSC_SUCCESS);
3176: }

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

3181:   Not Collective

3183:   Input Parameter:
3184: . sym - The section symmetry

3186:   Output Parameter:
3187: . type - The index set type name

3189:   Level: developer

3191: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3192: @*/
3193: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3194: {
3195:   PetscFunctionBegin;
3197:   PetscAssertPointer(type, 2);
3198:   *type = ((PetscObject)sym)->type_name;
3199:   PetscFunctionReturn(PETSC_SUCCESS);
3200: }

3202: /*@C
3203:   PetscSectionSymRegister - Registers a new section symmetry implementation

3205:   Not Collective

3207:   Input Parameters:
3208: + sname    - The name of a new user-defined creation routine
3209: - function - The creation routine itself

3211:   Level: developer

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

3216: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3217: @*/
3218: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3219: {
3220:   PetscFunctionBegin;
3221:   PetscCall(ISInitializePackage());
3222:   PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3223:   PetscFunctionReturn(PETSC_SUCCESS);
3224: }

3226: /*@
3227:   PetscSectionSymDestroy - Destroys a section symmetry.

3229:   Collective

3231:   Input Parameter:
3232: . sym - the section symmetry

3234:   Level: developer

3236: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3237: @*/
3238: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3239: {
3240:   SymWorkLink link, next;

3242:   PetscFunctionBegin;
3243:   if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3245:   if (--((PetscObject)(*sym))->refct > 0) {
3246:     *sym = NULL;
3247:     PetscFunctionReturn(PETSC_SUCCESS);
3248:   }
3249:   if ((*sym)->ops->destroy) PetscCall((*(*sym)->ops->destroy)(*sym));
3250:   PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3251:   for (link = (*sym)->workin; link; link = next) {
3252:     PetscInt    **perms = (PetscInt **)link->perms;
3253:     PetscScalar **rots  = (PetscScalar **)link->rots;
3254:     PetscCall(PetscFree2(perms, rots));
3255:     next = link->next;
3256:     PetscCall(PetscFree(link));
3257:   }
3258:   (*sym)->workin = NULL;
3259:   PetscCall(PetscHeaderDestroy(sym));
3260:   PetscFunctionReturn(PETSC_SUCCESS);
3261: }

3263: /*@C
3264:   PetscSectionSymView - Displays a section symmetry

3266:   Collective

3268:   Input Parameters:
3269: + sym    - the index set
3270: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.

3272:   Level: developer

3274: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3275: @*/
3276: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3277: {
3278:   PetscFunctionBegin;
3280:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3282:   PetscCheckSameComm(sym, 1, viewer, 2);
3283:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3284:   PetscTryTypeMethod(sym, view, viewer);
3285:   PetscFunctionReturn(PETSC_SUCCESS);
3286: }

3288: /*@
3289:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

3291:   Collective

3293:   Input Parameters:
3294: + section - the section describing data layout
3295: - sym     - the symmetry describing the affect of orientation on the access of the data

3297:   Level: developer

3299: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3300: @*/
3301: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3302: {
3303:   PetscFunctionBegin;
3305:   PetscCall(PetscSectionSymDestroy(&(section->sym)));
3306:   if (sym) {
3308:     PetscCheckSameComm(section, 1, sym, 2);
3309:     PetscCall(PetscObjectReference((PetscObject)sym));
3310:   }
3311:   section->sym = sym;
3312:   PetscFunctionReturn(PETSC_SUCCESS);
3313: }

3315: /*@
3316:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

3318:   Not Collective

3320:   Input Parameter:
3321: . section - the section describing data layout

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

3326:   Level: developer

3328: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3329: @*/
3330: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3331: {
3332:   PetscFunctionBegin;
3334:   *sym = section->sym;
3335:   PetscFunctionReturn(PETSC_SUCCESS);
3336: }

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

3341:   Collective

3343:   Input Parameters:
3344: + section - the section describing data layout
3345: . field   - the field number
3346: - sym     - the symmetry describing the affect of orientation on the access of the data

3348:   Level: developer

3350: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3351: @*/
3352: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3353: {
3354:   PetscFunctionBegin;
3356:   PetscSectionCheckValidField(field, section->numFields);
3357:   PetscCall(PetscSectionSetSym(section->field[field], sym));
3358:   PetscFunctionReturn(PETSC_SUCCESS);
3359: }

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

3364:   Collective

3366:   Input Parameters:
3367: + section - the section describing data layout
3368: - field   - the field number

3370:   Output Parameter:
3371: . sym - the symmetry describing the affect of orientation on the access of the data

3373:   Level: developer

3375: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3376: @*/
3377: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3378: {
3379:   PetscFunctionBegin;
3381:   PetscSectionCheckValidField(field, section->numFields);
3382:   *sym = section->field[field]->sym;
3383:   PetscFunctionReturn(PETSC_SUCCESS);
3384: }

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

3389:   Not Collective

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

3398:   Output Parameters:
3399: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3400: - rots  - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3401:     identity).

3403:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3404: .vb
3405:      const PetscInt    **perms;
3406:      const PetscScalar **rots;
3407:      PetscInt            lOffset;

3409:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3410:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3411:        PetscInt           point = points[2*i], dof, sOffset;
3412:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3413:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3415:        PetscSectionGetDof(section,point,&dof);
3416:        PetscSectionGetOffset(section,point,&sOffset);

3418:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
3419:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
3420:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
3421:        lOffset += dof;
3422:      }
3423:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3424: .ve

3426:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3427: .vb
3428:      const PetscInt    **perms;
3429:      const PetscScalar **rots;
3430:      PetscInt            lOffset;

3432:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3433:      for (i = 0, lOffset = 0; i < numPoints; i++) {
3434:        PetscInt           point = points[2*i], dof, sOffset;
3435:        const PetscInt    *perm  = perms ? perms[i] : NULL;
3436:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

3438:        PetscSectionGetDof(section,point,&dof);
3439:        PetscSectionGetOffset(section,point,&sOff);

3441:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3442:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
3443:        offset += dof;
3444:      }
3445:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3446: .ve

3448:   Level: developer

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

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

3455: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3456: @*/
3457: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3458: {
3459:   PetscSectionSym sym;

3461:   PetscFunctionBegin;
3463:   if (numPoints) PetscAssertPointer(points, 3);
3464:   if (perms) *perms = NULL;
3465:   if (rots) *rots = NULL;
3466:   sym = section->sym;
3467:   if (sym && (perms || rots)) {
3468:     SymWorkLink link;

3470:     if (sym->workin) {
3471:       link        = sym->workin;
3472:       sym->workin = sym->workin->next;
3473:     } else {
3474:       PetscCall(PetscNew(&link));
3475:     }
3476:     if (numPoints > link->numPoints) {
3477:       PetscInt    **perms = (PetscInt **)link->perms;
3478:       PetscScalar **rots  = (PetscScalar **)link->rots;
3479:       PetscCall(PetscFree2(perms, rots));
3480:       PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3481:       link->numPoints = numPoints;
3482:     }
3483:     link->next   = sym->workout;
3484:     sym->workout = link;
3485:     PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3486:     PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3487:     PetscCall((*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots));
3488:     if (perms) *perms = link->perms;
3489:     if (rots) *rots = link->rots;
3490:   }
3491:   PetscFunctionReturn(PETSC_SUCCESS);
3492: }

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

3497:   Not Collective

3499:   Input Parameters:
3500: + section   - the section
3501: . numPoints - the number of points
3502: . points    - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3503:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3504:     context, see `DMPlexGetConeOrientation()`).
3505: . perms     - The permutations for the given orientations: set to `NULL` at conclusion
3506: - rots      - The field rotations symmetries for the given orientations: set to `NULL` at conclusion

3508:   Level: developer

3510: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3511: @*/
3512: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3513: {
3514:   PetscSectionSym sym;

3516:   PetscFunctionBegin;
3518:   sym = section->sym;
3519:   if (sym && (perms || rots)) {
3520:     SymWorkLink *p, link;

3522:     for (p = &sym->workout; (link = *p); p = &link->next) {
3523:       if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3524:         *p          = link->next;
3525:         link->next  = sym->workin;
3526:         sym->workin = link;
3527:         if (perms) *perms = NULL;
3528:         if (rots) *rots = NULL;
3529:         PetscFunctionReturn(PETSC_SUCCESS);
3530:       }
3531:     }
3532:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3533:   }
3534:   PetscFunctionReturn(PETSC_SUCCESS);
3535: }

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

3540:   Not Collective

3542:   Input Parameters:
3543: + section   - the section
3544: . field     - the field of the section
3545: . numPoints - the number of points
3546: - points    - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3547:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3548:     context, see `DMPlexGetConeOrientation()`).

3550:   Output Parameters:
3551: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3552: - rots  - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3553:     identity).

3555:   Level: developer

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

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

3562: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3563: @*/
3564: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3565: {
3566:   PetscFunctionBegin;
3568:   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);
3569:   PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3570:   PetscFunctionReturn(PETSC_SUCCESS);
3571: }

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

3576:   Not Collective

3578:   Input Parameters:
3579: + section   - the section
3580: . field     - the field number
3581: . numPoints - the number of points
3582: . points    - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3583:     arbitrary integer: its interpretation is up to sym.  Orientations are used by `DM`: for their interpretation in that
3584:     context, see `DMPlexGetConeOrientation()`).
3585: . perms     - The permutations for the given orientations: set to NULL at conclusion
3586: - rots      - The field rotations symmetries for the given orientations: set to NULL at conclusion

3588:   Level: developer

3590: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3591: @*/
3592: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3593: {
3594:   PetscFunctionBegin;
3596:   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);
3597:   PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3598:   PetscFunctionReturn(PETSC_SUCCESS);
3599: }

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

3604:   Not Collective

3606:   Input Parameter:
3607: . sym - the `PetscSectionSym`

3609:   Output Parameter:
3610: . nsym - the equivalent symmetries

3612:   Level: developer

3614: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3615: @*/
3616: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3617: {
3618:   PetscFunctionBegin;
3621:   PetscTryTypeMethod(sym, copy, nsym);
3622:   PetscFunctionReturn(PETSC_SUCCESS);
3623: }

3625: /*@
3626:   PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`

3628:   Collective

3630:   Input Parameters:
3631: + sym         - the `PetscSectionSym`
3632: - migrationSF - the distribution map from roots to leaves

3634:   Output Parameter:
3635: . dsym - the redistributed symmetries

3637:   Level: developer

3639: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3640: @*/
3641: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3642: {
3643:   PetscFunctionBegin;
3646:   PetscAssertPointer(dsym, 3);
3647:   PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3648:   PetscFunctionReturn(PETSC_SUCCESS);
3649: }

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

3654:   Not Collective

3656:   Input Parameter:
3657: . s - the global `PetscSection`

3659:   Output Parameter:
3660: . flg - the flag

3662:   Level: developer

3664: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3665: @*/
3666: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3667: {
3668:   PetscFunctionBegin;
3670:   *flg = s->useFieldOff;
3671:   PetscFunctionReturn(PETSC_SUCCESS);
3672: }

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

3677:   Not Collective

3679:   Input Parameters:
3680: + s   - the global `PetscSection`
3681: - flg - the flag

3683:   Level: developer

3685: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3686: @*/
3687: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3688: {
3689:   PetscFunctionBegin;
3691:   s->useFieldOff = flg;
3692:   PetscFunctionReturn(PETSC_SUCCESS);
3693: }

3695: #define PetscSectionExpandPoints_Loop(TYPE) \
3696:   do { \
3697:     PetscInt i, n, o0, o1, size; \
3698:     TYPE    *a0 = (TYPE *)origArray, *a1; \
3699:     PetscCall(PetscSectionGetStorageSize(s, &size)); \
3700:     PetscCall(PetscMalloc1(size, &a1)); \
3701:     for (i = 0; i < npoints; i++) { \
3702:       PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3703:       PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3704:       PetscCall(PetscSectionGetDof(s, i, &n)); \
3705:       PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \
3706:     } \
3707:     *newArray = (void *)a1; \
3708:   } while (0)

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

3713:   Not Collective

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

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

3725:   Level: developer

3727: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3728: @*/
3729: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3730: {
3731:   PetscSection    s;
3732:   const PetscInt *points_;
3733:   PetscInt        i, n, npoints, pStart, pEnd;
3734:   PetscMPIInt     unitsize;

3736:   PetscFunctionBegin;
3738:   PetscAssertPointer(origArray, 3);
3740:   if (newSection) PetscAssertPointer(newSection, 5);
3741:   if (newArray) PetscAssertPointer(newArray, 6);
3742:   PetscCallMPI(MPI_Type_size(dataType, &unitsize));
3743:   PetscCall(ISGetLocalSize(points, &npoints));
3744:   PetscCall(ISGetIndices(points, &points_));
3745:   PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
3746:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
3747:   PetscCall(PetscSectionSetChart(s, 0, npoints));
3748:   for (i = 0; i < npoints; i++) {
3749:     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);
3750:     PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
3751:     PetscCall(PetscSectionSetDof(s, i, n));
3752:   }
3753:   PetscCall(PetscSectionSetUp(s));
3754:   if (newArray) {
3755:     if (dataType == MPIU_INT) {
3756:       PetscSectionExpandPoints_Loop(PetscInt);
3757:     } else if (dataType == MPIU_SCALAR) {
3758:       PetscSectionExpandPoints_Loop(PetscScalar);
3759:     } else if (dataType == MPIU_REAL) {
3760:       PetscSectionExpandPoints_Loop(PetscReal);
3761:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3762:   }
3763:   if (newSection) {
3764:     *newSection = s;
3765:   } else {
3766:     PetscCall(PetscSectionDestroy(&s));
3767:   }
3768:   PetscCall(ISRestoreIndices(points, &points_));
3769:   PetscFunctionReturn(PETSC_SUCCESS);
3770: }