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