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;
10: DMPlexGetDepth(dm, &depth);
11: DMPlexGetDepthLabel(dm,&depthLabel);
12: DMGetNumFields(dm, &Nf);
13: PetscCalloc1(Nf, &isFE);
14: for (f = 0; f < Nf; ++f) {
15: PetscObject obj;
16: PetscClassId id;
18: DMGetField(dm, f, NULL, &obj);
19: PetscObjectGetClassId(obj, &id);
20: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
21: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
22: }
24: PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
25: if (Nf > 0) {
26: PetscSectionSetNumFields(*section, Nf);
27: if (numComp) {
28: for (f = 0; f < Nf; ++f) {
29: PetscSectionSetFieldComponents(*section, f, numComp[f]);
30: if (isFE[f]) {
31: PetscFE fe;
32: PetscDualSpace dspace;
33: const PetscInt ***perms;
34: const PetscScalar ***flips;
35: const PetscInt *numDof;
37: DMGetField(dm,f,NULL,(PetscObject *) &fe);
38: PetscFEGetDualSpace(fe,&dspace);
39: PetscDualSpaceGetSymmetries(dspace,&perms,&flips);
40: PetscDualSpaceGetNumDof(dspace,&numDof);
41: if (perms || flips) {
42: DM K;
43: PetscInt sph, spdepth;
44: PetscSectionSym sym;
46: PetscDualSpaceGetDM(dspace,&K);
47: DMPlexGetDepth(K, &spdepth);
48: PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)*section),depthLabel,&sym);
49: for (sph = 0; sph <= spdepth; sph++) {
50: PetscDualSpace hspace;
51: PetscInt kStart, kEnd;
52: PetscInt kConeSize, h = sph + (depth - spdepth);
53: const PetscInt **perms0 = NULL;
54: const PetscScalar **flips0 = NULL;
56: PetscDualSpaceGetHeightSubspace(dspace, sph, &hspace);
57: DMPlexGetHeightStratum(K, h, &kStart, &kEnd);
58: if (!hspace) continue;
59: PetscDualSpaceGetSymmetries(hspace,&perms,&flips);
60: if (perms) perms0 = perms[0];
61: if (flips) flips0 = flips[0];
62: if (!(perms0 || flips0)) continue;
63: {
64: DMPolytopeType ct;
65: /* The number of arrangements is no longer based on the number of faces */
66: DMPlexGetCellType(K, kStart, &ct);
67: kConeSize = DMPolytopeTypeGetNumArrangments(ct) / 2;
68: }
69: PetscSectionSymLabelSetStratum(sym,depth - h,numDof[depth - h],-kConeSize,kConeSize,PETSC_USE_POINTER,perms0 ? &perms0[-kConeSize] : NULL,flips0 ? &flips0[-kConeSize] : NULL);
70: }
71: PetscSectionSetFieldSym(*section,f,sym);
72: PetscSectionSymDestroy(&sym);
73: }
74: }
75: }
76: }
77: }
78: DMPlexGetChart(dm, &pStart, &pEnd);
79: PetscSectionSetChart(*section, pStart, pEnd);
80: PetscFree(isFE);
81: return 0;
82: }
84: /* Set the number of dof on each point and separate by fields */
85: static PetscErrorCode DMPlexCreateSectionDof(DM dm, DMLabel label[],const PetscInt numDof[], PetscSection section)
86: {
87: DMLabel depthLabel;
88: DMPolytopeType ct;
89: PetscInt depth, cellHeight, pStart = 0, pEnd = 0;
90: PetscInt Nf, f, Nds, n, dim, d, dep, p;
91: PetscBool *isFE, hasCohesive = PETSC_FALSE;
93: DMGetDimension(dm, &dim);
94: DMPlexGetDepth(dm, &depth);
95: DMPlexGetDepthLabel(dm,&depthLabel);
96: DMGetNumFields(dm, &Nf);
97: DMGetNumDS(dm, &Nds);
98: for (n = 0; n < Nds; ++n) {
99: PetscDS ds;
100: PetscBool isCohesive;
102: DMGetRegionNumDS(dm, n, NULL, NULL, &ds);
103: PetscDSIsCohesive(ds, &isCohesive);
104: if (isCohesive) {hasCohesive = PETSC_TRUE; break;}
105: }
106: PetscMalloc1(Nf, &isFE);
107: for (f = 0; f < Nf; ++f) {
108: PetscObject obj;
109: PetscClassId id;
111: DMGetField(dm, f, NULL, &obj);
112: PetscObjectGetClassId(obj, &id);
113: /* User is allowed to put a "placeholder" field in (c.f. DMCreateDS) */
114: isFE[f] = id == PETSCFE_CLASSID ? PETSC_TRUE : PETSC_FALSE;
115: }
117: DMPlexGetVTKCellHeight(dm, &cellHeight);
118: for (f = 0; f < Nf; ++f) {
119: PetscBool avoidTensor;
121: DMGetFieldAvoidTensor(dm, f, &avoidTensor);
122: avoidTensor = (avoidTensor || hasCohesive) ? PETSC_TRUE : PETSC_FALSE;
123: if (label && label[f]) {
124: IS pointIS;
125: const PetscInt *points;
126: PetscInt n;
128: DMLabelGetStratumIS(label[f], 1, &pointIS);
129: if (!pointIS) continue;
130: ISGetLocalSize(pointIS, &n);
131: ISGetIndices(pointIS, &points);
132: for (p = 0; p < n; ++p) {
133: const PetscInt point = points[p];
134: PetscInt dof, d;
136: DMPlexGetCellType(dm, point, &ct);
137: DMLabelGetValue(depthLabel, point, &d);
138: /* If this is a tensor prism point, use dof for one dimension lower */
139: switch (ct) {
140: case DM_POLYTOPE_POINT_PRISM_TENSOR:
141: case DM_POLYTOPE_SEG_PRISM_TENSOR:
142: case DM_POLYTOPE_TRI_PRISM_TENSOR:
143: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
144: if (hasCohesive) {--d;} break;
145: default: break;
146: }
147: dof = d < 0 ? 0 : numDof[f*(dim+1)+d];
148: PetscSectionSetFieldDof(section, point, f, dof);
149: PetscSectionAddDof(section, point, dof);
150: }
151: ISRestoreIndices(pointIS, &points);
152: ISDestroy(&pointIS);
153: } else {
154: for (dep = 0; dep <= depth - cellHeight; ++dep) {
155: /* Cases: dim > depth (cell-vertex mesh), dim == depth (fully interpolated), dim < depth (interpolated submesh) */
156: d = dim <= depth ? dep : (!dep ? 0 : dim);
157: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
158: for (p = pStart; p < pEnd; ++p) {
159: const PetscInt dof = numDof[f*(dim+1)+d];
161: DMPlexGetCellType(dm, p, &ct);
162: switch (ct) {
163: case DM_POLYTOPE_POINT_PRISM_TENSOR:
164: case DM_POLYTOPE_SEG_PRISM_TENSOR:
165: case DM_POLYTOPE_TRI_PRISM_TENSOR:
166: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
167: if (avoidTensor && isFE[f]) continue;
168: default: break;
169: }
170: PetscSectionSetFieldDof(section, p, f, dof);
171: PetscSectionAddDof(section, p, dof);
172: }
173: }
174: }
175: }
176: PetscFree(isFE);
177: return 0;
178: }
180: /* Set the number of dof on each point and separate by fields
181: If bcComps is NULL or the IS is NULL, constrain every dof on the point
182: */
183: static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
184: {
185: PetscInt Nf;
186: PetscInt bc;
187: PetscSection aSec;
189: PetscSectionGetNumFields(section, &Nf);
190: for (bc = 0; bc < numBC; ++bc) {
191: PetscInt field = 0;
192: const PetscInt *comp;
193: const PetscInt *idx;
194: PetscInt Nc = 0, cNc = -1, n, i;
196: if (Nf) {
197: field = bcField[bc];
198: PetscSectionGetFieldComponents(section, field, &Nc);
199: }
200: if (bcComps && bcComps[bc]) ISGetLocalSize(bcComps[bc], &cNc);
201: if (bcComps && bcComps[bc]) ISGetIndices(bcComps[bc], &comp);
202: ISGetLocalSize(bcPoints[bc], &n);
203: ISGetIndices(bcPoints[bc], &idx);
204: for (i = 0; i < n; ++i) {
205: const PetscInt p = idx[i];
206: PetscInt numConst;
208: if (Nf) {
209: PetscSectionGetFieldDof(section, p, field, &numConst);
210: } else {
211: PetscSectionGetDof(section, p, &numConst);
212: }
213: /* If Nc <= 0, constrain every dof on the point */
214: if (cNc > 0) {
215: /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
216: and that those dofs are numbered n*Nc+c */
217: if (Nf) {
219: numConst = (numConst/Nc) * cNc;
220: } else {
221: numConst = PetscMin(numConst, cNc);
222: }
223: }
224: if (Nf) PetscSectionAddFieldConstraintDof(section, p, field, numConst);
225: PetscSectionAddConstraintDof(section, p, numConst);
226: }
227: ISRestoreIndices(bcPoints[bc], &idx);
228: if (bcComps && bcComps[bc]) ISRestoreIndices(bcComps[bc], &comp);
229: }
230: DMPlexGetAnchors(dm, &aSec, NULL);
231: if (aSec) {
232: PetscInt aStart, aEnd, a;
234: PetscSectionGetChart(aSec, &aStart, &aEnd);
235: for (a = aStart; a < aEnd; a++) {
236: PetscInt dof, f;
238: PetscSectionGetDof(aSec, a, &dof);
239: if (dof) {
240: /* if there are point-to-point constraints, then all dofs are constrained */
241: PetscSectionGetDof(section, a, &dof);
242: PetscSectionSetConstraintDof(section, a, dof);
243: for (f = 0; f < Nf; f++) {
244: PetscSectionGetFieldDof(section, a, f, &dof);
245: PetscSectionSetFieldConstraintDof(section, a, f, dof);
246: }
247: }
248: }
249: }
250: return 0;
251: }
253: /* Set the constrained field indices on each point
254: If bcComps is NULL or the IS is NULL, constrain every dof on the point
255: */
256: static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
257: {
258: PetscSection aSec;
259: PetscInt *indices;
260: PetscInt Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;
262: PetscSectionGetNumFields(section, &Nf);
263: if (!Nf) return 0;
264: /* Initialize all field indices to -1 */
265: PetscSectionGetChart(section, &pStart, &pEnd);
266: for (p = pStart; p < pEnd; ++p) {PetscSectionGetConstraintDof(section, p, &cdof)); maxDof = PetscMax(maxDof, cdof;}
267: PetscMalloc1(maxDof, &indices);
268: for (d = 0; d < maxDof; ++d) indices[d] = -1;
269: for (p = pStart; p < pEnd; ++p) for (f = 0; f < Nf; ++f) PetscSectionSetFieldConstraintIndices(section, p, f, indices);
270: /* Handle BC constraints */
271: for (bc = 0; bc < numBC; ++bc) {
272: const PetscInt field = bcField[bc];
273: const PetscInt *comp, *idx;
274: PetscInt Nc, cNc = -1, n, i;
276: PetscSectionGetFieldComponents(section, field, &Nc);
277: if (bcComps && bcComps[bc]) ISGetLocalSize(bcComps[bc], &cNc);
278: if (bcComps && bcComps[bc]) ISGetIndices(bcComps[bc], &comp);
279: ISGetLocalSize(bcPoints[bc], &n);
280: ISGetIndices(bcPoints[bc], &idx);
281: for (i = 0; i < n; ++i) {
282: const PetscInt p = idx[i];
283: const PetscInt *find;
284: PetscInt fdof, fcdof, c, j;
286: PetscSectionGetFieldDof(section, p, field, &fdof);
287: if (!fdof) continue;
288: if (cNc < 0) {
289: for (d = 0; d < fdof; ++d) indices[d] = d;
290: fcdof = fdof;
291: } else {
292: /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
293: and that those dofs are numbered n*Nc+c */
294: PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
295: PetscSectionGetFieldConstraintIndices(section, p, field, &find);
296: /* Get indices constrained by previous bcs */
297: for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
298: for (j = 0; j < fdof/Nc; ++j) for (c = 0; c < cNc; ++c) indices[d++] = j*Nc + comp[c];
299: PetscSortRemoveDupsInt(&d, indices);
300: for (c = d; c < fcdof; ++c) indices[c] = -1;
301: fcdof = d;
302: }
303: PetscSectionSetFieldConstraintDof(section, p, field, fcdof);
304: PetscSectionSetFieldConstraintIndices(section, p, field, indices);
305: }
306: if (bcComps && bcComps[bc]) ISRestoreIndices(bcComps[bc], &comp);
307: ISRestoreIndices(bcPoints[bc], &idx);
308: }
309: /* Handle anchors */
310: DMPlexGetAnchors(dm, &aSec, NULL);
311: if (aSec) {
312: PetscInt aStart, aEnd, a;
314: for (d = 0; d < maxDof; ++d) indices[d] = d;
315: PetscSectionGetChart(aSec, &aStart, &aEnd);
316: for (a = aStart; a < aEnd; a++) {
317: PetscInt dof, f;
319: PetscSectionGetDof(aSec, a, &dof);
320: if (dof) {
321: /* if there are point-to-point constraints, then all dofs are constrained */
322: for (f = 0; f < Nf; f++) {
323: PetscSectionSetFieldConstraintIndices(section, a, f, indices);
324: }
325: }
326: }
327: }
328: PetscFree(indices);
329: return 0;
330: }
332: /* Set the constrained indices on each point */
333: static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
334: {
335: PetscInt *indices;
336: PetscInt Nf, maxDof, pStart, pEnd, p, f, d;
338: PetscSectionGetNumFields(section, &Nf);
339: PetscSectionGetMaxDof(section, &maxDof);
340: PetscSectionGetChart(section, &pStart, &pEnd);
341: PetscMalloc1(maxDof, &indices);
342: for (d = 0; d < maxDof; ++d) indices[d] = -1;
343: for (p = pStart; p < pEnd; ++p) {
344: PetscInt cdof, d;
346: PetscSectionGetConstraintDof(section, p, &cdof);
347: if (cdof) {
348: if (Nf) {
349: PetscInt numConst = 0, foff = 0;
351: for (f = 0; f < Nf; ++f) {
352: const PetscInt *find;
353: PetscInt fcdof, fdof;
355: PetscSectionGetFieldDof(section, p, f, &fdof);
356: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
357: /* Change constraint numbering from field component to local dof number */
358: PetscSectionGetFieldConstraintIndices(section, p, f, &find);
359: for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
360: numConst += fcdof;
361: foff += fdof;
362: }
363: if (cdof != numConst) PetscSectionSetConstraintDof(section, p, numConst);
364: } else {
365: for (d = 0; d < cdof; ++d) indices[d] = d;
366: }
367: PetscSectionSetConstraintIndices(section, p, indices);
368: }
369: }
370: PetscFree(indices);
371: return 0;
372: }
374: /*@C
375: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
377: Not Collective
379: Input Parameters:
380: + dm - The DMPlex object
381: . label - The label indicating the mesh support of each field, or NULL for the whole mesh
382: . numComp - An array of size numFields that holds the number of components for each field
383: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
384: . numBC - The number of boundary conditions
385: . bcField - An array of size numBC giving the field number for each boundry condition
386: . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
387: . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
388: - perm - Optional permutation of the chart, or NULL
390: Output Parameter:
391: . section - The PetscSection object
393: Notes:
394: numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the
395: number of dof for field 0 on each edge.
397: The chart permutation is the same one set using PetscSectionSetPermutation()
399: Level: developer
401: TODO: How is this related to DMCreateLocalSection()
403: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
404: @*/
405: 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)
406: {
407: PetscSection aSec;
409: DMPlexCreateSectionFields(dm, numComp, section);
410: DMPlexCreateSectionDof(dm, label, numDof, *section);
411: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
412: if (perm) PetscSectionSetPermutation(*section, perm);
413: PetscSectionSetFromOptions(*section);
414: PetscSectionSetUp(*section);
415: DMPlexGetAnchors(dm,&aSec,NULL);
416: if (numBC || aSec) {
417: DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
418: DMPlexCreateSectionBCIndices(dm, *section);
419: }
420: PetscSectionViewFromOptions(*section,NULL,"-section_view");
421: return 0;
422: }
424: PetscErrorCode DMCreateLocalSection_Plex(DM dm)
425: {
426: PetscSection section;
427: DMLabel *labels;
428: IS *bcPoints, *bcComps;
429: PetscBool *isFE;
430: PetscInt *bcFields, *numComp, *numDof;
431: PetscInt depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f;
432: PetscInt cStart, cEnd, cEndInterior;
434: DMGetNumFields(dm, &Nf);
435: DMGetDimension(dm, &dim);
436: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
437: /* FE and FV boundary conditions are handled slightly differently */
438: PetscMalloc1(Nf, &isFE);
439: for (f = 0; f < Nf; ++f) {
440: PetscObject obj;
441: PetscClassId id;
443: DMGetField(dm, f, NULL, &obj);
444: PetscObjectGetClassId(obj, &id);
445: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
446: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
447: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
448: }
449: /* Allocate boundary point storage for FEM boundaries */
450: DMGetNumDS(dm, &Nds);
451: for (s = 0; s < Nds; ++s) {
452: PetscDS dsBC;
453: PetscInt numBd, bd;
455: DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
456: PetscDSGetNumBoundary(dsBC, &numBd);
458: for (bd = 0; bd < numBd; ++bd) {
459: PetscInt field;
460: DMBoundaryConditionType type;
461: DMLabel label;
463: PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL);
464: if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
465: }
466: }
467: /* Add ghost cell boundaries for FVM */
468: DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
469: for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
470: PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps);
471: /* Constrain ghost cells for FV */
472: for (f = 0; f < Nf; ++f) {
473: PetscInt *newidx, c;
475: if (isFE[f] || cEndInterior < 0) continue;
476: PetscMalloc1(cEnd-cEndInterior,&newidx);
477: for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
478: bcFields[bc] = f;
479: ISCreateGeneral(PETSC_COMM_SELF, cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
480: }
481: /* Handle FEM Dirichlet boundaries */
482: DMGetNumDS(dm, &Nds);
483: for (s = 0; s < Nds; ++s) {
484: PetscDS dsBC;
485: PetscInt numBd, bd;
487: DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC);
488: PetscDSGetNumBoundary(dsBC, &numBd);
489: for (bd = 0; bd < numBd; ++bd) {
490: DMLabel label;
491: const PetscInt *comps;
492: const PetscInt *values;
493: PetscInt bd2, field, numComps, numValues;
494: DMBoundaryConditionType type;
495: PetscBool duplicate = PETSC_FALSE;
497: PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, &numValues, &values, &field, &numComps, &comps, NULL, NULL, NULL);
498: if (!isFE[field] || !label) continue;
499: /* Only want to modify label once */
500: for (bd2 = 0; bd2 < bd; ++bd2) {
501: DMLabel l;
503: PetscDSGetBoundary(dsBC, bd2, NULL, NULL, NULL, &l, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
504: duplicate = l == label ? PETSC_TRUE : PETSC_FALSE;
505: if (duplicate) break;
506: }
507: if (!duplicate && (isFE[field])) {
508: /* don't complete cells, which are just present to give orientation to the boundary */
509: DMPlexLabelComplete(dm, label);
510: }
511: /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
512: if (type & DM_BC_ESSENTIAL) {
513: PetscInt *newidx;
514: PetscInt n, newn = 0, p, v;
516: bcFields[bc] = field;
517: if (numComps) ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);
518: for (v = 0; v < numValues; ++v) {
519: IS tmp;
520: const PetscInt *idx;
522: DMLabelGetStratumIS(label, values[v], &tmp);
523: if (!tmp) continue;
524: ISGetLocalSize(tmp, &n);
525: ISGetIndices(tmp, &idx);
526: if (isFE[field]) {
527: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
528: } else {
529: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
530: }
531: ISRestoreIndices(tmp, &idx);
532: ISDestroy(&tmp);
533: }
534: PetscMalloc1(newn, &newidx);
535: newn = 0;
536: for (v = 0; v < numValues; ++v) {
537: IS tmp;
538: const PetscInt *idx;
540: DMLabelGetStratumIS(label, values[v], &tmp);
541: if (!tmp) continue;
542: ISGetLocalSize(tmp, &n);
543: ISGetIndices(tmp, &idx);
544: if (isFE[field]) {
545: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
546: } else {
547: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
548: }
549: ISRestoreIndices(tmp, &idx);
550: ISDestroy(&tmp);
551: }
552: ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
553: }
554: }
555: }
556: /* Handle discretization */
557: PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);
558: for (f = 0; f < Nf; ++f) {
559: labels[f] = dm->fields[f].label;
560: if (isFE[f]) {
561: PetscFE fe = (PetscFE) dm->fields[f].disc;
562: const PetscInt *numFieldDof;
563: PetscInt fedim, d;
565: PetscFEGetNumComponents(fe, &numComp[f]);
566: PetscFEGetNumDof(fe, &numFieldDof);
567: PetscFEGetSpatialDimension(fe, &fedim);
568: for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
569: } else {
570: PetscFV fv = (PetscFV) dm->fields[f].disc;
572: PetscFVGetNumComponents(fv, &numComp[f]);
573: numDof[f*(dim+1)+dim] = numComp[f];
574: }
575: }
576: DMPlexGetDepth(dm, &depth);
577: for (f = 0; f < Nf; ++f) {
578: PetscInt d;
579: for (d = 1; d < dim; ++d) {
581: }
582: }
583: DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);
584: for (f = 0; f < Nf; ++f) {
585: PetscFE fe;
586: const char *name;
588: DMGetField(dm, f, NULL, (PetscObject *) &fe);
589: if (!((PetscObject) fe)->name) continue;
590: PetscObjectGetName((PetscObject) fe, &name);
591: PetscSectionSetFieldName(section, f, name);
592: }
593: DMSetLocalSection(dm, section);
594: PetscSectionDestroy(§ion);
595: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]));PetscCall(ISDestroy(&bcComps[bc]);}
596: PetscFree3(bcFields,bcPoints,bcComps);
597: PetscFree3(labels,numComp,numDof);
598: PetscFree(isFE);
599: return 0;
600: }