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:   DMPlexGetDepth(dm, &depth);
 11:   DMPlexGetDepthLabel(dm, &depthLabel);
 12:   DMGetNumFields(dm, &Nf);
 13:   PetscCalloc1(Nf, &isFE);
 14:   for (f = 0; f < Nf; ++f) {
 15:     PetscObject  obj;
 16:     PetscClassId id;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

245:     PetscSectionGetChart(aSec, &aStart, &aEnd);
246:     for (a = aStart; a < aEnd; a++) {
247:       PetscInt dof, f;

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

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

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

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

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

335:     for (d = 0; d < maxDof; ++d) indices[d] = d;
336:     PetscSectionGetChart(aSec, &aStart, &aEnd);
337:     for (a = aStart; a < aEnd; a++) {
338:       PetscInt dof, f;

340:       PetscSectionGetDof(aSec, a, &dof);
341:       if (dof) {
342:         /* if there are point-to-point constraints, then all dofs are constrained */
343:         for (f = 0; f < Nf; f++) PetscSectionSetFieldConstraintIndices(section, a, f, indices);
344:       }
345:     }
346:   }
347:   PetscFree(indices);
348:   return 0;
349: }

351: /* Set the constrained indices on each point */
352: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
353: {
354:   PetscInt *indices;
355:   PetscInt  Nf, maxDof, pStart, pEnd, p, f, d;

357:   PetscSectionGetNumFields(section, &Nf);
358:   PetscSectionGetMaxDof(section, &maxDof);
359:   PetscSectionGetChart(section, &pStart, &pEnd);
360:   PetscMalloc1(maxDof, &indices);
361:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
362:   for (p = pStart; p < pEnd; ++p) {
363:     PetscInt cdof, d;

365:     PetscSectionGetConstraintDof(section, p, &cdof);
366:     if (cdof) {
367:       if (Nf) {
368:         PetscInt numConst = 0, foff = 0;

370:         for (f = 0; f < Nf; ++f) {
371:           const PetscInt *find;
372:           PetscInt        fcdof, fdof;

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

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

396:   Not Collective

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

409:   Output Parameter:
410: . section - The `PetscSection` object

412:   Level: developer

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

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

420:   Developer Note:
421:   This is used by `DMCreateLocalSection()`?

423: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
424: @*/
425: 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)
426: {
427:   PetscSection aSec;

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

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

454:   DMGetNumFields(dm, &Nf);
455:   DMGetDimension(dm, &dim);
456:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
457:   /* FE and FV boundary conditions are handled slightly differently */
458:   PetscMalloc1(Nf, &isFE);
459:   for (f = 0; f < Nf; ++f) {
460:     PetscObject  obj;
461:     PetscClassId id;

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

477:     DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
478:     PetscDSGetNumBoundary(dsBC, &numBd);
480:     for (bd = 0; bd < numBd; ++bd) {
481:       PetscInt                field;
482:       DMBoundaryConditionType type;
483:       DMLabel                 label;

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

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

512:     DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
513:     PetscDSGetNumBoundary(dsBC, &numBd);
514:     for (bd = 0; bd < numBd; ++bd) {
515:       DMLabel                 label;
516:       const PetscInt         *comps;
517:       const PetscInt         *values;
518:       PetscInt                bd2, field, numComps, numValues;
519:       DMBoundaryConditionType type;
520:       PetscBool               duplicate = PETSC_FALSE;

522:       PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, &numValues, &values, &field, &numComps, &comps, NULL, NULL, NULL);
523:       if (!isFE[field] || !label) continue;
524:       /* Only want to modify label once */
525:       for (bd2 = 0; bd2 < bd; ++bd2) {
526:         DMLabel l;

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

537:         bcFields[bc] = field;
538:         if (numComps) ISCreateGeneral(PETSC_COMM_SELF, numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);
539:         for (v = 0; v < numValues; ++v) {
540:           IS              tmp;
541:           const PetscInt *idx;

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

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

590:       PetscFEGetNumComponents(fe, &numComp[f]);
591:       PetscFEGetNumDof(fe, &numFieldDof);
592:       PetscFEGetSpatialDimension(fe, &fedim);
593:       for (d = 0; d < PetscMin(dim, fedim) + 1; ++d) numDof[f * (dim + 1) + d] = numFieldDof[d];
594:     } else {
595:       PetscFV fv = (PetscFV)dm->fields[f].disc;

597:       PetscFVGetNumComponents(fv, &numComp[f]);
598:       numDof[f * (dim + 1) + dim] = numComp[f];
599:     }
600:   }
601:   DMPlexGetDepth(dm, &depth);
602:   for (f = 0; f < Nf; ++f) {
603:     PetscInt d;
605:   }
606:   DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
607:   for (f = 0; f < Nf; ++f) {
608:     PetscFE     fe;
609:     const char *name;

611:     DMGetField(dm, f, NULL, (PetscObject *)&fe);
612:     if (!((PetscObject)fe)->name) continue;
613:     PetscObjectGetName((PetscObject)fe, &name);
614:     PetscSectionSetFieldName(section, f, name);
615:   }
616:   DMSetLocalSection(dm, section);
617:   PetscSectionDestroy(&section);
618:   for (bc = 0; bc < numBC; ++bc) {
619:     ISDestroy(&bcPoints[bc]);
620:     ISDestroy(&bcComps[bc]);
621:   }
622:   PetscFree3(bcFields, bcPoints, bcComps);
623:   PetscFree3(labels, numComp, numDof);
624:   PetscFree(isFE);
625:   return 0;
626: }