Actual source code: plexsection.c

      1: #include <petsc/private/dmpleximpl.h>

      3: /* Set the number of dof on each point and separate by fields */
      4: static PetscErrorCode DMPlexCreateSectionFields(DM dm, const PetscInt numComp[], PetscSection *section)
      5: {
      6:   DMLabel    depthLabel;
      7:   PetscInt   depth, Nf, f, pStart, pEnd;
      8:   PetscBool *isFE;

     10:   PetscFunctionBegin;
     11:   PetscCall(DMPlexGetDepth(dm, &depth));
     12:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
     13:   PetscCall(DMGetNumFields(dm, &Nf));
     14:   PetscCall(PetscCalloc1(Nf, &isFE));
     15:   for (f = 0; f < Nf; ++f) {
     16:     PetscObject  obj;
     17:     PetscClassId id;

     19:     PetscCall(DMGetField(dm, f, NULL, &obj));
     20:     PetscCall(PetscObjectGetClassId(obj, &id));
     21:     if (id == PETSCFE_CLASSID) {
     22:       isFE[f] = PETSC_TRUE;
     23:     } else if (id == PETSCFV_CLASSID) {
     24:       isFE[f] = PETSC_FALSE;
     25:     }
     26:   }

     28:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), section));
     29:   if (Nf > 0) {
     30:     PetscCall(PetscSectionSetNumFields(*section, Nf));
     31:     if (numComp) {
     32:       for (f = 0; f < Nf; ++f) {
     33:         PetscCall(PetscSectionSetFieldComponents(*section, f, numComp[f]));
     34:         if (isFE[f]) {
     35:           PetscFE              fe;
     36:           PetscDualSpace       dspace;
     37:           const PetscInt    ***perms;
     38:           const PetscScalar ***flips;
     39:           const PetscInt      *numDof;

     41:           PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
     42:           PetscCall(PetscFEGetDualSpace(fe, &dspace));
     43:           PetscCall(PetscDualSpaceGetSymmetries(dspace, &perms, &flips));
     44:           PetscCall(PetscDualSpaceGetNumDof(dspace, &numDof));
     45:           if (perms || flips) {
     46:             DM              K;
     47:             PetscInt        sph, spdepth;
     48:             PetscSectionSym sym;

     50:             PetscCall(PetscDualSpaceGetDM(dspace, &K));
     51:             PetscCall(DMPlexGetDepth(K, &spdepth));
     52:             PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)*section), depthLabel, &sym));
     53:             for (sph = 0; sph <= spdepth; sph++) {
     54:               PetscDualSpace      hspace;
     55:               PetscInt            kStart, kEnd;
     56:               PetscInt            kConeSize, h = sph + (depth - spdepth);
     57:               const PetscInt    **perms0 = NULL;
     58:               const PetscScalar **flips0 = NULL;

     60:               PetscCall(PetscDualSpaceGetHeightSubspace(dspace, sph, &hspace));
     61:               PetscCall(DMPlexGetHeightStratum(K, h, &kStart, &kEnd));
     62:               if (!hspace) continue;
     63:               PetscCall(PetscDualSpaceGetSymmetries(hspace, &perms, &flips));
     64:               if (perms) perms0 = perms[0];
     65:               if (flips) flips0 = flips[0];
     66:               if (!(perms0 || flips0)) continue;
     67:               {
     68:                 DMPolytopeType ct;
     69:                 /* The number of arrangements is no longer based on the number of faces */
     70:                 PetscCall(DMPlexGetCellType(K, kStart, &ct));
     71:                 kConeSize = DMPolytopeTypeGetNumArrangements(ct) / 2;
     72:               }
     73:               PetscCall(PetscSectionSymLabelSetStratum(sym, depth - h, numDof[depth - h], -kConeSize, kConeSize, PETSC_USE_POINTER, PetscSafePointerPlusOffset(perms0, -kConeSize), PetscSafePointerPlusOffset(flips0, -kConeSize)));
     74:             }
     75:             PetscCall(PetscSectionSetFieldSym(*section, f, sym));
     76:             PetscCall(PetscSectionSymDestroy(&sym));
     77:           }
     78:         }
     79:       }
     80:     }
     81:   }
     82:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
     83:   PetscCall(PetscSectionSetChart(*section, pStart, pEnd));
     84:   PetscCall(PetscFree(isFE));
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: /* Set the number of dof on each point and separate by fields */
     89: static PetscErrorCode DMPlexCreateSectionDof(DM dm, DMLabel label[], const PetscInt numDof[], PetscSection section)
     90: {
     91:   DMLabel        depthLabel;
     92:   DMPolytopeType ct;
     93:   PetscInt       depth, cellHeight, pStart = 0, pEnd = 0;
     94:   PetscInt       Nf, f, Nds, n, dim, d, dep, p;
     95:   PetscBool     *isFE, hasCohesive = PETSC_FALSE;

     97:   PetscFunctionBegin;
     98:   PetscCall(DMGetDimension(dm, &dim));
     99:   PetscCall(DMPlexGetDepth(dm, &depth));
    100:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
    101:   PetscCall(DMGetNumFields(dm, &Nf));
    102:   PetscCall(DMGetNumDS(dm, &Nds));
    103:   for (n = 0; n < Nds; ++n) {
    104:     PetscDS   ds;
    105:     PetscBool isCohesive;

    107:     PetscCall(DMGetRegionNumDS(dm, n, NULL, NULL, &ds, NULL));
    108:     PetscCall(PetscDSIsCohesive(ds, &isCohesive));
    109:     if (isCohesive) {
    110:       hasCohesive = PETSC_TRUE;
    111:       break;
    112:     }
    113:   }
    114:   PetscCall(PetscMalloc1(Nf, &isFE));
    115:   for (f = 0; f < Nf; ++f) {
    116:     PetscObject  obj;
    117:     PetscClassId id;

    119:     PetscCall(DMGetField(dm, f, NULL, &obj));
    120:     PetscCall(PetscObjectGetClassId(obj, &id));
    121:     /* User is allowed to put a "placeholder" field in (c.f. DMCreateDS) */
    122:     isFE[f] = id == PETSCFE_CLASSID ? PETSC_TRUE : PETSC_FALSE;
    123:   }

    125:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
    126:   for (f = 0; f < Nf; ++f) {
    127:     PetscBool avoidTensor;

    129:     PetscCall(DMGetFieldAvoidTensor(dm, f, &avoidTensor));
    130:     avoidTensor = (avoidTensor || hasCohesive) ? PETSC_TRUE : PETSC_FALSE;
    131:     if (label && label[f]) {
    132:       IS              pointIS;
    133:       const PetscInt *points;
    134:       PetscInt        n;

    136:       PetscCall(DMLabelGetStratumIS(label[f], 1, &pointIS));
    137:       if (!pointIS) continue;
    138:       PetscCall(ISGetLocalSize(pointIS, &n));
    139:       PetscCall(ISGetIndices(pointIS, &points));
    140:       for (p = 0; p < n; ++p) {
    141:         const PetscInt point = points[p];
    142:         PetscInt       dof, d;

    144:         PetscCall(DMPlexGetCellType(dm, point, &ct));
    145:         PetscCall(DMLabelGetValue(depthLabel, point, &d));
    146:         /* If this is a tensor prism point, use dof for one dimension lower */
    147:         switch (ct) {
    148:         case DM_POLYTOPE_POINT_PRISM_TENSOR:
    149:         case DM_POLYTOPE_SEG_PRISM_TENSOR:
    150:         case DM_POLYTOPE_TRI_PRISM_TENSOR:
    151:         case DM_POLYTOPE_QUAD_PRISM_TENSOR:
    152:           if (hasCohesive) --d;
    153:           break;
    154:         default:
    155:           break;
    156:         }
    157:         dof = d < 0 ? 0 : numDof[f * (dim + 1) + d];
    158:         PetscCall(PetscSectionSetFieldDof(section, point, f, dof));
    159:         PetscCall(PetscSectionAddDof(section, point, dof));
    160:       }
    161:       PetscCall(ISRestoreIndices(pointIS, &points));
    162:       PetscCall(ISDestroy(&pointIS));
    163:     } else {
    164:       for (dep = 0; dep <= depth - cellHeight; ++dep) {
    165:         /* Cases: dim > depth (cell-vertex mesh), dim == depth (fully interpolated), dim < depth (interpolated submesh) */
    166:         d = dim <= depth ? dep : (!dep ? 0 : dim);
    167:         PetscCall(DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd));
    168:         for (p = pStart; p < pEnd; ++p) {
    169:           const PetscInt dof = numDof[f * (dim + 1) + d];

    171:           PetscCall(DMPlexGetCellType(dm, p, &ct));
    172:           switch (ct) {
    173:           case DM_POLYTOPE_POINT_PRISM_TENSOR:
    174:           case DM_POLYTOPE_SEG_PRISM_TENSOR:
    175:           case DM_POLYTOPE_TRI_PRISM_TENSOR:
    176:           case DM_POLYTOPE_QUAD_PRISM_TENSOR:
    177:             if (avoidTensor && isFE[f]) continue;
    178:           default:
    179:             break;
    180:           }
    181:           PetscCall(PetscSectionSetFieldDof(section, p, f, dof));
    182:           PetscCall(PetscSectionAddDof(section, p, dof));
    183:         }
    184:       }
    185:     }
    186:   }
    187:   PetscCall(PetscFree(isFE));
    188:   PetscFunctionReturn(PETSC_SUCCESS);
    189: }

    191: /* Set the number of dof on each point and separate by fields
    192:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
    193: */
    194: static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
    195: {
    196:   PetscInt     Nf;
    197:   PetscInt     bc;
    198:   PetscSection aSec;

    200:   PetscFunctionBegin;
    201:   PetscCall(PetscSectionGetNumFields(section, &Nf));
    202:   for (bc = 0; bc < numBC; ++bc) {
    203:     PetscInt        field = 0;
    204:     const PetscInt *comp;
    205:     const PetscInt *idx;
    206:     PetscInt        Nc = 0, cNc = -1, n, i;

    208:     if (Nf) {
    209:       field = bcField[bc];
    210:       PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
    211:     }
    212:     if (bcComps && bcComps[bc]) PetscCall(ISGetLocalSize(bcComps[bc], &cNc));
    213:     if (bcComps && bcComps[bc]) PetscCall(ISGetIndices(bcComps[bc], &comp));
    214:     if (bcPoints[bc]) {
    215:       PetscCall(ISGetLocalSize(bcPoints[bc], &n));
    216:       PetscCall(ISGetIndices(bcPoints[bc], &idx));
    217:       for (i = 0; i < n; ++i) {
    218:         const PetscInt p = idx[i];
    219:         PetscInt       numConst;

    221:         if (Nf) {
    222:           PetscCall(PetscSectionGetFieldDof(section, p, field, &numConst));
    223:         } else {
    224:           PetscCall(PetscSectionGetDof(section, p, &numConst));
    225:         }
    226:         /* If Nc <= 0, constrain every dof on the point */
    227:         if (cNc > 0) {
    228:           /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
    229:              and that those dofs are numbered n*Nc+c */
    230:           if (Nf) {
    231:             PetscCheck(numConst % Nc == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " has %" PetscInt_FMT " dof which is not divisible by %" PetscInt_FMT " field components", p, numConst, Nc);
    232:             numConst = (numConst / Nc) * cNc;
    233:           } else {
    234:             numConst = PetscMin(numConst, cNc);
    235:           }
    236:         }
    237:         if (Nf) PetscCall(PetscSectionAddFieldConstraintDof(section, p, field, numConst));
    238:         PetscCall(PetscSectionAddConstraintDof(section, p, numConst));
    239:       }
    240:       PetscCall(ISRestoreIndices(bcPoints[bc], &idx));
    241:     }
    242:     if (bcComps && bcComps[bc]) PetscCall(ISRestoreIndices(bcComps[bc], &comp));
    243:   }
    244:   PetscCall(DMPlexGetAnchors(dm, &aSec, NULL));
    245:   if (aSec) {
    246:     PetscInt aStart, aEnd, a;

    248:     PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
    249:     for (a = aStart; a < aEnd; a++) {
    250:       PetscInt dof, f;

    252:       PetscCall(PetscSectionGetDof(aSec, a, &dof));
    253:       if (dof) {
    254:         /* if there are point-to-point constraints, then all dofs are constrained */
    255:         PetscCall(PetscSectionGetDof(section, a, &dof));
    256:         PetscCall(PetscSectionSetConstraintDof(section, a, dof));
    257:         for (f = 0; f < Nf; f++) {
    258:           PetscCall(PetscSectionGetFieldDof(section, a, f, &dof));
    259:           PetscCall(PetscSectionSetFieldConstraintDof(section, a, f, dof));
    260:         }
    261:       }
    262:     }
    263:   }
    264:   PetscFunctionReturn(PETSC_SUCCESS);
    265: }

    267: /* Set the constrained field indices on each point
    268:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
    269: */
    270: static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
    271: {
    272:   PetscSection aSec;
    273:   PetscInt    *indices;
    274:   PetscInt     Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;

    276:   PetscFunctionBegin;
    277:   PetscCall(PetscSectionGetNumFields(section, &Nf));
    278:   if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
    279:   /* Initialize all field indices to -1 */
    280:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
    281:   for (p = pStart; p < pEnd; ++p) {
    282:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
    283:     maxDof = PetscMax(maxDof, cdof);
    284:   }
    285:   PetscCall(PetscMalloc1(maxDof, &indices));
    286:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
    287:   for (p = pStart; p < pEnd; ++p)
    288:     for (f = 0; f < Nf; ++f) PetscCall(PetscSectionSetFieldConstraintIndices(section, p, f, indices));
    289:   /* Handle BC constraints */
    290:   for (bc = 0; bc < numBC; ++bc) {
    291:     const PetscInt  field = bcField[bc];
    292:     const PetscInt *comp, *idx;
    293:     PetscInt        Nc, cNc = -1, n, i;

    295:     PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
    296:     if (bcComps && bcComps[bc]) PetscCall(ISGetLocalSize(bcComps[bc], &cNc));
    297:     if (bcComps && bcComps[bc]) PetscCall(ISGetIndices(bcComps[bc], &comp));
    298:     if (bcPoints[bc]) {
    299:       PetscCall(ISGetLocalSize(bcPoints[bc], &n));
    300:       PetscCall(ISGetIndices(bcPoints[bc], &idx));
    301:       for (i = 0; i < n; ++i) {
    302:         const PetscInt  p = idx[i];
    303:         const PetscInt *find;
    304:         PetscInt        fdof, fcdof, c, j;

    306:         PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
    307:         if (!fdof) continue;
    308:         if (cNc < 0) {
    309:           for (d = 0; d < fdof; ++d) indices[d] = d;
    310:           fcdof = fdof;
    311:         } else {
    312:           /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
    313:              and that those dofs are numbered n*Nc+c */
    314:           PetscCall(PetscSectionGetFieldConstraintDof(section, p, field, &fcdof));
    315:           PetscCall(PetscSectionGetFieldConstraintIndices(section, p, field, &find));
    316:           /* Get indices constrained by previous bcs */
    317:           for (d = 0; d < fcdof; ++d) {
    318:             if (find[d] < 0) break;
    319:             indices[d] = find[d];
    320:           }
    321:           for (j = 0; j < fdof / Nc; ++j)
    322:             for (c = 0; c < cNc; ++c) indices[d++] = j * Nc + comp[c];
    323:           PetscCall(PetscSortRemoveDupsInt(&d, indices));
    324:           for (c = d; c < fcdof; ++c) indices[c] = -1;
    325:           fcdof = d;
    326:         }
    327:         PetscCall(PetscSectionSetFieldConstraintDof(section, p, field, fcdof));
    328:         PetscCall(PetscSectionSetFieldConstraintIndices(section, p, field, indices));
    329:       }
    330:       PetscCall(ISRestoreIndices(bcPoints[bc], &idx));
    331:     }
    332:     if (bcComps && bcComps[bc]) PetscCall(ISRestoreIndices(bcComps[bc], &comp));
    333:   }
    334:   /* Handle anchors */
    335:   PetscCall(DMPlexGetAnchors(dm, &aSec, NULL));
    336:   if (aSec) {
    337:     PetscInt aStart, aEnd, a;

    339:     for (d = 0; d < maxDof; ++d) indices[d] = d;
    340:     PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
    341:     for (a = aStart; a < aEnd; a++) {
    342:       PetscInt dof, f;

    344:       PetscCall(PetscSectionGetDof(aSec, a, &dof));
    345:       if (dof) {
    346:         /* if there are point-to-point constraints, then all dofs are constrained */
    347:         for (f = 0; f < Nf; f++) PetscCall(PetscSectionSetFieldConstraintIndices(section, a, f, indices));
    348:       }
    349:     }
    350:   }
    351:   PetscCall(PetscFree(indices));
    352:   PetscFunctionReturn(PETSC_SUCCESS);
    353: }

    355: /* Set the constrained indices on each point */
    356: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
    357: {
    358:   PetscInt *indices;
    359:   PetscInt  Nf, maxDof, pStart, pEnd, p, f, d;

    361:   PetscFunctionBegin;
    362:   PetscCall(PetscSectionGetNumFields(section, &Nf));
    363:   PetscCall(PetscSectionGetMaxDof(section, &maxDof));
    364:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
    365:   PetscCall(PetscMalloc1(maxDof, &indices));
    366:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
    367:   for (p = pStart; p < pEnd; ++p) {
    368:     PetscInt cdof, d;

    370:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
    371:     if (cdof) {
    372:       if (Nf) {
    373:         PetscInt numConst = 0, foff = 0;

    375:         for (f = 0; f < Nf; ++f) {
    376:           const PetscInt *find;
    377:           PetscInt        fcdof, fdof;

    379:           PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
    380:           PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
    381:           /* Change constraint numbering from field component to local dof number */
    382:           PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &find));
    383:           for (d = 0; d < fcdof; ++d) indices[numConst + d] = find[d] + foff;
    384:           numConst += fcdof;
    385:           foff += fdof;
    386:         }
    387:         if (cdof != numConst) PetscCall(PetscSectionSetConstraintDof(section, p, numConst));
    388:       } else {
    389:         for (d = 0; d < cdof; ++d) indices[d] = d;
    390:       }
    391:       PetscCall(PetscSectionSetConstraintIndices(section, p, indices));
    392:     }
    393:   }
    394:   PetscCall(PetscFree(indices));
    395:   PetscFunctionReturn(PETSC_SUCCESS);
    396: }

    398: /*@C
    399:   DMPlexCreateSection - Create a `PetscSection` based upon the dof layout specification provided.

    401:   Not Collective

    403:   Input Parameters:
    404: + dm       - The `DMPLEX` object
    405: . label    - The label indicating the mesh support of each field, or NULL for the whole mesh
    406: . numComp  - An array of size numFields that holds the number of components for each field
    407: . numDof   - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
    408: . numBC    - The number of boundary conditions
    409: . bcField  - An array of size numBC giving the field number for each boundary condition
    410: . bcComps  - [Optional] An array of size numBC giving an `IS` holding the field components to which each boundary condition applies
    411: . bcPoints - An array of size numBC giving an `IS` holding the `DMPLEX` points to which each boundary condition applies
    412: - perm     - Optional permutation of the chart, or NULL

    414:   Output Parameter:
    415: . section - The `PetscSection` object

    417:   Level: developer

    419:   Notes:
    420:   numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the
    421:   number of dof for field 0 on each edge.

    423:   The chart permutation is the same one set using `PetscSectionSetPermutation()`

    425: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
    426: @*/
    427: PetscErrorCode DMPlexCreateSection(DM dm, DMLabel label[], const PetscInt numComp[], const PetscInt numDof[], PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], IS perm, PetscSection *section)
    428: {
    429:   PetscSection aSec;

    431:   PetscFunctionBegin;
    432:   PetscCall(DMPlexCreateSectionFields(dm, numComp, section));
    433:   PetscCall(DMPlexCreateSectionDof(dm, label, numDof, *section));
    434:   PetscCall(DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section));
    435:   if (perm) PetscCall(PetscSectionSetPermutation(*section, perm));
    436:   PetscCall(PetscSectionSetFromOptions(*section));
    437:   PetscCall(PetscSectionSetUp(*section));
    438:   PetscCall(DMPlexGetAnchors(dm, &aSec, NULL));
    439:   if (numBC || aSec) {
    440:     PetscCall(DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section));
    441:     PetscCall(DMPlexCreateSectionBCIndices(dm, *section));
    442:   }
    443:   PetscCall(PetscSectionViewFromOptions(*section, NULL, "-section_view"));
    444:   PetscFunctionReturn(PETSC_SUCCESS);
    445: }

    447: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
    448: {
    449:   PetscSection section;
    450:   DMLabel     *labels;
    451:   IS          *bcPoints, *bcComps, permIS;
    452:   PetscBT      blockStarts;
    453:   PetscBool   *isFE;
    454:   PetscInt    *bcFields, *numComp, *numDof;
    455:   PetscInt     depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f;
    456:   PetscInt     cStart, cEnd, cEndInterior;

    458:   PetscFunctionBegin;
    459:   PetscCall(DMGetNumFields(dm, &Nf));
    460:   PetscCall(DMGetDimension(dm, &dim));
    461:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    462:   /* FE and FV boundary conditions are handled slightly differently */
    463:   PetscCall(PetscMalloc1(Nf, &isFE));
    464:   for (f = 0; f < Nf; ++f) {
    465:     PetscObject  obj;
    466:     PetscClassId id;

    468:     PetscCall(DMGetField(dm, f, NULL, &obj));
    469:     PetscCall(PetscObjectGetClassId(obj, &id));
    470:     if (id == PETSCFE_CLASSID) {
    471:       isFE[f] = PETSC_TRUE;
    472:     } else if (id == PETSCFV_CLASSID) {
    473:       isFE[f] = PETSC_FALSE;
    474:     } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
    475:   }
    476:   /* Allocate boundary point storage for FEM boundaries */
    477:   PetscCall(DMGetNumDS(dm, &Nds));
    478:   for (s = 0; s < Nds; ++s) {
    479:     PetscDS  dsBC;
    480:     PetscInt numBd, bd;

    482:     PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC, NULL));
    483:     PetscCall(PetscDSGetNumBoundary(dsBC, &numBd));
    484:     PetscCheck(Nf || !numBd, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "number of fields is zero and number of boundary conditions is nonzero (this should never happen)");
    485:     for (bd = 0; bd < numBd; ++bd) {
    486:       PetscInt                field;
    487:       DMBoundaryConditionType type;
    488:       DMLabel                 label;

    490:       PetscCall(PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL));
    491:       if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
    492:     }
    493:   }
    494:   /* Add ghost cell boundaries for FVM */
    495:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
    496:   for (f = 0; f < Nf; ++f)
    497:     if (!isFE[f] && cEndInterior >= 0) ++numBC;
    498:   PetscCall(PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps));
    499:   /* Constrain ghost cells for FV */
    500:   for (f = 0; f < Nf; ++f) {
    501:     PetscInt *newidx, c;

    503:     if (isFE[f] || cEndInterior < 0) continue;
    504:     PetscCall(PetscMalloc1(cEnd - cEndInterior, &newidx));
    505:     for (c = cEndInterior; c < cEnd; ++c) newidx[c - cEndInterior] = c;
    506:     bcFields[bc] = f;
    507:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cEnd - cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]));
    508:   }
    509:   /* Complete labels for boundary conditions */
    510:   PetscCall(DMCompleteBCLabels_Internal(dm));
    511:   /* Handle FEM Dirichlet boundaries */
    512:   PetscCall(DMGetNumDS(dm, &Nds));
    513:   for (s = 0; s < Nds; ++s) {
    514:     PetscDS  dsBC;
    515:     PetscInt numBd, bd;

    517:     PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC, NULL));
    518:     PetscCall(PetscDSGetNumBoundary(dsBC, &numBd));
    519:     for (bd = 0; bd < numBd; ++bd) {
    520:       DMLabel                 label;
    521:       const PetscInt         *comps;
    522:       const PetscInt         *values;
    523:       PetscInt                bd2, field, numComps, numValues;
    524:       DMBoundaryConditionType type;
    525:       PetscBool               duplicate = PETSC_FALSE;

    527:       PetscCall(PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, &numValues, &values, &field, &numComps, &comps, NULL, NULL, NULL));
    528:       if (!isFE[field] || !label) continue;
    529:       /* Only want to modify label once */
    530:       for (bd2 = 0; bd2 < bd; ++bd2) {
    531:         DMLabel l;

    533:         PetscCall(PetscDSGetBoundary(dsBC, bd2, NULL, NULL, NULL, &l, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    534:         duplicate = l == label ? PETSC_TRUE : PETSC_FALSE;
    535:         if (duplicate) break;
    536:       }
    537:       /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
    538:       if (type & DM_BC_ESSENTIAL) {
    539:         PetscInt *newidx;
    540:         PetscInt  n, newn = 0, p, v;

    542:         bcFields[bc] = field;
    543:         if (numComps) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]));
    544:         for (v = 0; v < numValues; ++v) {
    545:           IS              tmp;
    546:           const PetscInt *idx;

    548:           PetscCall(DMLabelGetStratumIS(label, values[v], &tmp));
    549:           if (!tmp) continue;
    550:           PetscCall(ISGetLocalSize(tmp, &n));
    551:           PetscCall(ISGetIndices(tmp, &idx));
    552:           if (isFE[field]) {
    553:             for (p = 0; p < n; ++p)
    554:               if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
    555:           } else {
    556:             for (p = 0; p < n; ++p)
    557:               if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
    558:           }
    559:           PetscCall(ISRestoreIndices(tmp, &idx));
    560:           PetscCall(ISDestroy(&tmp));
    561:         }
    562:         PetscCall(PetscMalloc1(newn, &newidx));
    563:         newn = 0;
    564:         for (v = 0; v < numValues; ++v) {
    565:           IS              tmp;
    566:           const PetscInt *idx;

    568:           PetscCall(DMLabelGetStratumIS(label, values[v], &tmp));
    569:           if (!tmp) continue;
    570:           PetscCall(ISGetLocalSize(tmp, &n));
    571:           PetscCall(ISGetIndices(tmp, &idx));
    572:           if (isFE[field]) {
    573:             for (p = 0; p < n; ++p)
    574:               if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
    575:           } else {
    576:             for (p = 0; p < n; ++p)
    577:               if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
    578:           }
    579:           PetscCall(ISRestoreIndices(tmp, &idx));
    580:           PetscCall(ISDestroy(&tmp));
    581:         }
    582:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]));
    583:       }
    584:     }
    585:   }
    586:   /* Handle discretization */
    587:   PetscCall(PetscCalloc3(Nf, &labels, Nf, &numComp, Nf * (dim + 1), &numDof));
    588:   for (f = 0; f < Nf; ++f) {
    589:     labels[f] = dm->fields[f].label;
    590:     if (isFE[f]) {
    591:       PetscFE         fe = (PetscFE)dm->fields[f].disc;
    592:       const PetscInt *numFieldDof;
    593:       PetscInt        fedim, d;

    595:       PetscCall(PetscFEGetNumComponents(fe, &numComp[f]));
    596:       PetscCall(PetscFEGetNumDof(fe, &numFieldDof));
    597:       PetscCall(PetscFEGetSpatialDimension(fe, &fedim));
    598:       for (d = 0; d < PetscMin(dim, fedim) + 1; ++d) numDof[f * (dim + 1) + d] = numFieldDof[d];
    599:     } else {
    600:       PetscFV fv = (PetscFV)dm->fields[f].disc;

    602:       PetscCall(PetscFVGetNumComponents(fv, &numComp[f]));
    603:       numDof[f * (dim + 1) + dim] = numComp[f];
    604:     }
    605:   }
    606:   PetscCall(DMPlexGetDepth(dm, &depth));
    607:   for (f = 0; f < Nf; ++f) {
    608:     PetscInt d;
    609:     for (d = 1; d < dim; ++d) PetscCheck(numDof[f * (dim + 1) + d] <= 0 || depth >= dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
    610:   }
    611:   PetscCall(DMCreateSectionPermutation(dm, &permIS, &blockStarts));
    612:   PetscCall(DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, permIS, &section));
    613:   section->blockStarts = blockStarts;
    614:   PetscCall(ISDestroy(&permIS));
    615:   for (f = 0; f < Nf; ++f) {
    616:     PetscFE     fe;
    617:     const char *name;

    619:     PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
    620:     if (!((PetscObject)fe)->name) continue;
    621:     PetscCall(PetscObjectGetName((PetscObject)fe, &name));
    622:     PetscCall(PetscSectionSetFieldName(section, f, name));
    623:   }
    624:   PetscCall(DMSetLocalSection(dm, section));
    625:   PetscCall(PetscSectionDestroy(&section));
    626:   for (bc = 0; bc < numBC; ++bc) {
    627:     PetscCall(ISDestroy(&bcPoints[bc]));
    628:     PetscCall(ISDestroy(&bcComps[bc]));
    629:   }
    630:   PetscCall(PetscFree3(bcFields, bcPoints, bcComps));
    631:   PetscCall(PetscFree3(labels, numComp, numDof));
    632:   PetscCall(PetscFree(isFE));
    633:   /* Checking for CEED usage */
    634:   {
    635:     PetscBool useCeed;

    637:     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
    638:     if (useCeed) {
    639:       PetscCall(DMPlexSetUseMatClosurePermutation(dm, PETSC_FALSE));
    640:       PetscCall(DMUseTensorOrder(dm, PETSC_TRUE));
    641:     }
    642:   }
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }