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;

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

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

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

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

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

 58:               PetscDualSpaceGetHeightSubspace(dspace, sph, &hspace);
 59:               DMPlexGetHeightStratum(K, h, &kStart, &kEnd);
 60:               if (!hspace) continue;
 61:               PetscDualSpaceGetSymmetries(hspace,&perms,&flips);
 62:               if (perms) perms0 = perms[0];
 63:               if (flips) flips0 = flips[0];
 64:               if (!(perms0 || flips0)) continue;
 65:               DMPlexGetConeSize(K,kStart,&kConeSize);
 66:               PetscSectionSymLabelSetStratum(sym,depth - h,numDof[depth - h],-kConeSize,kConeSize,PETSC_USE_POINTER,perms0 ? &perms0[-kConeSize] : NULL,flips0 ? &flips0[-kConeSize] : NULL);
 67:             }
 68:             PetscSectionSetFieldSym(*section,f,sym);
 69:             PetscSectionSymDestroy(&sym);
 70:           }
 71:         }
 72:       }
 73:     }
 74:   }
 75:   DMPlexGetChart(dm, &pStart, &pEnd);
 76:   PetscSectionSetChart(*section, pStart, pEnd);
 77:   PetscFree(isFE);
 78:   return(0);
 79: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

343:   PetscSectionGetNumFields(section, &Nf);
344:   PetscSectionGetMaxDof(section, &maxDof);
345:   PetscSectionGetChart(section, &pStart, &pEnd);
346:   PetscMalloc1(maxDof, &indices);
347:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
348:   for (p = pStart; p < pEnd; ++p) {
349:     PetscInt cdof, d;

351:     PetscSectionGetConstraintDof(section, p, &cdof);
352:     if (cdof) {
353:       if (Nf) {
354:         PetscInt numConst = 0, foff = 0;

356:         for (f = 0; f < Nf; ++f) {
357:           const PetscInt *find;
358:           PetscInt        fcdof, fdof;

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

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

382:   Not Collective

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

395:   Output Parameter:
396: . section - The PetscSection object

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

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

404:   Level: developer

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

408: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
409: @*/
410: 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)
411: {
412:   PetscSection   aSec;

416:   DMPlexCreateSectionFields(dm, numComp, section);
417:   DMPlexCreateSectionDof(dm, label, numDof, *section);
418:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
419:   if (perm) {PetscSectionSetPermutation(*section, perm);}
420:   PetscSectionSetFromOptions(*section);
421:   PetscSectionSetUp(*section);
422:   DMPlexGetAnchors(dm,&aSec,NULL);
423:   if (numBC || aSec) {
424:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
425:     DMPlexCreateSectionBCIndices(dm, *section);
426:   }
427:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
428:   return(0);
429: }

431: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
432: {
433:   PetscSection   section;
434:   DMLabel       *labels;
435:   IS            *bcPoints, *bcComps;
436:   PetscBool     *isFE;
437:   PetscInt      *bcFields, *numComp, *numDof;
438:   PetscInt       depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f;
439:   PetscInt       cStart, cEnd, cEndInterior;

443:   DMGetNumFields(dm, &Nf);
444:   DMGetDimension(dm, &dim);
445:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
446:   /* FE and FV boundary conditions are handled slightly differently */
447:   PetscMalloc1(Nf, &isFE);
448:   for (f = 0; f < Nf; ++f) {
449:     PetscObject  obj;
450:     PetscClassId id;

452:     DMGetField(dm, f, NULL, &obj);
453:     PetscObjectGetClassId(obj, &id);
454:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
455:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
456:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
457:   }
458:   /* Allocate boundary point storage for FEM boundaries */
459:   DMGetNumDS(dm, &Nds);
460:   for (s = 0; s < Nds; ++s) {
461:     PetscDS  dsBC;
462:     PetscInt numBd, bd;

464:     DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
465:     PetscDSGetNumBoundary(dsBC, &numBd);
466:     if (!Nf && numBd) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "number of fields is zero and number of boundary conditions is nonzero (this should never happen)");
467:     for (bd = 0; bd < numBd; ++bd) {
468:       PetscInt                field;
469:       DMBoundaryConditionType type;
470:       const char             *labelName;
471:       DMLabel                 label;

473:       PetscDSGetBoundary(dsBC, bd, &type, NULL, &labelName, &field, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
474:       DMGetLabel(dm, labelName, &label);
475:       if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
476:     }
477:   }
478:   /* Add ghost cell boundaries for FVM */
479:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
480:   for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
481:   PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps);
482:   /* Constrain ghost cells for FV */
483:   for (f = 0; f < Nf; ++f) {
484:     PetscInt *newidx, c;

486:     if (isFE[f] || cEndInterior < 0) continue;
487:     PetscMalloc1(cEnd-cEndInterior,&newidx);
488:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
489:     bcFields[bc] = f;
490:     ISCreateGeneral(PETSC_COMM_SELF, cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
491:   }
492:   /* Handle FEM Dirichlet boundaries */
493:   DMGetNumDS(dm, &Nds);
494:   for (s = 0; s < Nds; ++s) {
495:     PetscDS  dsBC;
496:     PetscInt numBd, bd;

498:     DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
499:     PetscDSGetNumBoundary(dsBC, &numBd);
500:     for (bd = 0; bd < numBd; ++bd) {
501:       const char             *bdLabel;
502:       DMLabel                 label;
503:       const PetscInt         *comps;
504:       const PetscInt         *values;
505:       PetscInt                bd2, field, numComps, numValues;
506:       DMBoundaryConditionType type;
507:       PetscBool               duplicate = PETSC_FALSE;

509:       PetscDSGetBoundary(dsBC, bd, &type, NULL, &bdLabel, &field, &numComps, &comps, NULL, NULL, &numValues, &values, NULL);
510:       DMGetLabel(dm, bdLabel, &label);
511:       if (!isFE[field] || !label) continue;
512:       /* Only want to modify label once */
513:       for (bd2 = 0; bd2 < bd; ++bd2) {
514:         const char *bdname;
515:         PetscDSGetBoundary(dsBC, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
516:         PetscStrcmp(bdname, bdLabel, &duplicate);
517:         if (duplicate) break;
518:       }
519:       if (!duplicate && (isFE[field])) {
520:         /* don't complete cells, which are just present to give orientation to the boundary */
521:         DMPlexLabelComplete(dm, label);
522:       }
523:       /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
524:       if (type & DM_BC_ESSENTIAL) {
525:         PetscInt       *newidx;
526:         PetscInt        n, newn = 0, p, v;

528:         bcFields[bc] = field;
529:         if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
530:         for (v = 0; v < numValues; ++v) {
531:           IS              tmp;
532:           const PetscInt *idx;

534:           DMGetStratumIS(dm, bdLabel, values[v], &tmp);
535:           if (!tmp) continue;
536:           ISGetLocalSize(tmp, &n);
537:           ISGetIndices(tmp, &idx);
538:           if (isFE[field]) {
539:             for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
540:           } else {
541:             for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
542:           }
543:           ISRestoreIndices(tmp, &idx);
544:           ISDestroy(&tmp);
545:         }
546:         PetscMalloc1(newn, &newidx);
547:         newn = 0;
548:         for (v = 0; v < numValues; ++v) {
549:           IS              tmp;
550:           const PetscInt *idx;

552:           DMGetStratumIS(dm, bdLabel, values[v], &tmp);
553:           if (!tmp) continue;
554:           ISGetLocalSize(tmp, &n);
555:           ISGetIndices(tmp, &idx);
556:           if (isFE[field]) {
557:             for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
558:           } else {
559:             for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
560:           }
561:           ISRestoreIndices(tmp, &idx);
562:           ISDestroy(&tmp);
563:         }
564:         ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
565:       }
566:     }
567:   }
568:   /* Handle discretization */
569:   PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);
570:   for (f = 0; f < Nf; ++f) {
571:     labels[f] = dm->fields[f].label;
572:     if (isFE[f]) {
573:       PetscFE         fe = (PetscFE) dm->fields[f].disc;
574:       const PetscInt *numFieldDof;
575:       PetscInt        fedim, d;

577:       PetscFEGetNumComponents(fe, &numComp[f]);
578:       PetscFEGetNumDof(fe, &numFieldDof);
579:       PetscFEGetSpatialDimension(fe, &fedim);
580:       for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
581:     } else {
582:       PetscFV fv = (PetscFV) dm->fields[f].disc;

584:       PetscFVGetNumComponents(fv, &numComp[f]);
585:       numDof[f*(dim+1)+dim] = numComp[f];
586:     }
587:   }
588:   DMPlexGetDepth(dm, &depth);
589:   for (f = 0; f < Nf; ++f) {
590:     PetscInt d;
591:     for (d = 1; d < dim; ++d) {
592:       if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
593:     }
594:   }
595:   DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
596:   for (f = 0; f < Nf; ++f) {
597:     PetscFE     fe;
598:     const char *name;

600:     DMGetField(dm, f, NULL, (PetscObject *) &fe);
601:     if (!((PetscObject) fe)->name) continue;
602:     PetscObjectGetName((PetscObject) fe, &name);
603:     PetscSectionSetFieldName(section, f, name);
604:   }
605:   DMSetLocalSection(dm, section);
606:   PetscSectionDestroy(&section);
607:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
608:   PetscFree3(bcFields,bcPoints,bcComps);
609:   PetscFree3(labels,numComp,numDof);
610:   PetscFree(isFE);
611:   return(0);
612: }