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