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 = DMPolytopeTypeGetNumArrangments(ct) / 2;
 72:               }
 73:               PetscCall(PetscSectionSymLabelSetStratum(sym, depth - h, numDof[depth - h], -kConeSize, kConeSize, PETSC_USE_POINTER, perms0 ? &perms0[-kConeSize] : NULL, flips0 ? &flips0[-kConeSize] : NULL));
 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:   Developer Notes:
426:   This is used by `DMCreateLocalSection()`?

428: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
429: @*/
430: 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)
431: {
432:   PetscSection aSec;

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

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

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

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

484:     PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC, NULL));
485:     PetscCall(PetscDSGetNumBoundary(dsBC, &numBd));
486:     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)");
487:     for (bd = 0; bd < numBd; ++bd) {
488:       PetscInt                field;
489:       DMBoundaryConditionType type;
490:       DMLabel                 label;

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

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

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

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

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

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

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

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

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

604:       PetscCall(PetscFVGetNumComponents(fv, &numComp[f]));
605:       numDof[f * (dim + 1) + dim] = numComp[f];
606:     }
607:   }
608:   PetscCall(DMPlexGetDepth(dm, &depth));
609:   for (f = 0; f < Nf; ++f) {
610:     PetscInt d;
611:     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.");
612:   }
613:   PetscCall(DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section));
614:   for (f = 0; f < Nf; ++f) {
615:     PetscFE     fe;
616:     const char *name;

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

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