Actual source code: plexsection.c
petsc-3.13.6 2020-09-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: 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, §ion);
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(§ion);
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: }