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)      {isFE[f] = PETSC_TRUE;}
 21:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
 22:   }

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

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

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

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

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

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

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

111:     DMGetField(dm, f, NULL, &obj);
112:     PetscObjectGetClassId(obj, &id);
113:     /* User is allowed to put a "placeholder" field in (c.f. DMCreateDS) */
114:     isFE[f] = id == PETSCFE_CLASSID ? PETSC_TRUE : PETSC_FALSE;
115:   }

117:   DMPlexGetVTKCellHeight(dm, &cellHeight);
118:   for (f = 0; f < Nf; ++f) {
119:     PetscBool avoidTensor;

121:     DMGetFieldAvoidTensor(dm, f, &avoidTensor);
122:     avoidTensor = (avoidTensor || hasCohesive) ? PETSC_TRUE : PETSC_FALSE;
123:     if (label && label[f]) {
124:       IS              pointIS;
125:       const PetscInt *points;
126:       PetscInt        n;

128:       DMLabelGetStratumIS(label[f], 1, &pointIS);
129:       if (!pointIS) continue;
130:       ISGetLocalSize(pointIS, &n);
131:       ISGetIndices(pointIS, &points);
132:       for (p = 0; p < n; ++p) {
133:         const PetscInt point = points[p];
134:         PetscInt       dof, d;

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

161:           DMPlexGetCellType(dm, p, &ct);
162:           switch (ct) {
163:             case DM_POLYTOPE_POINT_PRISM_TENSOR:
164:             case DM_POLYTOPE_SEG_PRISM_TENSOR:
165:             case DM_POLYTOPE_TRI_PRISM_TENSOR:
166:             case DM_POLYTOPE_QUAD_PRISM_TENSOR:
167:               if (avoidTensor && isFE[f]) continue;
168:             default: break;
169:           }
170:           PetscSectionSetFieldDof(section, p, f, dof);
171:           PetscSectionAddDof(section, p, dof);
172:         }
173:       }
174:     }
175:   }
176:   PetscFree(isFE);
177:   return 0;
178: }

180: /* Set the number of dof on each point and separate by fields
181:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
182: */
183: static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
184: {
185:   PetscInt       Nf;
186:   PetscInt       bc;
187:   PetscSection   aSec;

189:   PetscSectionGetNumFields(section, &Nf);
190:   for (bc = 0; bc < numBC; ++bc) {
191:     PetscInt        field = 0;
192:     const PetscInt *comp;
193:     const PetscInt *idx;
194:     PetscInt        Nc = 0, cNc = -1, n, i;

196:     if (Nf) {
197:       field = bcField[bc];
198:       PetscSectionGetFieldComponents(section, field, &Nc);
199:     }
200:     if (bcComps && bcComps[bc]) ISGetLocalSize(bcComps[bc], &cNc);
201:     if (bcComps && bcComps[bc]) ISGetIndices(bcComps[bc], &comp);
202:     ISGetLocalSize(bcPoints[bc], &n);
203:     ISGetIndices(bcPoints[bc], &idx);
204:     for (i = 0; i < n; ++i) {
205:       const PetscInt p = idx[i];
206:       PetscInt       numConst;

208:       if (Nf) {
209:         PetscSectionGetFieldDof(section, p, field, &numConst);
210:       } else {
211:         PetscSectionGetDof(section, p, &numConst);
212:       }
213:       /* If Nc <= 0, constrain every dof on the point */
214:       if (cNc > 0) {
215:         /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
216:            and that those dofs are numbered n*Nc+c */
217:         if (Nf) {
219:           numConst = (numConst/Nc) * cNc;
220:         } else {
221:           numConst = PetscMin(numConst, cNc);
222:         }
223:       }
224:       if (Nf) PetscSectionAddFieldConstraintDof(section, p, field, numConst);
225:       PetscSectionAddConstraintDof(section, p, numConst);
226:     }
227:     ISRestoreIndices(bcPoints[bc], &idx);
228:     if (bcComps && bcComps[bc]) ISRestoreIndices(bcComps[bc], &comp);
229:   }
230:   DMPlexGetAnchors(dm, &aSec, NULL);
231:   if (aSec) {
232:     PetscInt aStart, aEnd, a;

234:     PetscSectionGetChart(aSec, &aStart, &aEnd);
235:     for (a = aStart; a < aEnd; a++) {
236:       PetscInt dof, f;

238:       PetscSectionGetDof(aSec, a, &dof);
239:       if (dof) {
240:         /* if there are point-to-point constraints, then all dofs are constrained */
241:         PetscSectionGetDof(section, a, &dof);
242:         PetscSectionSetConstraintDof(section, a, dof);
243:         for (f = 0; f < Nf; f++) {
244:           PetscSectionGetFieldDof(section, a, f, &dof);
245:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
246:         }
247:       }
248:     }
249:   }
250:   return 0;
251: }

253: /* Set the constrained field indices on each point
254:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
255: */
256: static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
257: {
258:   PetscSection   aSec;
259:   PetscInt      *indices;
260:   PetscInt       Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;

262:   PetscSectionGetNumFields(section, &Nf);
263:   if (!Nf) return 0;
264:   /* Initialize all field indices to -1 */
265:   PetscSectionGetChart(section, &pStart, &pEnd);
266:   for (p = pStart; p < pEnd; ++p) {PetscSectionGetConstraintDof(section, p, &cdof)); maxDof = PetscMax(maxDof, cdof;}
267:   PetscMalloc1(maxDof, &indices);
268:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
269:   for (p = pStart; p < pEnd; ++p) for (f = 0; f < Nf; ++f) PetscSectionSetFieldConstraintIndices(section, p, f, indices);
270:   /* Handle BC constraints */
271:   for (bc = 0; bc < numBC; ++bc) {
272:     const PetscInt  field = bcField[bc];
273:     const PetscInt *comp, *idx;
274:     PetscInt        Nc, cNc = -1, n, i;

276:     PetscSectionGetFieldComponents(section, field, &Nc);
277:     if (bcComps && bcComps[bc]) ISGetLocalSize(bcComps[bc], &cNc);
278:     if (bcComps && bcComps[bc]) ISGetIndices(bcComps[bc], &comp);
279:     ISGetLocalSize(bcPoints[bc], &n);
280:     ISGetIndices(bcPoints[bc], &idx);
281:     for (i = 0; i < n; ++i) {
282:       const PetscInt  p = idx[i];
283:       const PetscInt *find;
284:       PetscInt        fdof, fcdof, c, j;

286:       PetscSectionGetFieldDof(section, p, field, &fdof);
287:       if (!fdof) continue;
288:       if (cNc < 0) {
289:         for (d = 0; d < fdof; ++d) indices[d] = d;
290:         fcdof = fdof;
291:       } else {
292:         /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
293:            and that those dofs are numbered n*Nc+c */
294:         PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
295:         PetscSectionGetFieldConstraintIndices(section, p, field, &find);
296:         /* Get indices constrained by previous bcs */
297:         for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
298:         for (j = 0; j < fdof/Nc; ++j) for (c = 0; c < cNc; ++c) indices[d++] = j*Nc + comp[c];
299:         PetscSortRemoveDupsInt(&d, indices);
300:         for (c = d; c < fcdof; ++c) indices[c] = -1;
301:         fcdof = d;
302:       }
303:       PetscSectionSetFieldConstraintDof(section, p, field, fcdof);
304:       PetscSectionSetFieldConstraintIndices(section, p, field, indices);
305:     }
306:     if (bcComps && bcComps[bc]) ISRestoreIndices(bcComps[bc], &comp);
307:     ISRestoreIndices(bcPoints[bc], &idx);
308:   }
309:   /* Handle anchors */
310:   DMPlexGetAnchors(dm, &aSec, NULL);
311:   if (aSec) {
312:     PetscInt aStart, aEnd, a;

314:     for (d = 0; d < maxDof; ++d) indices[d] = d;
315:     PetscSectionGetChart(aSec, &aStart, &aEnd);
316:     for (a = aStart; a < aEnd; a++) {
317:       PetscInt dof, f;

319:       PetscSectionGetDof(aSec, a, &dof);
320:       if (dof) {
321:         /* if there are point-to-point constraints, then all dofs are constrained */
322:         for (f = 0; f < Nf; f++) {
323:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
324:         }
325:       }
326:     }
327:   }
328:   PetscFree(indices);
329:   return 0;
330: }

332: /* Set the constrained indices on each point */
333: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
334: {
335:   PetscInt      *indices;
336:   PetscInt       Nf, maxDof, pStart, pEnd, p, f, d;

338:   PetscSectionGetNumFields(section, &Nf);
339:   PetscSectionGetMaxDof(section, &maxDof);
340:   PetscSectionGetChart(section, &pStart, &pEnd);
341:   PetscMalloc1(maxDof, &indices);
342:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
343:   for (p = pStart; p < pEnd; ++p) {
344:     PetscInt cdof, d;

346:     PetscSectionGetConstraintDof(section, p, &cdof);
347:     if (cdof) {
348:       if (Nf) {
349:         PetscInt numConst = 0, foff = 0;

351:         for (f = 0; f < Nf; ++f) {
352:           const PetscInt *find;
353:           PetscInt        fcdof, fdof;

355:           PetscSectionGetFieldDof(section, p, f, &fdof);
356:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
357:           /* Change constraint numbering from field component to local dof number */
358:           PetscSectionGetFieldConstraintIndices(section, p, f, &find);
359:           for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
360:           numConst += fcdof;
361:           foff     += fdof;
362:         }
363:         if (cdof != numConst) PetscSectionSetConstraintDof(section, p, numConst);
364:       } else {
365:         for (d = 0; d < cdof; ++d) indices[d] = d;
366:       }
367:       PetscSectionSetConstraintIndices(section, p, indices);
368:     }
369:   }
370:   PetscFree(indices);
371:   return 0;
372: }

374: /*@C
375:   DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.

377:   Not Collective

379:   Input Parameters:
380: + dm        - The DMPlex object
381: . label     - The label indicating the mesh support of each field, or NULL for the whole mesh
382: . numComp   - An array of size numFields that holds the number of components for each field
383: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
384: . numBC     - The number of boundary conditions
385: . bcField   - An array of size numBC giving the field number for each boundry condition
386: . bcComps   - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
387: . bcPoints  - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
388: - perm      - Optional permutation of the chart, or NULL

390:   Output Parameter:
391: . section - The PetscSection object

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

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

399:   Level: developer

401:   TODO: How is this related to DMCreateLocalSection()

403: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
404: @*/
405: 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)
406: {
407:   PetscSection   aSec;

409:   DMPlexCreateSectionFields(dm, numComp, section);
410:   DMPlexCreateSectionDof(dm, label, numDof, *section);
411:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
412:   if (perm) PetscSectionSetPermutation(*section, perm);
413:   PetscSectionSetFromOptions(*section);
414:   PetscSectionSetUp(*section);
415:   DMPlexGetAnchors(dm,&aSec,NULL);
416:   if (numBC || aSec) {
417:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
418:     DMPlexCreateSectionBCIndices(dm, *section);
419:   }
420:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
421:   return 0;
422: }

424: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
425: {
426:   PetscSection   section;
427:   DMLabel       *labels;
428:   IS            *bcPoints, *bcComps;
429:   PetscBool     *isFE;
430:   PetscInt      *bcFields, *numComp, *numDof;
431:   PetscInt       depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f;
432:   PetscInt       cStart, cEnd, cEndInterior;

434:   DMGetNumFields(dm, &Nf);
435:   DMGetDimension(dm, &dim);
436:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
437:   /* FE and FV boundary conditions are handled slightly differently */
438:   PetscMalloc1(Nf, &isFE);
439:   for (f = 0; f < Nf; ++f) {
440:     PetscObject  obj;
441:     PetscClassId id;

443:     DMGetField(dm, f, NULL, &obj);
444:     PetscObjectGetClassId(obj, &id);
445:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
446:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
447:     else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
448:   }
449:   /* Allocate boundary point storage for FEM boundaries */
450:   DMGetNumDS(dm, &Nds);
451:   for (s = 0; s < Nds; ++s) {
452:     PetscDS  dsBC;
453:     PetscInt numBd, bd;

455:     DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
456:     PetscDSGetNumBoundary(dsBC, &numBd);
458:     for (bd = 0; bd < numBd; ++bd) {
459:       PetscInt                field;
460:       DMBoundaryConditionType type;
461:       DMLabel                 label;

463:       PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL);
464:       if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
465:     }
466:   }
467:   /* Add ghost cell boundaries for FVM */
468:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
469:   for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
470:   PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps);
471:   /* Constrain ghost cells for FV */
472:   for (f = 0; f < Nf; ++f) {
473:     PetscInt *newidx, c;

475:     if (isFE[f] || cEndInterior < 0) continue;
476:     PetscMalloc1(cEnd-cEndInterior,&newidx);
477:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
478:     bcFields[bc] = f;
479:     ISCreateGeneral(PETSC_COMM_SELF, cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
480:   }
481:   /* Handle FEM Dirichlet boundaries */
482:   DMGetNumDS(dm, &Nds);
483:   for (s = 0; s < Nds; ++s) {
484:     PetscDS  dsBC;
485:     PetscInt numBd, bd;

487:     DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
488:     PetscDSGetNumBoundary(dsBC, &numBd);
489:     for (bd = 0; bd < numBd; ++bd) {
490:       DMLabel                 label;
491:       const PetscInt         *comps;
492:       const PetscInt         *values;
493:       PetscInt                bd2, field, numComps, numValues;
494:       DMBoundaryConditionType type;
495:       PetscBool               duplicate = PETSC_FALSE;

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

503:         PetscDSGetBoundary(dsBC, bd2, NULL, NULL, NULL, &l, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
504:         duplicate = l == label ? PETSC_TRUE : PETSC_FALSE;
505:         if (duplicate) break;
506:       }
507:       if (!duplicate && (isFE[field])) {
508:         /* don't complete cells, which are just present to give orientation to the boundary */
509:         DMPlexLabelComplete(dm, label);
510:       }
511:       /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
512:       if (type & DM_BC_ESSENTIAL) {
513:         PetscInt       *newidx;
514:         PetscInt        n, newn = 0, p, v;

516:         bcFields[bc] = field;
517:         if (numComps) ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);
518:         for (v = 0; v < numValues; ++v) {
519:           IS              tmp;
520:           const PetscInt *idx;

522:           DMLabelGetStratumIS(label, values[v], &tmp);
523:           if (!tmp) continue;
524:           ISGetLocalSize(tmp, &n);
525:           ISGetIndices(tmp, &idx);
526:           if (isFE[field]) {
527:             for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
528:           } else {
529:             for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
530:           }
531:           ISRestoreIndices(tmp, &idx);
532:           ISDestroy(&tmp);
533:         }
534:         PetscMalloc1(newn, &newidx);
535:         newn = 0;
536:         for (v = 0; v < numValues; ++v) {
537:           IS              tmp;
538:           const PetscInt *idx;

540:           DMLabelGetStratumIS(label, values[v], &tmp);
541:           if (!tmp) continue;
542:           ISGetLocalSize(tmp, &n);
543:           ISGetIndices(tmp, &idx);
544:           if (isFE[field]) {
545:             for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
546:           } else {
547:             for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
548:           }
549:           ISRestoreIndices(tmp, &idx);
550:           ISDestroy(&tmp);
551:         }
552:         ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
553:       }
554:     }
555:   }
556:   /* Handle discretization */
557:   PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);
558:   for (f = 0; f < Nf; ++f) {
559:     labels[f] = dm->fields[f].label;
560:     if (isFE[f]) {
561:       PetscFE         fe = (PetscFE) dm->fields[f].disc;
562:       const PetscInt *numFieldDof;
563:       PetscInt        fedim, d;

565:       PetscFEGetNumComponents(fe, &numComp[f]);
566:       PetscFEGetNumDof(fe, &numFieldDof);
567:       PetscFEGetSpatialDimension(fe, &fedim);
568:       for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
569:     } else {
570:       PetscFV fv = (PetscFV) dm->fields[f].disc;

572:       PetscFVGetNumComponents(fv, &numComp[f]);
573:       numDof[f*(dim+1)+dim] = numComp[f];
574:     }
575:   }
576:   DMPlexGetDepth(dm, &depth);
577:   for (f = 0; f < Nf; ++f) {
578:     PetscInt d;
579:     for (d = 1; d < dim; ++d) {
581:     }
582:   }
583:   DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
584:   for (f = 0; f < Nf; ++f) {
585:     PetscFE     fe;
586:     const char *name;

588:     DMGetField(dm, f, NULL, (PetscObject *) &fe);
589:     if (!((PetscObject) fe)->name) continue;
590:     PetscObjectGetName((PetscObject) fe, &name);
591:     PetscSectionSetFieldName(section, f, name);
592:   }
593:   DMSetLocalSection(dm, section);
594:   PetscSectionDestroy(&section);
595:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]));PetscCall(ISDestroy(&bcComps[bc]);}
596:   PetscFree3(bcFields,bcPoints,bcComps);
597:   PetscFree3(labels,numComp,numDof);
598:   PetscFree(isFE);
599:   return 0;
600: }