1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3 /* Set the number of dof on each point and separate by fields */
DMPlexCreateSectionFields(DM dm,const PetscInt numComp[],PetscSection * section)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;
9
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;
18
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 }
27
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;
40
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;
49
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;
59
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 }
87
88 /* Set the number of dof on each point and separate by fields */
DMPlexCreateSectionDof(DM dm,DMLabel label[],const PetscInt numDof[],PetscSection section)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;
96
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;
106
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;
118
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 }
124
125 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
126 for (f = 0; f < Nf; ++f) {
127 PetscBool avoidTensor;
128
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;
135
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;
143
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];
170
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 }
190
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 */
DMPlexCreateSectionBCDof(DM dm,PetscInt numBC,const PetscInt bcField[],const IS bcComps[],const IS bcPoints[],PetscSection section)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;
199
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;
207
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;
220
221 if (Nf) {
222 PetscCall(PetscSectionGetFieldDof(section, p, field, &numConst));
223 } else {
224 PetscCall(PetscSectionGetDof(section, p, &numConst));
225 }
226 /* If Nc <= 0, constrain every dof on the point */
227 if (cNc > 0) {
228 /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
229 and that those dofs are numbered n*Nc+c */
230 if (Nf) {
231 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);
232 numConst = (numConst / Nc) * cNc;
233 } else {
234 numConst = PetscMin(numConst, cNc);
235 }
236 }
237 if (Nf) PetscCall(PetscSectionAddFieldConstraintDof(section, p, field, numConst));
238 PetscCall(PetscSectionAddConstraintDof(section, p, numConst));
239 }
240 PetscCall(ISRestoreIndices(bcPoints[bc], &idx));
241 }
242 if (bcComps && bcComps[bc]) PetscCall(ISRestoreIndices(bcComps[bc], &comp));
243 }
244 PetscCall(DMPlexGetAnchors(dm, &aSec, NULL));
245 if (aSec) {
246 PetscInt aStart, aEnd, a;
247
248 PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
249 for (a = aStart; a < aEnd; a++) {
250 PetscInt dof, f;
251
252 PetscCall(PetscSectionGetDof(aSec, a, &dof));
253 if (dof) {
254 /* if there are point-to-point constraints, then all dofs are constrained */
255 PetscCall(PetscSectionGetDof(section, a, &dof));
256 PetscCall(PetscSectionSetConstraintDof(section, a, dof));
257 for (f = 0; f < Nf; f++) {
258 PetscCall(PetscSectionGetFieldDof(section, a, f, &dof));
259 PetscCall(PetscSectionSetFieldConstraintDof(section, a, f, dof));
260 }
261 }
262 }
263 }
264 PetscFunctionReturn(PETSC_SUCCESS);
265 }
266
267 /* Set the constrained field indices on each point
268 If bcComps is NULL or the IS is NULL, constrain every dof on the point
269 */
DMPlexCreateSectionBCIndicesField(DM dm,PetscInt numBC,const PetscInt bcField[],const IS bcComps[],const IS bcPoints[],PetscSection section)270 static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
271 {
272 PetscSection aSec;
273 PetscInt *indices;
274 PetscInt Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d;
275
276 PetscFunctionBegin;
277 PetscCall(PetscSectionGetNumFields(section, &Nf));
278 if (!Nf) PetscFunctionReturn(PETSC_SUCCESS);
279 /* Initialize all field indices to -1 */
280 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
281 for (p = pStart; p < pEnd; ++p) {
282 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
283 maxDof = PetscMax(maxDof, cdof);
284 }
285 PetscCall(PetscMalloc1(maxDof, &indices));
286 for (d = 0; d < maxDof; ++d) indices[d] = -1;
287 for (p = pStart; p < pEnd; ++p)
288 for (f = 0; f < Nf; ++f) PetscCall(PetscSectionSetFieldConstraintIndices(section, p, f, indices));
289 /* Handle BC constraints */
290 for (bc = 0; bc < numBC; ++bc) {
291 const PetscInt field = bcField[bc];
292 const PetscInt *comp, *idx;
293 PetscInt Nc, cNc = -1, n, i;
294
295 PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
296 if (bcComps && bcComps[bc]) PetscCall(ISGetLocalSize(bcComps[bc], &cNc));
297 if (bcComps && bcComps[bc]) PetscCall(ISGetIndices(bcComps[bc], &comp));
298 if (bcPoints[bc]) {
299 PetscCall(ISGetLocalSize(bcPoints[bc], &n));
300 PetscCall(ISGetIndices(bcPoints[bc], &idx));
301 for (i = 0; i < n; ++i) {
302 const PetscInt p = idx[i];
303 const PetscInt *find;
304 PetscInt fdof, fcdof, c, j;
305
306 PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
307 if (!fdof) continue;
308 if (cNc < 0) {
309 for (d = 0; d < fdof; ++d) indices[d] = d;
310 fcdof = fdof;
311 } else {
312 /* We assume that a point may have multiple "nodes", which are collections of Nc dofs,
313 and that those dofs are numbered n*Nc+c */
314 PetscCall(PetscSectionGetFieldConstraintDof(section, p, field, &fcdof));
315 PetscCall(PetscSectionGetFieldConstraintIndices(section, p, field, &find));
316 /* Get indices constrained by previous bcs */
317 for (d = 0; d < fcdof; ++d) {
318 if (find[d] < 0) break;
319 indices[d] = find[d];
320 }
321 for (j = 0; j < fdof / Nc; ++j)
322 for (c = 0; c < cNc; ++c) indices[d++] = j * Nc + comp[c];
323 PetscCall(PetscSortRemoveDupsInt(&d, indices));
324 for (c = d; c < fcdof; ++c) indices[c] = -1;
325 fcdof = d;
326 }
327 PetscCall(PetscSectionSetFieldConstraintDof(section, p, field, fcdof));
328 PetscCall(PetscSectionSetFieldConstraintIndices(section, p, field, indices));
329 }
330 PetscCall(ISRestoreIndices(bcPoints[bc], &idx));
331 }
332 if (bcComps && bcComps[bc]) PetscCall(ISRestoreIndices(bcComps[bc], &comp));
333 }
334 /* Handle anchors */
335 PetscCall(DMPlexGetAnchors(dm, &aSec, NULL));
336 if (aSec) {
337 PetscInt aStart, aEnd, a;
338
339 for (d = 0; d < maxDof; ++d) indices[d] = d;
340 PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
341 for (a = aStart; a < aEnd; a++) {
342 PetscInt dof, f;
343
344 PetscCall(PetscSectionGetDof(aSec, a, &dof));
345 if (dof) {
346 /* if there are point-to-point constraints, then all dofs are constrained */
347 for (f = 0; f < Nf; f++) PetscCall(PetscSectionSetFieldConstraintIndices(section, a, f, indices));
348 }
349 }
350 }
351 PetscCall(PetscFree(indices));
352 PetscFunctionReturn(PETSC_SUCCESS);
353 }
354
355 /* Set the constrained indices on each point */
DMPlexCreateSectionBCIndices(DM dm,PetscSection section)356 static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
357 {
358 PetscInt *indices;
359 PetscInt Nf, maxDof, pStart, pEnd, p, f, d;
360
361 PetscFunctionBegin;
362 PetscCall(PetscSectionGetNumFields(section, &Nf));
363 PetscCall(PetscSectionGetMaxDof(section, &maxDof));
364 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
365 PetscCall(PetscMalloc1(maxDof, &indices));
366 for (d = 0; d < maxDof; ++d) indices[d] = -1;
367 for (p = pStart; p < pEnd; ++p) {
368 PetscInt cdof, d;
369
370 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
371 if (cdof) {
372 if (Nf) {
373 PetscInt numConst = 0, foff = 0;
374
375 for (f = 0; f < Nf; ++f) {
376 const PetscInt *find;
377 PetscInt fcdof, fdof;
378
379 PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
380 PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
381 /* Change constraint numbering from field component to local dof number */
382 PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &find));
383 for (d = 0; d < fcdof; ++d) indices[numConst + d] = find[d] + foff;
384 numConst += fcdof;
385 foff += fdof;
386 }
387 if (cdof != numConst) PetscCall(PetscSectionSetConstraintDof(section, p, numConst));
388 } else {
389 for (d = 0; d < cdof; ++d) indices[d] = d;
390 }
391 PetscCall(PetscSectionSetConstraintIndices(section, p, indices));
392 }
393 }
394 PetscCall(PetscFree(indices));
395 PetscFunctionReturn(PETSC_SUCCESS);
396 }
397
398 /*@C
399 DMPlexCreateSection - Create a `PetscSection` based upon the dof layout specification provided.
400
401 Not Collective
402
403 Input Parameters:
404 + dm - The `DMPLEX` object
405 . label - An array of `DMLabel` of length `numFields` indicating the mesh support of each field, or `NULL` for the whole mesh
406 . numComp - An array of size `numFields` that holds the number of components for each field
407 . 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
408 . numBC - The number of boundary conditions
409 . bcField - An array of size `numBC` giving the field number for each boundary condition
410 . bcComps - [Optional] An array of size `numBC` of `IS` holding the field components to which each boundary condition applies
411 . bcPoints - An array of size `numBC` of `IS` holding the `DMPLEX` points to which each boundary condition applies
412 - perm - Optional permutation of the chart, or `NULL`
413
414 Output Parameter:
415 . section - The `PetscSection` object
416
417 Level: developer
418
419 Notes:
420 numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the
421 number of dof for field 0 on each edge.
422
423 The chart permutation is the same one set using `PetscSectionSetPermutation()`
424
425 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
426 @*/
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)427 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)
428 {
429 PetscSection aSec;
430
431 PetscFunctionBegin;
432 PetscCall(DMPlexCreateSectionFields(dm, numComp, section));
433 PetscCall(DMPlexCreateSectionDof(dm, label, numDof, *section));
434 PetscCall(DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section));
435 if (perm) PetscCall(PetscSectionSetPermutation(*section, perm));
436 PetscCall(PetscSectionSetFromOptions(*section));
437 PetscCall(PetscSectionSetUp(*section));
438 PetscCall(DMPlexGetAnchors(dm, &aSec, NULL));
439 if (numBC || aSec) {
440 PetscCall(DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section));
441 PetscCall(DMPlexCreateSectionBCIndices(dm, *section));
442 }
443 PetscCall(PetscSectionViewFromOptions(*section, NULL, "-section_view"));
444 PetscFunctionReturn(PETSC_SUCCESS);
445 }
446
DMCreateLocalSection_Plex(DM dm)447 PetscErrorCode DMCreateLocalSection_Plex(DM dm)
448 {
449 PetscSection section;
450 DMLabel *labels;
451 IS *bcPoints, *bcComps, permIS;
452 PetscBT blockStarts;
453 PetscBool *isFE;
454 PetscInt *bcFields, *numComp, *numDof;
455 PetscInt depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f;
456 PetscInt cStart, cEnd, cEndInterior;
457
458 PetscFunctionBegin;
459 PetscCall(DMGetNumFields(dm, &Nf));
460 PetscCall(DMGetDimension(dm, &dim));
461 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
462 /* FE and FV boundary conditions are handled slightly differently */
463 PetscCall(PetscMalloc1(Nf, &isFE));
464 for (f = 0; f < Nf; ++f) {
465 PetscObject obj;
466 PetscClassId id;
467
468 PetscCall(DMGetField(dm, f, NULL, &obj));
469 PetscCall(PetscObjectGetClassId(obj, &id));
470 if (id == PETSCFE_CLASSID) {
471 isFE[f] = PETSC_TRUE;
472 } else if (id == PETSCFV_CLASSID) {
473 isFE[f] = PETSC_FALSE;
474 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
475 }
476 /* Allocate boundary point storage for FEM boundaries */
477 PetscCall(DMGetNumDS(dm, &Nds));
478 for (s = 0; s < Nds; ++s) {
479 PetscDS dsBC;
480 PetscInt numBd, bd;
481
482 PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC, NULL));
483 PetscCall(PetscDSGetNumBoundary(dsBC, &numBd));
484 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)");
485 for (bd = 0; bd < numBd; ++bd) {
486 PetscInt field;
487 DMBoundaryConditionType type;
488 DMLabel label;
489
490 PetscCall(PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL));
491 if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC;
492 }
493 }
494 /* Add ghost cell boundaries for FVM */
495 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
496 for (f = 0; f < Nf; ++f)
497 if (!isFE[f] && cEndInterior >= 0) ++numBC;
498 PetscCall(PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps));
499 /* Constrain ghost cells for FV */
500 for (f = 0; f < Nf; ++f) {
501 PetscInt *newidx, c;
502
503 if (isFE[f] || cEndInterior < 0) continue;
504 PetscCall(PetscMalloc1(cEnd - cEndInterior, &newidx));
505 for (c = cEndInterior; c < cEnd; ++c) newidx[c - cEndInterior] = c;
506 bcFields[bc] = f;
507 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cEnd - cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]));
508 }
509 /* Complete labels for boundary conditions */
510 PetscCall(DMCompleteBCLabels_Internal(dm));
511 /* Handle FEM Dirichlet boundaries */
512 PetscCall(DMGetNumDS(dm, &Nds));
513 for (s = 0; s < Nds; ++s) {
514 PetscDS dsBC;
515 PetscInt numBd, bd;
516
517 PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC, NULL));
518 PetscCall(PetscDSGetNumBoundary(dsBC, &numBd));
519 for (bd = 0; bd < numBd; ++bd) {
520 DMLabel label;
521 const PetscInt *comps;
522 const PetscInt *values;
523 PetscInt bd2, field, numComps, numValues;
524 DMBoundaryConditionType type;
525 PetscBool duplicate = PETSC_FALSE;
526
527 PetscCall(PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, &numValues, &values, &field, &numComps, &comps, NULL, NULL, NULL));
528 if (!isFE[field] || !label) continue;
529 /* Only want to modify label once */
530 for (bd2 = 0; bd2 < bd; ++bd2) {
531 DMLabel l;
532
533 PetscCall(PetscDSGetBoundary(dsBC, bd2, NULL, NULL, NULL, &l, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
534 duplicate = l == label ? PETSC_TRUE : PETSC_FALSE;
535 if (duplicate) break;
536 }
537 /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
538 if (type & DM_BC_ESSENTIAL) {
539 PetscInt *newidx;
540 PetscInt n, newn = 0, p, v;
541
542 bcFields[bc] = field;
543 if (numComps) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]));
544 for (v = 0; v < numValues; ++v) {
545 IS tmp;
546 const PetscInt *idx;
547
548 PetscCall(DMLabelGetStratumIS(label, values[v], &tmp));
549 if (!tmp) continue;
550 PetscCall(ISGetLocalSize(tmp, &n));
551 PetscCall(ISGetIndices(tmp, &idx));
552 if (isFE[field]) {
553 for (p = 0; p < n; ++p)
554 if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
555 } else {
556 for (p = 0; p < n; ++p)
557 if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
558 }
559 PetscCall(ISRestoreIndices(tmp, &idx));
560 PetscCall(ISDestroy(&tmp));
561 }
562 PetscCall(PetscMalloc1(newn, &newidx));
563 newn = 0;
564 for (v = 0; v < numValues; ++v) {
565 IS tmp;
566 const PetscInt *idx;
567
568 PetscCall(DMLabelGetStratumIS(label, values[v], &tmp));
569 if (!tmp) continue;
570 PetscCall(ISGetLocalSize(tmp, &n));
571 PetscCall(ISGetIndices(tmp, &idx));
572 if (isFE[field]) {
573 for (p = 0; p < n; ++p)
574 if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
575 } else {
576 for (p = 0; p < n; ++p)
577 if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
578 }
579 PetscCall(ISRestoreIndices(tmp, &idx));
580 PetscCall(ISDestroy(&tmp));
581 }
582 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]));
583 }
584 }
585 }
586 /* Handle discretization */
587 PetscCall(PetscCalloc3(Nf, &labels, Nf, &numComp, Nf * (dim + 1), &numDof));
588 for (f = 0; f < Nf; ++f) {
589 labels[f] = dm->fields[f].label;
590 if (isFE[f]) {
591 PetscFE fe = (PetscFE)dm->fields[f].disc;
592 const PetscInt *numFieldDof;
593 PetscInt fedim, d;
594
595 PetscCall(PetscFEGetNumComponents(fe, &numComp[f]));
596 PetscCall(PetscFEGetNumDof(fe, &numFieldDof));
597 PetscCall(PetscFEGetSpatialDimension(fe, &fedim));
598 for (d = 0; d < PetscMin(dim, fedim) + 1; ++d) numDof[f * (dim + 1) + d] = numFieldDof[d];
599 } else {
600 PetscFV fv = (PetscFV)dm->fields[f].disc;
601
602 PetscCall(PetscFVGetNumComponents(fv, &numComp[f]));
603 numDof[f * (dim + 1) + dim] = numComp[f];
604 }
605 }
606 PetscCall(DMPlexGetDepth(dm, &depth));
607 for (f = 0; f < Nf; ++f) {
608 PetscInt d;
609 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.");
610 }
611 PetscCall(DMCreateSectionPermutation(dm, &permIS, &blockStarts));
612 PetscCall(DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, permIS, §ion));
613 section->blockStarts = blockStarts;
614 PetscCall(ISDestroy(&permIS));
615 for (f = 0; f < Nf; ++f) {
616 PetscFE fe;
617 const char *name;
618
619 PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
620 if (!((PetscObject)fe)->name) continue;
621 PetscCall(PetscObjectGetName((PetscObject)fe, &name));
622 PetscCall(PetscSectionSetFieldName(section, f, name));
623 }
624 PetscCall(DMSetLocalSection(dm, section));
625 PetscCall(PetscSectionDestroy(§ion));
626 for (bc = 0; bc < numBC; ++bc) {
627 PetscCall(ISDestroy(&bcPoints[bc]));
628 PetscCall(ISDestroy(&bcComps[bc]));
629 }
630 PetscCall(PetscFree3(bcFields, bcPoints, bcComps));
631 PetscCall(PetscFree3(labels, numComp, numDof));
632 PetscCall(PetscFree(isFE));
633 /* Checking for CEED usage */
634 {
635 PetscBool useCeed;
636
637 PetscCall(DMPlexGetUseCeed(dm, &useCeed));
638 if (useCeed) {
639 PetscCall(DMPlexSetUseMatClosurePermutation(dm, PETSC_FALSE));
640 PetscCall(DMUseTensorOrder(dm, PETSC_TRUE));
641 }
642 }
643 PetscFunctionReturn(PETSC_SUCCESS);
644 }
645