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 <petscsection.h>
      7: #include <petscsf.h>
      8: #include <petscviewer.h>

     10: PetscClassId PETSC_SECTION_CLASSID;

     12: /*@
     13:   PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.

     15:   Collective

     17:   Input Parameters:
     18: + comm - the MPI communicator
     19: - s    - pointer to the section

     21:   Level: beginner

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

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

     36: .seealso: PetscSection, PetscSectionDestroy()
     37: @*/
     38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
     39: {

     44:   ISInitializePackage();

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

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

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

     73:   Collective

     75:   Input Parameter:
     76: . section - the PetscSection

     78:   Output Parameter:
     79: . newSection - the copy

     81:   Level: intermediate

     83: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
     84: @*/
     85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
     86: {
     87:   PetscSectionSym sym;
     88:   IS              perm;
     89:   PetscInt        numFields, f, c, pStart, pEnd, p;
     90:   PetscErrorCode  ierr;

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

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

    122:     PetscSectionGetDof(section, p, &dof);
    123:     PetscSectionSetDof(newSection, p, dof);
    124:     PetscSectionGetConstraintDof(section, p, &cdof);
    125:     if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
    126:     for (f = 0; f < numFields; ++f) {
    127:       PetscSectionGetFieldDof(section, p, f, &dof);
    128:       PetscSectionSetFieldDof(newSection, p, f, dof);
    129:       if (cdof) {
    130:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
    131:         if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
    132:       }
    133:     }
    134:   }
    135:   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:     PetscSectionGetOffset(section, p, &off);
    142:     PetscSectionSetOffset(newSection, p, off);
    143:     PetscSectionGetConstraintDof(section, p, &cdof);
    144:     if (cdof) {
    145:       PetscSectionGetConstraintIndices(section, p, &cInd);
    146:       PetscSectionSetConstraintIndices(newSection, p, cInd);
    147:       for (f = 0; f < numFields; ++f) {
    148:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
    149:         if (fcdof) {
    150:           PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
    151:           PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
    152:         }
    153:       }
    154:     }
    155:   }
    156:   return(0);
    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: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
    173: @*/
    174: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
    175: {

    181:   PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
    182:   PetscSectionCopy(section, *newSection);
    183:   return(0);
    184: }

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

    189:   Collective on PetscSection

    191:   Input Parameter:
    192: . section - the PetscSection

    194:   Options Database:
    195: . -petscsection_point_major the dof order

    197:   Level: intermediate

    199: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
    200: @*/
    201: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
    202: {

    207:   PetscObjectOptionsBegin((PetscObject) s);
    208:   PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
    209:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
    210:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
    211:   PetscOptionsEnd();
    212:   PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
    213:   return(0);
    214: }

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

    219:   Collective on PetscSection

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

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

    228:   Level: intermediate

    230:   Notes:
    231:   Field names are disregarded.

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

    248:   flg = PETSC_FALSE;

    250:   MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
    251:   if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
    252:     *congruent = PETSC_FALSE;
    253:     return(0);
    254:   }

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

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

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

    272:     PetscSectionGetDof(s1, p, &n1);
    273:     PetscSectionGetDof(s2, p, &n2);
    274:     if (n1 != n2) goto not_congruent;

    276:     PetscSectionGetConstraintDof(s1, p, &ncdof);
    277:     PetscSectionGetConstraintDof(s2, p, &n2);
    278:     if (ncdof != n2) goto not_congruent;

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

    286:   PetscSectionGetNumFields(s1, &nfields);
    287:   PetscSectionGetNumFields(s2, &n2);
    288:   if (nfields != n2) goto not_congruent;

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

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

    300:       PetscSectionGetFieldDof(s1, p, f, &n1);
    301:       PetscSectionGetFieldDof(s2, p, f, &n2);
    302:       if (n1 != n2) goto not_congruent;

    304:       PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
    305:       PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
    306:       if (nfcdof != n2) goto not_congruent;

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

    315:   flg = PETSC_TRUE;
    316: not_congruent:
    317:   MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
    318:   return(0);
    319: }

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

    324:   Not collective

    326:   Input Parameter:
    327: . s - the PetscSection

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

    332:   Level: intermediate

    334: .seealso: PetscSectionSetNumFields()
    335: @*/
    336: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
    337: {
    341:   *numFields = s->numFields;
    342:   return(0);
    343: }

    345: /*@
    346:   PetscSectionSetNumFields - Sets the number of fields.

    348:   Not collective

    350:   Input Parameters:
    351: + s - the PetscSection
    352: - numFields - the number of fields

    354:   Level: intermediate

    356: .seealso: PetscSectionGetNumFields()
    357: @*/
    358: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
    359: {
    360:   PetscInt       f;

    365:   if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
    366:   PetscSectionReset(s);

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

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

    378:     PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
    379:     PetscSNPrintf(name, 64, "Field_%D", f);
    380:     PetscStrallocpy(name, (char **) &s->fieldNames[f]);
    381:     PetscSNPrintf(name, 64, "Component_0");
    382:     PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
    383:     PetscStrallocpy(name, (char **) &s->compNames[f][0]);
    384:   }
    385:   return(0);
    386: }

    388: /*@C
    389:   PetscSectionGetFieldName - Returns the name of a field in the PetscSection

    391:   Not Collective

    393:   Input Parameters:
    394: + s     - the PetscSection
    395: - field - the field number

    397:   Output Parameter:
    398: . fieldName - the field name

    400:   Level: intermediate

    402: .seealso: PetscSectionSetFieldName()
    403: @*/
    404: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
    405: {
    409:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    410:   *fieldName = s->fieldNames[field];
    411:   return(0);
    412: }

    414: /*@C
    415:   PetscSectionSetFieldName - Sets the name of a field in the PetscSection

    417:   Not Collective

    419:   Input Parameters:
    420: + s     - the PetscSection
    421: . field - the field number
    422: - fieldName - the field name

    424:   Level: intermediate

    426: .seealso: PetscSectionGetFieldName()
    427: @*/
    428: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
    429: {

    435:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    436:   PetscFree(s->fieldNames[field]);
    437:   PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
    438:   return(0);
    439: }

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

    444:   Not Collective

    446:   Input Parameters:
    447: + s     - the PetscSection
    448: . field - the field number
    449: . comp  - the component number
    450: - compName - the component name

    452:   Level: intermediate

    454: .seealso: PetscSectionSetComponentName()
    455: @*/
    456: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
    457: {
    461:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    462:   if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
    463:   *compName = s->compNames[field][comp];
    464:   return(0);
    465: }

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

    470:   Not Collective

    472:   Input Parameters:
    473: + s     - the PetscSection
    474: . field - the field number
    475: . comp  - the component number
    476: - compName - the component name

    478:   Level: intermediate

    480: .seealso: PetscSectionGetComponentName()
    481: @*/
    482: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
    483: {

    489:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    490:   if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
    491:   PetscFree(s->compNames[field][comp]);
    492:   PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
    493:   return(0);
    494: }

    496: /*@
    497:   PetscSectionGetFieldComponents - Returns the number of field components for the given field.

    499:   Not collective

    501:   Input Parameters:
    502: + s - the PetscSection
    503: - field - the field number

    505:   Output Parameter:
    506: . numComp - the number of field components

    508:   Level: intermediate

    510: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
    511: @*/
    512: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
    513: {
    517:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    518:   *numComp = s->numFieldComponents[field];
    519:   return(0);
    520: }

    522: /*@
    523:   PetscSectionSetFieldComponents - Sets the number of field components for the given field.

    525:   Not collective

    527:   Input Parameters:
    528: + s - the PetscSection
    529: . field - the field number
    530: - numComp - the number of field components

    532:   Level: intermediate

    534: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
    535: @*/
    536: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
    537: {
    539:   PetscInt c;
    540:   char name[64];

    544:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    545:   if (s->compNames) {
    546:     for (c = 0; c < s->numFieldComponents[field]; ++c) {
    547:       PetscFree(s->compNames[field][c]);
    548:     }
    549:     PetscFree(s->compNames[field]);
    550:   }

    552:   s->numFieldComponents[field] = numComp;
    553:   if (numComp) {
    554:     PetscMalloc1(numComp, (char ***) &s->compNames[field]);
    555:     for (c = 0; c < numComp; ++c) {
    556:       PetscSNPrintf(name, 64, "%D", c);
    557:       PetscStrallocpy(name, (char **) &s->compNames[field][c]);
    558:     }
    559:   }
    560:   return(0);
    561: }

    563: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
    564: {

    568:   if (!s->bc) {
    569:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
    570:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
    571:   }
    572:   return(0);
    573: }

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

    578:   Not collective

    580:   Input Parameter:
    581: . s - the PetscSection

    583:   Output Parameters:
    584: + pStart - the first point
    585: - pEnd - one past the last point

    587:   Level: intermediate

    589: .seealso: PetscSectionSetChart(), PetscSectionCreate()
    590: @*/
    591: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
    592: {
    595:   if (pStart) *pStart = s->pStart;
    596:   if (pEnd)   *pEnd   = s->pEnd;
    597:   return(0);
    598: }

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

    603:   Not collective

    605:   Input Parameters:
    606: + s - the PetscSection
    607: . pStart - the first point
    608: - pEnd - one past the last point

    610:   Level: intermediate

    612: .seealso: PetscSectionGetChart(), PetscSectionCreate()
    613: @*/
    614: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
    615: {
    616:   PetscInt       f;

    621:   /* Cannot Reset() because it destroys field information */
    622:   s->setup = PETSC_FALSE;
    623:   PetscSectionDestroy(&s->bc);
    624:   PetscFree(s->bcIndices);
    625:   PetscFree2(s->atlasDof, s->atlasOff);

    627:   s->pStart = pStart;
    628:   s->pEnd   = pEnd;
    629:   PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
    630:   PetscArrayzero(s->atlasDof, pEnd - pStart);
    631:   for (f = 0; f < s->numFields; ++f) {
    632:     PetscSectionSetChart(s->field[f], pStart, pEnd);
    633:   }
    634:   return(0);
    635: }

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

    640:   Not collective

    642:   Input Parameter:
    643: . s - the PetscSection

    645:   Output Parameters:
    646: . perm - The permutation as an IS

    648:   Level: intermediate

    650: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
    651: @*/
    652: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
    653: {
    657:   return(0);
    658: }

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

    663:   Not collective

    665:   Input Parameters:
    666: + s - the PetscSection
    667: - perm - the permutation of points

    669:   Level: intermediate

    671: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
    672: @*/
    673: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
    674: {

    680:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
    681:   if (s->perm != perm) {
    682:     ISDestroy(&s->perm);
    683:     if (perm) {
    684:       s->perm = perm;
    685:       PetscObjectReference((PetscObject) s->perm);
    686:     }
    687:   }
    688:   return(0);
    689: }

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

    694:   Not collective

    696:   Input Parameter:
    697: . s - the PetscSection

    699:   Output Parameter:
    700: . pm - the flag for point major ordering

    702:   Level: intermediate

    704: .seealso: PetscSectionSetPointMajor()
    705: @*/
    706: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
    707: {
    711:   *pm = s->pointMajor;
    712:   return(0);
    713: }

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

    718:   Not collective

    720:   Input Parameters:
    721: + s  - the PetscSection
    722: - pm - the flag for point major ordering

    724:   Not collective

    726:   Level: intermediate

    728: .seealso: PetscSectionGetPointMajor()
    729: @*/
    730: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
    731: {
    734:   if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
    735:   s->pointMajor = pm;
    736:   return(0);
    737: }

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

    742:   Not collective

    744:   Input Parameters:
    745: + s - the PetscSection
    746: - point - the point

    748:   Output Parameter:
    749: . numDof - the number of dof

    751:   Level: intermediate

    753: .seealso: PetscSectionSetDof(), PetscSectionCreate()
    754: @*/
    755: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
    756: {
    760:   if (PetscDefined(USE_DEBUG)) {
    761:     if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
    762:   }
    763:   *numDof = s->atlasDof[point - s->pStart];
    764:   return(0);
    765: }

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

    770:   Not collective

    772:   Input Parameters:
    773: + s - the PetscSection
    774: . point - the point
    775: - numDof - the number of dof

    777:   Level: intermediate

    779: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
    780: @*/
    781: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
    782: {
    785:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
    786:   s->atlasDof[point - s->pStart] = numDof;
    787:   return(0);
    788: }

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

    793:   Not collective

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

    800:   Level: intermediate

    802: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
    803: @*/
    804: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
    805: {
    808:   if (PetscDefined(USE_DEBUG)) {
    809:     if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
    810:   }
    811:   s->atlasDof[point - s->pStart] += numDof;
    812:   return(0);
    813: }

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

    818:   Not collective

    820:   Input Parameters:
    821: + s - the PetscSection
    822: . point - the point
    823: - field - the field

    825:   Output Parameter:
    826: . numDof - the number of dof

    828:   Level: intermediate

    830: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
    831: @*/
    832: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
    833: {

    838:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    839:   PetscSectionGetDof(s->field[field], point, numDof);
    840:   return(0);
    841: }

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

    846:   Not collective

    848:   Input Parameters:
    849: + s - the PetscSection
    850: . point - the point
    851: . field - the field
    852: - numDof - the number of dof

    854:   Level: intermediate

    856: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
    857: @*/
    858: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
    859: {

    864:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    865:   PetscSectionSetDof(s->field[field], point, numDof);
    866:   return(0);
    867: }

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

    872:   Not collective

    874:   Input Parameters:
    875: + s - the PetscSection
    876: . point - the point
    877: . field - the field
    878: - numDof - the number of dof

    880:   Level: intermediate

    882: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
    883: @*/
    884: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
    885: {

    890:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    891:   PetscSectionAddDof(s->field[field], point, numDof);
    892:   return(0);
    893: }

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

    898:   Not collective

    900:   Input Parameters:
    901: + s - the PetscSection
    902: - point - the point

    904:   Output Parameter:
    905: . numDof - the number of dof which are fixed by constraints

    907:   Level: intermediate

    909: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
    910: @*/
    911: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
    912: {

    918:   if (s->bc) {
    919:     PetscSectionGetDof(s->bc, point, numDof);
    920:   } else *numDof = 0;
    921:   return(0);
    922: }

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

    927:   Not collective

    929:   Input Parameters:
    930: + s - the PetscSection
    931: . point - the point
    932: - numDof - the number of dof which are fixed by constraints

    934:   Level: intermediate

    936: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
    937: @*/
    938: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
    939: {

    944:   if (numDof) {
    945:     PetscSectionCheckConstraints_Static(s);
    946:     PetscSectionSetDof(s->bc, point, numDof);
    947:   }
    948:   return(0);
    949: }

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

    954:   Not collective

    956:   Input Parameters:
    957: + s - the PetscSection
    958: . point - the point
    959: - numDof - the number of additional dof which are fixed by constraints

    961:   Level: intermediate

    963: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
    964: @*/
    965: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
    966: {

    971:   if (numDof) {
    972:     PetscSectionCheckConstraints_Static(s);
    973:     PetscSectionAddDof(s->bc, point, numDof);
    974:   }
    975:   return(0);
    976: }

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

    981:   Not collective

    983:   Input Parameters:
    984: + s - the PetscSection
    985: . point - the point
    986: - field - the field

    988:   Output Parameter:
    989: . numDof - the number of dof which are fixed by constraints

    991:   Level: intermediate

    993: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
    994: @*/
    995: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
    996: {

    1002:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    1003:   PetscSectionGetConstraintDof(s->field[field], point, numDof);
    1004:   return(0);
    1005: }

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

    1010:   Not collective

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

    1018:   Level: intermediate

    1020: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
    1021: @*/
    1022: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
    1023: {

    1028:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    1029:   PetscSectionSetConstraintDof(s->field[field], point, numDof);
    1030:   return(0);
    1031: }

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

    1036:   Not collective

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

    1044:   Level: intermediate

    1046: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
    1047: @*/
    1048: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
    1049: {

    1054:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    1055:   PetscSectionAddConstraintDof(s->field[field], point, numDof);
    1056:   return(0);
    1057: }

    1059: /*@
    1060:   PetscSectionSetUpBC - Setup the subsections describing boundary conditions.

    1062:   Not collective

    1064:   Input Parameter:
    1065: . s - the PetscSection

    1067:   Level: advanced

    1069: .seealso: PetscSectionSetUp(), PetscSectionCreate()
    1070: @*/
    1071: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
    1072: {

    1077:   if (s->bc) {
    1078:     const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;

    1080:     PetscSectionSetUp(s->bc);
    1081:     PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
    1082:   }
    1083:   return(0);
    1084: }

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

    1089:   Not collective

    1091:   Input Parameter:
    1092: . s - the PetscSection

    1094:   Level: intermediate

    1096: .seealso: PetscSectionCreate()
    1097: @*/
    1098: PetscErrorCode PetscSectionSetUp(PetscSection s)
    1099: {
    1100:   const PetscInt *pind   = NULL;
    1101:   PetscInt        offset = 0, foff, p, f;
    1102:   PetscErrorCode  ierr;

    1106:   if (s->setup) return(0);
    1107:   s->setup = PETSC_TRUE;
    1108:   /* Set offsets and field offsets for all points */
    1109:   /*   Assume that all fields have the same chart */
    1110:   if (s->perm) {ISGetIndices(s->perm, &pind);}
    1111:   if (s->pointMajor) {
    1112:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
    1113:       const PetscInt q = pind ? pind[p] : p;

    1115:       /* Set point offset */
    1116:       s->atlasOff[q] = offset;
    1117:       offset        += s->atlasDof[q];
    1118:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[q]);
    1119:       /* Set field offset */
    1120:       for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
    1121:         PetscSection sf = s->field[f];

    1123:         sf->atlasOff[q] = foff;
    1124:         foff += sf->atlasDof[q];
    1125:       }
    1126:     }
    1127:   } else {
    1128:     /* Set field offsets for all points */
    1129:     for (f = 0; f < s->numFields; ++f) {
    1130:       PetscSection sf = s->field[f];

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

    1135:         sf->atlasOff[q] = offset;
    1136:         offset += sf->atlasDof[q];
    1137:       }
    1138:     }
    1139:     /* Disable point offsets since these are unused */
    1140:     for (p = 0; p < s->pEnd - s->pStart; ++p) {
    1141:       s->atlasOff[p] = -1;
    1142:       s->maxDof      = PetscMax(s->maxDof, s->atlasDof[p]);
    1143:     }
    1144:   }
    1145:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
    1146:   /* Setup BC sections */
    1147:   PetscSectionSetUpBC(s);
    1148:   for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
    1149:   return(0);
    1150: }

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

    1155:   Not collective

    1157:   Input Parameters:
    1158: . s - the PetscSection

    1160:   Output Parameter:
    1161: . maxDof - the maximum dof

    1163:   Level: intermediate

    1165: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
    1166: @*/
    1167: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
    1168: {
    1172:   *maxDof = s->maxDof;
    1173:   return(0);
    1174: }

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

    1179:   Not collective

    1181:   Input Parameter:
    1182: . s - the PetscSection

    1184:   Output Parameter:
    1185: . size - the size of an array which can hold all the dofs

    1187:   Level: intermediate

    1189: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
    1190: @*/
    1191: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
    1192: {
    1193:   PetscInt p, n = 0;

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

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

    1206:   Not collective

    1208:   Input Parameter:
    1209: . s - the PetscSection

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

    1214:   Level: intermediate

    1216: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
    1217: @*/
    1218: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
    1219: {
    1220:   PetscInt p, n = 0;

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

    1233: /*@
    1234:   PetscSectionCreateGlobalSection - Create a section describing the global field layout using
    1235:   the local section and an SF describing the section point overlap.

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

    1243:   Output Parameter:
    1244: . gsection - The PetscSection for the global field layout

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

    1248:   Level: intermediate

    1250: .seealso: PetscSectionCreate()
    1251: @*/
    1252: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
    1253: {
    1254:   PetscSection    gs;
    1255:   const PetscInt *pind = NULL;
    1256:   PetscInt       *recv = NULL, *neg = NULL;
    1257:   PetscInt        pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
    1258:   PetscErrorCode  ierr;

    1266:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
    1267:   PetscSectionGetChart(s, &pStart, &pEnd);
    1268:   PetscSectionSetChart(gs, pStart, pEnd);
    1269:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
    1270:   nlocal = nroots;              /* The local/leaf space matches global/root space */
    1271:   /* Must allocate for all points visible to SF, which may be more than this section */
    1272:   if (nroots >= 0) {             /* nroots < 0 means that the graph has not been set, only happens in serial */
    1273:     PetscSFGetLeafRange(sf, NULL, &maxleaf);
    1274:     if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
    1275:     if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
    1276:     PetscMalloc2(nroots,&neg,nlocal,&recv);
    1277:     PetscArrayzero(neg,nroots);
    1278:   }
    1279:   /* Mark all local points with negative dof */
    1280:   for (p = pStart; p < pEnd; ++p) {
    1281:     PetscSectionGetDof(s, p, &dof);
    1282:     PetscSectionSetDof(gs, p, dof);
    1283:     PetscSectionGetConstraintDof(s, p, &cdof);
    1284:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
    1285:     if (neg) neg[p] = -(dof+1);
    1286:   }
    1287:   PetscSectionSetUpBC(gs);
    1288:   if (gs->bcIndices) {PetscArraycpy(gs->bcIndices, s->bcIndices,gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]);}
    1289:   if (nroots >= 0) {
    1290:     PetscArrayzero(recv,nlocal);
    1291:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
    1292:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
    1293:     for (p = pStart; p < pEnd; ++p) {
    1294:       if (recv[p] < 0) {
    1295:         gs->atlasDof[p-pStart] = recv[p];
    1296:         PetscSectionGetDof(s, p, &dof);
    1297:         if (-(recv[p]+1) != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %D for point %D is not the unconstrained %D", -(recv[p]+1), p, dof);
    1298:       }
    1299:     }
    1300:   }
    1301:   /* Calculate new sizes, get process offset, and calculate point offsets */
    1302:   if (s->perm) {ISGetIndices(s->perm, &pind);}
    1303:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
    1304:     const PetscInt q = pind ? pind[p] : p;

    1306:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
    1307:     gs->atlasOff[q] = off;
    1308:     off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
    1309:   }
    1310:   if (!localOffsets) {
    1311:     MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
    1312:     globalOff -= off;
    1313:   }
    1314:   for (p = pStart, off = 0; p < pEnd; ++p) {
    1315:     gs->atlasOff[p-pStart] += globalOff;
    1316:     if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
    1317:   }
    1318:   if (s->perm) {ISRestoreIndices(s->perm, &pind);}
    1319:   /* Put in negative offsets for ghost points */
    1320:   if (nroots >= 0) {
    1321:     PetscArrayzero(recv,nlocal);
    1322:     PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
    1323:     PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
    1324:     for (p = pStart; p < pEnd; ++p) {
    1325:       if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
    1326:     }
    1327:   }
    1328:   PetscFree2(neg,recv);
    1329:   PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
    1330:   *gsection = gs;
    1331:   return(0);
    1332: }

    1334: /*@
    1335:   PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
    1336:   the local section and an SF describing the section point overlap.

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

    1345:   Output Parameter:
    1346:   . gsection - The PetscSection for the global field layout

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

    1350:   Level: advanced

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

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

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

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

    1434:   Collective on comm

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

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

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

    1445:   Level: advanced

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

    1455:   PetscSectionGetChart(s, &pStart, &pEnd);
    1456:   for (p = pStart; p < pEnd; ++p) {
    1457:     PetscInt dof;

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

    1469: /*@
    1470:   PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.

    1472:   Collective on comm

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

    1478:   Output Parameter:
    1479: . layout - The dof layout for the section

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

    1483:   Level: advanced

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

    1495:   PetscSectionGetChart(s, &pStart, &pEnd);
    1496:   for (p = pStart; p < pEnd; ++p) {
    1497:     PetscInt dof,cdof;

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

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

    1513:   Not collective

    1515:   Input Parameters:
    1516: + s - the PetscSection
    1517: - point - the point

    1519:   Output Parameter:
    1520: . offset - the offset

    1522:   Level: intermediate

    1524: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
    1525: @*/
    1526: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
    1527: {
    1531:   if (PetscDefined(USE_DEBUG)) {
    1532:     if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
    1533:   }
    1534:   *offset = s->atlasOff[point - s->pStart];
    1535:   return(0);
    1536: }

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

    1541:   Not collective

    1543:   Input Parameters:
    1544: + s - the PetscSection
    1545: . point - the point
    1546: - offset - the offset

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

    1550:   Level: intermediate

    1552: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
    1553: @*/
    1554: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
    1555: {
    1558:   if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %D should be in [%D, %D)", point, s->pStart, s->pEnd);
    1559:   s->atlasOff[point - s->pStart] = offset;
    1560:   return(0);
    1561: }

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

    1566:   Not collective

    1568:   Input Parameters:
    1569: + s - the PetscSection
    1570: . point - the point
    1571: - field - the field

    1573:   Output Parameter:
    1574: . offset - the offset

    1576:   Level: intermediate

    1578: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
    1579: @*/
    1580: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
    1581: {

    1587:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    1588:   PetscSectionGetOffset(s->field[field], point, offset);
    1589:   return(0);
    1590: }

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

    1595:   Not collective

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

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

    1605:   Level: intermediate

    1607: .seealso: PetscSectionGetFieldOffset(), PetscSectionSetOffset(), PetscSectionCreate(), PetscSectionSetUp()
    1608: @*/
    1609: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
    1610: {

    1615:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    1616:   PetscSectionSetOffset(s->field[field], point, offset);
    1617:   return(0);
    1618: }

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

    1623:   Not collective

    1625:   Input Parameters:
    1626: + s - the PetscSection
    1627: . point - the point
    1628: - field - the field

    1630:   Output Parameter:
    1631: . offset - the offset

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

    1636:   Level: advanced

    1638: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
    1639: @*/
    1640: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
    1641: {
    1642:   PetscInt       off, foff;

    1648:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    1649:   PetscSectionGetOffset(s, point, &off);
    1650:   PetscSectionGetOffset(s->field[field], point, &foff);
    1651:   *offset = foff - off;
    1652:   return(0);
    1653: }

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

    1658:   Not collective

    1660:   Input Parameter:
    1661: . s - the PetscSection

    1663:   Output Parameters:
    1664: + start - the minimum offset
    1665: - end   - one more than the maximum offset

    1667:   Level: intermediate

    1669: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
    1670: @*/
    1671: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
    1672: {
    1673:   PetscInt       os = 0, oe = 0, pStart, pEnd, p;

    1678:   if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
    1679:   PetscSectionGetChart(s, &pStart, &pEnd);
    1680:   for (p = 0; p < pEnd-pStart; ++p) {
    1681:     PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];

    1683:     if (off >= 0) {
    1684:       os = PetscMin(os, off);
    1685:       oe = PetscMax(oe, off+dof);
    1686:     }
    1687:   }
    1688:   if (start) *start = os;
    1689:   if (end)   *end   = oe;
    1690:   return(0);
    1691: }

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

    1696:   Collective on s

    1698:   Input Parameter:
    1699: + s      - the PetscSection
    1700: . len    - the number of subfields
    1701: - fields - the subfield numbers

    1703:   Output Parameter:
    1704: . subs   - the subsection

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

    1708:   Level: advanced

    1710: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
    1711: @*/
    1712: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
    1713: {
    1714:   PetscInt       nF, f, c, pStart, pEnd, p, maxCdof = 0;

    1718:   if (!len) return(0);
    1722:   PetscSectionGetNumFields(s, &nF);
    1723:   if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
    1724:   PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
    1725:   PetscSectionSetNumFields(*subs, len);
    1726:   for (f = 0; f < len; ++f) {
    1727:     const char *name   = NULL;
    1728:     PetscInt   numComp = 0;

    1730:     PetscSectionGetFieldName(s, fields[f], &name);
    1731:     PetscSectionSetFieldName(*subs, f, name);
    1732:     PetscSectionGetFieldComponents(s, fields[f], &numComp);
    1733:     PetscSectionSetFieldComponents(*subs, f, numComp);
    1734:     for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
    1735:       PetscSectionGetComponentName(s, fields[f], c, &name);
    1736:       PetscSectionSetComponentName(*subs, f, c, name);
    1737:     }
    1738:   }
    1739:   PetscSectionGetChart(s, &pStart, &pEnd);
    1740:   PetscSectionSetChart(*subs, pStart, pEnd);
    1741:   for (p = pStart; p < pEnd; ++p) {
    1742:     PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;

    1744:     for (f = 0; f < len; ++f) {
    1745:       PetscSectionGetFieldDof(s, p, fields[f], &fdof);
    1746:       PetscSectionSetFieldDof(*subs, p, f, fdof);
    1747:       PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
    1748:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
    1749:       dof  += fdof;
    1750:       cdof += cfdof;
    1751:     }
    1752:     PetscSectionSetDof(*subs, p, dof);
    1753:     if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
    1754:     maxCdof = PetscMax(cdof, maxCdof);
    1755:   }
    1756:   PetscSectionSetUp(*subs);
    1757:   if (maxCdof) {
    1758:     PetscInt *indices;

    1760:     PetscMalloc1(maxCdof, &indices);
    1761:     for (p = pStart; p < pEnd; ++p) {
    1762:       PetscInt cdof;

    1764:       PetscSectionGetConstraintDof(*subs, p, &cdof);
    1765:       if (cdof) {
    1766:         const PetscInt *oldIndices = NULL;
    1767:         PetscInt       fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;

    1769:         for (f = 0; f < len; ++f) {
    1770:           PetscInt oldFoff = 0, oldf;

    1772:           PetscSectionGetFieldDof(s, p, fields[f], &fdof);
    1773:           PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
    1774:           PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
    1775:           /* This can be sped up if we assume sorted fields */
    1776:           for (oldf = 0; oldf < fields[f]; ++oldf) {
    1777:             PetscInt oldfdof = 0;
    1778:             PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
    1779:             oldFoff += oldfdof;
    1780:           }
    1781:           for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
    1782:           PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
    1783:           numConst += cfdof;
    1784:           fOff     += fdof;
    1785:         }
    1786:         PetscSectionSetConstraintIndices(*subs, p, indices);
    1787:       }
    1788:     }
    1789:     PetscFree(indices);
    1790:   }
    1791:   return(0);
    1792: }

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

    1797:   Collective on s

    1799:   Input Parameters:
    1800: + s     - the input sections
    1801: - len   - the number of input sections

    1803:   Output Parameter:
    1804: . supers - the supersection

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

    1808:   Level: advanced

    1810: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
    1811: @*/
    1812: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
    1813: {
    1814:   PetscInt       Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;

    1818:   if (!len) return(0);
    1819:   for (i = 0; i < len; ++i) {
    1820:     PetscInt nf, pStarti, pEndi;

    1822:     PetscSectionGetNumFields(s[i], &nf);
    1823:     PetscSectionGetChart(s[i], &pStarti, &pEndi);
    1824:     pStart = PetscMin(pStart, pStarti);
    1825:     pEnd   = PetscMax(pEnd,   pEndi);
    1826:     Nf += nf;
    1827:   }
    1828:   PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
    1829:   PetscSectionSetNumFields(*supers, Nf);
    1830:   for (i = 0, f = 0; i < len; ++i) {
    1831:     PetscInt nf, fi, ci;

    1833:     PetscSectionGetNumFields(s[i], &nf);
    1834:     for (fi = 0; fi < nf; ++fi, ++f) {
    1835:       const char *name   = NULL;
    1836:       PetscInt   numComp = 0;

    1838:       PetscSectionGetFieldName(s[i], fi, &name);
    1839:       PetscSectionSetFieldName(*supers, f, name);
    1840:       PetscSectionGetFieldComponents(s[i], fi, &numComp);
    1841:       PetscSectionSetFieldComponents(*supers, f, numComp);
    1842:       for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
    1843:         PetscSectionGetComponentName(s[i], fi, ci, &name);
    1844:         PetscSectionSetComponentName(*supers, f, ci, name);
    1845:       }
    1846:     }
    1847:   }
    1848:   PetscSectionSetChart(*supers, pStart, pEnd);
    1849:   for (p = pStart; p < pEnd; ++p) {
    1850:     PetscInt dof = 0, cdof = 0;

    1852:     for (i = 0, f = 0; i < len; ++i) {
    1853:       PetscInt nf, fi, pStarti, pEndi;
    1854:       PetscInt fdof = 0, cfdof = 0;

    1856:       PetscSectionGetNumFields(s[i], &nf);
    1857:       PetscSectionGetChart(s[i], &pStarti, &pEndi);
    1858:       if ((p < pStarti) || (p >= pEndi)) continue;
    1859:       for (fi = 0; fi < nf; ++fi, ++f) {
    1860:         PetscSectionGetFieldDof(s[i], p, fi, &fdof);
    1861:         PetscSectionAddFieldDof(*supers, p, f, fdof);
    1862:         PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
    1863:         if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
    1864:         dof  += fdof;
    1865:         cdof += cfdof;
    1866:       }
    1867:     }
    1868:     PetscSectionSetDof(*supers, p, dof);
    1869:     if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
    1870:     maxCdof = PetscMax(cdof, maxCdof);
    1871:   }
    1872:   PetscSectionSetUp(*supers);
    1873:   if (maxCdof) {
    1874:     PetscInt *indices;

    1876:     PetscMalloc1(maxCdof, &indices);
    1877:     for (p = pStart; p < pEnd; ++p) {
    1878:       PetscInt cdof;

    1880:       PetscSectionGetConstraintDof(*supers, p, &cdof);
    1881:       if (cdof) {
    1882:         PetscInt dof, numConst = 0, fOff = 0;

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

    1888:           PetscSectionGetNumFields(s[i], &nf);
    1889:           PetscSectionGetChart(s[i], &pStarti, &pEndi);
    1890:           if ((p < pStarti) || (p >= pEndi)) continue;
    1891:           for (fi = 0; fi < nf; ++fi, ++f) {
    1892:             PetscSectionGetFieldDof(s[i], p, fi, &fdof);
    1893:             PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
    1894:             PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
    1895:             for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
    1896:             PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
    1897:             numConst += cfdof;
    1898:           }
    1899:           PetscSectionGetDof(s[i], p, &dof);
    1900:           fOff += dof;
    1901:         }
    1902:         PetscSectionSetConstraintIndices(*supers, p, indices);
    1903:       }
    1904:     }
    1905:     PetscFree(indices);
    1906:   }
    1907:   return(0);
    1908: }

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

    1913:   Collective on s

    1915:   Input Parameters:
    1916: + s           - the PetscSection
    1917: - subpointMap - a sorted list of points in the original mesh which are in the submesh

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

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

    1924:   Level: advanced

    1926: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
    1927: @*/
    1928: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
    1929: {
    1930:   const PetscInt *points = NULL, *indices = NULL;
    1931:   PetscInt       numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;

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

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

    1964:     PetscFindInt(p, numSubpoints, points, &subp);
    1965:     if (subp < 0) continue;
    1966:     for (f = 0; f < numFields; ++f) {
    1967:       PetscSectionGetFieldDof(s, p, f, &fdof);
    1968:       PetscSectionSetFieldDof(*subs, subp, f, fdof);
    1969:       PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
    1970:       if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
    1971:     }
    1972:     PetscSectionGetDof(s, p, &dof);
    1973:     PetscSectionSetDof(*subs, subp, dof);
    1974:     PetscSectionGetConstraintDof(s, p, &cdof);
    1975:     if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
    1976:   }
    1977:   PetscSectionSetUp(*subs);
    1978:   /* Change offsets to original offsets */
    1979:   for (p = pStart; p < pEnd; ++p) {
    1980:     PetscInt off, foff = 0;

    1982:     PetscFindInt(p, numSubpoints, points, &subp);
    1983:     if (subp < 0) continue;
    1984:     for (f = 0; f < numFields; ++f) {
    1985:       PetscSectionGetFieldOffset(s, p, f, &foff);
    1986:       PetscSectionSetFieldOffset(*subs, subp, f, foff);
    1987:     }
    1988:     PetscSectionGetOffset(s, p, &off);
    1989:     PetscSectionSetOffset(*subs, subp, off);
    1990:   }
    1991:   /* Copy constraint indices */
    1992:   for (subp = 0; subp < numSubpoints; ++subp) {
    1993:     PetscInt cdof;

    1995:     PetscSectionGetConstraintDof(*subs, subp, &cdof);
    1996:     if (cdof) {
    1997:       for (f = 0; f < numFields; ++f) {
    1998:         PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
    1999:         PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
    2000:       }
    2001:       PetscSectionGetConstraintIndices(s, points[subp], &indices);
    2002:       PetscSectionSetConstraintIndices(*subs, subp, indices);
    2003:     }
    2004:   }
    2005:   if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
    2006:   return(0);
    2007: }

    2009: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
    2010: {
    2011:   PetscInt       p;
    2012:   PetscMPIInt    rank;

    2016:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
    2017:   PetscViewerASCIIPushSynchronized(viewer);
    2018:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
    2019:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
    2020:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
    2021:       PetscInt b;

    2023:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
    2024:       if (s->bcIndices) {
    2025:         for (b = 0; b < s->bc->atlasDof[p]; ++b) {
    2026:           PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
    2027:         }
    2028:       }
    2029:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
    2030:     } else {
    2031:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
    2032:     }
    2033:   }
    2034:   PetscViewerFlush(viewer);
    2035:   PetscViewerASCIIPopSynchronized(viewer);
    2036:   if (s->sym) {
    2037:     PetscViewerASCIIPushTab(viewer);
    2038:     PetscSectionSymView(s->sym,viewer);
    2039:     PetscViewerASCIIPopTab(viewer);
    2040:   }
    2041:   return(0);
    2042: }

    2044: /*@C
    2045:    PetscSectionViewFromOptions - View from Options

    2047:    Collective on PetscSection

    2049:    Input Parameters:
    2050: +  A - the PetscSection object to view
    2051: .  obj - Optional object
    2052: -  name - command line option

    2054:    Level: intermediate
    2055: .seealso:  PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
    2056: @*/
    2057: PetscErrorCode  PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
    2058: {

    2063:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
    2064:   return(0);
    2065: }

    2067: /*@C
    2068:   PetscSectionView - Views a PetscSection

    2070:   Collective on PetscSection

    2072:   Input Parameters:
    2073: + s - the PetscSection object to view
    2074: - v - the viewer

    2076:   Level: beginner

    2078: .seealso PetscSectionCreate(), PetscSectionDestroy()
    2079: @*/
    2080: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
    2081: {
    2082:   PetscBool      isascii;
    2083:   PetscInt       f;

    2088:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
    2090:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
    2091:   if (isascii) {
    2092:     PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
    2093:     if (s->numFields) {
    2094:       PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
    2095:       for (f = 0; f < s->numFields; ++f) {
    2096:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
    2097:         PetscSectionView_ASCII(s->field[f], viewer);
    2098:       }
    2099:     } else {
    2100:       PetscSectionView_ASCII(s, viewer);
    2101:     }
    2102:   }
    2103:   return(0);
    2104: }

    2106: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
    2107: {
    2109:   PetscSectionClosurePermVal clVal;

    2112:   if (!section->clHash) return(0);
    2113:   kh_foreach_value(section->clHash, clVal, {
    2114:       PetscFree(clVal.perm);
    2115:       PetscFree(clVal.invPerm);
    2116:     });
    2117:   kh_destroy(ClPerm, section->clHash);
    2118:   section->clHash = NULL;
    2119:   return(0);
    2120: }

    2122: /*@
    2123:   PetscSectionReset - Frees all section data.

    2125:   Not collective

    2127:   Input Parameters:
    2128: . s - the PetscSection

    2130:   Level: beginner

    2132: .seealso: PetscSection, PetscSectionCreate()
    2133: @*/
    2134: PetscErrorCode PetscSectionReset(PetscSection s)
    2135: {
    2136:   PetscInt       f, c;

    2141:   for (f = 0; f < s->numFields; ++f) {
    2142:     PetscSectionDestroy(&s->field[f]);
    2143:     PetscFree(s->fieldNames[f]);
    2144:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
    2145:       PetscFree(s->compNames[f][c]);
    2146:     }
    2147:     PetscFree(s->compNames[f]);
    2148:   }
    2149:   PetscFree(s->numFieldComponents);
    2150:   PetscFree(s->fieldNames);
    2151:   PetscFree(s->compNames);
    2152:   PetscFree(s->field);
    2153:   PetscSectionDestroy(&s->bc);
    2154:   PetscFree(s->bcIndices);
    2155:   PetscFree2(s->atlasDof, s->atlasOff);
    2156:   PetscSectionDestroy(&s->clSection);
    2157:   ISDestroy(&s->clPoints);
    2158:   ISDestroy(&s->perm);
    2159:   PetscSectionResetClosurePermutation(s);
    2160:   PetscSectionSymDestroy(&s->sym);
    2161:   PetscSectionDestroy(&s->clSection);
    2162:   ISDestroy(&s->clPoints);

    2164:   s->pStart    = -1;
    2165:   s->pEnd      = -1;
    2166:   s->maxDof    = 0;
    2167:   s->setup     = PETSC_FALSE;
    2168:   s->numFields = 0;
    2169:   s->clObj     = NULL;
    2170:   return(0);
    2171: }

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

    2176:   Not collective

    2178:   Input Parameters:
    2179: . s - the PetscSection

    2181:   Level: beginner

    2183: .seealso: PetscSection, PetscSectionCreate()
    2184: @*/
    2185: PetscErrorCode PetscSectionDestroy(PetscSection *s)
    2186: {

    2190:   if (!*s) return(0);
    2192:   if (--((PetscObject)(*s))->refct > 0) {
    2193:     *s = NULL;
    2194:     return(0);
    2195:   }
    2196:   PetscSectionReset(*s);
    2197:   PetscHeaderDestroy(s);
    2198:   return(0);
    2199: }

    2201: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
    2202: {
    2203:   const PetscInt p = point - s->pStart;

    2207:   *values = &baseArray[s->atlasOff[p]];
    2208:   return(0);
    2209: }

    2211: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
    2212: {
    2213:   PetscInt       *array;
    2214:   const PetscInt p           = point - s->pStart;
    2215:   const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
    2216:   PetscInt       cDim        = 0;

    2221:   PetscSectionGetConstraintDof(s, p, &cDim);
    2222:   array = &baseArray[s->atlasOff[p]];
    2223:   if (!cDim) {
    2224:     if (orientation >= 0) {
    2225:       const PetscInt dim = s->atlasDof[p];
    2226:       PetscInt       i;

    2228:       if (mode == INSERT_VALUES) {
    2229:         for (i = 0; i < dim; ++i) array[i] = values[i];
    2230:       } else {
    2231:         for (i = 0; i < dim; ++i) array[i] += values[i];
    2232:       }
    2233:     } else {
    2234:       PetscInt offset = 0;
    2235:       PetscInt j      = -1, field, i;

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

    2240:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
    2241:         offset += dim;
    2242:       }
    2243:     }
    2244:   } else {
    2245:     if (orientation >= 0) {
    2246:       const PetscInt dim  = s->atlasDof[p];
    2247:       PetscInt       cInd = 0, i;
    2248:       const PetscInt *cDof;

    2250:       PetscSectionGetConstraintIndices(s, point, &cDof);
    2251:       if (mode == INSERT_VALUES) {
    2252:         for (i = 0; i < dim; ++i) {
    2253:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
    2254:           array[i] = values[i];
    2255:         }
    2256:       } else {
    2257:         for (i = 0; i < dim; ++i) {
    2258:           if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
    2259:           array[i] += values[i];
    2260:         }
    2261:       }
    2262:     } else {
    2263:       const PetscInt *cDof;
    2264:       PetscInt       offset  = 0;
    2265:       PetscInt       cOffset = 0;
    2266:       PetscInt       j       = 0, field;

    2268:       PetscSectionGetConstraintIndices(s, point, &cDof);
    2269:       for (field = 0; field < s->numFields; ++field) {
    2270:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
    2271:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
    2272:         const PetscInt sDim = dim - tDim;
    2273:         PetscInt       cInd = 0, i,k;

    2275:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
    2276:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
    2277:           array[j] = values[k];
    2278:         }
    2279:         offset  += dim;
    2280:         cOffset += dim - tDim;
    2281:       }
    2282:     }
    2283:   }
    2284:   return(0);
    2285: }

    2287: /*@C
    2288:   PetscSectionHasConstraints - Determine whether a section has constrained dofs

    2290:   Not collective

    2292:   Input Parameter:
    2293: . s - The PetscSection

    2295:   Output Parameter:
    2296: . hasConstraints - flag indicating that the section has constrained dofs

    2298:   Level: intermediate

    2300: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
    2301: @*/
    2302: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
    2303: {
    2307:   *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
    2308:   return(0);
    2309: }

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

    2314:   Not collective

    2316:   Input Parameters:
    2317: + s     - The PetscSection
    2318: - point - The point

    2320:   Output Parameter:
    2321: . indices - The constrained dofs

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

    2325:   Level: intermediate

    2327: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
    2328: @*/
    2329: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
    2330: {

    2335:   if (s->bc) {
    2336:     VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
    2337:   } else *indices = NULL;
    2338:   return(0);
    2339: }

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

    2344:   Not collective

    2346:   Input Parameters:
    2347: + s     - The PetscSection
    2348: . point - The point
    2349: - indices - The constrained dofs

    2351:   Note: The Fortran is PetscSectionSetConstraintIndicesF90()

    2353:   Level: intermediate

    2355: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
    2356: @*/
    2357: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
    2358: {

    2363:   if (s->bc) {
    2364:     VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
    2365:   }
    2366:   return(0);
    2367: }

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

    2372:   Not collective

    2374:   Input Parameters:
    2375: + s     - The PetscSection
    2376: . field  - The field number
    2377: - point - The point

    2379:   Output Parameter:
    2380: . indices - The constrained dofs sorted in ascending order

    2382:   Notes:
    2383:   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().

    2385:   Fortran Note:
    2386:   In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()

    2388:   Level: intermediate

    2390: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
    2391: @*/
    2392: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
    2393: {

    2398:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    2399:   PetscSectionGetConstraintIndices(s->field[field], point, indices);
    2400:   return(0);
    2401: }

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

    2406:   Not collective

    2408:   Input Parameters:
    2409: + s       - The PetscSection
    2410: . point   - The point
    2411: . field   - The field number
    2412: - indices - The constrained dofs

    2414:   Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()

    2416:   Level: intermediate

    2418: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
    2419: @*/
    2420: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
    2421: {

    2426:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    2427:   PetscSectionSetConstraintIndices(s->field[field], point, indices);
    2428:   return(0);
    2429: }

    2431: /*@
    2432:   PetscSectionPermute - Reorder the section according to the input point permutation

    2434:   Collective on PetscSection

    2436:   Input Parameter:
    2437: + section - The PetscSection object
    2438: - perm - The point permutation, old point p becomes new point perm[p]

    2440:   Output Parameter:
    2441: . sectionNew - The permuted PetscSection

    2443:   Level: intermediate

    2445: .seealso: MatPermute()
    2446: @*/
    2447: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
    2448: {
    2449:   PetscSection    s = section, sNew;
    2450:   const PetscInt *perm;
    2451:   PetscInt        numFields, f, c, numPoints, pStart, pEnd, p;
    2452:   PetscErrorCode  ierr;

    2458:   PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
    2459:   PetscSectionGetNumFields(s, &numFields);
    2460:   if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
    2461:   for (f = 0; f < numFields; ++f) {
    2462:     const char *name;
    2463:     PetscInt    numComp;

    2465:     PetscSectionGetFieldName(s, f, &name);
    2466:     PetscSectionSetFieldName(sNew, f, name);
    2467:     PetscSectionGetFieldComponents(s, f, &numComp);
    2468:     PetscSectionSetFieldComponents(sNew, f, numComp);
    2469:     for (c = 0; c < s->numFieldComponents[f]; ++c) {
    2470:       PetscSectionGetComponentName(s, f, c, &name);
    2471:       PetscSectionSetComponentName(sNew, f, c, name);
    2472:     }
    2473:   }
    2474:   ISGetLocalSize(permutation, &numPoints);
    2475:   ISGetIndices(permutation, &perm);
    2476:   PetscSectionGetChart(s, &pStart, &pEnd);
    2477:   PetscSectionSetChart(sNew, pStart, pEnd);
    2478:   if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
    2479:   for (p = pStart; p < pEnd; ++p) {
    2480:     PetscInt dof, cdof;

    2482:     PetscSectionGetDof(s, p, &dof);
    2483:     PetscSectionSetDof(sNew, perm[p], dof);
    2484:     PetscSectionGetConstraintDof(s, p, &cdof);
    2485:     if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
    2486:     for (f = 0; f < numFields; ++f) {
    2487:       PetscSectionGetFieldDof(s, p, f, &dof);
    2488:       PetscSectionSetFieldDof(sNew, perm[p], f, dof);
    2489:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
    2490:       if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
    2491:     }
    2492:   }
    2493:   PetscSectionSetUp(sNew);
    2494:   for (p = pStart; p < pEnd; ++p) {
    2495:     const PetscInt *cind;
    2496:     PetscInt        cdof;

    2498:     PetscSectionGetConstraintDof(s, p, &cdof);
    2499:     if (cdof) {
    2500:       PetscSectionGetConstraintIndices(s, p, &cind);
    2501:       PetscSectionSetConstraintIndices(sNew, perm[p], cind);
    2502:     }
    2503:     for (f = 0; f < numFields; ++f) {
    2504:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
    2505:       if (cdof) {
    2506:         PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
    2507:         PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
    2508:       }
    2509:     }
    2510:   }
    2511:   ISRestoreIndices(permutation, &perm);
    2512:   *sectionNew = sNew;
    2513:   return(0);
    2514: }

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

    2519:   Collective on section

    2521:   Input Parameters:
    2522: + section   - The PetscSection
    2523: . obj       - A PetscObject which serves as the key for this index
    2524: . clSection - Section giving the size of the closure of each point
    2525: - clPoints  - IS giving the points in each closure

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

    2529:   Level: advanced

    2531: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
    2532: @*/
    2533: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
    2534: {

    2541:   if (section->clObj != obj) {PetscSectionResetClosurePermutation(section);}
    2542:   section->clObj     = obj;
    2543:   PetscObjectReference((PetscObject)clSection);
    2544:   PetscObjectReference((PetscObject)clPoints);
    2545:   PetscSectionDestroy(&section->clSection);
    2546:   ISDestroy(&section->clPoints);
    2547:   section->clSection = clSection;
    2548:   section->clPoints  = clPoints;
    2549:   return(0);
    2550: }

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

    2555:   Collective on section

    2557:   Input Parameters:
    2558: + section   - The PetscSection
    2559: - obj       - A PetscObject which serves as the key for this index

    2561:   Output Parameters:
    2562: + clSection - Section giving the size of the closure of each point
    2563: - clPoints  - IS giving the points in each closure

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

    2567:   Level: advanced

    2569: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
    2570: @*/
    2571: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
    2572: {
    2574:   if (section->clObj == obj) {
    2575:     if (clSection) *clSection = section->clSection;
    2576:     if (clPoints)  *clPoints  = section->clPoints;
    2577:   } else {
    2578:     if (clSection) *clSection = NULL;
    2579:     if (clPoints)  *clPoints  = NULL;
    2580:   }
    2581:   return(0);
    2582: }

    2584: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
    2585: {
    2586:   PetscInt       i;
    2587:   khiter_t iter;
    2588:   int new_entry;
    2589:   PetscSectionClosurePermKey key = {depth, clSize};
    2590:   PetscSectionClosurePermVal *val;

    2594:   if (section->clObj != obj) {
    2595:     PetscSectionDestroy(&section->clSection);
    2596:     ISDestroy(&section->clPoints);
    2597:   }
    2598:   section->clObj = obj;
    2599:   if (!section->clHash) {PetscClPermCreate(&section->clHash);}
    2600:   iter = kh_put(ClPerm, section->clHash, key, &new_entry);
    2601:   val = &kh_val(section->clHash, iter);
    2602:   if (!new_entry) {
    2603:     PetscFree(val->perm);
    2604:     PetscFree(val->invPerm);
    2605:   }
    2606:   if (mode == PETSC_COPY_VALUES) {
    2607:     PetscMalloc1(clSize, &val->perm);
    2608:     PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
    2609:     PetscArraycpy(val->perm, clPerm, clSize);
    2610:   } else if (mode == PETSC_OWN_POINTER) {
    2611:     val->perm = clPerm;
    2612:   } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
    2613:   PetscMalloc1(clSize, &val->invPerm);
    2614:   for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
    2615:   return(0);
    2616: }

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

    2621:   Not Collective

    2623:   Input Parameters:
    2624: + section - The PetscSection
    2625: . obj     - A PetscObject which serves as the key for this index (usually a DM)
    2626: . depth   - Depth of points on which to apply the given permutation
    2627: - perm    - Permutation of the cell dof closure

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

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

    2637:   Level: intermediate

    2639: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
    2640: @*/
    2641: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
    2642: {
    2643:   const PetscInt *clPerm = NULL;
    2644:   PetscInt        clSize = 0;
    2645:   PetscErrorCode  ierr;

    2648:   if (perm) {
    2649:     ISGetLocalSize(perm, &clSize);
    2650:     ISGetIndices(perm, &clPerm);
    2651:   }
    2652:   PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
    2653:   if (perm) {ISRestoreIndices(perm, &clPerm);}
    2654:   return(0);
    2655: }

    2657: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
    2658: {

    2662:   if (section->clObj == obj) {
    2663:     PetscSectionClosurePermKey k = {depth, size};
    2664:     PetscSectionClosurePermVal v;
    2665:     PetscClPermGet(section->clHash, k, &v);
    2666:     if (perm) *perm = v.perm;
    2667:   } else {
    2668:     if (perm) *perm = NULL;
    2669:   }
    2670:   return(0);
    2671: }

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

    2676:   Not collective

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

    2684:   Output Parameter:
    2685: . perm - The dof closure permutation

    2687:   Note:
    2688:   The user must destroy the IS that is returned.

    2690:   Level: intermediate

    2692: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
    2693: @*/
    2694: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
    2695: {
    2696:   const PetscInt *clPerm;
    2697:   PetscErrorCode  ierr;

    2700:   PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
    2701:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
    2702:   return(0);
    2703: }

    2705: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
    2706: {

    2710:   if (section->clObj == obj && section->clHash) {
    2711:     PetscSectionClosurePermKey k = {depth, size};
    2712:     PetscSectionClosurePermVal v;
    2713:     PetscClPermGet(section->clHash, k, &v);
    2714:     if (perm) *perm = v.invPerm;
    2715:   } else {
    2716:     if (perm) *perm = NULL;
    2717:   }
    2718:   return(0);
    2719: }

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

    2724:   Not collective

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

    2732:   Output Parameters:
    2733: . perm - The dof closure permutation

    2735:   Note:
    2736:   The user must destroy the IS that is returned.

    2738:   Level: intermediate

    2740: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
    2741: @*/
    2742: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
    2743: {
    2744:   const PetscInt *clPerm;
    2745:   PetscErrorCode  ierr;

    2748:   PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
    2749:   ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
    2750:   return(0);
    2751: }

    2753: /*@
    2754:   PetscSectionGetField - Get the subsection associated with a single field

    2756:   Input Parameters:
    2757: + s     - The PetscSection
    2758: - field - The field number

    2760:   Output Parameter:
    2761: . subs  - The subsection for the given field

    2763:   Level: intermediate

    2765: .seealso: PetscSectionSetNumFields()
    2766: @*/
    2767: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
    2768: {
    2772:   if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %D should be in [%D, %D)", field, 0, s->numFields);
    2773:   *subs = s->field[field];
    2774:   return(0);
    2775: }

    2777: PetscClassId      PETSC_SECTION_SYM_CLASSID;
    2778: PetscFunctionList PetscSectionSymList = NULL;

    2780: /*@
    2781:   PetscSectionSymCreate - Creates an empty PetscSectionSym object.

    2783:   Collective

    2785:   Input Parameter:
    2786: . comm - the MPI communicator

    2788:   Output Parameter:
    2789: . sym - pointer to the new set of symmetries

    2791:   Level: developer

    2793: .seealso: PetscSectionSym, PetscSectionSymDestroy()
    2794: @*/
    2795: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
    2796: {

    2801:   ISInitializePackage();
    2802:   PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
    2803:   return(0);
    2804: }

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

    2809:   Collective on PetscSectionSym

    2811:   Input Parameters:
    2812: + sym    - The section symmetry object
    2813: - method - The name of the section symmetry type

    2815:   Level: developer

    2817: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
    2818: @*/
    2819: PetscErrorCode  PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
    2820: {
    2821:   PetscErrorCode (*r)(PetscSectionSym);
    2822:   PetscBool      match;

    2827:   PetscObjectTypeCompare((PetscObject) sym, method, &match);
    2828:   if (match) return(0);

    2830:   PetscFunctionListFind(PetscSectionSymList,method,&r);
    2831:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
    2832:   if (sym->ops->destroy) {
    2833:     (*sym->ops->destroy)(sym);
    2834:     sym->ops->destroy = NULL;
    2835:   }
    2836:   (*r)(sym);
    2837:   PetscObjectChangeTypeName((PetscObject)sym,method);
    2838:   return(0);
    2839: }

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

    2845:   Not Collective

    2847:   Input Parameter:
    2848: . sym  - The section symmetry

    2850:   Output Parameter:
    2851: . type - The index set type name

    2853:   Level: developer

    2855: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
    2856: @*/
    2857: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
    2858: {
    2862:   *type = ((PetscObject)sym)->type_name;
    2863:   return(0);
    2864: }

    2866: /*@C
    2867:   PetscSectionSymRegister - Adds a new section symmetry implementation

    2869:   Not Collective

    2871:   Input Parameters:
    2872: + name        - The name of a new user-defined creation routine
    2873: - create_func - The creation routine itself

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

    2878:   Level: developer

    2880: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
    2881: @*/
    2882: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
    2883: {

    2887:   ISInitializePackage();
    2888:   PetscFunctionListAdd(&PetscSectionSymList,sname,function);
    2889:   return(0);
    2890: }

    2892: /*@
    2893:    PetscSectionSymDestroy - Destroys a section symmetry.

    2895:    Collective on PetscSectionSym

    2897:    Input Parameters:
    2898: .  sym - the section symmetry

    2900:    Level: developer

    2902: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
    2903: @*/
    2904: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
    2905: {
    2906:   SymWorkLink    link,next;

    2910:   if (!*sym) return(0);
    2912:   if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return(0);}
    2913:   if ((*sym)->ops->destroy) {
    2914:     (*(*sym)->ops->destroy)(*sym);
    2915:   }
    2916:   if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
    2917:   for (link=(*sym)->workin; link; link=next) {
    2918:     next = link->next;
    2919:     PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
    2920:     PetscFree(link);
    2921:   }
    2922:   (*sym)->workin = NULL;
    2923:   PetscHeaderDestroy(sym);
    2924:   return(0);
    2925: }

    2927: /*@C
    2928:    PetscSectionSymView - Displays a section symmetry

    2930:    Collective on PetscSectionSym

    2932:    Input Parameters:
    2933: +  sym - the index set
    2934: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

    2936:    Level: developer

    2938: .seealso: PetscViewerASCIIOpen()
    2939: @*/
    2940: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
    2941: {

    2946:   if (!viewer) {
    2947:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
    2948:   }
    2951:   PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
    2952:   if (sym->ops->view) {
    2953:     (*sym->ops->view)(sym,viewer);
    2954:   }
    2955:   return(0);
    2956: }

    2958: /*@
    2959:   PetscSectionSetSym - Set the symmetries for the data referred to by the section

    2961:   Collective on PetscSection

    2963:   Input Parameters:
    2964: + section - the section describing data layout
    2965: - sym - the symmetry describing the affect of orientation on the access of the data

    2967:   Level: developer

    2969: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
    2970: @*/
    2971: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
    2972: {

    2977:   PetscSectionSymDestroy(&(section->sym));
    2978:   if (sym) {
    2981:     PetscObjectReference((PetscObject) sym);
    2982:   }
    2983:   section->sym = sym;
    2984:   return(0);
    2985: }

    2987: /*@
    2988:   PetscSectionGetSym - Get the symmetries for the data referred to by the section

    2990:   Not collective

    2992:   Input Parameters:
    2993: . section - the section describing data layout

    2995:   Output Parameters:
    2996: . sym - the symmetry describing the affect of orientation on the access of the data

    2998:   Level: developer

    3000: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
    3001: @*/
    3002: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
    3003: {
    3006:   *sym = section->sym;
    3007:   return(0);
    3008: }

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

    3013:   Collective on PetscSection

    3015:   Input Parameters:
    3016: + section - the section describing data layout
    3017: . field - the field number
    3018: - sym - the symmetry describing the affect of orientation on the access of the data

    3020:   Level: developer

    3022: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
    3023: @*/
    3024: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
    3025: {

    3030:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
    3031:   PetscSectionSetSym(section->field[field],sym);
    3032:   return(0);
    3033: }

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

    3038:   Collective on PetscSection

    3040:   Input Parameters:
    3041: + section - the section describing data layout
    3042: - field - the field number

    3044:   Output Parameters:
    3045: . sym - the symmetry describing the affect of orientation on the access of the data

    3047:   Level: developer

    3049: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
    3050: @*/
    3051: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
    3052: {
    3055:   if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
    3056:   *sym = section->field[field]->sym;
    3057:   return(0);
    3058: }

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

    3063:   Not collective

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

    3072:   Output Parameter:
    3073: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
    3074: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
    3075:     identity).

    3077:   Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
    3078: .vb
    3079:      const PetscInt    **perms;
    3080:      const PetscScalar **rots;
    3081:      PetscInt            lOffset;

    3083:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
    3084:      for (i = 0, lOffset = 0; i < numPoints; i++) {
    3085:        PetscInt           point = points[2*i], dof, sOffset;
    3086:        const PetscInt    *perm  = perms ? perms[i] : NULL;
    3087:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

    3089:        PetscSectionGetDof(section,point,&dof);
    3090:        PetscSectionGetOffset(section,point,&sOffset);

    3092:        if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]]  = sArray[sOffset + j];}}
    3093:        else      {for (j = 0; j < dof; j++) {lArray[lOffset +      j ]  = sArray[sOffset + j];}}
    3094:        if (rot)  {for (j = 0; j < dof; j++) {lArray[lOffset +      j ] *= rot[j];             }}
    3095:        lOffset += dof;
    3096:      }
    3097:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
    3098: .ve

    3100:   Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
    3101: .vb
    3102:      const PetscInt    **perms;
    3103:      const PetscScalar **rots;
    3104:      PetscInt            lOffset;

    3106:      PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
    3107:      for (i = 0, lOffset = 0; i < numPoints; i++) {
    3108:        PetscInt           point = points[2*i], dof, sOffset;
    3109:        const PetscInt    *perm  = perms ? perms[i] : NULL;
    3110:        const PetscScalar *rot   = rots  ? rots[i]  : NULL;

    3112:        PetscSectionGetDof(section,point,&dof);
    3113:        PetscSectionGetOffset(section,point,&sOff);

    3115:        if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
    3116:        else      {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset +      j ] * (rot ? PetscConj(rot[     j ]) : 1.);}}
    3117:        offset += dof;
    3118:      }
    3119:      PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
    3120: .ve

    3122:   Level: developer

    3124: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
    3125: @*/
    3126: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
    3127: {
    3128:   PetscSectionSym sym;
    3129:   PetscErrorCode  ierr;

    3134:   if (perms) *perms = NULL;
    3135:   if (rots)  *rots  = NULL;
    3136:   sym = section->sym;
    3137:   if (sym && (perms || rots)) {
    3138:     SymWorkLink link;

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

    3162: /*@C
    3163:   PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()

    3165:   Not collective

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

    3174:   Output Parameter:
    3175: + perms - The permutations for the given orientations: set to NULL at conclusion
    3176: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion

    3178:   Level: developer

    3180: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
    3181: @*/
    3182: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
    3183: {
    3184:   PetscSectionSym sym;

    3188:   sym = section->sym;
    3189:   if (sym && (perms || rots)) {
    3190:     SymWorkLink *p,link;

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

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

    3210:   Not collective

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

    3220:   Output Parameter:
    3221: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
    3222: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
    3223:     identity).

    3225:   Level: developer

    3227: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
    3228: @*/
    3229: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
    3230: {

    3235:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
    3236:   PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
    3237:   return(0);
    3238: }

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

    3243:   Not collective

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

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

    3257:   Level: developer

    3259: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
    3260: @*/
    3261: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
    3262: {

    3267:   if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
    3268:   PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
    3269:   return(0);
    3270: }

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

    3275:   Not collective

    3277:   Input Parameter:
    3278: . s - the global PetscSection

    3280:   Output Parameters:
    3281: . flg - the flag

    3283:   Level: developer

    3285: .seealso: PetscSectionSetChart(), PetscSectionCreate()
    3286: @*/
    3287: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
    3288: {
    3291:   *flg = s->useFieldOff;
    3292:   return(0);
    3293: }

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

    3298:   Not collective

    3300:   Input Parameters:
    3301: + s   - the global PetscSection
    3302: - flg - the flag

    3304:   Level: developer

    3306: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
    3307: @*/
    3308: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
    3309: {
    3312:   s->useFieldOff = flg;
    3313:   return(0);
    3314: }

    3316: #define PetscSectionExpandPoints_Loop(TYPE) \
    3317: { \
    3318:   PetscInt i, n, o0, o1, size; \
    3319:   TYPE *a0 = (TYPE*)origArray, *a1; \
    3320:   PetscSectionGetStorageSize(s, &size); \
    3321:   PetscMalloc1(size, &a1); \
    3322:   for (i=0; i<npoints; i++) { \
    3323:     PetscSectionGetOffset(origSection, points_[i], &o0); \
    3324:     PetscSectionGetOffset(s, i, &o1); \
    3325:     PetscSectionGetDof(s, i, &n); \
    3326:     PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
    3327:   } \
    3328:   *newArray = (void*)a1; \
    3329: }

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

    3334:   Not collective

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

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

    3346:   Level: developer

    3348: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
    3349: @*/
    3350: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
    3351: {
    3352:   PetscSection        s;
    3353:   const PetscInt      *points_;
    3354:   PetscInt            i, n, npoints, pStart, pEnd;
    3355:   PetscMPIInt         unitsize;
    3356:   PetscErrorCode      ierr;

    3364:   MPI_Type_size(dataType, &unitsize);
    3365:   ISGetLocalSize(points, &npoints);
    3366:   ISGetIndices(points, &points_);
    3367:   PetscSectionGetChart(origSection, &pStart, &pEnd);
    3368:   PetscSectionCreate(PETSC_COMM_SELF, &s);
    3369:   PetscSectionSetChart(s, 0, npoints);
    3370:   for (i=0; i<npoints; i++) {
    3371:     if (PetscUnlikely(points_[i] < pStart || points_[i] >= pEnd)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %d (index %d) in input IS out of input section's chart", points_[i], i);
    3372:     PetscSectionGetDof(origSection, points_[i], &n);
    3373:     PetscSectionSetDof(s, i, n);
    3374:   }
    3375:   PetscSectionSetUp(s);
    3376:   if (newArray) {
    3377:     if (dataType == MPIU_INT)           {PetscSectionExpandPoints_Loop(PetscInt);}
    3378:     else if (dataType == MPIU_SCALAR)   {PetscSectionExpandPoints_Loop(PetscScalar);}
    3379:     else if (dataType == MPIU_REAL)     {PetscSectionExpandPoints_Loop(PetscReal);}
    3380:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
    3381:   }
    3382:   if (newSection) {
    3383:     *newSection = s;
    3384:   } else {
    3385:     PetscSectionDestroy(&s);
    3386:   }
    3387:   return(0);
    3388: }