Actual source code: plexsection.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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        h;
 46:             PetscSectionSym sym;

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

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

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

 91:   DMGetDimension(dm, &dim);
 92:   DMPlexGetDepth(dm, &depth);
 93:   DMPlexGetDepthLabel(dm,&depthLabel);
 94:   DMGetNumFields(dm, &Nf);
 95:   PetscMalloc1(Nf, &isFE);
 96:   for (f = 0; f < Nf; ++f) {
 97:     PetscObject  obj;
 98:     PetscClassId id;

100:     DMGetField(dm, f, NULL, &obj);
101:     PetscObjectGetClassId(obj, &id);
102:     /* User is allowed to put a "placeholder" field in (c.f. DMCreateDS) */
103:     isFE[f] = id == PETSCFE_CLASSID ? PETSC_TRUE : PETSC_FALSE;
104:   }

106:   DMPlexGetVTKCellHeight(dm, &cellHeight);
107:   for (f = 0; f < Nf; ++f) {
108:     if (label && label[f]) {
109:       IS              pointIS;
110:       const PetscInt *points;
111:       PetscInt        n;

113:       DMLabelGetStratumIS(label[f], 1, &pointIS);
114:       if (!pointIS) continue;
115:       ISGetLocalSize(pointIS, &n);
116:       ISGetIndices(pointIS, &points);
117:       for (p = 0; p < n; ++p) {
118:         const PetscInt point = points[p];
119:         PetscInt       dof, d;

121:         DMPlexGetCellType(dm, point, &ct);
122:         DMLabelGetValue(depthLabel, point, &d);
123:         /* If this is a tensor prism point, use dof for one dimension lower */
124:         switch (ct) {
125:           case DM_POLYTOPE_POINT_PRISM_TENSOR:
126:           case DM_POLYTOPE_SEG_PRISM_TENSOR:
127:           case DM_POLYTOPE_TRI_PRISM_TENSOR:
128:           case DM_POLYTOPE_QUAD_PRISM_TENSOR:
129:             --d;break;
130:           default: break;
131:         }
132:         dof  = d < 0 ? 0 : numDof[f*(dim+1)+d];
133:         PetscSectionSetFieldDof(section, point, f, dof);
134:         PetscSectionAddDof(section, point, dof);
135:       }
136:       ISRestoreIndices(pointIS, &points);
137:       ISDestroy(&pointIS);
138:     } else {
139:       for (dep = 0; dep <= depth - cellHeight; ++dep) {
140:         /* Cases: dim > depth (cell-vertex mesh), dim == depth (fully interpolated), dim < depth (interpolated submesh) */
141:         d    = dim <= depth ? dep : (!dep ? 0 : dim);
142:         DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
143:         for (p = pStart; p < pEnd; ++p) {
144:           const PetscInt dof = numDof[f*(dim+1)+d];

146:           DMPlexGetCellType(dm, p, &ct);
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 (isFE[f]) continue;
153:             default: break;
154:           }
155:           PetscSectionSetFieldDof(section, p, f, dof);
156:           PetscSectionAddDof(section, p, dof);
157:         }
158:       }
159:     }
160:   }
161:   PetscFree(isFE);
162:   return(0);
163: }

165: /* Set the number of dof on each point and separate by fields
166:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
167: */
168: static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
169: {
170:   PetscInt       Nf;
171:   PetscInt       bc;
172:   PetscSection   aSec;

176:   PetscSectionGetNumFields(section, &Nf);
177:   for (bc = 0; bc < numBC; ++bc) {
178:     PetscInt        field = 0;
179:     const PetscInt *comp;
180:     const PetscInt *idx;
181:     PetscInt        Nc = 0, cNc = -1, n, i;

183:     if (Nf) {
184:       field = bcField[bc];
185:       PetscSectionGetFieldComponents(section, field, &Nc);
186:     }
187:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &cNc);}
188:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
189:     ISGetLocalSize(bcPoints[bc], &n);
190:     ISGetIndices(bcPoints[bc], &idx);
191:     for (i = 0; i < n; ++i) {
192:       const PetscInt p = idx[i];
193:       PetscInt       numConst;

195:       if (Nf) {
196:         PetscSectionGetFieldDof(section, p, field, &numConst);
197:       } else {
198:         PetscSectionGetDof(section, p, &numConst);
199:       }
200:       /* If Nc <= 0, constrain every dof on the point */
201:       if (cNc > 0) {
202:         /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
203:            and that those dofs are numbered n*Nc+c */
204:         if (Nf) {
205:           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);
206:           numConst = (numConst/Nc) * cNc;
207:         } else {
208:           numConst = PetscMin(numConst, cNc);
209:         }
210:       }
211:       if (Nf) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
212:       PetscSectionAddConstraintDof(section, p, numConst);
213:     }
214:     ISRestoreIndices(bcPoints[bc], &idx);
215:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
216:   }
217:   DMPlexGetAnchors(dm, &aSec, NULL);
218:   if (aSec) {
219:     PetscInt aStart, aEnd, a;

221:     PetscSectionGetChart(aSec, &aStart, &aEnd);
222:     for (a = aStart; a < aEnd; a++) {
223:       PetscInt dof, f;

225:       PetscSectionGetDof(aSec, a, &dof);
226:       if (dof) {
227:         /* if there are point-to-point constraints, then all dofs are constrained */
228:         PetscSectionGetDof(section, a, &dof);
229:         PetscSectionSetConstraintDof(section, a, dof);
230:         for (f = 0; f < Nf; f++) {
231:           PetscSectionGetFieldDof(section, a, f, &dof);
232:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
233:         }
234:       }
235:     }
236:   }
237:   return(0);
238: }

240: /* Set the constrained field indices on each point
241:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
242: */
243: static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
244: {
245:   PetscSection   aSec;
246:   PetscInt      *indices;
247:   PetscInt       Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;

251:   PetscSectionGetNumFields(section, &Nf);
252:   if (!Nf) return(0);
253:   /* Initialize all field indices to -1 */
254:   PetscSectionGetChart(section, &pStart, &pEnd);
255:   for (p = pStart; p < pEnd; ++p) {PetscSectionGetConstraintDof(section, p, &cdof); maxDof = PetscMax(maxDof, cdof);}
256:   PetscMalloc1(maxDof, &indices);
257:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
258:   for (p = pStart; p < pEnd; ++p) for (f = 0; f < Nf; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
259:   /* Handle BC constraints */
260:   for (bc = 0; bc < numBC; ++bc) {
261:     const PetscInt  field = bcField[bc];
262:     const PetscInt *comp, *idx;
263:     PetscInt        Nc, cNc = -1, n, i;

265:     PetscSectionGetFieldComponents(section, field, &Nc);
266:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &cNc);}
267:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
268:     ISGetLocalSize(bcPoints[bc], &n);
269:     ISGetIndices(bcPoints[bc], &idx);
270:     for (i = 0; i < n; ++i) {
271:       const PetscInt  p = idx[i];
272:       const PetscInt *find;
273:       PetscInt        fdof, fcdof, c, j;

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

303:     for (d = 0; d < maxDof; ++d) indices[d] = d;
304:     PetscSectionGetChart(aSec, &aStart, &aEnd);
305:     for (a = aStart; a < aEnd; a++) {
306:       PetscInt dof, f;

308:       PetscSectionGetDof(aSec, a, &dof);
309:       if (dof) {
310:         /* if there are point-to-point constraints, then all dofs are constrained */
311:         for (f = 0; f < Nf; f++) {
312:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
313:         }
314:       }
315:     }
316:   }
317:   PetscFree(indices);
318:   return(0);
319: }

321: /* Set the constrained indices on each point */
322: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
323: {
324:   PetscInt      *indices;
325:   PetscInt       Nf, maxDof, pStart, pEnd, p, f, d;

329:   PetscSectionGetNumFields(section, &Nf);
330:   PetscSectionGetMaxDof(section, &maxDof);
331:   PetscSectionGetChart(section, &pStart, &pEnd);
332:   PetscMalloc1(maxDof, &indices);
333:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
334:   for (p = pStart; p < pEnd; ++p) {
335:     PetscInt cdof, d;

337:     PetscSectionGetConstraintDof(section, p, &cdof);
338:     if (cdof) {
339:       if (Nf) {
340:         PetscInt numConst = 0, foff = 0;

342:         for (f = 0; f < Nf; ++f) {
343:           const PetscInt *find;
344:           PetscInt        fcdof, fdof;

346:           PetscSectionGetFieldDof(section, p, f, &fdof);
347:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
348:           /* Change constraint numbering from field component to local dof number */
349:           PetscSectionGetFieldConstraintIndices(section, p, f, &find);
350:           for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
351:           numConst += fcdof;
352:           foff     += fdof;
353:         }
354:         if (cdof != numConst) {PetscSectionSetConstraintDof(section, p, numConst);}
355:       } else {
356:         for (d = 0; d < cdof; ++d) indices[d] = d;
357:       }
358:       PetscSectionSetConstraintIndices(section, p, indices);
359:     }
360:   }
361:   PetscFree(indices);
362:   return(0);
363: }

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

368:   Not Collective

370:   Input Parameters:
371: + dm        - The DMPlex object
372: . label     - The label indicating the mesh support of each field, or NULL for the whole mesh
373: . numComp   - An array of size numFields that holds the number of components for each field
374: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
375: . numBC     - The number of boundary conditions
376: . bcField   - An array of size numBC giving the field number for each boundry condition
377: . bcComps   - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
378: . bcPoints  - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
379: - perm      - Optional permutation of the chart, or NULL

381:   Output Parameter:
382: . section - The PetscSection object

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

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

390:   Level: developer

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

394: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
395: @*/
396: 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)
397: {
398:   PetscSection   aSec;

402:   DMPlexCreateSectionFields(dm, numComp, section);
403:   DMPlexCreateSectionDof(dm, label, numDof, *section);
404:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
405:   if (perm) {PetscSectionSetPermutation(*section, perm);}
406:   PetscSectionSetFromOptions(*section);
407:   PetscSectionSetUp(*section);
408:   DMPlexGetAnchors(dm,&aSec,NULL);
409:   if (numBC || aSec) {
410:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
411:     DMPlexCreateSectionBCIndices(dm, *section);
412:   }
413:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
414:   return(0);
415: }

417: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
418: {
419:   PetscDS        probBC;
420:   PetscSection   section;
421:   DMLabel       *labels;
422:   IS            *bcPoints, *bcComps;
423:   PetscBool     *isFE;
424:   PetscInt      *bcFields, *numComp, *numDof;
425:   PetscInt       depth, dim, numBd, numBC = 0, Nf, bd, bc = 0, f;
426:   PetscInt       cStart, cEnd, cEndInterior;

430:   DMGetNumFields(dm, &Nf);
431:   /* FE and FV boundary conditions are handled slightly differently */
432:   PetscMalloc1(Nf, &isFE);
433:   for (f = 0; f < Nf; ++f) {
434:     PetscObject  obj;
435:     PetscClassId id;

437:     DMGetField(dm, f, NULL, &obj);
438:     PetscObjectGetClassId(obj, &id);
439:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
440:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
441:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
442:   }
443:   /* Allocate boundary point storage for FEM boundaries */
444:   DMPlexGetDepth(dm, &depth);
445:   DMGetDimension(dm, &dim);
446:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
447:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
448:   DMGetDS(dm, &probBC);
449:   PetscDSGetNumBoundary(probBC, &numBd);
450:   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)");
451:   for (bd = 0; bd < numBd; ++bd) {
452:     PetscInt                field;
453:     DMBoundaryConditionType type;
454:     const char             *labelName;
455:     DMLabel                 label;

457:     PetscDSGetBoundary(probBC, bd, &type, NULL, &labelName, &field, NULL, NULL, NULL, NULL, NULL, NULL);
458:     DMGetLabel(dm,labelName,&label);
459:     if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
460:   }
461:   /* Add ghost cell boundaries for FVM */
462:   for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
463:   PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
464:   /* Constrain ghost cells for FV */
465:   for (f = 0; f < Nf; ++f) {
466:     PetscInt *newidx, c;

468:     if (isFE[f] || cEndInterior < 0) continue;
469:     PetscMalloc1(cEnd-cEndInterior,&newidx);
470:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
471:     bcFields[bc] = f;
472:     ISCreateGeneral(PETSC_COMM_SELF, cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
473:   }
474:   /* Handle FEM Dirichlet boundaries */
475:   for (bd = 0; bd < numBd; ++bd) {
476:     const char             *bdLabel;
477:     DMLabel                 label;
478:     const PetscInt         *comps;
479:     const PetscInt         *values;
480:     PetscInt                bd2, field, numComps, numValues;
481:     DMBoundaryConditionType type;
482:     PetscBool               duplicate = PETSC_FALSE;

484:     PetscDSGetBoundary(probBC, bd, &type, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
485:     DMGetLabel(dm, bdLabel, &label);
486:     if (!isFE[field] || !label) continue;
487:     /* Only want to modify label once */
488:     for (bd2 = 0; bd2 < bd; ++bd2) {
489:       const char *bdname;
490:       PetscDSGetBoundary(probBC, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
491:       PetscStrcmp(bdname, bdLabel, &duplicate);
492:       if (duplicate) break;
493:     }
494:     if (!duplicate && (isFE[field])) {
495:       /* don't complete cells, which are just present to give orientation to the boundary */
496:       DMPlexLabelComplete(dm, label);
497:     }
498:     /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
499:     if (type & DM_BC_ESSENTIAL) {
500:       PetscInt       *newidx;
501:       PetscInt        n, newn = 0, p, v;

503:       bcFields[bc] = field;
504:       if (numComps) {ISCreateGeneral(PETSC_COMM_SELF, numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
505:       for (v = 0; v < numValues; ++v) {
506:         IS              tmp;
507:         const PetscInt *idx;

509:         DMGetStratumIS(dm, bdLabel, values[v], &tmp);
510:         if (!tmp) continue;
511:         ISGetLocalSize(tmp, &n);
512:         ISGetIndices(tmp, &idx);
513:         if (isFE[field]) {
514:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
515:         } else {
516:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
517:         }
518:         ISRestoreIndices(tmp, &idx);
519:         ISDestroy(&tmp);
520:       }
521:       PetscMalloc1(newn,&newidx);
522:       newn = 0;
523:       for (v = 0; v < numValues; ++v) {
524:         IS              tmp;
525:         const PetscInt *idx;

527:         DMGetStratumIS(dm, bdLabel, values[v], &tmp);
528:         if (!tmp) continue;
529:         ISGetLocalSize(tmp, &n);
530:         ISGetIndices(tmp, &idx);
531:         if (isFE[field]) {
532:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
533:         } else {
534:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
535:         }
536:         ISRestoreIndices(tmp, &idx);
537:         ISDestroy(&tmp);
538:       }
539:       ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
540:     }
541:   }
542:   /* Handle discretization */
543:   PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);
544:   for (f = 0; f < Nf; ++f) {
545:     labels[f] = dm->fields[f].label;
546:     if (isFE[f]) {
547:       PetscFE         fe = (PetscFE) dm->fields[f].disc;
548:       const PetscInt *numFieldDof;
549:       PetscInt        fedim, d;

551:       PetscFEGetNumComponents(fe, &numComp[f]);
552:       PetscFEGetNumDof(fe, &numFieldDof);
553:       PetscFEGetSpatialDimension(fe, &fedim);
554:       for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
555:     } else {
556:       PetscFV fv = (PetscFV) dm->fields[f].disc;

558:       PetscFVGetNumComponents(fv, &numComp[f]);
559:       numDof[f*(dim+1)+dim] = numComp[f];
560:     }
561:   }
562:   for (f = 0; f < Nf; ++f) {
563:     PetscInt d;
564:     for (d = 1; d < dim; ++d) {
565:       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.");
566:     }
567:   }
568:   DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
569:   for (f = 0; f < Nf; ++f) {
570:     PetscFE     fe;
571:     const char *name;

573:     DMGetField(dm, f, NULL, (PetscObject *) &fe);
574:     if (!((PetscObject) fe)->name) continue;
575:     PetscObjectGetName((PetscObject) fe, &name);
576:     PetscSectionSetFieldName(section, f, name);
577:   }
578:   DMSetLocalSection(dm, section);
579:   PetscSectionDestroy(&section);
580:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
581:   PetscFree3(bcFields,bcPoints,bcComps);
582:   PetscFree3(labels,numComp,numDof);
583:   PetscFree(isFE);
584:   return(0);
585: }