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