Actual source code: plexsection.c
petsc-3.14.6 2021-03-30
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, dim, d, dep, p;
88: PetscBool *isFE;
92: DMGetDimension(dm, &dim);
93: DMPlexGetDepth(dm, &depth);
94: DMPlexGetDepthLabel(dm,&depthLabel);
95: DMGetNumFields(dm, &Nf);
96: PetscMalloc1(Nf, &isFE);
97: for (f = 0; f < Nf; ++f) {
98: PetscObject obj;
99: PetscClassId id;
101: DMGetField(dm, f, NULL, &obj);
102: PetscObjectGetClassId(obj, &id);
103: /* User is allowed to put a "placeholder" field in (c.f. DMCreateDS) */
104: isFE[f] = id == PETSCFE_CLASSID ? PETSC_TRUE : PETSC_FALSE;
105: }
107: DMPlexGetVTKCellHeight(dm, &cellHeight);
108: for (f = 0; f < Nf; ++f) {
109: if (label && label[f]) {
110: IS pointIS;
111: const PetscInt *points;
112: PetscInt n;
114: DMLabelGetStratumIS(label[f], 1, &pointIS);
115: if (!pointIS) continue;
116: ISGetLocalSize(pointIS, &n);
117: ISGetIndices(pointIS, &points);
118: for (p = 0; p < n; ++p) {
119: const PetscInt point = points[p];
120: PetscInt dof, d;
122: DMPlexGetCellType(dm, point, &ct);
123: DMLabelGetValue(depthLabel, point, &d);
124: /* If this is a tensor prism point, use dof for one dimension lower */
125: switch (ct) {
126: case DM_POLYTOPE_POINT_PRISM_TENSOR:
127: case DM_POLYTOPE_SEG_PRISM_TENSOR:
128: case DM_POLYTOPE_TRI_PRISM_TENSOR:
129: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
130: --d;break;
131: default: break;
132: }
133: dof = d < 0 ? 0 : numDof[f*(dim+1)+d];
134: PetscSectionSetFieldDof(section, point, f, dof);
135: PetscSectionAddDof(section, point, dof);
136: }
137: ISRestoreIndices(pointIS, &points);
138: ISDestroy(&pointIS);
139: } else {
140: for (dep = 0; dep <= depth - cellHeight; ++dep) {
141: /* Cases: dim > depth (cell-vertex mesh), dim == depth (fully interpolated), dim < depth (interpolated submesh) */
142: d = dim <= depth ? dep : (!dep ? 0 : dim);
143: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
144: for (p = pStart; p < pEnd; ++p) {
145: const PetscInt dof = numDof[f*(dim+1)+d];
147: DMPlexGetCellType(dm, p, &ct);
148: switch (ct) {
149: case DM_POLYTOPE_POINT_PRISM_TENSOR:
150: case DM_POLYTOPE_SEG_PRISM_TENSOR:
151: case DM_POLYTOPE_TRI_PRISM_TENSOR:
152: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
153: if (isFE[f]) continue;
154: default: break;
155: }
156: PetscSectionSetFieldDof(section, p, f, dof);
157: PetscSectionAddDof(section, p, dof);
158: }
159: }
160: }
161: }
162: PetscFree(isFE);
163: return(0);
164: }
166: /* Set the number of dof on each point and separate by fields
167: If bcComps is NULL or the IS is NULL, constrain every dof on the point
168: */
169: static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
170: {
171: PetscInt Nf;
172: PetscInt bc;
173: PetscSection aSec;
177: PetscSectionGetNumFields(section, &Nf);
178: for (bc = 0; bc < numBC; ++bc) {
179: PetscInt field = 0;
180: const PetscInt *comp;
181: const PetscInt *idx;
182: PetscInt Nc = 0, cNc = -1, n, i;
184: if (Nf) {
185: field = bcField[bc];
186: PetscSectionGetFieldComponents(section, field, &Nc);
187: }
188: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &cNc);}
189: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
190: ISGetLocalSize(bcPoints[bc], &n);
191: ISGetIndices(bcPoints[bc], &idx);
192: for (i = 0; i < n; ++i) {
193: const PetscInt p = idx[i];
194: PetscInt numConst;
196: if (Nf) {
197: PetscSectionGetFieldDof(section, p, field, &numConst);
198: } else {
199: PetscSectionGetDof(section, p, &numConst);
200: }
201: /* If Nc <= 0, constrain every dof on the point */
202: if (cNc > 0) {
203: /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
204: and that those dofs are numbered n*Nc+c */
205: if (Nf) {
206: 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);
207: numConst = (numConst/Nc) * cNc;
208: } else {
209: numConst = PetscMin(numConst, cNc);
210: }
211: }
212: if (Nf) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
213: PetscSectionAddConstraintDof(section, p, numConst);
214: }
215: ISRestoreIndices(bcPoints[bc], &idx);
216: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
217: }
218: DMPlexGetAnchors(dm, &aSec, NULL);
219: if (aSec) {
220: PetscInt aStart, aEnd, a;
222: PetscSectionGetChart(aSec, &aStart, &aEnd);
223: for (a = aStart; a < aEnd; a++) {
224: PetscInt dof, f;
226: PetscSectionGetDof(aSec, a, &dof);
227: if (dof) {
228: /* if there are point-to-point constraints, then all dofs are constrained */
229: PetscSectionGetDof(section, a, &dof);
230: PetscSectionSetConstraintDof(section, a, dof);
231: for (f = 0; f < Nf; f++) {
232: PetscSectionGetFieldDof(section, a, f, &dof);
233: PetscSectionSetFieldConstraintDof(section, a, f, dof);
234: }
235: }
236: }
237: }
238: return(0);
239: }
241: /* Set the constrained field indices on each point
242: If bcComps is NULL or the IS is NULL, constrain every dof on the point
243: */
244: static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
245: {
246: PetscSection aSec;
247: PetscInt *indices;
248: PetscInt Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;
252: PetscSectionGetNumFields(section, &Nf);
253: if (!Nf) return(0);
254: /* Initialize all field indices to -1 */
255: PetscSectionGetChart(section, &pStart, &pEnd);
256: for (p = pStart; p < pEnd; ++p) {PetscSectionGetConstraintDof(section, p, &cdof); maxDof = PetscMax(maxDof, cdof);}
257: PetscMalloc1(maxDof, &indices);
258: for (d = 0; d < maxDof; ++d) indices[d] = -1;
259: for (p = pStart; p < pEnd; ++p) for (f = 0; f < Nf; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
260: /* Handle BC constraints */
261: for (bc = 0; bc < numBC; ++bc) {
262: const PetscInt field = bcField[bc];
263: const PetscInt *comp, *idx;
264: PetscInt Nc, cNc = -1, n, i;
266: PetscSectionGetFieldComponents(section, field, &Nc);
267: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &cNc);}
268: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
269: ISGetLocalSize(bcPoints[bc], &n);
270: ISGetIndices(bcPoints[bc], &idx);
271: for (i = 0; i < n; ++i) {
272: const PetscInt p = idx[i];
273: const PetscInt *find;
274: PetscInt fdof, fcdof, c, j;
276: PetscSectionGetFieldDof(section, p, field, &fdof);
277: if (!fdof) continue;
278: if (cNc < 0) {
279: for (d = 0; d < fdof; ++d) indices[d] = d;
280: fcdof = fdof;
281: } else {
282: /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
283: and that those dofs are numbered n*Nc+c */
284: PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
285: PetscSectionGetFieldConstraintIndices(section, p, field, &find);
286: /* Get indices constrained by previous bcs */
287: for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
288: for (j = 0; j < fdof/Nc; ++j) for (c = 0; c < cNc; ++c) indices[d++] = j*Nc + comp[c];
289: PetscSortRemoveDupsInt(&d, indices);
290: for (c = d; c < fcdof; ++c) indices[c] = -1;
291: fcdof = d;
292: }
293: PetscSectionSetFieldConstraintDof(section, p, field, fcdof);
294: PetscSectionSetFieldConstraintIndices(section, p, field, indices);
295: }
296: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
297: ISRestoreIndices(bcPoints[bc], &idx);
298: }
299: /* Handle anchors */
300: DMPlexGetAnchors(dm, &aSec, NULL);
301: if (aSec) {
302: PetscInt aStart, aEnd, a;
304: for (d = 0; d < maxDof; ++d) indices[d] = d;
305: PetscSectionGetChart(aSec, &aStart, &aEnd);
306: for (a = aStart; a < aEnd; a++) {
307: PetscInt dof, f;
309: PetscSectionGetDof(aSec, a, &dof);
310: if (dof) {
311: /* if there are point-to-point constraints, then all dofs are constrained */
312: for (f = 0; f < Nf; f++) {
313: PetscSectionSetFieldConstraintIndices(section, a, f, indices);
314: }
315: }
316: }
317: }
318: PetscFree(indices);
319: return(0);
320: }
322: /* Set the constrained indices on each point */
323: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
324: {
325: PetscInt *indices;
326: PetscInt Nf, maxDof, pStart, pEnd, p, f, d;
330: PetscSectionGetNumFields(section, &Nf);
331: PetscSectionGetMaxDof(section, &maxDof);
332: PetscSectionGetChart(section, &pStart, &pEnd);
333: PetscMalloc1(maxDof, &indices);
334: for (d = 0; d < maxDof; ++d) indices[d] = -1;
335: for (p = pStart; p < pEnd; ++p) {
336: PetscInt cdof, d;
338: PetscSectionGetConstraintDof(section, p, &cdof);
339: if (cdof) {
340: if (Nf) {
341: PetscInt numConst = 0, foff = 0;
343: for (f = 0; f < Nf; ++f) {
344: const PetscInt *find;
345: PetscInt fcdof, fdof;
347: PetscSectionGetFieldDof(section, p, f, &fdof);
348: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
349: /* Change constraint numbering from field component to local dof number */
350: PetscSectionGetFieldConstraintIndices(section, p, f, &find);
351: for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
352: numConst += fcdof;
353: foff += fdof;
354: }
355: if (cdof != numConst) {PetscSectionSetConstraintDof(section, p, numConst);}
356: } else {
357: for (d = 0; d < cdof; ++d) indices[d] = d;
358: }
359: PetscSectionSetConstraintIndices(section, p, indices);
360: }
361: }
362: PetscFree(indices);
363: return(0);
364: }
366: /*@C
367: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
369: Not Collective
371: Input Parameters:
372: + dm - The DMPlex object
373: . label - The label indicating the mesh support of each field, or NULL for the whole mesh
374: . numComp - An array of size numFields that holds the number of components for each field
375: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
376: . numBC - The number of boundary conditions
377: . bcField - An array of size numBC giving the field number for each boundry condition
378: . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
379: . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
380: - perm - Optional permutation of the chart, or NULL
382: Output Parameter:
383: . section - The PetscSection object
385: Notes:
386: numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the
387: number of dof for field 0 on each edge.
389: The chart permutation is the same one set using PetscSectionSetPermutation()
391: Level: developer
393: TODO: How is this related to DMCreateLocalSection()
395: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
396: @*/
397: 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)
398: {
399: PetscSection aSec;
403: DMPlexCreateSectionFields(dm, numComp, section);
404: DMPlexCreateSectionDof(dm, label, numDof, *section);
405: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
406: if (perm) {PetscSectionSetPermutation(*section, perm);}
407: PetscSectionSetFromOptions(*section);
408: PetscSectionSetUp(*section);
409: DMPlexGetAnchors(dm,&aSec,NULL);
410: if (numBC || aSec) {
411: DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
412: DMPlexCreateSectionBCIndices(dm, *section);
413: }
414: PetscSectionViewFromOptions(*section,NULL,"-section_view");
415: return(0);
416: }
418: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
419: {
420: PetscSection section;
421: DMLabel *labels;
422: IS *bcPoints, *bcComps;
423: PetscBool *isFE;
424: PetscInt *bcFields, *numComp, *numDof;
425: PetscInt depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f;
426: PetscInt cStart, cEnd, cEndInterior;
430: DMGetNumFields(dm, &Nf);
431: DMGetDimension(dm, &dim);
432: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
433: /* FE and FV boundary conditions are handled slightly differently */
434: PetscMalloc1(Nf, &isFE);
435: for (f = 0; f < Nf; ++f) {
436: PetscObject obj;
437: PetscClassId id;
439: DMGetField(dm, f, NULL, &obj);
440: PetscObjectGetClassId(obj, &id);
441: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
442: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
443: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
444: }
445: /* Allocate boundary point storage for FEM boundaries */
446: DMGetNumDS(dm, &Nds);
447: for (s = 0; s < Nds; ++s) {
448: PetscDS dsBC;
449: PetscInt numBd, bd;
451: DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
452: PetscDSGetNumBoundary(dsBC, &numBd);
453: 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)");
454: for (bd = 0; bd < numBd; ++bd) {
455: PetscInt field;
456: DMBoundaryConditionType type;
457: const char *labelName;
458: DMLabel label;
460: PetscDSGetBoundary(dsBC, bd, &type, NULL, &labelName, &field, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
461: DMGetLabel(dm, labelName, &label);
462: if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
463: }
464: }
465: /* Add ghost cell boundaries for FVM */
466: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
467: for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
468: PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps);
469: /* Constrain ghost cells for FV */
470: for (f = 0; f < Nf; ++f) {
471: PetscInt *newidx, c;
473: if (isFE[f] || cEndInterior < 0) continue;
474: PetscMalloc1(cEnd-cEndInterior,&newidx);
475: for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
476: bcFields[bc] = f;
477: ISCreateGeneral(PETSC_COMM_SELF, cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
478: }
479: /* Handle FEM Dirichlet boundaries */
480: DMGetNumDS(dm, &Nds);
481: for (s = 0; s < Nds; ++s) {
482: PetscDS dsBC;
483: PetscInt numBd, bd;
485: DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
486: PetscDSGetNumBoundary(dsBC, &numBd);
487: for (bd = 0; bd < numBd; ++bd) {
488: const char *bdLabel;
489: DMLabel label;
490: const PetscInt *comps;
491: const PetscInt *values;
492: PetscInt bd2, field, numComps, numValues;
493: DMBoundaryConditionType type;
494: PetscBool duplicate = PETSC_FALSE;
496: PetscDSGetBoundary(dsBC, bd, &type, NULL, &bdLabel, &field, &numComps, &comps, NULL, NULL, &numValues, &values, NULL);
497: DMGetLabel(dm, bdLabel, &label);
498: if (!isFE[field] || !label) continue;
499: /* Only want to modify label once */
500: for (bd2 = 0; bd2 < bd; ++bd2) {
501: const char *bdname;
502: PetscDSGetBoundary(dsBC, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
503: PetscStrcmp(bdname, bdLabel, &duplicate);
504: if (duplicate) break;
505: }
506: if (!duplicate && (isFE[field])) {
507: /* don't complete cells, which are just present to give orientation to the boundary */
508: DMPlexLabelComplete(dm, label);
509: }
510: /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
511: if (type & DM_BC_ESSENTIAL) {
512: PetscInt *newidx;
513: PetscInt n, newn = 0, p, v;
515: bcFields[bc] = field;
516: if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
517: for (v = 0; v < numValues; ++v) {
518: IS tmp;
519: const PetscInt *idx;
521: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
522: if (!tmp) continue;
523: ISGetLocalSize(tmp, &n);
524: ISGetIndices(tmp, &idx);
525: if (isFE[field]) {
526: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
527: } else {
528: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
529: }
530: ISRestoreIndices(tmp, &idx);
531: ISDestroy(&tmp);
532: }
533: PetscMalloc1(newn, &newidx);
534: newn = 0;
535: for (v = 0; v < numValues; ++v) {
536: IS tmp;
537: const PetscInt *idx;
539: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
540: if (!tmp) continue;
541: ISGetLocalSize(tmp, &n);
542: ISGetIndices(tmp, &idx);
543: if (isFE[field]) {
544: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
545: } else {
546: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
547: }
548: ISRestoreIndices(tmp, &idx);
549: ISDestroy(&tmp);
550: }
551: ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
552: }
553: }
554: }
555: /* Handle discretization */
556: PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);
557: for (f = 0; f < Nf; ++f) {
558: labels[f] = dm->fields[f].label;
559: if (isFE[f]) {
560: PetscFE fe = (PetscFE) dm->fields[f].disc;
561: const PetscInt *numFieldDof;
562: PetscInt fedim, d;
564: PetscFEGetNumComponents(fe, &numComp[f]);
565: PetscFEGetNumDof(fe, &numFieldDof);
566: PetscFEGetSpatialDimension(fe, &fedim);
567: for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
568: } else {
569: PetscFV fv = (PetscFV) dm->fields[f].disc;
571: PetscFVGetNumComponents(fv, &numComp[f]);
572: numDof[f*(dim+1)+dim] = numComp[f];
573: }
574: }
575: DMPlexGetDepth(dm, &depth);
576: for (f = 0; f < Nf; ++f) {
577: PetscInt d;
578: for (d = 1; d < dim; ++d) {
579: 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.");
580: }
581: }
582: DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);
583: for (f = 0; f < Nf; ++f) {
584: PetscFE fe;
585: const char *name;
587: DMGetField(dm, f, NULL, (PetscObject *) &fe);
588: if (!((PetscObject) fe)->name) continue;
589: PetscObjectGetName((PetscObject) fe, &name);
590: PetscSectionSetFieldName(section, f, name);
591: }
592: DMSetLocalSection(dm, section);
593: PetscSectionDestroy(§ion);
594: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
595: PetscFree3(bcFields,bcPoints,bcComps);
596: PetscFree3(labels,numComp,numDof);
597: PetscFree(isFE);
598: return(0);
599: }