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