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 { 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 = DMPolytopeTypeGetNumArrangments(ct) / 2; 72 } 73 PetscCall(PetscSectionSymLabelSetStratum(sym, depth - h, numDof[depth - h], -kConeSize, kConeSize, PETSC_USE_POINTER, perms0 ? &perms0[-kConeSize] : NULL, flips0 ? &flips0[-kConeSize] : NULL)); 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(0); 86 } 87 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; 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)); 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(0); 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 */ 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(0); 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 */ 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(0); 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(0); 353 } 354 355 /* Set the constrained indices on each point */ 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(0); 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 - The label 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*(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 boundry condition 410 . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies 411 . bcPoints - An array of size numBC giving an IS holding the Plex 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 Notes: 418 numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the 419 number of dof for field 0 on each edge. 420 421 The chart permutation is the same one set using PetscSectionSetPermutation() 422 423 Level: developer 424 425 TODO: How is this related to DMCreateLocalSection() 426 427 .seealso: `DMPlexCreate()`, `PetscSectionCreate()`, `PetscSectionSetPermutation()` 428 @*/ 429 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) 430 { 431 PetscSection aSec; 432 433 PetscFunctionBegin; 434 PetscCall(DMPlexCreateSectionFields(dm, numComp, section)); 435 PetscCall(DMPlexCreateSectionDof(dm, label, numDof, *section)); 436 PetscCall(DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section)); 437 if (perm) PetscCall(PetscSectionSetPermutation(*section, perm)); 438 PetscCall(PetscSectionSetFromOptions(*section)); 439 PetscCall(PetscSectionSetUp(*section)); 440 PetscCall(DMPlexGetAnchors(dm, &aSec, NULL)); 441 if (numBC || aSec) { 442 PetscCall(DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section)); 443 PetscCall(DMPlexCreateSectionBCIndices(dm, *section)); 444 } 445 PetscCall(PetscSectionViewFromOptions(*section, NULL, "-section_view")); 446 PetscFunctionReturn(0); 447 } 448 449 PetscErrorCode DMCreateLocalSection_Plex(DM dm) 450 { 451 PetscSection section; 452 DMLabel *labels; 453 IS *bcPoints, *bcComps; 454 PetscBool *isFE; 455 PetscInt *bcFields, *numComp, *numDof; 456 PetscInt depth, dim, numBC = 0, Nf, Nds, s, bc = 0, f; 457 PetscInt cStart, cEnd, cEndInterior; 458 459 PetscFunctionBegin; 460 PetscCall(DMGetNumFields(dm, &Nf)); 461 PetscCall(DMGetDimension(dm, &dim)); 462 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 463 /* FE and FV boundary conditions are handled slightly differently */ 464 PetscCall(PetscMalloc1(Nf, &isFE)); 465 for (f = 0; f < Nf; ++f) { 466 PetscObject obj; 467 PetscClassId id; 468 469 PetscCall(DMGetField(dm, f, NULL, &obj)); 470 PetscCall(PetscObjectGetClassId(obj, &id)); 471 if (id == PETSCFE_CLASSID) { 472 isFE[f] = PETSC_TRUE; 473 } else if (id == PETSCFV_CLASSID) { 474 isFE[f] = PETSC_FALSE; 475 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f); 476 } 477 /* Allocate boundary point storage for FEM boundaries */ 478 PetscCall(DMGetNumDS(dm, &Nds)); 479 for (s = 0; s < Nds; ++s) { 480 PetscDS dsBC; 481 PetscInt numBd, bd; 482 483 PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC)); 484 PetscCall(PetscDSGetNumBoundary(dsBC, &numBd)); 485 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)"); 486 for (bd = 0; bd < numBd; ++bd) { 487 PetscInt field; 488 DMBoundaryConditionType type; 489 DMLabel label; 490 491 PetscCall(PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL)); 492 if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC; 493 } 494 } 495 /* Add ghost cell boundaries for FVM */ 496 PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL)); 497 for (f = 0; f < Nf; ++f) 498 if (!isFE[f] && cEndInterior >= 0) ++numBC; 499 PetscCall(PetscCalloc3(numBC, &bcFields, numBC, &bcPoints, numBC, &bcComps)); 500 /* Constrain ghost cells for FV */ 501 for (f = 0; f < Nf; ++f) { 502 PetscInt *newidx, c; 503 504 if (isFE[f] || cEndInterior < 0) continue; 505 PetscCall(PetscMalloc1(cEnd - cEndInterior, &newidx)); 506 for (c = cEndInterior; c < cEnd; ++c) newidx[c - cEndInterior] = c; 507 bcFields[bc] = f; 508 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cEnd - cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++])); 509 } 510 /* Complete labels for boundary conditions */ 511 PetscCall(DMCompleteBCLabels_Internal(dm)); 512 /* Handle FEM Dirichlet boundaries */ 513 PetscCall(DMGetNumDS(dm, &Nds)); 514 for (s = 0; s < Nds; ++s) { 515 PetscDS dsBC; 516 PetscInt numBd, bd; 517 518 PetscCall(DMGetRegionNumDS(dm, s, NULL, NULL, &dsBC)); 519 PetscCall(PetscDSGetNumBoundary(dsBC, &numBd)); 520 for (bd = 0; bd < numBd; ++bd) { 521 DMLabel label; 522 const PetscInt *comps; 523 const PetscInt *values; 524 PetscInt bd2, field, numComps, numValues; 525 DMBoundaryConditionType type; 526 PetscBool duplicate = PETSC_FALSE; 527 528 PetscCall(PetscDSGetBoundary(dsBC, bd, NULL, &type, NULL, &label, &numValues, &values, &field, &numComps, &comps, NULL, NULL, NULL)); 529 if (!isFE[field] || !label) continue; 530 /* Only want to modify label once */ 531 for (bd2 = 0; bd2 < bd; ++bd2) { 532 DMLabel l; 533 534 PetscCall(PetscDSGetBoundary(dsBC, bd2, NULL, NULL, NULL, &l, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 535 duplicate = l == label ? PETSC_TRUE : PETSC_FALSE; 536 if (duplicate) break; 537 } 538 /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */ 539 if (type & DM_BC_ESSENTIAL) { 540 PetscInt *newidx; 541 PetscInt n, newn = 0, p, v; 542 543 bcFields[bc] = field; 544 if (numComps) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numComps, comps, PETSC_COPY_VALUES, &bcComps[bc])); 545 for (v = 0; v < numValues; ++v) { 546 IS tmp; 547 const PetscInt *idx; 548 549 PetscCall(DMLabelGetStratumIS(label, values[v], &tmp)); 550 if (!tmp) continue; 551 PetscCall(ISGetLocalSize(tmp, &n)); 552 PetscCall(ISGetIndices(tmp, &idx)); 553 if (isFE[field]) { 554 for (p = 0; p < n; ++p) 555 if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn; 556 } else { 557 for (p = 0; p < n; ++p) 558 if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn; 559 } 560 PetscCall(ISRestoreIndices(tmp, &idx)); 561 PetscCall(ISDestroy(&tmp)); 562 } 563 PetscCall(PetscMalloc1(newn, &newidx)); 564 newn = 0; 565 for (v = 0; v < numValues; ++v) { 566 IS tmp; 567 const PetscInt *idx; 568 569 PetscCall(DMLabelGetStratumIS(label, values[v], &tmp)); 570 if (!tmp) continue; 571 PetscCall(ISGetLocalSize(tmp, &n)); 572 PetscCall(ISGetIndices(tmp, &idx)); 573 if (isFE[field]) { 574 for (p = 0; p < n; ++p) 575 if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p]; 576 } else { 577 for (p = 0; p < n; ++p) 578 if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p]; 579 } 580 PetscCall(ISRestoreIndices(tmp, &idx)); 581 PetscCall(ISDestroy(&tmp)); 582 } 583 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++])); 584 } 585 } 586 } 587 /* Handle discretization */ 588 PetscCall(PetscCalloc3(Nf, &labels, Nf, &numComp, Nf * (dim + 1), &numDof)); 589 for (f = 0; f < Nf; ++f) { 590 labels[f] = dm->fields[f].label; 591 if (isFE[f]) { 592 PetscFE fe = (PetscFE)dm->fields[f].disc; 593 const PetscInt *numFieldDof; 594 PetscInt fedim, d; 595 596 PetscCall(PetscFEGetNumComponents(fe, &numComp[f])); 597 PetscCall(PetscFEGetNumDof(fe, &numFieldDof)); 598 PetscCall(PetscFEGetSpatialDimension(fe, &fedim)); 599 for (d = 0; d < PetscMin(dim, fedim) + 1; ++d) numDof[f * (dim + 1) + d] = numFieldDof[d]; 600 } else { 601 PetscFV fv = (PetscFV)dm->fields[f].disc; 602 603 PetscCall(PetscFVGetNumComponents(fv, &numComp[f])); 604 numDof[f * (dim + 1) + dim] = numComp[f]; 605 } 606 } 607 PetscCall(DMPlexGetDepth(dm, &depth)); 608 for (f = 0; f < Nf; ++f) { 609 PetscInt d; 610 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."); 611 } 612 PetscCall(DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion)); 613 for (f = 0; f < Nf; ++f) { 614 PetscFE fe; 615 const char *name; 616 617 PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe)); 618 if (!((PetscObject)fe)->name) continue; 619 PetscCall(PetscObjectGetName((PetscObject)fe, &name)); 620 PetscCall(PetscSectionSetFieldName(section, f, name)); 621 } 622 PetscCall(DMSetLocalSection(dm, section)); 623 PetscCall(PetscSectionDestroy(§ion)); 624 for (bc = 0; bc < numBC; ++bc) { 625 PetscCall(ISDestroy(&bcPoints[bc])); 626 PetscCall(ISDestroy(&bcComps[bc])); 627 } 628 PetscCall(PetscFree3(bcFields, bcPoints, bcComps)); 629 PetscCall(PetscFree3(labels, numComp, numDof)); 630 PetscCall(PetscFree(isFE)); 631 PetscFunctionReturn(0); 632 } 633