Actual source code: plexsection.c
petsc-3.12.5 2020-03-29
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: PetscInt *pMax;
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: PetscMalloc1(depth+1, &pMax);
107: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
108: DMPlexGetVTKCellHeight(dm, &cellHeight);
109: for (f = 0; f < Nf; ++f) {
110: if (label && label[f]) {
111: IS pointIS;
112: const PetscInt *points;
113: PetscInt n;
115: DMLabelGetStratumIS(label[f], 1, &pointIS);
116: if (!pointIS) continue;
117: ISGetLocalSize(pointIS, &n);
118: ISGetIndices(pointIS, &points);
119: for (p = 0; p < n; ++p) {
120: const PetscInt point = points[p];
121: PetscInt dof, d;
123: DMLabelGetValue(depthLabel, point, &d);
124: DMPlexGetDepthStratum(dm, d, NULL, &pEnd);
125: /* If this is a hybrid point, use dof for one dimension lower */
126: if ((point >= pMax[d]) && (point < pEnd)) --d;
127: dof = d < 0 ? 0 : numDof[f*(dim+1)+d];
128: PetscSectionSetFieldDof(section, point, f, dof);
129: PetscSectionAddDof(section, point, dof);
130: }
131: ISRestoreIndices(pointIS, &points);
132: ISDestroy(&pointIS);
133: } else {
134: for (dep = 0; dep <= depth - cellHeight; ++dep) {
135: /* Cases: dim > depth (cell-vertex mesh), dim == depth (fully interpolated), dim < depth (interpolated submesh) */
136: d = dim <= depth ? dep : (!dep ? 0 : dim);
137: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
138: pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
139: for (p = pStart; p < pEnd; ++p) {
140: const PetscInt dof = numDof[f*(dim+1)+d];
142: if (isFE[f] && p >= pMax[dep]) continue;
143: PetscSectionSetFieldDof(section, p, f, dof);
144: PetscSectionAddDof(section, p, dof);
145: }
146: }
147: }
148: }
149: PetscFree(pMax);
150: PetscFree(isFE);
151: return(0);
152: }
154: /* Set the number of dof on each point and separate by fields
155: If bcComps is NULL or the IS is NULL, constrain every dof on the point
156: */
157: static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
158: {
159: PetscInt Nf;
160: PetscInt bc;
161: PetscSection aSec;
165: PetscSectionGetNumFields(section, &Nf);
166: for (bc = 0; bc < numBC; ++bc) {
167: PetscInt field = 0;
168: const PetscInt *comp;
169: const PetscInt *idx;
170: PetscInt Nc = -1, n, i;
172: if (Nf) field = bcField[bc];
173: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
174: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
175: ISGetLocalSize(bcPoints[bc], &n);
176: ISGetIndices(bcPoints[bc], &idx);
177: for (i = 0; i < n; ++i) {
178: const PetscInt p = idx[i];
179: PetscInt numConst;
181: if (Nf) {
182: PetscSectionGetFieldDof(section, p, field, &numConst);
183: } else {
184: PetscSectionGetDof(section, p, &numConst);
185: }
186: /* If Nc < 0, constrain every dof on the point */
187: /* TODO: Matt, this only works if there is one node on the point. We need to handle numDofs > NumComponents */
188: if (Nc > 0) numConst = PetscMin(numConst, Nc);
189: if (Nf) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
190: PetscSectionAddConstraintDof(section, p, numConst);
191: }
192: ISRestoreIndices(bcPoints[bc], &idx);
193: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
194: }
195: DMPlexGetAnchors(dm, &aSec, NULL);
196: if (aSec) {
197: PetscInt aStart, aEnd, a;
199: PetscSectionGetChart(aSec, &aStart, &aEnd);
200: for (a = aStart; a < aEnd; a++) {
201: PetscInt dof, f;
203: PetscSectionGetDof(aSec, a, &dof);
204: if (dof) {
205: /* if there are point-to-point constraints, then all dofs are constrained */
206: PetscSectionGetDof(section, a, &dof);
207: PetscSectionSetConstraintDof(section, a, dof);
208: for (f = 0; f < Nf; f++) {
209: PetscSectionGetFieldDof(section, a, f, &dof);
210: PetscSectionSetFieldConstraintDof(section, a, f, dof);
211: }
212: }
213: }
214: }
215: return(0);
216: }
218: /* Set the constrained field indices on each point
219: If bcComps is NULL or the IS is NULL, constrain every dof on the point
220: */
221: static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
222: {
223: PetscSection aSec;
224: PetscInt *indices;
225: PetscInt Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;
229: PetscSectionGetNumFields(section, &Nf);
230: if (!Nf) return(0);
231: /* Initialize all field indices to -1 */
232: PetscSectionGetChart(section, &pStart, &pEnd);
233: for (p = pStart; p < pEnd; ++p) {PetscSectionGetConstraintDof(section, p, &cdof); maxDof = PetscMax(maxDof, cdof);}
234: PetscMalloc1(maxDof, &indices);
235: for (d = 0; d < maxDof; ++d) indices[d] = -1;
236: for (p = pStart; p < pEnd; ++p) for (f = 0; f < Nf; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
237: /* Handle BC constraints */
238: for (bc = 0; bc < numBC; ++bc) {
239: const PetscInt field = bcField[bc];
240: const PetscInt *comp, *idx;
241: PetscInt Nc = -1, n, i;
243: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
244: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
245: ISGetLocalSize(bcPoints[bc], &n);
246: ISGetIndices(bcPoints[bc], &idx);
247: for (i = 0; i < n; ++i) {
248: const PetscInt p = idx[i];
249: const PetscInt *find;
250: PetscInt fdof, fcdof, c;
252: PetscSectionGetFieldDof(section, p, field, &fdof);
253: if (!fdof) continue;
254: if (Nc < 0) {
255: for (d = 0; d < fdof; ++d) indices[d] = d;
256: fcdof = fdof;
257: } else {
258: PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
259: PetscSectionGetFieldConstraintIndices(section, p, field, &find);
260: for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
261: for (c = 0; c < Nc; ++c) indices[d++] = comp[c];
262: PetscSortRemoveDupsInt(&d, indices);
263: for (c = d; c < fcdof; ++c) indices[c] = -1;
264: fcdof = d;
265: }
266: PetscSectionSetFieldConstraintDof(section, p, field, fcdof);
267: PetscSectionSetFieldConstraintIndices(section, p, field, indices);
268: }
269: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
270: ISRestoreIndices(bcPoints[bc], &idx);
271: }
272: /* Handle anchors */
273: DMPlexGetAnchors(dm, &aSec, NULL);
274: if (aSec) {
275: PetscInt aStart, aEnd, a;
277: for (d = 0; d < maxDof; ++d) indices[d] = d;
278: PetscSectionGetChart(aSec, &aStart, &aEnd);
279: for (a = aStart; a < aEnd; a++) {
280: PetscInt dof, f;
282: PetscSectionGetDof(aSec, a, &dof);
283: if (dof) {
284: /* if there are point-to-point constraints, then all dofs are constrained */
285: for (f = 0; f < Nf; f++) {
286: PetscSectionSetFieldConstraintIndices(section, a, f, indices);
287: }
288: }
289: }
290: }
291: PetscFree(indices);
292: return(0);
293: }
295: /* Set the constrained indices on each point */
296: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
297: {
298: PetscInt *indices;
299: PetscInt Nf, maxDof, pStart, pEnd, p, f, d;
303: PetscSectionGetNumFields(section, &Nf);
304: PetscSectionGetMaxDof(section, &maxDof);
305: PetscSectionGetChart(section, &pStart, &pEnd);
306: PetscMalloc1(maxDof, &indices);
307: for (d = 0; d < maxDof; ++d) indices[d] = -1;
308: for (p = pStart; p < pEnd; ++p) {
309: PetscInt cdof, d;
311: PetscSectionGetConstraintDof(section, p, &cdof);
312: if (cdof) {
313: if (Nf) {
314: PetscInt numConst = 0, foff = 0;
316: for (f = 0; f < Nf; ++f) {
317: const PetscInt *find;
318: PetscInt fcdof, fdof;
320: PetscSectionGetFieldDof(section, p, f, &fdof);
321: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
322: /* Change constraint numbering from field component to local dof number */
323: PetscSectionGetFieldConstraintIndices(section, p, f, &find);
324: for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
325: numConst += fcdof;
326: foff += fdof;
327: }
328: if (cdof != numConst) {PetscSectionSetConstraintDof(section, p, numConst);}
329: } else {
330: for (d = 0; d < cdof; ++d) indices[d] = d;
331: }
332: PetscSectionSetConstraintIndices(section, p, indices);
333: }
334: }
335: PetscFree(indices);
336: return(0);
337: }
339: /*@C
340: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
342: Not Collective
344: Input Parameters:
345: + dm - The DMPlex object
346: . numComp - An array of size numFields that holds the number of components for each field
347: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
348: . label - The label indicating the mesh support of each field, or NULL for the whole mesh
349: . numBC - The number of boundary conditions
350: . bcField - An array of size numBC giving the field number for each boundry condition
351: . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
352: . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
353: - perm - Optional permutation of the chart, or NULL
355: Output Parameter:
356: . section - The PetscSection object
358: Notes:
359: numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the
360: number of dof for field 0 on each edge.
362: The chart permutation is the same one set using PetscSectionSetPermutation()
364: Level: developer
366: Fortran Notes:
367: A Fortran 90 version is available as DMPlexCreateSectionF90()
369: TODO: How is this related to DMCreateLocalSection()
371: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
372: @*/
373: 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)
374: {
375: PetscSection aSec;
379: DMPlexCreateSectionFields(dm, numComp, section);
380: DMPlexCreateSectionDof(dm, label, numDof, *section);
381: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
382: if (perm) {PetscSectionSetPermutation(*section, perm);}
383: PetscSectionSetFromOptions(*section);
384: PetscSectionSetUp(*section);
385: DMPlexGetAnchors(dm,&aSec,NULL);
386: if (numBC || aSec) {
387: DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
388: DMPlexCreateSectionBCIndices(dm, *section);
389: }
390: PetscSectionViewFromOptions(*section,NULL,"-section_view");
391: return(0);
392: }
394: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
395: {
396: PetscDS probBC;
397: PetscSection section;
398: DMLabel *labels;
399: IS *bcPoints, *bcComps;
400: PetscBool *isFE;
401: PetscInt *bcFields, *numComp, *numDof;
402: PetscInt depth, dim, numBd, numBC = 0, Nf, bd, bc = 0, f;
403: PetscInt cStart, cEnd, cEndInterior;
407: DMGetNumFields(dm, &Nf);
408: /* FE and FV boundary conditions are handled slightly differently */
409: PetscMalloc1(Nf, &isFE);
410: for (f = 0; f < Nf; ++f) {
411: PetscObject obj;
412: PetscClassId id;
414: DMGetField(dm, f, NULL, &obj);
415: PetscObjectGetClassId(obj, &id);
416: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
417: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
418: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
419: }
420: /* Allocate boundary point storage for FEM boundaries */
421: DMPlexGetDepth(dm, &depth);
422: DMGetDimension(dm, &dim);
423: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
424: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
425: DMGetDS(dm, &probBC);
426: PetscDSGetNumBoundary(probBC, &numBd);
427: 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)");
428: for (bd = 0; bd < numBd; ++bd) {
429: PetscInt field;
430: DMBoundaryConditionType type;
431: const char *labelName;
432: DMLabel label;
434: PetscDSGetBoundary(probBC, bd, &type, NULL, &labelName, &field, NULL, NULL, NULL, NULL, NULL, NULL);
435: DMGetLabel(dm,labelName,&label);
436: if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
437: }
438: /* Add ghost cell boundaries for FVM */
439: for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
440: PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
441: /* Constrain ghost cells for FV */
442: for (f = 0; f < Nf; ++f) {
443: PetscInt *newidx, c;
445: if (isFE[f] || cEndInterior < 0) continue;
446: PetscMalloc1(cEnd-cEndInterior,&newidx);
447: for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
448: bcFields[bc] = f;
449: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
450: }
451: /* Handle FEM Dirichlet boundaries */
452: for (bd = 0; bd < numBd; ++bd) {
453: const char *bdLabel;
454: DMLabel label;
455: const PetscInt *comps;
456: const PetscInt *values;
457: PetscInt bd2, field, numComps, numValues;
458: DMBoundaryConditionType type;
459: PetscBool duplicate = PETSC_FALSE;
461: PetscDSGetBoundary(probBC, bd, &type, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
462: DMGetLabel(dm, bdLabel, &label);
463: if (!isFE[field] || !label) continue;
464: /* Only want to modify label once */
465: for (bd2 = 0; bd2 < bd; ++bd2) {
466: const char *bdname;
467: PetscDSGetBoundary(probBC, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
468: PetscStrcmp(bdname, bdLabel, &duplicate);
469: if (duplicate) break;
470: }
471: if (!duplicate && (isFE[field])) {
472: /* don't complete cells, which are just present to give orientation to the boundary */
473: DMPlexLabelComplete(dm, label);
474: }
475: /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
476: if (type & DM_BC_ESSENTIAL) {
477: PetscInt *newidx;
478: PetscInt n, newn = 0, p, v;
480: bcFields[bc] = field;
481: if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
482: for (v = 0; v < numValues; ++v) {
483: IS tmp;
484: const PetscInt *idx;
486: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
487: if (!tmp) continue;
488: ISGetLocalSize(tmp, &n);
489: ISGetIndices(tmp, &idx);
490: if (isFE[field]) {
491: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
492: } else {
493: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
494: }
495: ISRestoreIndices(tmp, &idx);
496: ISDestroy(&tmp);
497: }
498: PetscMalloc1(newn,&newidx);
499: newn = 0;
500: for (v = 0; v < numValues; ++v) {
501: IS tmp;
502: const PetscInt *idx;
504: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
505: if (!tmp) continue;
506: ISGetLocalSize(tmp, &n);
507: ISGetIndices(tmp, &idx);
508: if (isFE[field]) {
509: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
510: } else {
511: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
512: }
513: ISRestoreIndices(tmp, &idx);
514: ISDestroy(&tmp);
515: }
516: ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
517: }
518: }
519: /* Handle discretization */
520: PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);
521: for (f = 0; f < Nf; ++f) {
522: labels[f] = dm->fields[f].label;
523: if (isFE[f]) {
524: PetscFE fe = (PetscFE) dm->fields[f].disc;
525: const PetscInt *numFieldDof;
526: PetscInt fedim, d;
528: PetscFEGetNumComponents(fe, &numComp[f]);
529: PetscFEGetNumDof(fe, &numFieldDof);
530: PetscFEGetSpatialDimension(fe, &fedim);
531: for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
532: } else {
533: PetscFV fv = (PetscFV) dm->fields[f].disc;
535: PetscFVGetNumComponents(fv, &numComp[f]);
536: numDof[f*(dim+1)+dim] = numComp[f];
537: }
538: }
539: for (f = 0; f < Nf; ++f) {
540: PetscInt d;
541: for (d = 1; d < dim; ++d) {
542: 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.");
543: }
544: }
545: DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);
546: for (f = 0; f < Nf; ++f) {
547: PetscFE fe;
548: const char *name;
550: if (isFE[f]) {
551: DMGetField(dm, f, NULL, (PetscObject *) &fe);
552: PetscObjectGetName((PetscObject) fe, &name);
553: PetscSectionSetFieldName(section, f, name);
554: }
555: }
556: DMSetLocalSection(dm, section);
557: PetscSectionDestroy(§ion);
558: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
559: PetscFree3(bcFields,bcPoints,bcComps);
560: PetscFree3(labels,numComp,numDof);
561: PetscFree(isFE);
562: return(0);
563: }