1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/ 3 #include <petscsf.h> 4 5 static PetscErrorCode DMPlexCellIsHybrid_Internal(DM dm, PetscInt p, PetscBool *isHybrid) 6 { 7 DMPolytopeType ct; 8 9 PetscFunctionBegin; 10 PetscCall(DMPlexGetCellType(dm, p, &ct)); 11 switch (ct) { 12 case DM_POLYTOPE_POINT_PRISM_TENSOR: 13 case DM_POLYTOPE_SEG_PRISM_TENSOR: 14 case DM_POLYTOPE_TRI_PRISM_TENSOR: 15 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 16 *isHybrid = PETSC_TRUE; 17 break; 18 default: 19 *isHybrid = PETSC_FALSE; 20 break; 21 } 22 PetscFunctionReturn(PETSC_SUCCESS); 23 } 24 25 static PetscErrorCode DMPlexGetTensorPrismBounds_Internal(DM dm, PetscInt dim, PetscInt *cStart, PetscInt *cEnd) 26 { 27 DMLabel ctLabel; 28 29 PetscFunctionBegin; 30 if (cStart) *cStart = -1; 31 if (cEnd) *cEnd = -1; 32 PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel)); 33 switch (dim) { 34 case 1: 35 PetscCall(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd)); 36 break; 37 case 2: 38 PetscCall(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd)); 39 break; 40 case 3: 41 PetscCall(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd)); 42 if (*cStart < 0) PetscCall(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_QUAD_PRISM_TENSOR, cStart, cEnd)); 43 break; 44 default: 45 PetscFunctionReturn(PETSC_SUCCESS); 46 } 47 PetscFunctionReturn(PETSC_SUCCESS); 48 } 49 50 PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM dm, PetscInt val, PetscInt cellHeight, DMLabel label, PetscBool missing_only) 51 { 52 PetscInt depth, pStart, pEnd, fStart, fEnd, f, supportSize, nroots = -1, nleaves = -1, defval; 53 PetscSF sf; 54 const PetscSFNode *iremote = NULL; 55 const PetscInt *ilocal = NULL; 56 PetscInt *leafData; 57 58 PetscFunctionBegin; 59 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 60 PetscCall(DMPlexGetDepth(dm, &depth)); 61 if (depth >= cellHeight + 1) { 62 PetscCall(DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd)); 63 } else { 64 /* Note DMPlexGetHeightStratum() returns fStart, fEnd = pStart, pEnd */ 65 /* if height > depth, which is not what we want here. */ 66 fStart = 0; 67 fEnd = 0; 68 } 69 PetscCall(DMLabelGetDefaultValue(label, &defval)); 70 PetscCall(PetscCalloc1(pEnd - pStart, &leafData)); 71 leafData = PetscSafePointerPlusOffset(leafData, -pStart); 72 PetscCall(DMGetPointSF(dm, &sf)); 73 if (sf) PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &ilocal, &iremote)); 74 if (sf && nroots >= 0) { 75 PetscInt cStart, cEnd, c, i; 76 PetscInt *rootData, *rootData1, *cellOwners, hasTwoSupportCells = -2; 77 const PetscInt *support; 78 PetscMPIInt rank; 79 80 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 81 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 82 PetscCall(PetscCalloc3(pEnd - pStart, &rootData, pEnd - pStart, &rootData1, cEnd - cStart, &cellOwners)); 83 rootData -= pStart; 84 rootData1 -= pStart; 85 for (c = cStart; c < cEnd; ++c) cellOwners[c - cStart] = rank; 86 for (i = 0; i < nleaves; ++i) { 87 c = ilocal ? ilocal[i] : i; 88 if (c >= cStart && c < cEnd) cellOwners[c - cStart] = iremote[i].rank; 89 } 90 for (f = fStart; f < fEnd; ++f) { 91 PetscCall(DMPlexGetSupportSize(dm, f, &supportSize)); 92 if (supportSize == 1) { 93 PetscCall(DMPlexGetSupport(dm, f, &support)); 94 leafData[f] = cellOwners[support[0] - cStart]; 95 } else { 96 /* TODO: When using DMForest, we could have a parent facet (a coarse facet) */ 97 /* supportSize of which alone does not tell us if it is an interior */ 98 /* facet or an exterior facet. Those facets can be identified by checking */ 99 /* if they are in the parent tree. We should probably skip those parent */ 100 /* facets here, which will allow for including the following check, and */ 101 /* see if they are exterior facets or not afterwards by checking if the */ 102 /* children are exterior or not. */ 103 /* PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid facet support size (%" PetscInt_FMT ") on facet (%" PetscInt_FMT ")", supportSize, f); */ 104 leafData[f] = hasTwoSupportCells; /* some negative PetscInt */ 105 } 106 rootData[f] = leafData[f]; 107 rootData1[f] = leafData[f]; 108 } 109 PetscCall(PetscSFReduceBegin(sf, MPIU_INT, leafData, rootData, MPI_MIN)); 110 PetscCall(PetscSFReduceBegin(sf, MPIU_INT, leafData, rootData1, MPI_MAX)); 111 PetscCall(PetscSFReduceEnd(sf, MPIU_INT, leafData, rootData, MPI_MIN)); 112 PetscCall(PetscSFReduceEnd(sf, MPIU_INT, leafData, rootData1, MPI_MAX)); 113 for (f = fStart; f < fEnd; ++f) { 114 /* Store global support size of f. */ 115 /* Facet f is an interior facet if and only if one of the following two is satisfied: */ 116 /* 1. supportSize is 2 on some rank. */ 117 /* 2. supportSize is 1 on any rank that can see f, but f is on a partition boundary; */ 118 /* i.e., rootData[f] < rootData1[f]. */ 119 rootData[f] = (rootData[f] == hasTwoSupportCells || (rootData[f] < rootData1[f])) ? 2 : 1; 120 leafData[f] = rootData[f]; 121 } 122 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, rootData, leafData, MPI_REPLACE)); 123 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, rootData, leafData, MPI_REPLACE)); 124 rootData += pStart; 125 rootData1 += pStart; 126 PetscCall(PetscFree3(rootData, rootData1, cellOwners)); 127 } else { 128 for (f = fStart; f < fEnd; ++f) { 129 PetscCall(DMPlexGetSupportSize(dm, f, &supportSize)); 130 leafData[f] = supportSize; 131 } 132 } 133 for (f = fStart; f < fEnd; ++f) { 134 if (leafData[f] == 1) { 135 if (val < 0) { 136 PetscInt *closure = NULL; 137 PetscInt clSize, cl, cval; 138 139 PetscAssert(!missing_only, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented"); 140 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure)); 141 for (cl = 0; cl < clSize * 2; cl += 2) { 142 PetscCall(DMLabelGetValue(label, closure[cl], &cval)); 143 if (cval < 0) continue; 144 PetscCall(DMLabelSetValue(label, f, cval)); 145 break; 146 } 147 if (cl == clSize * 2) PetscCall(DMLabelSetValue(label, f, 1)); 148 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure)); 149 } else { 150 if (missing_only) { 151 PetscInt fval; 152 PetscCall(DMLabelGetValue(label, f, &fval)); 153 if (fval != defval) PetscCall(DMLabelClearValue(label, f, fval)); 154 else PetscCall(DMLabelSetValue(label, f, val)); 155 } else { 156 PetscCall(DMLabelSetValue(label, f, val)); 157 } 158 } 159 } else { 160 /* TODO: See the above comment on DMForest */ 161 /* PetscCheck(leafData[f] == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid facet support size (%" PetscInt_FMT ") on facet (%" PetscInt_FMT ")", leafData[f], f); */ 162 } 163 } 164 leafData = PetscSafePointerPlusOffset(leafData, pStart); 165 PetscCall(PetscFree(leafData)); 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 /*@ 170 DMPlexMarkBoundaryFaces - Mark all faces on the boundary 171 172 Collective 173 174 Input Parameters: 175 + dm - The original `DM` 176 - val - The marker value, or `PETSC_DETERMINE` to use some value in the closure (or 1 if none are found) 177 178 Output Parameter: 179 . label - The `DMLabel` marking boundary faces with the given value 180 181 Level: developer 182 183 Note: 184 This function will use the point `PetscSF` from the input `DM` and the ownership of the support cells to exclude points on the partition boundary from being marked. If you also wish to mark the partition boundary, you can use `DMSetPointSF()` to temporarily set it to `NULL`, and then reset it to the original object after the call. 185 186 In DMForest there can be facets support sizes of which alone can not determine whether they are on the boundary. Currently, this function is not guaranteed to produce the correct result in such case. 187 188 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabelCreate()`, `DMCreateLabel()` 189 @*/ 190 PetscErrorCode DMPlexMarkBoundaryFaces(DM dm, PetscInt val, DMLabel label) 191 { 192 DM plex; 193 DMPlexInterpolatedFlag flg; 194 195 PetscFunctionBegin; 196 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 197 PetscCall(DMConvert(dm, DMPLEX, &plex)); 198 PetscCall(DMPlexIsInterpolated(plex, &flg)); 199 PetscCheck(flg == DMPLEX_INTERPOLATED_FULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not fully interpolated on this rank"); 200 PetscCall(DMPlexMarkBoundaryFaces_Internal(plex, val, 0, label, PETSC_FALSE)); 201 PetscCall(DMDestroy(&plex)); 202 PetscFunctionReturn(PETSC_SUCCESS); 203 } 204 205 static PetscErrorCode DMPlexLabelComplete_Internal(DM dm, DMLabel label, PetscBool completeCells) 206 { 207 IS valueIS; 208 PetscSF sfPoint; 209 const PetscInt *values; 210 PetscInt numValues, v, cStart, cEnd, nroots; 211 212 PetscFunctionBegin; 213 PetscCall(DMLabelGetNumValues(label, &numValues)); 214 PetscCall(DMLabelGetValueIS(label, &valueIS)); 215 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 216 PetscCall(ISGetIndices(valueIS, &values)); 217 for (v = 0; v < numValues; ++v) { 218 IS pointIS; 219 const PetscInt *points; 220 PetscInt numPoints, p; 221 222 PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints)); 223 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 224 PetscCall(ISGetIndices(pointIS, &points)); 225 for (p = 0; p < numPoints; ++p) { 226 PetscInt q = points[p]; 227 PetscInt *closure = NULL; 228 PetscInt closureSize, c; 229 230 if (cStart <= q && q < cEnd && !completeCells) { /* skip cells */ 231 continue; 232 } 233 PetscCall(DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure)); 234 for (c = 0; c < closureSize * 2; c += 2) PetscCall(DMLabelSetValue(label, closure[c], values[v])); 235 PetscCall(DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure)); 236 } 237 PetscCall(ISRestoreIndices(pointIS, &points)); 238 PetscCall(ISDestroy(&pointIS)); 239 } 240 PetscCall(ISRestoreIndices(valueIS, &values)); 241 PetscCall(ISDestroy(&valueIS)); 242 PetscCall(DMGetPointSF(dm, &sfPoint)); 243 PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 244 if (nroots >= 0) { 245 DMLabel lblRoots, lblLeaves; 246 IS valueIS, pointIS; 247 const PetscInt *values; 248 PetscInt numValues, v; 249 250 /* Pull point contributions from remote leaves into local roots */ 251 PetscCall(DMLabelGather(label, sfPoint, &lblLeaves)); 252 PetscCall(DMLabelGetValueIS(lblLeaves, &valueIS)); 253 PetscCall(ISGetLocalSize(valueIS, &numValues)); 254 PetscCall(ISGetIndices(valueIS, &values)); 255 for (v = 0; v < numValues; ++v) { 256 const PetscInt value = values[v]; 257 258 PetscCall(DMLabelGetStratumIS(lblLeaves, value, &pointIS)); 259 PetscCall(DMLabelInsertIS(label, pointIS, value)); 260 PetscCall(ISDestroy(&pointIS)); 261 } 262 PetscCall(ISRestoreIndices(valueIS, &values)); 263 PetscCall(ISDestroy(&valueIS)); 264 PetscCall(DMLabelDestroy(&lblLeaves)); 265 /* Push point contributions from roots into remote leaves */ 266 PetscCall(DMLabelDistribute(label, sfPoint, &lblRoots)); 267 PetscCall(DMLabelGetValueIS(lblRoots, &valueIS)); 268 PetscCall(ISGetLocalSize(valueIS, &numValues)); 269 PetscCall(ISGetIndices(valueIS, &values)); 270 for (v = 0; v < numValues; ++v) { 271 const PetscInt value = values[v]; 272 273 PetscCall(DMLabelGetStratumIS(lblRoots, value, &pointIS)); 274 PetscCall(DMLabelInsertIS(label, pointIS, value)); 275 PetscCall(ISDestroy(&pointIS)); 276 } 277 PetscCall(ISRestoreIndices(valueIS, &values)); 278 PetscCall(ISDestroy(&valueIS)); 279 PetscCall(DMLabelDestroy(&lblRoots)); 280 } 281 PetscFunctionReturn(PETSC_SUCCESS); 282 } 283 284 /*@ 285 DMPlexLabelComplete - Starting with a label marking points on a surface, we add the transitive closure to the surface 286 287 Input Parameters: 288 + dm - The `DM` 289 - label - A `DMLabel` marking the surface points 290 291 Output Parameter: 292 . label - A `DMLabel` marking all surface points in the transitive closure 293 294 Level: developer 295 296 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelCohesiveComplete()` 297 @*/ 298 PetscErrorCode DMPlexLabelComplete(DM dm, DMLabel label) 299 { 300 PetscFunctionBegin; 301 PetscCall(DMPlexLabelComplete_Internal(dm, label, PETSC_TRUE)); 302 PetscFunctionReturn(PETSC_SUCCESS); 303 } 304 305 /*@ 306 DMPlexLabelAddCells - Starting with a label marking points on a surface, we add a cell for each point 307 308 Input Parameters: 309 + dm - The `DM` 310 - label - A `DMLabel` marking the surface points 311 312 Output Parameter: 313 . label - A `DMLabel` incorporating cells 314 315 Level: developer 316 317 Note: 318 The cells allow FEM boundary conditions to be applied using the cell geometry 319 320 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelAddFaceCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()` 321 @*/ 322 PetscErrorCode DMPlexLabelAddCells(DM dm, DMLabel label) 323 { 324 IS valueIS; 325 const PetscInt *values; 326 PetscInt numValues, v, csStart, csEnd, chStart, chEnd; 327 328 PetscFunctionBegin; 329 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &csStart, &csEnd)); 330 PetscCall(DMPlexGetHeightStratum(dm, 0, &chStart, &chEnd)); 331 PetscCall(DMLabelGetNumValues(label, &numValues)); 332 PetscCall(DMLabelGetValueIS(label, &valueIS)); 333 PetscCall(ISGetIndices(valueIS, &values)); 334 for (v = 0; v < numValues; ++v) { 335 IS pointIS; 336 const PetscInt *points; 337 PetscInt numPoints, p; 338 339 PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints)); 340 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 341 PetscCall(ISGetIndices(pointIS, &points)); 342 for (p = 0; p < numPoints; ++p) { 343 const PetscInt point = points[p]; 344 PetscInt *closure = NULL; 345 PetscInt closureSize, cl, h, cStart, cEnd; 346 DMPolytopeType ct; 347 348 // If the point is a hybrid, allow hybrid cells 349 PetscCall(DMPlexGetCellType(dm, point, &ct)); 350 PetscCall(DMPlexGetPointHeight(dm, point, &h)); 351 if (DMPolytopeTypeIsHybrid(ct)) { 352 cStart = chStart; 353 cEnd = chEnd; 354 } else { 355 cStart = csStart; 356 cEnd = csEnd; 357 } 358 359 PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure)); 360 for (cl = closureSize - 1; cl > 0; --cl) { 361 const PetscInt cell = closure[cl * 2]; 362 if (cell >= cStart && cell < cEnd) { 363 PetscCall(DMLabelSetValue(label, cell, values[v])); 364 break; 365 } 366 } 367 PetscCall(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure)); 368 } 369 PetscCall(ISRestoreIndices(pointIS, &points)); 370 PetscCall(ISDestroy(&pointIS)); 371 } 372 PetscCall(ISRestoreIndices(valueIS, &values)); 373 PetscCall(ISDestroy(&valueIS)); 374 PetscFunctionReturn(PETSC_SUCCESS); 375 } 376 377 /*@ 378 DMPlexLabelAddFaceCells - Starting with a label marking faces on a surface, we add a cell for each face 379 380 Input Parameters: 381 + dm - The `DM` 382 - label - A `DMLabel` marking the surface points 383 384 Output Parameter: 385 . label - A `DMLabel` incorporating cells 386 387 Level: developer 388 389 Note: 390 The cells allow FEM boundary conditions to be applied using the cell geometry 391 392 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelAddCells()`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()` 393 @*/ 394 PetscErrorCode DMPlexLabelAddFaceCells(DM dm, DMLabel label) 395 { 396 IS valueIS; 397 const PetscInt *values; 398 PetscInt numValues, v, cStart, cEnd, fStart, fEnd; 399 400 PetscFunctionBegin; 401 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 402 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 403 PetscCall(DMLabelGetNumValues(label, &numValues)); 404 PetscCall(DMLabelGetValueIS(label, &valueIS)); 405 PetscCall(ISGetIndices(valueIS, &values)); 406 for (v = 0; v < numValues; ++v) { 407 IS pointIS; 408 const PetscInt *points; 409 PetscInt numPoints, p; 410 411 PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints)); 412 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 413 PetscCall(ISGetIndices(pointIS, &points)); 414 for (p = 0; p < numPoints; ++p) { 415 const PetscInt face = points[p]; 416 PetscInt *closure = NULL; 417 PetscInt closureSize, cl; 418 419 if ((face < fStart) || (face >= fEnd)) continue; 420 PetscCall(DMPlexGetTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure)); 421 for (cl = closureSize - 1; cl > 0; --cl) { 422 const PetscInt cell = closure[cl * 2]; 423 if (cell >= cStart && cell < cEnd) { 424 PetscCall(DMLabelSetValue(label, cell, values[v])); 425 break; 426 } 427 } 428 PetscCall(DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure)); 429 } 430 PetscCall(ISRestoreIndices(pointIS, &points)); 431 PetscCall(ISDestroy(&pointIS)); 432 } 433 PetscCall(ISRestoreIndices(valueIS, &values)); 434 PetscCall(ISDestroy(&valueIS)); 435 PetscFunctionReturn(PETSC_SUCCESS); 436 } 437 438 /*@ 439 DMPlexLabelClearCells - Remove cells from a label 440 441 Input Parameters: 442 + dm - The `DM` 443 - label - A `DMLabel` marking surface points and their adjacent cells 444 445 Output Parameter: 446 . label - A `DMLabel` without cells 447 448 Level: developer 449 450 Note: 451 This undoes `DMPlexLabelAddCells()` or `DMPlexLabelAddFaceCells()` 452 453 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexLabelComplete()`, `DMPlexLabelCohesiveComplete()`, `DMPlexLabelAddCells()` 454 @*/ 455 PetscErrorCode DMPlexLabelClearCells(DM dm, DMLabel label) 456 { 457 IS valueIS; 458 const PetscInt *values; 459 PetscInt numValues, v, cStart, cEnd; 460 461 PetscFunctionBegin; 462 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 463 PetscCall(DMLabelGetNumValues(label, &numValues)); 464 PetscCall(DMLabelGetValueIS(label, &valueIS)); 465 PetscCall(ISGetIndices(valueIS, &values)); 466 for (v = 0; v < numValues; ++v) { 467 IS pointIS; 468 const PetscInt *points; 469 PetscInt numPoints, p; 470 471 PetscCall(DMLabelGetStratumSize(label, values[v], &numPoints)); 472 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 473 PetscCall(ISGetIndices(pointIS, &points)); 474 for (p = 0; p < numPoints; ++p) { 475 PetscInt point = points[p]; 476 477 if (point >= cStart && point < cEnd) PetscCall(DMLabelClearValue(label, point, values[v])); 478 } 479 PetscCall(ISRestoreIndices(pointIS, &points)); 480 PetscCall(ISDestroy(&pointIS)); 481 } 482 PetscCall(ISRestoreIndices(valueIS, &values)); 483 PetscCall(ISDestroy(&valueIS)); 484 PetscFunctionReturn(PETSC_SUCCESS); 485 } 486 487 /* take (oldEnd, added) pairs, ordered by height and convert them to (oldstart, newstart) pairs, ordered by ascending 488 * index (skipping first, which is (0,0)) */ 489 static inline PetscErrorCode DMPlexShiftPointSetUp_Internal(PetscInt depth, PetscInt depthShift[]) 490 { 491 PetscInt d, off = 0; 492 493 PetscFunctionBegin; 494 /* sort by (oldend): yes this is an O(n^2) sort, we expect depth <= 3 */ 495 for (d = 0; d < depth; d++) { 496 PetscInt firstd = d; 497 PetscInt firstStart = depthShift[2 * d]; 498 PetscInt e; 499 500 for (e = d + 1; e <= depth; e++) { 501 if (depthShift[2 * e] < firstStart) { 502 firstd = e; 503 firstStart = depthShift[2 * d]; 504 } 505 } 506 if (firstd != d) { 507 PetscInt swap[2]; 508 509 e = firstd; 510 swap[0] = depthShift[2 * d]; 511 swap[1] = depthShift[2 * d + 1]; 512 depthShift[2 * d] = depthShift[2 * e]; 513 depthShift[2 * d + 1] = depthShift[2 * e + 1]; 514 depthShift[2 * e] = swap[0]; 515 depthShift[2 * e + 1] = swap[1]; 516 } 517 } 518 /* convert (oldstart, added) to (oldstart, newstart) */ 519 for (d = 0; d <= depth; d++) { 520 off += depthShift[2 * d + 1]; 521 depthShift[2 * d + 1] = depthShift[2 * d] + off; 522 } 523 PetscFunctionReturn(PETSC_SUCCESS); 524 } 525 526 /* depthShift is a list of (old, new) pairs */ 527 static inline PetscInt DMPlexShiftPoint_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[]) 528 { 529 PetscInt d; 530 PetscInt newOff = 0; 531 532 for (d = 0; d <= depth; d++) { 533 if (p < depthShift[2 * d]) return p + newOff; 534 else newOff = depthShift[2 * d + 1] - depthShift[2 * d]; 535 } 536 return p + newOff; 537 } 538 539 /* depthShift is a list of (old, new) pairs */ 540 static inline PetscInt DMPlexShiftPointInverse_Internal(PetscInt p, PetscInt depth, PetscInt depthShift[]) 541 { 542 PetscInt d; 543 PetscInt newOff = 0; 544 545 for (d = 0; d <= depth; d++) { 546 if (p < depthShift[2 * d + 1]) return p + newOff; 547 else newOff = depthShift[2 * d] - depthShift[2 * d + 1]; 548 } 549 return p + newOff; 550 } 551 552 static PetscErrorCode DMPlexShiftSizes_Internal(DM dm, PetscInt depthShift[], DM dmNew) 553 { 554 PetscInt depth = 0, d, pStart, pEnd, p; 555 DMLabel depthLabel; 556 557 PetscFunctionBegin; 558 PetscCall(DMPlexGetDepth(dm, &depth)); 559 if (depth < 0) PetscFunctionReturn(PETSC_SUCCESS); 560 /* Step 1: Expand chart */ 561 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 562 pEnd = DMPlexShiftPoint_Internal(pEnd, depth, depthShift); 563 PetscCall(DMPlexSetChart(dmNew, pStart, pEnd)); 564 PetscCall(DMCreateLabel(dmNew, "depth")); 565 PetscCall(DMPlexGetDepthLabel(dmNew, &depthLabel)); 566 PetscCall(DMCreateLabel(dmNew, "celltype")); 567 /* Step 2: Set cone and support sizes */ 568 for (d = 0; d <= depth; ++d) { 569 PetscInt pStartNew, pEndNew; 570 IS pIS; 571 572 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 573 pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift); 574 pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift); 575 PetscCall(ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS)); 576 PetscCall(DMLabelSetStratumIS(depthLabel, d, pIS)); 577 PetscCall(ISDestroy(&pIS)); 578 for (p = pStart; p < pEnd; ++p) { 579 PetscInt newp = DMPlexShiftPoint_Internal(p, depth, depthShift); 580 PetscInt size; 581 DMPolytopeType ct; 582 583 PetscCall(DMPlexGetConeSize(dm, p, &size)); 584 PetscCall(DMPlexSetConeSize(dmNew, newp, size)); 585 PetscCall(DMPlexGetSupportSize(dm, p, &size)); 586 PetscCall(DMPlexSetSupportSize(dmNew, newp, size)); 587 PetscCall(DMPlexGetCellType(dm, p, &ct)); 588 PetscCall(DMPlexSetCellType(dmNew, newp, ct)); 589 } 590 } 591 PetscFunctionReturn(PETSC_SUCCESS); 592 } 593 594 static PetscErrorCode DMPlexShiftPoints_Internal(DM dm, PetscInt depthShift[], DM dmNew) 595 { 596 PetscInt *newpoints; 597 PetscInt depth = 0, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, pStart, pEnd, p; 598 599 PetscFunctionBegin; 600 PetscCall(DMPlexGetDepth(dm, &depth)); 601 if (depth < 0) PetscFunctionReturn(PETSC_SUCCESS); 602 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 603 PetscCall(DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew)); 604 PetscCall(PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)), &newpoints)); 605 /* Step 5: Set cones and supports */ 606 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 607 for (p = pStart; p < pEnd; ++p) { 608 const PetscInt *points = NULL, *orientations = NULL; 609 PetscInt size, sizeNew, i, newp = DMPlexShiftPoint_Internal(p, depth, depthShift); 610 611 PetscCall(DMPlexGetConeSize(dm, p, &size)); 612 PetscCall(DMPlexGetCone(dm, p, &points)); 613 PetscCall(DMPlexGetConeOrientation(dm, p, &orientations)); 614 for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift); 615 PetscCall(DMPlexSetCone(dmNew, newp, newpoints)); 616 PetscCall(DMPlexSetConeOrientation(dmNew, newp, orientations)); 617 PetscCall(DMPlexGetSupportSize(dm, p, &size)); 618 PetscCall(DMPlexGetSupportSize(dmNew, newp, &sizeNew)); 619 PetscCall(DMPlexGetSupport(dm, p, &points)); 620 for (i = 0; i < size; ++i) newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift); 621 for (i = size; i < sizeNew; ++i) newpoints[i] = 0; 622 PetscCall(DMPlexSetSupport(dmNew, newp, newpoints)); 623 } 624 PetscCall(PetscFree(newpoints)); 625 PetscFunctionReturn(PETSC_SUCCESS); 626 } 627 628 static PetscErrorCode DMPlexShiftCoordinates_Internal(DM dm, PetscInt depthShift[], DM dmNew) 629 { 630 PetscSection coordSection, newCoordSection; 631 Vec coordinates, newCoordinates; 632 PetscScalar *coords, *newCoords; 633 PetscInt coordSize, sStart, sEnd; 634 PetscInt dim, depth = 0, cStart, cEnd, cStartNew, cEndNew, c, vStart, vEnd, vStartNew, vEndNew, v; 635 PetscBool hasCells; 636 637 PetscFunctionBegin; 638 PetscCall(DMGetCoordinateDim(dm, &dim)); 639 PetscCall(DMSetCoordinateDim(dmNew, dim)); 640 PetscCall(DMPlexGetDepth(dm, &depth)); 641 /* Step 8: Convert coordinates */ 642 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 643 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 644 PetscCall(DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew)); 645 PetscCall(DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew)); 646 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 647 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection)); 648 PetscCall(PetscSectionSetNumFields(newCoordSection, 1)); 649 PetscCall(PetscSectionSetFieldComponents(newCoordSection, 0, dim)); 650 PetscCall(PetscSectionGetChart(coordSection, &sStart, &sEnd)); 651 hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE; 652 PetscCall(PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew)); 653 if (hasCells) { 654 for (c = cStart; c < cEnd; ++c) { 655 PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof; 656 657 PetscCall(PetscSectionGetDof(coordSection, c, &dof)); 658 PetscCall(PetscSectionSetDof(newCoordSection, cNew, dof)); 659 PetscCall(PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof)); 660 } 661 } 662 for (v = vStartNew; v < vEndNew; ++v) { 663 PetscCall(PetscSectionSetDof(newCoordSection, v, dim)); 664 PetscCall(PetscSectionSetFieldDof(newCoordSection, v, 0, dim)); 665 } 666 PetscCall(PetscSectionSetUp(newCoordSection)); 667 PetscCall(DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection)); 668 PetscCall(PetscSectionGetStorageSize(newCoordSection, &coordSize)); 669 PetscCall(VecCreate(PETSC_COMM_SELF, &newCoordinates)); 670 PetscCall(PetscObjectSetName((PetscObject)newCoordinates, "coordinates")); 671 PetscCall(VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE)); 672 PetscCall(VecSetBlockSize(newCoordinates, dim)); 673 PetscCall(VecSetType(newCoordinates, VECSTANDARD)); 674 PetscCall(DMSetCoordinatesLocal(dmNew, newCoordinates)); 675 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 676 PetscCall(VecGetArray(coordinates, &coords)); 677 PetscCall(VecGetArray(newCoordinates, &newCoords)); 678 if (hasCells) { 679 for (c = cStart; c < cEnd; ++c) { 680 PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d; 681 682 PetscCall(PetscSectionGetDof(coordSection, c, &dof)); 683 PetscCall(PetscSectionGetOffset(coordSection, c, &off)); 684 PetscCall(PetscSectionGetOffset(newCoordSection, cNew, &noff)); 685 for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d]; 686 } 687 } 688 for (v = vStart; v < vEnd; ++v) { 689 PetscInt dof, off, noff, d; 690 691 PetscCall(PetscSectionGetDof(coordSection, v, &dof)); 692 PetscCall(PetscSectionGetOffset(coordSection, v, &off)); 693 PetscCall(PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff)); 694 for (d = 0; d < dof; ++d) newCoords[noff + d] = coords[off + d]; 695 } 696 PetscCall(VecRestoreArray(coordinates, &coords)); 697 PetscCall(VecRestoreArray(newCoordinates, &newCoords)); 698 PetscCall(VecDestroy(&newCoordinates)); 699 PetscCall(PetscSectionDestroy(&newCoordSection)); 700 PetscFunctionReturn(PETSC_SUCCESS); 701 } 702 703 static PetscErrorCode DMPlexShiftSF_Single(DM dm, PetscInt depthShift[], PetscSF sf, PetscSF sfNew) 704 { 705 const PetscSFNode *remotePoints; 706 PetscSFNode *gremotePoints; 707 const PetscInt *localPoints; 708 PetscInt *glocalPoints, *newLocation, *newRemoteLocation; 709 PetscInt numRoots, numLeaves, l, pStart, pEnd, depth = 0, totShift = 0; 710 711 PetscFunctionBegin; 712 PetscCall(DMPlexGetDepth(dm, &depth)); 713 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 714 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 715 totShift = DMPlexShiftPoint_Internal(pEnd, depth, depthShift) - pEnd; 716 if (numRoots >= 0) { 717 PetscCall(PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation)); 718 for (l = 0; l < numRoots; ++l) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift); 719 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE)); 720 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE)); 721 PetscCall(PetscMalloc1(numLeaves, &glocalPoints)); 722 PetscCall(PetscMalloc1(numLeaves, &gremotePoints)); 723 for (l = 0; l < numLeaves; ++l) { 724 glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift); 725 gremotePoints[l].rank = remotePoints[l].rank; 726 gremotePoints[l].index = newRemoteLocation[localPoints[l]]; 727 } 728 PetscCall(PetscFree2(newLocation, newRemoteLocation)); 729 PetscCall(PetscSFSetGraph(sfNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER)); 730 } 731 PetscFunctionReturn(PETSC_SUCCESS); 732 } 733 734 static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew) 735 { 736 PetscSF sfPoint, sfPointNew; 737 PetscBool useNatural; 738 739 PetscFunctionBegin; 740 /* Step 9: Convert pointSF */ 741 PetscCall(DMGetPointSF(dm, &sfPoint)); 742 PetscCall(DMGetPointSF(dmNew, &sfPointNew)); 743 PetscCall(DMPlexShiftSF_Single(dm, depthShift, sfPoint, sfPointNew)); 744 /* Step 9b: Convert naturalSF */ 745 PetscCall(DMGetUseNatural(dm, &useNatural)); 746 if (useNatural) { 747 PetscSF sfNat, sfNatNew; 748 749 PetscCall(DMSetUseNatural(dmNew, useNatural)); 750 PetscCall(DMGetNaturalSF(dm, &sfNat)); 751 PetscCall(DMGetNaturalSF(dmNew, &sfNatNew)); 752 PetscCall(DMPlexShiftSF_Single(dm, depthShift, sfNat, sfNatNew)); 753 } 754 PetscFunctionReturn(PETSC_SUCCESS); 755 } 756 757 static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew) 758 { 759 PetscInt depth = 0, numLabels, l; 760 761 PetscFunctionBegin; 762 PetscCall(DMPlexGetDepth(dm, &depth)); 763 /* Step 10: Convert labels */ 764 PetscCall(DMGetNumLabels(dm, &numLabels)); 765 for (l = 0; l < numLabels; ++l) { 766 DMLabel label, newlabel; 767 const char *lname; 768 PetscBool isDepth, isDim; 769 IS valueIS; 770 const PetscInt *values; 771 PetscInt numValues, val; 772 773 PetscCall(DMGetLabelName(dm, l, &lname)); 774 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 775 if (isDepth) continue; 776 PetscCall(PetscStrcmp(lname, "dim", &isDim)); 777 if (isDim) continue; 778 PetscCall(DMCreateLabel(dmNew, lname)); 779 PetscCall(DMGetLabel(dm, lname, &label)); 780 PetscCall(DMGetLabel(dmNew, lname, &newlabel)); 781 PetscCall(DMLabelGetDefaultValue(label, &val)); 782 PetscCall(DMLabelSetDefaultValue(newlabel, val)); 783 PetscCall(DMLabelGetValueIS(label, &valueIS)); 784 PetscCall(ISGetLocalSize(valueIS, &numValues)); 785 PetscCall(ISGetIndices(valueIS, &values)); 786 for (val = 0; val < numValues; ++val) { 787 IS pointIS; 788 const PetscInt *points; 789 PetscInt numPoints, p; 790 791 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 792 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 793 PetscCall(ISGetIndices(pointIS, &points)); 794 for (p = 0; p < numPoints; ++p) { 795 const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift); 796 797 PetscCall(DMLabelSetValue(newlabel, newpoint, values[val])); 798 } 799 PetscCall(ISRestoreIndices(pointIS, &points)); 800 PetscCall(ISDestroy(&pointIS)); 801 } 802 PetscCall(ISRestoreIndices(valueIS, &values)); 803 PetscCall(ISDestroy(&valueIS)); 804 } 805 PetscFunctionReturn(PETSC_SUCCESS); 806 } 807 808 static PetscErrorCode DMPlexCreateVTKLabel_Internal(DM dm, PetscBool createGhostLabel, DM dmNew) 809 { 810 PetscSF sfPoint; 811 DMLabel vtkLabel, ghostLabel = NULL; 812 const PetscSFNode *leafRemote; 813 const PetscInt *leafLocal; 814 PetscInt cellHeight, cStart, cEnd, c, fStart, fEnd, f, numLeaves, l; 815 PetscMPIInt rank; 816 817 PetscFunctionBegin; 818 /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */ 819 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 820 PetscCall(DMGetPointSF(dm, &sfPoint)); 821 PetscCall(DMPlexGetVTKCellHeight(dmNew, &cellHeight)); 822 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 823 PetscCall(PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote)); 824 PetscCall(DMCreateLabel(dmNew, "vtk")); 825 PetscCall(DMGetLabel(dmNew, "vtk", &vtkLabel)); 826 if (createGhostLabel) { 827 PetscCall(DMCreateLabel(dmNew, "ghost")); 828 PetscCall(DMGetLabel(dmNew, "ghost", &ghostLabel)); 829 } 830 for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) { 831 for (; c < leafLocal[l] && c < cEnd; ++c) PetscCall(DMLabelSetValue(vtkLabel, c, 1)); 832 if (leafLocal[l] >= cEnd) break; 833 if (leafRemote[l].rank == rank) { 834 PetscCall(DMLabelSetValue(vtkLabel, c, 1)); 835 } else if (ghostLabel) PetscCall(DMLabelSetValue(ghostLabel, c, 2)); 836 } 837 for (; c < cEnd; ++c) PetscCall(DMLabelSetValue(vtkLabel, c, 1)); 838 if (ghostLabel) { 839 PetscCall(DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd)); 840 for (f = fStart; f < fEnd; ++f) { 841 PetscInt numCells; 842 843 PetscCall(DMPlexGetSupportSize(dmNew, f, &numCells)); 844 if (numCells < 2) { 845 PetscCall(DMLabelSetValue(ghostLabel, f, 1)); 846 } else { 847 const PetscInt *cells = NULL; 848 PetscInt vA, vB; 849 850 PetscCall(DMPlexGetSupport(dmNew, f, &cells)); 851 PetscCall(DMLabelGetValue(vtkLabel, cells[0], &vA)); 852 PetscCall(DMLabelGetValue(vtkLabel, cells[1], &vB)); 853 if (vA != 1 && vB != 1) PetscCall(DMLabelSetValue(ghostLabel, f, 1)); 854 } 855 } 856 } 857 PetscFunctionReturn(PETSC_SUCCESS); 858 } 859 860 static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew) 861 { 862 DM refTree; 863 PetscSection pSec; 864 PetscInt *parents, *childIDs; 865 866 PetscFunctionBegin; 867 PetscCall(DMPlexGetReferenceTree(dm, &refTree)); 868 PetscCall(DMPlexSetReferenceTree(dmNew, refTree)); 869 PetscCall(DMPlexGetTree(dm, &pSec, &parents, &childIDs, NULL, NULL)); 870 if (pSec) { 871 PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth; 872 PetscInt *childIDsShifted; 873 PetscSection pSecShifted; 874 875 PetscCall(PetscSectionGetChart(pSec, &pStart, &pEnd)); 876 PetscCall(DMPlexGetDepth(dm, &depth)); 877 pStartShifted = DMPlexShiftPoint_Internal(pStart, depth, depthShift); 878 pEndShifted = DMPlexShiftPoint_Internal(pEnd, depth, depthShift); 879 PetscCall(PetscMalloc2(pEndShifted - pStartShifted, &parentsShifted, pEndShifted - pStartShifted, &childIDsShifted)); 880 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dmNew), &pSecShifted)); 881 PetscCall(PetscSectionSetChart(pSecShifted, pStartShifted, pEndShifted)); 882 for (p = pStartShifted; p < pEndShifted; p++) { 883 /* start off assuming no children */ 884 PetscCall(PetscSectionSetDof(pSecShifted, p, 0)); 885 } 886 for (p = pStart; p < pEnd; p++) { 887 PetscInt dof; 888 PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift); 889 890 PetscCall(PetscSectionGetDof(pSec, p, &dof)); 891 PetscCall(PetscSectionSetDof(pSecShifted, pNew, dof)); 892 } 893 PetscCall(PetscSectionSetUp(pSecShifted)); 894 for (p = pStart; p < pEnd; p++) { 895 PetscInt dof; 896 PetscInt pNew = DMPlexShiftPoint_Internal(p, depth, depthShift); 897 898 PetscCall(PetscSectionGetDof(pSec, p, &dof)); 899 if (dof) { 900 PetscInt off, offNew; 901 902 PetscCall(PetscSectionGetOffset(pSec, p, &off)); 903 PetscCall(PetscSectionGetOffset(pSecShifted, pNew, &offNew)); 904 parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off], depth, depthShift); 905 childIDsShifted[offNew] = childIDs[off]; 906 } 907 } 908 PetscCall(DMPlexSetTree(dmNew, pSecShifted, parentsShifted, childIDsShifted)); 909 PetscCall(PetscFree2(parentsShifted, childIDsShifted)); 910 PetscCall(PetscSectionDestroy(&pSecShifted)); 911 } 912 PetscFunctionReturn(PETSC_SUCCESS); 913 } 914 915 static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm) 916 { 917 PetscSF sf; 918 IS valueIS; 919 const PetscInt *values, *leaves; 920 PetscInt *depthShift; 921 PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c; 922 const PetscReal *maxCell, *Lstart, *L; 923 924 PetscFunctionBegin; 925 PetscCall(DMGetPointSF(dm, &sf)); 926 PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); 927 nleaves = PetscMax(0, nleaves); 928 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 929 /* Count ghost cells */ 930 PetscCall(DMLabelGetValueIS(label, &valueIS)); 931 PetscCall(ISGetLocalSize(valueIS, &numFS)); 932 PetscCall(ISGetIndices(valueIS, &values)); 933 Ng = 0; 934 for (fs = 0; fs < numFS; ++fs) { 935 IS faceIS; 936 const PetscInt *faces; 937 PetscInt numFaces, f, numBdFaces = 0; 938 939 PetscCall(DMLabelGetStratumIS(label, values[fs], &faceIS)); 940 PetscCall(ISGetLocalSize(faceIS, &numFaces)); 941 PetscCall(ISGetIndices(faceIS, &faces)); 942 for (f = 0; f < numFaces; ++f) { 943 PetscInt numChildren; 944 945 PetscCall(PetscFindInt(faces[f], nleaves, leaves, &loc)); 946 PetscCall(DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL)); 947 /* non-local and ancestors points don't get to register ghosts */ 948 if (loc >= 0 || numChildren) continue; 949 if (faces[f] >= fStart && faces[f] < fEnd) ++numBdFaces; 950 } 951 Ng += numBdFaces; 952 PetscCall(ISRestoreIndices(faceIS, &faces)); 953 PetscCall(ISDestroy(&faceIS)); 954 } 955 PetscCall(DMPlexGetDepth(dm, &depth)); 956 PetscCall(PetscMalloc1(2 * (depth + 1), &depthShift)); 957 for (d = 0; d <= depth; d++) { 958 PetscInt dEnd; 959 960 PetscCall(DMPlexGetDepthStratum(dm, d, NULL, &dEnd)); 961 depthShift[2 * d] = dEnd; 962 depthShift[2 * d + 1] = 0; 963 } 964 if (depth >= 0) depthShift[2 * depth + 1] = Ng; 965 PetscCall(DMPlexShiftPointSetUp_Internal(depth, depthShift)); 966 PetscCall(DMPlexShiftSizes_Internal(dm, depthShift, gdm)); 967 /* Step 3: Set cone/support sizes for new points */ 968 PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 969 for (c = cEnd; c < cEnd + Ng; ++c) PetscCall(DMPlexSetConeSize(gdm, c, 1)); 970 for (fs = 0; fs < numFS; ++fs) { 971 IS faceIS; 972 const PetscInt *faces; 973 PetscInt numFaces, f; 974 975 PetscCall(DMLabelGetStratumIS(label, values[fs], &faceIS)); 976 PetscCall(ISGetLocalSize(faceIS, &numFaces)); 977 PetscCall(ISGetIndices(faceIS, &faces)); 978 for (f = 0; f < numFaces; ++f) { 979 PetscInt size, numChildren; 980 981 PetscCall(PetscFindInt(faces[f], nleaves, leaves, &loc)); 982 PetscCall(DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL)); 983 if (loc >= 0 || numChildren) continue; 984 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 985 PetscCall(DMPlexGetSupportSize(dm, faces[f], &size)); 986 PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %" PetscInt_FMT " with %" PetscInt_FMT " support cells", faces[f], size); 987 PetscCall(DMPlexSetSupportSize(gdm, faces[f] + Ng, 2)); 988 } 989 PetscCall(ISRestoreIndices(faceIS, &faces)); 990 PetscCall(ISDestroy(&faceIS)); 991 } 992 /* Step 4: Setup ghosted DM */ 993 PetscCall(DMSetUp(gdm)); 994 PetscCall(DMPlexShiftPoints_Internal(dm, depthShift, gdm)); 995 /* Step 6: Set cones and supports for new points */ 996 ghostCell = cEnd; 997 for (fs = 0; fs < numFS; ++fs) { 998 IS faceIS; 999 const PetscInt *faces; 1000 PetscInt numFaces, f; 1001 1002 PetscCall(DMLabelGetStratumIS(label, values[fs], &faceIS)); 1003 PetscCall(ISGetLocalSize(faceIS, &numFaces)); 1004 PetscCall(ISGetIndices(faceIS, &faces)); 1005 for (f = 0; f < numFaces; ++f) { 1006 PetscInt newFace = faces[f] + Ng, numChildren; 1007 1008 PetscCall(PetscFindInt(faces[f], nleaves, leaves, &loc)); 1009 PetscCall(DMPlexGetTreeChildren(dm, faces[f], &numChildren, NULL)); 1010 if (loc >= 0 || numChildren) continue; 1011 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 1012 PetscCall(DMPlexSetCone(gdm, ghostCell, &newFace)); 1013 PetscCall(DMPlexInsertSupport(gdm, newFace, 1, ghostCell)); 1014 ++ghostCell; 1015 } 1016 PetscCall(ISRestoreIndices(faceIS, &faces)); 1017 PetscCall(ISDestroy(&faceIS)); 1018 } 1019 PetscCall(ISRestoreIndices(valueIS, &values)); 1020 PetscCall(ISDestroy(&valueIS)); 1021 PetscCall(DMPlexShiftCoordinates_Internal(dm, depthShift, gdm)); 1022 PetscCall(DMPlexShiftSF_Internal(dm, depthShift, gdm)); 1023 PetscCall(DMPlexShiftLabels_Internal(dm, depthShift, gdm)); 1024 PetscCall(DMPlexCreateVTKLabel_Internal(dm, PETSC_TRUE, gdm)); 1025 PetscCall(DMPlexShiftTree_Internal(dm, depthShift, gdm)); 1026 PetscCall(PetscFree(depthShift)); 1027 for (c = cEnd; c < cEnd + Ng; ++c) PetscCall(DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST)); 1028 /* Step 7: Periodicity */ 1029 PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); 1030 PetscCall(DMSetPeriodicity(gdm, maxCell, Lstart, L)); 1031 if (numGhostCells) *numGhostCells = Ng; 1032 PetscFunctionReturn(PETSC_SUCCESS); 1033 } 1034 1035 /*@ 1036 DMPlexConstructGhostCells - Construct ghost cells which connect to every boundary face 1037 1038 Collective 1039 1040 Input Parameters: 1041 + dm - The original `DM` 1042 - labelName - The label specifying the boundary faces, or "Face Sets" if this is `NULL` 1043 1044 Output Parameters: 1045 + numGhostCells - The number of ghost cells added to the `DM` 1046 - dmGhosted - The new `DM` 1047 1048 Level: developer 1049 1050 Note: 1051 If no label exists of that name, one will be created marking all boundary faces 1052 1053 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()` 1054 @*/ 1055 PetscErrorCode DMPlexConstructGhostCells(DM dm, const char labelName[], PetscInt *numGhostCells, DM *dmGhosted) 1056 { 1057 DM gdm; 1058 DMLabel label; 1059 const char *name = labelName ? labelName : "Face Sets"; 1060 PetscInt dim, Ng = 0; 1061 PetscBool useCone, useClosure; 1062 1063 PetscFunctionBegin; 1064 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1065 if (numGhostCells) PetscAssertPointer(numGhostCells, 3); 1066 PetscAssertPointer(dmGhosted, 4); 1067 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &gdm)); 1068 PetscCall(DMSetType(gdm, DMPLEX)); 1069 PetscCall(DMGetDimension(dm, &dim)); 1070 PetscCall(DMSetDimension(gdm, dim)); 1071 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 1072 PetscCall(DMSetBasicAdjacency(gdm, useCone, useClosure)); 1073 PetscCall(DMGetLabel(dm, name, &label)); 1074 if (!label) { 1075 /* Get label for boundary faces */ 1076 PetscCall(DMCreateLabel(dm, name)); 1077 PetscCall(DMGetLabel(dm, name, &label)); 1078 PetscCall(DMPlexMarkBoundaryFaces(dm, 1, label)); 1079 } 1080 PetscCall(DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm)); 1081 PetscCall(DMCopyDisc(dm, gdm)); 1082 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, gdm)); 1083 gdm->setfromoptionscalled = dm->setfromoptionscalled; 1084 if (numGhostCells) *numGhostCells = Ng; 1085 *dmGhosted = gdm; 1086 PetscFunctionReturn(PETSC_SUCCESS); 1087 } 1088 1089 static PetscErrorCode DivideCells_Private(DM dm, DMLabel label, DMPlexPointQueue queue) 1090 { 1091 PetscInt dim, d, shift = 100, *pStart, *pEnd; 1092 1093 PetscFunctionBegin; 1094 PetscCall(DMGetDimension(dm, &dim)); 1095 PetscCall(PetscMalloc2(dim, &pStart, dim, &pEnd)); 1096 for (d = 0; d < dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d])); 1097 while (!DMPlexPointQueueEmpty(queue)) { 1098 PetscInt cell = -1; 1099 PetscInt *closure = NULL; 1100 PetscInt closureSize, cl, cval; 1101 1102 PetscCall(DMPlexPointQueueDequeue(queue, &cell)); 1103 PetscCall(DMLabelGetValue(label, cell, &cval)); 1104 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 1105 /* Mark points in the cell closure that touch the fault */ 1106 for (d = 0; d < dim; ++d) { 1107 for (cl = 0; cl < closureSize * 2; cl += 2) { 1108 const PetscInt clp = closure[cl]; 1109 PetscInt clval; 1110 1111 if ((clp < pStart[d]) || (clp >= pEnd[d])) continue; 1112 PetscCall(DMLabelGetValue(label, clp, &clval)); 1113 if (clval == -1) { 1114 const PetscInt *cone; 1115 PetscInt coneSize, c; 1116 1117 /* If a cone point touches the fault, then this point touches the fault */ 1118 PetscCall(DMPlexGetCone(dm, clp, &cone)); 1119 PetscCall(DMPlexGetConeSize(dm, clp, &coneSize)); 1120 for (c = 0; c < coneSize; ++c) { 1121 PetscInt cpval; 1122 1123 PetscCall(DMLabelGetValue(label, cone[c], &cpval)); 1124 if (cpval != -1) { 1125 PetscInt dep; 1126 1127 PetscCall(DMPlexGetPointDepth(dm, clp, &dep)); 1128 clval = cval < 0 ? -(shift + dep) : shift + dep; 1129 PetscCall(DMLabelSetValue(label, clp, clval)); 1130 break; 1131 } 1132 } 1133 } 1134 /* Mark neighbor cells through marked faces (these cells must also touch the fault) */ 1135 if (d == dim - 1 && clval != -1) { 1136 const PetscInt *support; 1137 PetscInt supportSize, s, nval; 1138 1139 PetscCall(DMPlexGetSupport(dm, clp, &support)); 1140 PetscCall(DMPlexGetSupportSize(dm, clp, &supportSize)); 1141 for (s = 0; s < supportSize; ++s) { 1142 PetscCall(DMLabelGetValue(label, support[s], &nval)); 1143 if (nval == -1) { 1144 PetscCall(DMLabelSetValue(label, support[s], clval < 0 ? clval - 1 : clval + 1)); 1145 PetscCall(DMPlexPointQueueEnqueue(queue, support[s])); 1146 } 1147 } 1148 } 1149 } 1150 } 1151 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 1152 } 1153 PetscCall(PetscFree2(pStart, pEnd)); 1154 PetscFunctionReturn(PETSC_SUCCESS); 1155 } 1156 1157 typedef struct { 1158 DM dm; 1159 DMPlexPointQueue queue; 1160 } PointDivision; 1161 1162 static PetscErrorCode divideCell(DMLabel label, PetscInt p, PetscInt val, void *ctx) 1163 { 1164 PointDivision *div = (PointDivision *)ctx; 1165 PetscInt cval = val < 0 ? val - 1 : val + 1; 1166 const PetscInt *support; 1167 PetscInt supportSize, s; 1168 1169 PetscFunctionBegin; 1170 PetscCall(DMPlexGetSupport(div->dm, p, &support)); 1171 PetscCall(DMPlexGetSupportSize(div->dm, p, &supportSize)); 1172 for (s = 0; s < supportSize; ++s) { 1173 PetscCall(DMLabelSetValue(label, support[s], cval)); 1174 PetscCall(DMPlexPointQueueEnqueue(div->queue, support[s])); 1175 } 1176 PetscFunctionReturn(PETSC_SUCCESS); 1177 } 1178 1179 /* Mark cells by label propagation */ 1180 static PetscErrorCode DMPlexLabelFaultHalo(DM dm, DMLabel faultLabel) 1181 { 1182 DMPlexPointQueue queue = NULL; 1183 PointDivision div; 1184 PetscSF pointSF; 1185 IS pointIS; 1186 const PetscInt *points; 1187 PetscBool empty; 1188 PetscInt dim, shift = 100, n, i; 1189 1190 PetscFunctionBegin; 1191 PetscCall(DMGetDimension(dm, &dim)); 1192 PetscCall(DMPlexPointQueueCreate(1024, &queue)); 1193 div.dm = dm; 1194 div.queue = queue; 1195 /* Enqueue cells on fault */ 1196 PetscCall(DMLabelGetStratumIS(faultLabel, shift + dim, &pointIS)); 1197 if (pointIS) { 1198 PetscCall(ISGetLocalSize(pointIS, &n)); 1199 PetscCall(ISGetIndices(pointIS, &points)); 1200 for (i = 0; i < n; ++i) PetscCall(DMPlexPointQueueEnqueue(queue, points[i])); 1201 PetscCall(ISRestoreIndices(pointIS, &points)); 1202 PetscCall(ISDestroy(&pointIS)); 1203 } 1204 PetscCall(DMLabelGetStratumIS(faultLabel, -(shift + dim), &pointIS)); 1205 if (pointIS) { 1206 PetscCall(ISGetLocalSize(pointIS, &n)); 1207 PetscCall(ISGetIndices(pointIS, &points)); 1208 for (i = 0; i < n; ++i) PetscCall(DMPlexPointQueueEnqueue(queue, points[i])); 1209 PetscCall(ISRestoreIndices(pointIS, &points)); 1210 PetscCall(ISDestroy(&pointIS)); 1211 } 1212 1213 PetscCall(DMGetPointSF(dm, &pointSF)); 1214 PetscCall(DMLabelPropagateBegin(faultLabel, pointSF)); 1215 /* While edge queue is not empty: */ 1216 PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty)); 1217 while (!empty) { 1218 PetscCall(DivideCells_Private(dm, faultLabel, queue)); 1219 PetscCall(DMLabelPropagatePush(faultLabel, pointSF, divideCell, &div)); 1220 PetscCall(DMPlexPointQueueEmptyCollective((PetscObject)dm, queue, &empty)); 1221 } 1222 PetscCall(DMLabelPropagateEnd(faultLabel, pointSF)); 1223 PetscCall(DMPlexPointQueueDestroy(&queue)); 1224 PetscFunctionReturn(PETSC_SUCCESS); 1225 } 1226 1227 /* 1228 We are adding three kinds of points here: 1229 Replicated: Copies of points which exist in the mesh, such as vertices identified across a fault 1230 Non-replicated: Points which exist on the fault, but are not replicated 1231 Ghost: These are shared fault faces which are not owned by this process. These do not produce hybrid cells and do not replicate 1232 Hybrid: Entirely new points, such as cohesive cells 1233 1234 When creating subsequent cohesive cells, we shift the old hybrid cells to the end of the numbering at 1235 each depth so that the new split/hybrid points can be inserted as a block. 1236 */ 1237 static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DMLabel splitLabel, DM sdm) 1238 { 1239 MPI_Comm comm; 1240 IS valueIS; 1241 PetscInt numSP = 0; /* The number of depths for which we have replicated points */ 1242 const PetscInt *values; /* List of depths for which we have replicated points */ 1243 IS *splitIS; 1244 IS *unsplitIS; 1245 IS ghostIS; 1246 PetscInt *numSplitPoints; /* The number of replicated points at each depth */ 1247 PetscInt *numUnsplitPoints; /* The number of non-replicated points at each depth which still give rise to hybrid points */ 1248 PetscInt *numHybridPoints; /* The number of new hybrid points at each depth */ 1249 PetscInt *numHybridPointsOld; /* The number of existing hybrid points at each depth */ 1250 PetscInt numGhostPoints; /* The number of unowned, shared fault faces */ 1251 const PetscInt **splitPoints; /* Replicated points for each depth */ 1252 const PetscInt **unsplitPoints; /* Non-replicated points for each depth */ 1253 const PetscInt *ghostPoints; /* Ghost fault faces */ 1254 PetscSection coordSection; 1255 Vec coordinates; 1256 PetscScalar *coords; 1257 PetscInt *depthMax; /* The first hybrid point at each depth in the original mesh */ 1258 PetscInt *depthEnd; /* The point limit at each depth in the original mesh */ 1259 PetscInt *depthShift; /* Number of replicated+hybrid points at each depth */ 1260 PetscInt *pMaxNew; /* The first replicated point at each depth in the new mesh, hybrids come after this */ 1261 PetscInt *coneNew, *coneONew, *supportNew; 1262 PetscInt shift = 100, shift2 = 200, depth = 0, dep, dim, d, sp, maxConeSize, maxSupportSize, maxConeSizeNew, maxSupportSizeNew, numLabels, vStart, vEnd, pEnd, p, v; 1263 1264 PetscFunctionBegin; 1265 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1266 PetscCall(DMGetDimension(dm, &dim)); 1267 PetscCall(DMPlexGetDepth(dm, &depth)); 1268 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1269 /* We do not want this label automatically computed, instead we compute it here */ 1270 PetscCall(DMCreateLabel(sdm, "celltype")); 1271 /* Count split points and add cohesive cells */ 1272 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 1273 PetscCall(PetscMalloc5(depth + 1, &depthMax, depth + 1, &depthEnd, 2 * (depth + 1), &depthShift, depth + 1, &pMaxNew, depth + 1, &numHybridPointsOld)); 1274 PetscCall(PetscMalloc7(depth + 1, &splitIS, depth + 1, &unsplitIS, depth + 1, &numSplitPoints, depth + 1, &numUnsplitPoints, depth + 1, &numHybridPoints, depth + 1, &splitPoints, depth + 1, &unsplitPoints)); 1275 for (d = 0; d <= depth; ++d) { 1276 PetscCall(DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d])); 1277 PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, d, &depthMax[d], NULL)); 1278 depthEnd[d] = pMaxNew[d]; 1279 depthMax[d] = depthMax[d] < 0 ? depthEnd[d] : depthMax[d]; 1280 numSplitPoints[d] = 0; 1281 numUnsplitPoints[d] = 0; 1282 numHybridPoints[d] = 0; 1283 numHybridPointsOld[d] = depthMax[d] < 0 ? 0 : depthEnd[d] - depthMax[d]; 1284 splitPoints[d] = NULL; 1285 unsplitPoints[d] = NULL; 1286 splitIS[d] = NULL; 1287 unsplitIS[d] = NULL; 1288 /* we are shifting the existing hybrid points with the stratum behind them, so 1289 * the split comes at the end of the normal points, i.e., at depthMax[d] */ 1290 depthShift[2 * d] = depthMax[d]; 1291 depthShift[2 * d + 1] = 0; 1292 } 1293 numGhostPoints = 0; 1294 ghostPoints = NULL; 1295 if (label) { 1296 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1297 PetscCall(ISGetLocalSize(valueIS, &numSP)); 1298 PetscCall(ISGetIndices(valueIS, &values)); 1299 } 1300 for (sp = 0; sp < numSP; ++sp) { 1301 const PetscInt dep = values[sp]; 1302 1303 if ((dep < 0) || (dep > depth)) continue; 1304 PetscCall(DMLabelGetStratumIS(label, dep, &splitIS[dep])); 1305 if (splitIS[dep]) { 1306 PetscCall(ISGetLocalSize(splitIS[dep], &numSplitPoints[dep])); 1307 PetscCall(ISGetIndices(splitIS[dep], &splitPoints[dep])); 1308 } 1309 PetscCall(DMLabelGetStratumIS(label, shift2 + dep, &unsplitIS[dep])); 1310 if (unsplitIS[dep]) { 1311 PetscCall(ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep])); 1312 PetscCall(ISGetIndices(unsplitIS[dep], &unsplitPoints[dep])); 1313 } 1314 } 1315 PetscCall(DMLabelGetStratumIS(label, shift2 + dim - 1, &ghostIS)); 1316 if (ghostIS) { 1317 PetscCall(ISGetLocalSize(ghostIS, &numGhostPoints)); 1318 PetscCall(ISGetIndices(ghostIS, &ghostPoints)); 1319 } 1320 /* Calculate number of hybrid points */ 1321 for (d = 1; d <= depth; ++d) numHybridPoints[d] = numSplitPoints[d - 1] + numUnsplitPoints[d - 1]; /* There is a hybrid cell/face/edge for every split face/edge/vertex */ 1322 for (d = 0; d <= depth; ++d) depthShift[2 * d + 1] = numSplitPoints[d] + numHybridPoints[d]; 1323 PetscCall(DMPlexShiftPointSetUp_Internal(depth, depthShift)); 1324 /* the end of the points in this stratum that come before the new points: 1325 * shifting pMaxNew[d] gets the new start of the next stratum, then count back the old hybrid points and the newly 1326 * added points */ 1327 for (d = 0; d <= depth; ++d) pMaxNew[d] = DMPlexShiftPoint_Internal(pMaxNew[d], depth, depthShift) - (numHybridPointsOld[d] + numSplitPoints[d] + numHybridPoints[d]); 1328 PetscCall(DMPlexShiftSizes_Internal(dm, depthShift, sdm)); 1329 /* Step 3: Set cone/support sizes for new points */ 1330 for (dep = 0; dep <= depth; ++dep) { 1331 for (p = 0; p < numSplitPoints[dep]; ++p) { 1332 const PetscInt oldp = splitPoints[dep][p]; 1333 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/; 1334 const PetscInt splitp = p + pMaxNew[dep]; 1335 const PetscInt *support; 1336 DMPolytopeType ct; 1337 PetscInt coneSize, supportSize, qf, qn, qp, e; 1338 1339 PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize)); 1340 PetscCall(DMPlexSetConeSize(sdm, splitp, coneSize)); 1341 PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1342 PetscCall(DMPlexSetSupportSize(sdm, splitp, supportSize)); 1343 PetscCall(DMPlexGetCellType(dm, oldp, &ct)); 1344 PetscCall(DMPlexSetCellType(sdm, splitp, ct)); 1345 if (dep == depth - 1) { 1346 const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1347 1348 /* Add cohesive cells, they are prisms */ 1349 PetscCall(DMPlexSetConeSize(sdm, hybcell, 2 + coneSize)); 1350 switch (coneSize) { 1351 case 2: 1352 PetscCall(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR)); 1353 break; 1354 case 3: 1355 PetscCall(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR)); 1356 break; 1357 case 4: 1358 PetscCall(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_QUAD_PRISM_TENSOR)); 1359 break; 1360 } 1361 /* Shared fault faces with only one support cell now have two with the cohesive cell */ 1362 /* TODO Check thaat oldp has rootdegree == 1 */ 1363 if (supportSize == 1) { 1364 const PetscInt *support; 1365 PetscInt val; 1366 1367 PetscCall(DMPlexGetSupport(dm, oldp, &support)); 1368 PetscCall(DMLabelGetValue(label, support[0], &val)); 1369 if (val < 0) PetscCall(DMPlexSetSupportSize(sdm, splitp, 2)); 1370 else PetscCall(DMPlexSetSupportSize(sdm, newp, 2)); 1371 } 1372 } else if (dep == 0) { 1373 const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1374 1375 PetscCall(DMPlexGetSupport(dm, oldp, &support)); 1376 for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) { 1377 PetscInt val; 1378 1379 PetscCall(DMLabelGetValue(label, support[e], &val)); 1380 if (val == 1) ++qf; 1381 if ((val == 1) || (val == (shift + 1))) ++qn; 1382 if ((val == 1) || (val == -(shift + 1))) ++qp; 1383 } 1384 /* Split old vertex: Edges into original vertex and new cohesive edge */ 1385 PetscCall(DMPlexSetSupportSize(sdm, newp, qn + 1)); 1386 /* Split new vertex: Edges into split vertex and new cohesive edge */ 1387 PetscCall(DMPlexSetSupportSize(sdm, splitp, qp + 1)); 1388 /* Add hybrid edge */ 1389 PetscCall(DMPlexSetConeSize(sdm, hybedge, 2)); 1390 PetscCall(DMPlexSetSupportSize(sdm, hybedge, qf)); 1391 PetscCall(DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR)); 1392 } else if (dep == dim - 2) { 1393 const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1394 1395 PetscCall(DMPlexGetSupport(dm, oldp, &support)); 1396 for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) { 1397 PetscInt val; 1398 1399 PetscCall(DMLabelGetValue(label, support[e], &val)); 1400 if (val == dim - 1) ++qf; 1401 if ((val == dim - 1) || (val == (shift + dim - 1))) ++qn; 1402 if ((val == dim - 1) || (val == -(shift + dim - 1))) ++qp; 1403 } 1404 /* Split old edge: Faces into original edge and cohesive face (positive side?) */ 1405 PetscCall(DMPlexSetSupportSize(sdm, newp, qn + 1)); 1406 /* Split new edge: Faces into split edge and cohesive face (negative side?) */ 1407 PetscCall(DMPlexSetSupportSize(sdm, splitp, qp + 1)); 1408 /* Add hybrid face */ 1409 PetscCall(DMPlexSetConeSize(sdm, hybface, 4)); 1410 PetscCall(DMPlexSetSupportSize(sdm, hybface, qf)); 1411 PetscCall(DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR)); 1412 } 1413 } 1414 } 1415 for (dep = 0; dep <= depth; ++dep) { 1416 for (p = 0; p < numUnsplitPoints[dep]; ++p) { 1417 const PetscInt oldp = unsplitPoints[dep][p]; 1418 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/; 1419 const PetscInt *support; 1420 PetscInt coneSize, supportSize, qf, e, s; 1421 1422 PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize)); 1423 PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1424 PetscCall(DMPlexGetSupport(dm, oldp, &support)); 1425 if (dep == 0) { 1426 const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep]; 1427 1428 /* Unsplit vertex: Edges into original vertex, split edges, and new cohesive edge twice */ 1429 for (s = 0, qf = 0; s < supportSize; ++s, ++qf) { 1430 PetscCall(PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e)); 1431 if (e >= 0) ++qf; 1432 } 1433 PetscCall(DMPlexSetSupportSize(sdm, newp, qf + 2)); 1434 /* Add hybrid edge */ 1435 PetscCall(DMPlexSetConeSize(sdm, hybedge, 2)); 1436 for (e = 0, qf = 0; e < supportSize; ++e) { 1437 PetscInt val; 1438 1439 PetscCall(DMLabelGetValue(label, support[e], &val)); 1440 /* Split and unsplit edges produce hybrid faces */ 1441 if (val == 1) ++qf; 1442 if (val == (shift2 + 1)) ++qf; 1443 } 1444 PetscCall(DMPlexSetSupportSize(sdm, hybedge, qf)); 1445 PetscCall(DMPlexSetCellType(sdm, hybedge, DM_POLYTOPE_POINT_PRISM_TENSOR)); 1446 } else if (dep == dim - 2) { 1447 const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep]; 1448 PetscInt val; 1449 1450 for (e = 0, qf = 0; e < supportSize; ++e) { 1451 PetscCall(DMLabelGetValue(label, support[e], &val)); 1452 if (val == dim - 1) qf += 2; 1453 else ++qf; 1454 } 1455 /* Unsplit edge: Faces into original edge, split face, and cohesive face twice */ 1456 PetscCall(DMPlexSetSupportSize(sdm, newp, qf + 2)); 1457 /* Add hybrid face */ 1458 for (e = 0, qf = 0; e < supportSize; ++e) { 1459 PetscCall(DMLabelGetValue(label, support[e], &val)); 1460 if (val == dim - 1) ++qf; 1461 } 1462 PetscCall(DMPlexSetConeSize(sdm, hybface, 4)); 1463 PetscCall(DMPlexSetSupportSize(sdm, hybface, qf)); 1464 PetscCall(DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR)); 1465 } 1466 } 1467 } 1468 /* Step 4: Setup split DM */ 1469 PetscCall(DMSetUp(sdm)); 1470 PetscCall(DMPlexShiftPoints_Internal(dm, depthShift, sdm)); 1471 PetscCall(DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew)); 1472 PetscCall(PetscMalloc3(PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneNew, PetscMax(maxConeSize, maxConeSizeNew) * 3, &coneONew, PetscMax(maxSupportSize, maxSupportSizeNew), &supportNew)); 1473 /* Step 6: Set cones and supports for new points */ 1474 for (dep = 0; dep <= depth; ++dep) { 1475 for (p = 0; p < numSplitPoints[dep]; ++p) { 1476 const PetscInt oldp = splitPoints[dep][p]; 1477 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/; 1478 const PetscInt splitp = p + pMaxNew[dep]; 1479 const PetscInt *cone, *support, *ornt; 1480 DMPolytopeType ct; 1481 PetscInt coneSize, supportSize, q, qf, qn, qp, v, e, s; 1482 1483 PetscCall(DMPlexGetCellType(dm, oldp, &ct)); 1484 PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize)); 1485 PetscCall(DMPlexGetCone(dm, oldp, &cone)); 1486 PetscCall(DMPlexGetConeOrientation(dm, oldp, &ornt)); 1487 PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1488 PetscCall(DMPlexGetSupport(dm, oldp, &support)); 1489 if (dep == depth - 1) { 1490 PetscBool hasUnsplit = PETSC_FALSE; 1491 const PetscInt hybcell = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1492 const PetscInt *supportF; 1493 1494 coneONew[0] = coneONew[1] = -1000; 1495 /* Split face: copy in old face to new face to start */ 1496 PetscCall(DMPlexGetSupport(sdm, newp, &supportF)); 1497 PetscCall(DMPlexSetSupport(sdm, splitp, supportF)); 1498 /* Split old face: old vertices/edges in cone so no change */ 1499 /* Split new face: new vertices/edges in cone */ 1500 for (q = 0; q < coneSize; ++q) { 1501 PetscCall(PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v)); 1502 if (v < 0) { 1503 PetscCall(PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v)); 1504 PetscCheck(v >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %" PetscInt_FMT " in split or unsplit points of depth %" PetscInt_FMT, cone[q], dep - 1); 1505 coneNew[2 + q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/; 1506 hasUnsplit = PETSC_TRUE; 1507 } else { 1508 coneNew[2 + q] = v + pMaxNew[dep - 1]; 1509 if (dep > 1) { 1510 const PetscInt *econe; 1511 PetscInt econeSize, r, vs, vu; 1512 1513 PetscCall(DMPlexGetConeSize(dm, cone[q], &econeSize)); 1514 PetscCall(DMPlexGetCone(dm, cone[q], &econe)); 1515 for (r = 0; r < econeSize; ++r) { 1516 PetscCall(PetscFindInt(econe[r], numSplitPoints[dep - 2], splitPoints[dep - 2], &vs)); 1517 PetscCall(PetscFindInt(econe[r], numUnsplitPoints[dep - 2], unsplitPoints[dep - 2], &vu)); 1518 if (vs >= 0) continue; 1519 PetscCheck(vu >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %" PetscInt_FMT " in split or unsplit points of depth %" PetscInt_FMT, econe[r], dep - 2); 1520 hasUnsplit = PETSC_TRUE; 1521 } 1522 } 1523 } 1524 } 1525 PetscCall(DMPlexSetCone(sdm, splitp, &coneNew[2])); 1526 PetscCall(DMPlexSetConeOrientation(sdm, splitp, ornt)); 1527 /* Face support */ 1528 PetscInt vals[2]; 1529 1530 PetscCall(DMLabelGetValue(label, support[0], &vals[0])); 1531 if (supportSize > 1) PetscCall(DMLabelGetValue(label, support[1], &vals[1])); 1532 else vals[1] = -vals[0]; 1533 PetscCheck(vals[0] * vals[1] < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid support labels %" PetscInt_FMT " %" PetscInt_FMT, vals[0], vals[1]); 1534 1535 for (s = 0; s < 2; ++s) { 1536 if (s >= supportSize) { 1537 if (vals[s] < 0) { 1538 /* Ghost old face: Replace negative side cell with cohesive cell */ 1539 PetscCall(DMPlexInsertSupport(sdm, newp, s, hybcell)); 1540 } else { 1541 /* Ghost new face: Replace positive side cell with cohesive cell */ 1542 PetscCall(DMPlexInsertSupport(sdm, splitp, s, hybcell)); 1543 } 1544 } else { 1545 if (vals[s] < 0) { 1546 /* Split old face: Replace negative side cell with cohesive cell */ 1547 PetscCall(DMPlexInsertSupport(sdm, newp, s, hybcell)); 1548 } else { 1549 /* Split new face: Replace positive side cell with cohesive cell */ 1550 PetscCall(DMPlexInsertSupport(sdm, splitp, s, hybcell)); 1551 } 1552 } 1553 } 1554 /* Get orientation for cohesive face using the positive side cell */ 1555 { 1556 const PetscInt *ncone, *nconeO; 1557 PetscInt nconeSize, nc, ocell; 1558 PetscBool flip = PETSC_FALSE; 1559 1560 if (supportSize > 1) { 1561 ocell = vals[0] < 0 ? support[1] : support[0]; 1562 } else { 1563 ocell = support[0]; 1564 flip = vals[0] < 0 ? PETSC_TRUE : PETSC_FALSE; 1565 } 1566 PetscCall(DMPlexGetConeSize(dm, ocell, &nconeSize)); 1567 PetscCall(DMPlexGetCone(dm, ocell, &ncone)); 1568 PetscCall(DMPlexGetConeOrientation(dm, ocell, &nconeO)); 1569 for (nc = 0; nc < nconeSize; ++nc) { 1570 if (ncone[nc] == oldp) { 1571 coneONew[0] = flip ? -(nconeO[nc] + 1) : nconeO[nc]; 1572 break; 1573 } 1574 } 1575 PetscCheck(nc < nconeSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate face %" PetscInt_FMT " in neighboring cell %" PetscInt_FMT, oldp, ocell); 1576 } 1577 /* Cohesive cell: Old and new split face, then new cohesive faces */ 1578 { 1579 const PetscInt No = DMPolytopeTypeGetNumArrangements(ct) / 2; 1580 PetscCheck(coneONew[0] >= -No && coneONew[0] < No, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid %s orientation %" PetscInt_FMT, DMPolytopeTypes[ct], coneONew[0]); 1581 } 1582 const PetscInt *arr = DMPolytopeTypeGetArrangement(ct, coneONew[0]); 1583 1584 coneNew[0] = newp; /* Extracted negative side orientation above */ 1585 coneNew[1] = splitp; 1586 coneONew[1] = coneONew[0]; 1587 for (q = 0; q < coneSize; ++q) { 1588 /* Hybrid faces must follow order from oriented end face */ 1589 const PetscInt qa = arr[q * 2 + 0]; 1590 const PetscInt qo = arr[q * 2 + 1]; 1591 DMPolytopeType ft = dep == 2 ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT; 1592 1593 PetscCall(PetscFindInt(cone[qa], numSplitPoints[dep - 1], splitPoints[dep - 1], &v)); 1594 if (v < 0) { 1595 PetscCall(PetscFindInt(cone[qa], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v)); 1596 coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1]; 1597 } else { 1598 coneNew[2 + q] = v + pMaxNew[dep] + numSplitPoints[dep]; 1599 } 1600 coneONew[2 + q] = DMPolytopeTypeComposeOrientation(ft, qo, ornt[qa]); 1601 } 1602 PetscCall(DMPlexSetCone(sdm, hybcell, coneNew)); 1603 PetscCall(DMPlexSetConeOrientation(sdm, hybcell, coneONew)); 1604 /* Label the hybrid cells on the boundary of the split */ 1605 if (hasUnsplit) PetscCall(DMLabelSetValue(label, -hybcell, dim)); 1606 } else if (dep == 0) { 1607 const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1608 1609 /* Split old vertex: Edges in old split faces and new cohesive edge */ 1610 for (e = 0, qn = 0; e < supportSize; ++e) { 1611 PetscInt val; 1612 1613 PetscCall(DMLabelGetValue(label, support[e], &val)); 1614 if ((val == 1) || (val == (shift + 1))) supportNew[qn++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/; 1615 } 1616 supportNew[qn] = hybedge; 1617 PetscCall(DMPlexSetSupport(sdm, newp, supportNew)); 1618 /* Split new vertex: Edges in new split faces and new cohesive edge */ 1619 for (e = 0, qp = 0; e < supportSize; ++e) { 1620 PetscInt val, edge; 1621 1622 PetscCall(DMLabelGetValue(label, support[e], &val)); 1623 if (val == 1) { 1624 PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge)); 1625 PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a split edge", support[e]); 1626 supportNew[qp++] = edge + pMaxNew[dep + 1]; 1627 } else if (val == -(shift + 1)) { 1628 supportNew[qp++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/; 1629 } 1630 } 1631 supportNew[qp] = hybedge; 1632 PetscCall(DMPlexSetSupport(sdm, splitp, supportNew)); 1633 /* Hybrid edge: Old and new split vertex */ 1634 coneNew[0] = newp; 1635 coneNew[1] = splitp; 1636 PetscCall(DMPlexSetCone(sdm, hybedge, coneNew)); 1637 for (e = 0, qf = 0; e < supportSize; ++e) { 1638 PetscInt val, edge; 1639 1640 PetscCall(DMLabelGetValue(label, support[e], &val)); 1641 if (val == 1) { 1642 PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge)); 1643 PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a split edge", support[e]); 1644 supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2]; 1645 } 1646 } 1647 PetscCall(DMPlexSetSupport(sdm, hybedge, supportNew)); 1648 } else if (dep == dim - 2) { 1649 const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1650 1651 /* Split old edge: old vertices in cone so no change */ 1652 /* Split new edge: new vertices in cone */ 1653 for (q = 0; q < coneSize; ++q) { 1654 PetscCall(PetscFindInt(cone[q], numSplitPoints[dep - 1], splitPoints[dep - 1], &v)); 1655 if (v < 0) { 1656 PetscCall(PetscFindInt(cone[q], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v)); 1657 PetscCheck(v >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %" PetscInt_FMT " in split or unsplit points of depth %" PetscInt_FMT, cone[q], dep - 1); 1658 coneNew[q] = DMPlexShiftPoint_Internal(cone[q], depth, depthShift) /*cone[q] + depthOffset[dep-1]*/; 1659 } else { 1660 coneNew[q] = v + pMaxNew[dep - 1]; 1661 } 1662 } 1663 PetscCall(DMPlexSetCone(sdm, splitp, coneNew)); 1664 /* Split old edge: Faces in positive side cells and old split faces */ 1665 for (e = 0, q = 0; e < supportSize; ++e) { 1666 PetscInt val; 1667 1668 PetscCall(DMLabelGetValue(label, support[e], &val)); 1669 if (val == dim - 1) { 1670 supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/; 1671 } else if (val == (shift + dim - 1)) { 1672 supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/; 1673 } 1674 } 1675 supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1676 PetscCall(DMPlexSetSupport(sdm, newp, supportNew)); 1677 /* Split new edge: Faces in negative side cells and new split faces */ 1678 for (e = 0, q = 0; e < supportSize; ++e) { 1679 PetscInt val, face; 1680 1681 PetscCall(DMLabelGetValue(label, support[e], &val)); 1682 if (val == dim - 1) { 1683 PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face)); 1684 PetscCheck(face >= 0, comm, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " is not a split face", support[e]); 1685 supportNew[q++] = face + pMaxNew[dep + 1]; 1686 } else if (val == -(shift + dim - 1)) { 1687 supportNew[q++] = DMPlexShiftPoint_Internal(support[e], depth, depthShift) /*support[e] + depthOffset[dep+1]*/; 1688 } 1689 } 1690 supportNew[q++] = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1]; 1691 PetscCall(DMPlexSetSupport(sdm, splitp, supportNew)); 1692 /* Hybrid face */ 1693 coneNew[0] = newp; 1694 coneNew[1] = splitp; 1695 for (v = 0; v < coneSize; ++v) { 1696 PetscInt vertex; 1697 PetscCall(PetscFindInt(cone[v], numSplitPoints[dep - 1], splitPoints[dep - 1], &vertex)); 1698 if (vertex < 0) { 1699 PetscCall(PetscFindInt(cone[v], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &vertex)); 1700 PetscCheck(vertex >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate point %" PetscInt_FMT " in split or unsplit points of depth %" PetscInt_FMT, cone[v], dep - 1); 1701 coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1]; 1702 } else { 1703 coneNew[2 + v] = vertex + pMaxNew[dep] + numSplitPoints[dep]; 1704 } 1705 } 1706 PetscCall(DMPlexSetCone(sdm, hybface, coneNew)); 1707 for (e = 0, qf = 0; e < supportSize; ++e) { 1708 PetscInt val, face; 1709 1710 PetscCall(DMLabelGetValue(label, support[e], &val)); 1711 if (val == dim - 1) { 1712 PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &face)); 1713 PetscCheck(face >= 0, comm, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " is not a split face", support[e]); 1714 supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2]; 1715 } 1716 } 1717 PetscCall(DMPlexSetSupport(sdm, hybface, supportNew)); 1718 } 1719 } 1720 } 1721 for (dep = 0; dep <= depth; ++dep) { 1722 for (p = 0; p < numUnsplitPoints[dep]; ++p) { 1723 const PetscInt oldp = unsplitPoints[dep][p]; 1724 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*oldp + depthOffset[dep]*/; 1725 const PetscInt *cone, *support; 1726 PetscInt coneSize, supportSize, supportSizeNew, q, qf, e, f, s; 1727 1728 PetscCall(DMPlexGetConeSize(dm, oldp, &coneSize)); 1729 PetscCall(DMPlexGetCone(dm, oldp, &cone)); 1730 PetscCall(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1731 PetscCall(DMPlexGetSupport(dm, oldp, &support)); 1732 if (dep == 0) { 1733 const PetscInt hybedge = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep]; 1734 1735 /* Unsplit vertex */ 1736 PetscCall(DMPlexGetSupportSize(sdm, newp, &supportSizeNew)); 1737 for (s = 0, q = 0; s < supportSize; ++s) { 1738 supportNew[q++] = DMPlexShiftPoint_Internal(support[s], depth, depthShift) /*support[s] + depthOffset[dep+1]*/; 1739 PetscCall(PetscFindInt(support[s], numSplitPoints[dep + 1], splitPoints[dep + 1], &e)); 1740 if (e >= 0) supportNew[q++] = e + pMaxNew[dep + 1]; 1741 } 1742 supportNew[q++] = hybedge; 1743 supportNew[q++] = hybedge; 1744 PetscCheck(q == supportSizeNew, comm, PETSC_ERR_ARG_WRONG, "Support size %" PetscInt_FMT " != %" PetscInt_FMT " for vertex %" PetscInt_FMT, q, supportSizeNew, newp); 1745 PetscCall(DMPlexSetSupport(sdm, newp, supportNew)); 1746 /* Hybrid edge */ 1747 coneNew[0] = newp; 1748 coneNew[1] = newp; 1749 PetscCall(DMPlexSetCone(sdm, hybedge, coneNew)); 1750 for (e = 0, qf = 0; e < supportSize; ++e) { 1751 PetscInt val, edge; 1752 1753 PetscCall(DMLabelGetValue(label, support[e], &val)); 1754 if (val == 1) { 1755 PetscCall(PetscFindInt(support[e], numSplitPoints[dep + 1], splitPoints[dep + 1], &edge)); 1756 PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a split edge", support[e]); 1757 supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2]; 1758 } else if (val == (shift2 + 1)) { 1759 PetscCall(PetscFindInt(support[e], numUnsplitPoints[dep + 1], unsplitPoints[dep + 1], &edge)); 1760 PetscCheck(edge >= 0, comm, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " is not a unsplit edge", support[e]); 1761 supportNew[qf++] = edge + pMaxNew[dep + 2] + numSplitPoints[dep + 2] + numSplitPoints[dep + 1]; 1762 } 1763 } 1764 PetscCall(DMPlexSetSupport(sdm, hybedge, supportNew)); 1765 } else if (dep == dim - 2) { 1766 const PetscInt hybface = p + pMaxNew[dep + 1] + numSplitPoints[dep + 1] + numSplitPoints[dep]; 1767 1768 /* Unsplit edge: Faces into original edge, split face, and hybrid face twice */ 1769 for (f = 0, qf = 0; f < supportSize; ++f) { 1770 PetscInt val, face; 1771 1772 PetscCall(DMLabelGetValue(label, support[f], &val)); 1773 if (val == dim - 1) { 1774 PetscCall(PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face)); 1775 PetscCheck(face >= 0, comm, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " is not a split face", support[f]); 1776 supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/; 1777 supportNew[qf++] = face + pMaxNew[dep + 1]; 1778 } else { 1779 supportNew[qf++] = DMPlexShiftPoint_Internal(support[f], depth, depthShift) /*support[f] + depthOffset[dep+1]*/; 1780 } 1781 } 1782 supportNew[qf++] = hybface; 1783 supportNew[qf++] = hybface; 1784 PetscCall(DMPlexGetSupportSize(sdm, newp, &supportSizeNew)); 1785 PetscCheck(qf == supportSizeNew, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %" PetscInt_FMT " is %" PetscInt_FMT " != %" PetscInt_FMT, newp, qf, supportSizeNew); 1786 PetscCall(DMPlexSetSupport(sdm, newp, supportNew)); 1787 /* Add hybrid face */ 1788 coneNew[0] = newp; 1789 coneNew[1] = newp; 1790 PetscCall(PetscFindInt(cone[0], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v)); 1791 PetscCheck(v >= 0, comm, PETSC_ERR_ARG_WRONG, "Vertex %" PetscInt_FMT " is not an unsplit vertex", cone[0]); 1792 coneNew[2] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1]; 1793 PetscCall(PetscFindInt(cone[1], numUnsplitPoints[dep - 1], unsplitPoints[dep - 1], &v)); 1794 PetscCheck(v >= 0, comm, PETSC_ERR_ARG_WRONG, "Vertex %" PetscInt_FMT " is not an unsplit vertex", cone[1]); 1795 coneNew[3] = v + pMaxNew[dep] + numSplitPoints[dep] + numSplitPoints[dep - 1]; 1796 PetscCall(DMPlexSetCone(sdm, hybface, coneNew)); 1797 for (f = 0, qf = 0; f < supportSize; ++f) { 1798 PetscInt val, face; 1799 1800 PetscCall(DMLabelGetValue(label, support[f], &val)); 1801 if (val == dim - 1) { 1802 PetscCall(PetscFindInt(support[f], numSplitPoints[dep + 1], splitPoints[dep + 1], &face)); 1803 supportNew[qf++] = face + pMaxNew[dep + 2] + numSplitPoints[dep + 2]; 1804 } 1805 } 1806 PetscCall(DMPlexGetSupportSize(sdm, hybface, &supportSizeNew)); 1807 PetscCheck(qf == supportSizeNew, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %" PetscInt_FMT " is %" PetscInt_FMT " != %" PetscInt_FMT, hybface, qf, supportSizeNew); 1808 PetscCall(DMPlexSetSupport(sdm, hybface, supportNew)); 1809 } 1810 } 1811 } 1812 /* Step 6b: Replace split points in negative side cones */ 1813 for (sp = 0; sp < numSP; ++sp) { 1814 PetscInt dep = values[sp]; 1815 IS pIS; 1816 PetscInt numPoints; 1817 const PetscInt *points; 1818 1819 if (dep >= 0) continue; 1820 PetscCall(DMLabelGetStratumIS(label, dep, &pIS)); 1821 if (!pIS) continue; 1822 dep = -dep - shift; 1823 PetscCall(ISGetLocalSize(pIS, &numPoints)); 1824 PetscCall(ISGetIndices(pIS, &points)); 1825 for (p = 0; p < numPoints; ++p) { 1826 const PetscInt oldp = points[p]; 1827 const PetscInt newp = DMPlexShiftPoint_Internal(oldp, depth, depthShift) /*depthOffset[dep] + oldp*/; 1828 const PetscInt *cone; 1829 PetscInt coneSize, c; 1830 /* PetscBool replaced = PETSC_FALSE; */ 1831 1832 /* Negative edge: replace split vertex */ 1833 /* Negative cell: replace split face */ 1834 PetscCall(DMPlexGetConeSize(sdm, newp, &coneSize)); 1835 PetscCall(DMPlexGetCone(sdm, newp, &cone)); 1836 for (c = 0; c < coneSize; ++c) { 1837 const PetscInt coldp = DMPlexShiftPointInverse_Internal(cone[c], depth, depthShift); 1838 PetscInt csplitp, cp, val; 1839 1840 PetscCall(DMLabelGetValue(label, coldp, &val)); 1841 if (val == dep - 1) { 1842 PetscCall(PetscFindInt(coldp, numSplitPoints[dep - 1], splitPoints[dep - 1], &cp)); 1843 PetscCheck(cp >= 0, comm, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " is not a split point of dimension %" PetscInt_FMT, oldp, dep - 1); 1844 csplitp = pMaxNew[dep - 1] + cp; 1845 PetscCall(DMPlexInsertCone(sdm, newp, c, csplitp)); 1846 /* replaced = PETSC_TRUE; */ 1847 } 1848 } 1849 /* Cells with only a vertex or edge on the submesh have no replacement */ 1850 /* PetscCheck(replaced,comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); */ 1851 } 1852 PetscCall(ISRestoreIndices(pIS, &points)); 1853 PetscCall(ISDestroy(&pIS)); 1854 } 1855 PetscCall(DMPlexReorderCohesiveSupports(sdm)); 1856 /* Step 7: Coordinates */ 1857 PetscCall(DMPlexShiftCoordinates_Internal(dm, depthShift, sdm)); 1858 PetscCall(DMGetCoordinateSection(sdm, &coordSection)); 1859 PetscCall(DMGetCoordinatesLocal(sdm, &coordinates)); 1860 PetscCall(VecGetArray(coordinates, &coords)); 1861 for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) { 1862 const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[0][v], depth, depthShift) /*depthOffset[0] + splitPoints[0][v]*/; 1863 const PetscInt splitp = pMaxNew[0] + v; 1864 PetscInt dof, off, soff, d; 1865 1866 PetscCall(PetscSectionGetDof(coordSection, newp, &dof)); 1867 PetscCall(PetscSectionGetOffset(coordSection, newp, &off)); 1868 PetscCall(PetscSectionGetOffset(coordSection, splitp, &soff)); 1869 for (d = 0; d < dof; ++d) coords[soff + d] = coords[off + d]; 1870 } 1871 PetscCall(VecRestoreArray(coordinates, &coords)); 1872 /* Step 8: SF, if I can figure this out we can split the mesh in parallel */ 1873 PetscCall(DMPlexShiftSF_Internal(dm, depthShift, sdm)); 1874 /* TODO We need to associate the ghost points with the correct replica */ 1875 /* Step 9: Labels */ 1876 PetscCall(DMPlexShiftLabels_Internal(dm, depthShift, sdm)); 1877 PetscCall(DMPlexCreateVTKLabel_Internal(dm, PETSC_FALSE, sdm)); 1878 PetscCall(DMGetNumLabels(sdm, &numLabels)); 1879 for (dep = 0; dep <= depth; ++dep) { 1880 for (p = 0; p < numSplitPoints[dep]; ++p) { 1881 const PetscInt newp = DMPlexShiftPoint_Internal(splitPoints[dep][p], depth, depthShift) /*depthOffset[dep] + splitPoints[dep][p]*/; 1882 const PetscInt splitp = pMaxNew[dep] + p; 1883 PetscInt l; 1884 1885 if (splitLabel) { 1886 const PetscInt val = 100 + dep; 1887 1888 PetscCall(DMLabelSetValue(splitLabel, newp, val)); 1889 PetscCall(DMLabelSetValue(splitLabel, splitp, -val)); 1890 } 1891 for (l = 0; l < numLabels; ++l) { 1892 DMLabel mlabel; 1893 const char *lname; 1894 PetscInt val; 1895 PetscBool isDepth; 1896 1897 PetscCall(DMGetLabelName(sdm, l, &lname)); 1898 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 1899 if (isDepth) continue; 1900 PetscCall(DMGetLabel(sdm, lname, &mlabel)); 1901 PetscCall(DMLabelGetValue(mlabel, newp, &val)); 1902 if (val >= 0) PetscCall(DMLabelSetValue(mlabel, splitp, val)); 1903 } 1904 } 1905 } 1906 for (sp = 0; sp < numSP; ++sp) { 1907 const PetscInt dep = values[sp]; 1908 1909 if ((dep < 0) || (dep > depth)) continue; 1910 if (splitIS[dep]) PetscCall(ISRestoreIndices(splitIS[dep], &splitPoints[dep])); 1911 PetscCall(ISDestroy(&splitIS[dep])); 1912 if (unsplitIS[dep]) PetscCall(ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep])); 1913 PetscCall(ISDestroy(&unsplitIS[dep])); 1914 } 1915 if (ghostIS) PetscCall(ISRestoreIndices(ghostIS, &ghostPoints)); 1916 PetscCall(ISDestroy(&ghostIS)); 1917 if (label) { 1918 PetscCall(ISRestoreIndices(valueIS, &values)); 1919 PetscCall(ISDestroy(&valueIS)); 1920 } 1921 for (d = 0; d <= depth; ++d) { 1922 PetscCall(DMPlexGetDepthStratum(sdm, d, NULL, &pEnd)); 1923 pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d]; 1924 } 1925 PetscCall(PetscFree3(coneNew, coneONew, supportNew)); 1926 PetscCall(PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld)); 1927 PetscCall(PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints)); 1928 PetscFunctionReturn(PETSC_SUCCESS); 1929 } 1930 1931 /*@C 1932 DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface 1933 1934 Collective 1935 1936 Input Parameters: 1937 + dm - The original `DM` 1938 - label - The `DMLabel` specifying the boundary faces (this could be auto-generated) 1939 1940 Output Parameters: 1941 + splitLabel - The `DMLabel` containing the split points, or `NULL` if no output is desired 1942 - dmSplit - The new `DM` 1943 1944 Level: developer 1945 1946 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexLabelCohesiveComplete()` 1947 @*/ 1948 PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DMLabel splitLabel, DM *dmSplit) 1949 { 1950 DM sdm; 1951 PetscInt dim; 1952 1953 PetscFunctionBegin; 1954 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1955 PetscAssertPointer(dmSplit, 4); 1956 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &sdm)); 1957 PetscCall(DMSetType(sdm, DMPLEX)); 1958 PetscCall(DMGetDimension(dm, &dim)); 1959 PetscCall(DMSetDimension(sdm, dim)); 1960 switch (dim) { 1961 case 2: 1962 case 3: 1963 PetscCall(DMPlexConstructCohesiveCells_Internal(dm, label, splitLabel, sdm)); 1964 break; 1965 default: 1966 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %" PetscInt_FMT, dim); 1967 } 1968 *dmSplit = sdm; 1969 PetscFunctionReturn(PETSC_SUCCESS); 1970 } 1971 1972 /* Returns the side of the surface for a given cell with a face on the surface */ 1973 static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos) 1974 { 1975 const PetscInt *cone, *ornt; 1976 PetscInt dim, coneSize, c; 1977 1978 PetscFunctionBegin; 1979 *pos = PETSC_TRUE; 1980 PetscCall(DMGetDimension(dm, &dim)); 1981 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 1982 PetscCall(DMPlexGetCone(dm, cell, &cone)); 1983 PetscCall(DMPlexGetConeOrientation(dm, cell, &ornt)); 1984 for (c = 0; c < coneSize; ++c) { 1985 if (cone[c] == face) { 1986 PetscInt o = ornt[c]; 1987 1988 if (subdm) { 1989 const PetscInt *subcone, *subornt; 1990 PetscInt subpoint, subface, subconeSize, sc; 1991 1992 PetscCall(PetscFindInt(cell, numSubpoints, subpoints, &subpoint)); 1993 PetscCall(PetscFindInt(face, numSubpoints, subpoints, &subface)); 1994 PetscCall(DMPlexGetConeSize(subdm, subpoint, &subconeSize)); 1995 PetscCall(DMPlexGetCone(subdm, subpoint, &subcone)); 1996 PetscCall(DMPlexGetConeOrientation(subdm, subpoint, &subornt)); 1997 for (sc = 0; sc < subconeSize; ++sc) { 1998 if (subcone[sc] == subface) { 1999 o = subornt[0]; 2000 break; 2001 } 2002 } 2003 PetscCheck(sc < subconeSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find subpoint %" PetscInt_FMT " (%" PetscInt_FMT ") in cone for subpoint %" PetscInt_FMT " (%" PetscInt_FMT ")", subface, face, subpoint, cell); 2004 } 2005 if (o >= 0) *pos = PETSC_TRUE; 2006 else *pos = PETSC_FALSE; 2007 break; 2008 } 2009 } 2010 PetscCheck(c != coneSize, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " in split face %" PetscInt_FMT " support does not have it in the cone", cell, face); 2011 PetscFunctionReturn(PETSC_SUCCESS); 2012 } 2013 2014 static PetscErrorCode CheckFaultEdge_Private(DM dm, DMLabel label, PetscBool split) 2015 { 2016 IS facePosIS, faceNegIS, dimIS; 2017 const PetscInt *points; 2018 PetscInt *closure = NULL, *inclosure = NULL; 2019 PetscInt dim, numPoints, shift = 100, shift2 = 200, debug = 0; 2020 2021 PetscFunctionBegin; 2022 PetscCall(DMGetDimension(dm, &dim)); 2023 // If any faces touching the fault divide cells on either side, 2024 // either split them, or unsplit the connection 2025 PetscCall(DMLabelGetStratumIS(label, shift + dim - 1, &facePosIS)); 2026 PetscCall(DMLabelGetStratumIS(label, -(shift + dim - 1), &faceNegIS)); 2027 if (!facePosIS || !faceNegIS) { 2028 PetscCall(ISDestroy(&facePosIS)); 2029 PetscCall(ISDestroy(&faceNegIS)); 2030 PetscFunctionReturn(PETSC_SUCCESS); 2031 } 2032 PetscCall(ISExpand(facePosIS, faceNegIS, &dimIS)); 2033 PetscCall(ISDestroy(&facePosIS)); 2034 PetscCall(ISDestroy(&faceNegIS)); 2035 PetscCall(ISGetLocalSize(dimIS, &numPoints)); 2036 PetscCall(ISGetIndices(dimIS, &points)); 2037 for (PetscInt p = 0; p < numPoints; ++p) { 2038 const PetscInt point = points[p]; 2039 const PetscInt *support; 2040 PetscInt supportSize, valA, valB; 2041 2042 PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); 2043 if (supportSize != 2) continue; 2044 PetscCall(DMPlexGetSupport(dm, point, &support)); 2045 PetscCall(DMLabelGetValue(label, support[0], &valA)); 2046 PetscCall(DMLabelGetValue(label, support[1], &valB)); 2047 if ((valA == -1) || (valB == -1)) continue; 2048 if (valA * valB > 0) continue; 2049 // Check that this face is not incident on only unsplit faces, 2050 // meaning has at least one split face 2051 { 2052 PetscBool split = PETSC_FALSE; 2053 PetscInt Ncl, val; 2054 2055 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &Ncl, &closure)); 2056 for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) { 2057 PetscCall(DMLabelGetValue(label, closure[cl], &val)); 2058 if (val >= 0 && val <= dim) { 2059 split = PETSC_TRUE; 2060 break; 2061 } 2062 } 2063 if (!split) continue; 2064 } 2065 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT " is impinging (%" PetscInt_FMT ":%" PetscInt_FMT ", %" PetscInt_FMT ":%" PetscInt_FMT ")\n", point, support[0], valA, support[1], valB)); 2066 if (split) { 2067 // Split the face 2068 PetscCall(DMLabelGetValue(label, point, &valA)); 2069 PetscCall(DMLabelClearValue(label, point, valA)); 2070 PetscCall(DMLabelSetValue(label, point, dim - 1)); 2071 /* Label its closure: 2072 unmarked: label as unsplit 2073 incident: relabel as split 2074 split: do nothing */ 2075 { 2076 PetscInt closureSize, cl, dep; 2077 2078 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2079 for (cl = 0; cl < closureSize * 2; cl += 2) { 2080 PetscCall(DMLabelGetValue(label, closure[cl], &valA)); 2081 if (valA == -1) { /* Mark as unsplit */ 2082 PetscCall(DMPlexGetPointDepth(dm, closure[cl], &dep)); 2083 PetscCall(DMLabelSetValue(label, closure[cl], shift2 + dep)); 2084 } else if ((valA >= shift && valA < shift2) || (valA <= -shift && valA > -shift2)) { 2085 PetscCall(DMPlexGetPointDepth(dm, closure[cl], &dep)); 2086 PetscCall(DMLabelClearValue(label, closure[cl], valA)); 2087 PetscCall(DMLabelSetValue(label, closure[cl], dep)); 2088 } 2089 } 2090 } 2091 } else { 2092 // Unsplit the incident faces and their closures 2093 PetscInt Ncl, dep, val; 2094 2095 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &Ncl, &closure)); 2096 for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) { 2097 PetscCall(DMLabelGetValue(label, closure[cl], &val)); 2098 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Point %" PetscInt_FMT ":%" PetscInt_FMT "\n", closure[cl], val)); 2099 if (val >= 0 && val <= dim) { 2100 PetscInt Nincl, inval, indep; 2101 2102 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Point %" PetscInt_FMT " is being unsplit\n", closure[cl])); 2103 PetscCall(DMPlexGetPointDepth(dm, closure[cl], &dep)); 2104 PetscCall(DMLabelClearValue(label, closure[cl], val)); 2105 PetscCall(DMLabelSetValue(label, closure[cl], shift2 + dep)); 2106 2107 PetscCall(DMPlexGetTransitiveClosure(dm, closure[cl], PETSC_TRUE, &Nincl, &inclosure)); 2108 for (PetscInt incl = 0; incl < Nincl * 2; incl += 2) { 2109 PetscCall(DMLabelGetValue(label, inclosure[cl], &inval)); 2110 if (inval >= 0 && inval <= dim) { 2111 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Point %" PetscInt_FMT " is being unsplit\n", inclosure[incl])); 2112 PetscCall(DMPlexGetPointDepth(dm, inclosure[incl], &indep)); 2113 PetscCall(DMLabelClearValue(label, inclosure[incl], inval)); 2114 PetscCall(DMLabelSetValue(label, inclosure[incl], shift2 + indep)); 2115 } 2116 } 2117 } 2118 } 2119 } 2120 } 2121 PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &inclosure)); 2122 PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure)); 2123 PetscCall(ISRestoreIndices(dimIS, &points)); 2124 PetscCall(ISDestroy(&dimIS)); 2125 PetscFunctionReturn(PETSC_SUCCESS); 2126 } 2127 2128 /*@ 2129 DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces 2130 to complete the surface 2131 2132 Input Parameters: 2133 + dm - The `DM` 2134 . label - A `DMLabel` marking the surface 2135 . blabel - A `DMLabel` marking the vertices on the boundary which will not be duplicated, or `NULL` to find them automatically 2136 . bvalue - Value of `DMLabel` marking the vertices on the boundary 2137 . flip - Flag to flip the submesh normal and replace points on the other side 2138 . split - Split faces impinging on the surface, rather than clamping the surface boundary 2139 - subdm - The `DM` associated with the label, or `NULL` 2140 2141 Output Parameter: 2142 . label - A `DMLabel` marking all surface points 2143 2144 Level: developer 2145 2146 Note: 2147 The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation. 2148 2149 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexConstructCohesiveCells()`, `DMPlexLabelComplete()` 2150 @*/ 2151 PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscInt bvalue, PetscBool flip, PetscBool split, DM subdm) 2152 { 2153 DMLabel depthLabel; 2154 IS dimIS, subpointIS = NULL; 2155 const PetscInt *points, *subpoints; 2156 const PetscInt rev = flip ? -1 : 1; 2157 PetscInt shift = 100, shift2 = 200, shift3 = split ? 300 : 0, dim, depth, numPoints, numSubpoints, p, val; 2158 2159 PetscFunctionBegin; 2160 PetscCall(DMPlexGetDepth(dm, &depth)); 2161 PetscCall(DMGetDimension(dm, &dim)); 2162 PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); 2163 if (subdm) { 2164 PetscCall(DMPlexGetSubpointIS(subdm, &subpointIS)); 2165 if (subpointIS) { 2166 PetscCall(ISGetLocalSize(subpointIS, &numSubpoints)); 2167 PetscCall(ISGetIndices(subpointIS, &subpoints)); 2168 } 2169 } 2170 /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */ 2171 PetscCall(DMLabelGetStratumIS(label, dim - 1, &dimIS)); 2172 if (!dimIS) goto divide; 2173 PetscCall(ISGetLocalSize(dimIS, &numPoints)); 2174 PetscCall(ISGetIndices(dimIS, &points)); 2175 for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */ 2176 const PetscInt *support; 2177 PetscInt supportSize, s; 2178 2179 PetscCall(DMPlexGetSupportSize(dm, points[p], &supportSize)); 2180 #if 0 2181 if (supportSize != 2) { 2182 const PetscInt *lp; 2183 PetscInt Nlp, pind; 2184 2185 /* Check that for a cell with a single support face, that face is in the SF */ 2186 /* THis check only works for the remote side. We would need root side information */ 2187 PetscCall(PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL)); 2188 PetscCall(PetscFindInt(points[p], Nlp, lp, &pind)); 2189 PetscCheck(pind >= 0,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Split face %" PetscInt_FMT " has %" PetscInt_FMT " != 2 supports, and the face is not shared with another process", points[p], supportSize); 2190 } 2191 #endif 2192 PetscCall(DMPlexGetSupport(dm, points[p], &support)); 2193 for (s = 0; s < supportSize; ++s) { 2194 const PetscInt *cone; 2195 PetscInt coneSize, c; 2196 PetscBool pos; 2197 2198 PetscCall(GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos)); 2199 if (pos) PetscCall(DMLabelSetValue(label, support[s], rev * (shift + dim))); 2200 else PetscCall(DMLabelSetValue(label, support[s], -rev * (shift + dim))); 2201 if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE; 2202 /* Put faces touching the fault in the label */ 2203 PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize)); 2204 PetscCall(DMPlexGetCone(dm, support[s], &cone)); 2205 for (c = 0; c < coneSize; ++c) { 2206 const PetscInt point = cone[c]; 2207 2208 PetscCall(DMLabelGetValue(label, point, &val)); 2209 if (val == -1) { 2210 PetscInt *closure = NULL; 2211 PetscInt closureSize, cl; 2212 2213 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2214 for (cl = 0; cl < closureSize * 2; cl += 2) { 2215 const PetscInt clp = closure[cl]; 2216 PetscInt bval = -1; 2217 2218 PetscCall(DMLabelGetValue(label, clp, &val)); 2219 if (blabel) PetscCall(DMLabelGetValue(blabel, clp, &bval)); 2220 if (val >= 0 && val < dim - 1 && bval < 0) { 2221 PetscCall(DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift + dim - 1 : -(shift + dim - 1))); 2222 break; 2223 } 2224 } 2225 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2226 } 2227 } 2228 } 2229 } 2230 PetscCall(ISRestoreIndices(dimIS, &points)); 2231 PetscCall(ISDestroy(&dimIS)); 2232 /* Mark boundary points as unsplit */ 2233 if (blabel) { 2234 IS bdIS; 2235 2236 PetscCall(DMLabelGetStratumIS(blabel, bvalue, &bdIS)); 2237 PetscCall(ISGetLocalSize(bdIS, &numPoints)); 2238 PetscCall(ISGetIndices(bdIS, &points)); 2239 for (p = 0; p < numPoints; ++p) { 2240 const PetscInt point = points[p]; 2241 PetscInt val, bval; 2242 2243 PetscCall(DMLabelGetValue(blabel, point, &bval)); 2244 if (bval >= 0) { 2245 PetscCall(DMLabelGetValue(label, point, &val)); 2246 if ((val < 0) || (val > dim)) { 2247 /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */ 2248 PetscCall(DMLabelClearValue(blabel, point, bval)); 2249 } 2250 } 2251 } 2252 for (p = 0; p < numPoints; ++p) { 2253 const PetscInt point = points[p]; 2254 PetscInt val, bval; 2255 2256 PetscCall(DMLabelGetValue(blabel, point, &bval)); 2257 if (bval >= 0) { 2258 const PetscInt *cone, *support; 2259 PetscInt coneSize, supportSize, s, valA, valB, valE; 2260 2261 /* Mark as unsplit */ 2262 PetscCall(DMLabelGetValue(label, point, &val)); 2263 PetscCheck(val >= 0 && val <= dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has label value %" PetscInt_FMT ", should be part of the fault", point, val); 2264 PetscCall(DMLabelClearValue(label, point, val)); 2265 PetscCall(DMLabelSetValue(label, point, shift2 + val)); 2266 /* Check for cross-edge 2267 A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */ 2268 if (val != 0) continue; 2269 PetscCall(DMPlexGetSupport(dm, point, &support)); 2270 PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); 2271 for (s = 0; s < supportSize; ++s) { 2272 PetscCall(DMPlexGetCone(dm, support[s], &cone)); 2273 PetscCall(DMPlexGetConeSize(dm, support[s], &coneSize)); 2274 PetscCheck(coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %" PetscInt_FMT " has %" PetscInt_FMT " vertices != 2", support[s], coneSize); 2275 PetscCall(DMLabelGetValue(blabel, cone[0], &valA)); 2276 PetscCall(DMLabelGetValue(blabel, cone[1], &valB)); 2277 PetscCall(DMLabelGetValue(blabel, support[s], &valE)); 2278 if (valE < 0 && valA >= 0 && valB >= 0 && cone[0] != cone[1]) PetscCall(DMLabelSetValue(blabel, support[s], 2)); 2279 } 2280 } 2281 } 2282 PetscCall(ISRestoreIndices(bdIS, &points)); 2283 PetscCall(ISDestroy(&bdIS)); 2284 } 2285 /* Mark ghost fault cells */ 2286 { 2287 PetscSF sf; 2288 const PetscInt *leaves; 2289 PetscInt Nl, l; 2290 2291 PetscCall(DMGetPointSF(dm, &sf)); 2292 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL)); 2293 PetscCall(DMLabelGetStratumIS(label, dim - 1, &dimIS)); 2294 if (!dimIS) goto divide; 2295 PetscCall(ISGetLocalSize(dimIS, &numPoints)); 2296 PetscCall(ISGetIndices(dimIS, &points)); 2297 if (Nl > 0) { 2298 for (p = 0; p < numPoints; ++p) { 2299 const PetscInt point = points[p]; 2300 PetscInt val; 2301 2302 PetscCall(PetscFindInt(point, Nl, leaves, &l)); 2303 if (l >= 0) { 2304 PetscInt *closure = NULL; 2305 PetscInt closureSize, cl; 2306 2307 PetscCall(DMLabelGetValue(label, point, &val)); 2308 PetscCheck((val == dim - 1) || (val == shift2 + dim - 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " has label value %" PetscInt_FMT ", should be a fault face", point, val); 2309 PetscCall(DMLabelClearValue(label, point, val)); 2310 PetscCall(DMLabelSetValue(label, point, shift3 + val)); 2311 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2312 for (cl = 2; cl < closureSize * 2; cl += 2) { 2313 const PetscInt clp = closure[cl]; 2314 2315 PetscCall(DMLabelGetValue(label, clp, &val)); 2316 PetscCheck(val != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " is missing from label, but is in the closure of a fault face", point); 2317 PetscCall(DMLabelClearValue(label, clp, val)); 2318 PetscCall(DMLabelSetValue(label, clp, shift3 + val)); 2319 } 2320 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2321 } 2322 } 2323 } 2324 PetscCall(ISRestoreIndices(dimIS, &points)); 2325 PetscCall(ISDestroy(&dimIS)); 2326 } 2327 divide: 2328 if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 2329 PetscCall(DMPlexLabelFaultHalo(dm, label)); 2330 PetscCall(CheckFaultEdge_Private(dm, label, split)); 2331 PetscFunctionReturn(PETSC_SUCCESS); 2332 } 2333 2334 /* Check that no cell have all vertices on the fault */ 2335 static PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm) 2336 { 2337 IS subpointIS; 2338 const PetscInt *dmpoints; 2339 PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd; 2340 2341 PetscFunctionBegin; 2342 if (!label) PetscFunctionReturn(PETSC_SUCCESS); 2343 PetscCall(DMLabelGetDefaultValue(label, &defaultValue)); 2344 PetscCall(DMPlexGetSubpointIS(subdm, &subpointIS)); 2345 if (!subpointIS) PetscFunctionReturn(PETSC_SUCCESS); 2346 PetscCall(DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd)); 2347 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 2348 PetscCall(ISGetIndices(subpointIS, &dmpoints)); 2349 for (c = cStart; c < cEnd; ++c) { 2350 PetscBool invalidCell = PETSC_TRUE; 2351 PetscInt *closure = NULL; 2352 PetscInt closureSize, cl; 2353 2354 PetscCall(DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure)); 2355 for (cl = 0; cl < closureSize * 2; cl += 2) { 2356 PetscInt value = 0; 2357 2358 if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue; 2359 PetscCall(DMLabelGetValue(label, closure[cl], &value)); 2360 if (value == defaultValue) { 2361 invalidCell = PETSC_FALSE; 2362 break; 2363 } 2364 } 2365 PetscCall(DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure)); 2366 if (invalidCell) { 2367 PetscCall(ISRestoreIndices(subpointIS, &dmpoints)); 2368 PetscCall(ISDestroy(&subpointIS)); 2369 PetscCall(DMDestroy(&subdm)); 2370 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %" PetscInt_FMT " has all of its vertices on the submesh.", dmpoints[c]); 2371 } 2372 } 2373 PetscCall(ISRestoreIndices(subpointIS, &dmpoints)); 2374 PetscFunctionReturn(PETSC_SUCCESS); 2375 } 2376 2377 /*@ 2378 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 2379 2380 Collective 2381 2382 Input Parameters: 2383 + dm - The original `DM` 2384 . label - The label specifying the interface vertices 2385 . bdlabel - The optional label specifying the interface boundary vertices 2386 - bdvalue - Value of optional label specifying the interface boundary vertices 2387 2388 Output Parameters: 2389 + hybridLabel - The label fully marking the interface, or `NULL` if no output is desired 2390 . splitLabel - The label containing the split points, or `NULL` if no output is desired 2391 . dmInterface - The new interface `DM`, or `NULL` 2392 - dmHybrid - The new `DM` with cohesive cells 2393 2394 Level: developer 2395 2396 Note: 2397 The hybridLabel indicates what parts of the original mesh impinged on the division surface. For points 2398 directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be 2399 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with 2400 one vertex 3 on the surface would be 6 (101) and 3 (0) in hybridLabel. If an edge 9 from the negative side of the 2401 surface also hits vertex 3, it would be 9 (-101) in hybridLabel. 2402 2403 The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original 2404 mesh. The label value is $\pm 100+dim$ for each point. For example, if two edges 10 and 14 in the hybrid resulting from 2405 splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel. 2406 2407 The dmInterface is a `DM` built from the original division surface. It has a label which can be retrieved using 2408 `DMPlexGetSubpointMap()` which maps each point back to the point in the surface of the original mesh. 2409 2410 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexConstructCohesiveCells()`, `DMPlexLabelCohesiveComplete()`, `DMPlexGetSubpointMap()`, `DMCreate()` 2411 @*/ 2412 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, PetscInt bdvalue, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid) 2413 { 2414 DM idm; 2415 DMLabel subpointMap, hlabel, slabel = NULL; 2416 PetscInt dim; 2417 2418 PetscFunctionBegin; 2419 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2420 if (label) PetscAssertPointer(label, 2); 2421 if (bdlabel) PetscAssertPointer(bdlabel, 3); 2422 if (hybridLabel) PetscAssertPointer(hybridLabel, 5); 2423 if (splitLabel) PetscAssertPointer(splitLabel, 6); 2424 if (dmInterface) PetscAssertPointer(dmInterface, 7); 2425 PetscAssertPointer(dmHybrid, 8); 2426 PetscCall(DMGetDimension(dm, &dim)); 2427 PetscCall(DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm)); 2428 PetscCall(DMPlexCheckValidSubmesh_Private(dm, label, idm)); 2429 PetscCall(DMPlexOrient(idm)); 2430 PetscCall(DMPlexGetSubpointMap(idm, &subpointMap)); 2431 PetscCall(DMLabelDuplicate(subpointMap, &hlabel)); 2432 PetscCall(DMLabelClearStratum(hlabel, dim)); 2433 if (splitLabel) { 2434 const char *name; 2435 char sname[PETSC_MAX_PATH_LEN]; 2436 2437 PetscCall(PetscObjectGetName((PetscObject)hlabel, &name)); 2438 PetscCall(PetscStrncpy(sname, name, sizeof(sname))); 2439 PetscCall(PetscStrlcat(sname, " split", sizeof(sname))); 2440 PetscCall(DMLabelCreate(PETSC_COMM_SELF, sname, &slabel)); 2441 } 2442 PetscCall(DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, bdvalue, PETSC_FALSE, PETSC_TRUE, idm)); 2443 if (dmInterface) { 2444 *dmInterface = idm; 2445 } else PetscCall(DMDestroy(&idm)); 2446 PetscCall(DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid)); 2447 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *dmHybrid)); 2448 if (hybridLabel) *hybridLabel = hlabel; 2449 else PetscCall(DMLabelDestroy(&hlabel)); 2450 if (splitLabel) *splitLabel = slabel; 2451 { 2452 DM cdm; 2453 DMLabel ctLabel; 2454 2455 /* We need to somehow share the celltype label with the coordinate dm */ 2456 PetscCall(DMGetCoordinateDM(*dmHybrid, &cdm)); 2457 PetscCall(DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel)); 2458 PetscCall(DMSetLabel(cdm, ctLabel)); 2459 } 2460 PetscFunctionReturn(PETSC_SUCCESS); 2461 } 2462 2463 /* Here we need the explicit assumption that: 2464 2465 For any marked cell, the marked vertices constitute a single face 2466 */ 2467 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 2468 { 2469 IS subvertexIS = NULL; 2470 const PetscInt *subvertices; 2471 PetscInt *pStart, *pEnd, pSize; 2472 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 2473 2474 PetscFunctionBegin; 2475 *numFaces = 0; 2476 *nFV = 0; 2477 PetscCall(DMPlexGetDepth(dm, &depth)); 2478 PetscCall(DMGetDimension(dm, &dim)); 2479 pSize = PetscMax(depth, dim) + 1; 2480 PetscCall(PetscMalloc2(pSize, &pStart, pSize, &pEnd)); 2481 for (d = 0; d <= depth; ++d) PetscCall(DMPlexGetSimplexOrBoxCells(dm, depth - d, &pStart[d], &pEnd[d])); 2482 /* Loop over initial vertices and mark all faces in the collective star() */ 2483 if (vertexLabel) PetscCall(DMLabelGetStratumIS(vertexLabel, value, &subvertexIS)); 2484 if (subvertexIS) { 2485 PetscCall(ISGetSize(subvertexIS, &numSubVerticesInitial)); 2486 PetscCall(ISGetIndices(subvertexIS, &subvertices)); 2487 } 2488 for (v = 0; v < numSubVerticesInitial; ++v) { 2489 const PetscInt vertex = subvertices[v]; 2490 PetscInt *star = NULL; 2491 PetscInt starSize, s, numCells = 0, c; 2492 2493 PetscCall(DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star)); 2494 for (s = 0; s < starSize * 2; s += 2) { 2495 const PetscInt point = star[s]; 2496 if (point >= pStart[depth] && point < pEnd[depth]) star[numCells++] = point; 2497 } 2498 for (c = 0; c < numCells; ++c) { 2499 const PetscInt cell = star[c]; 2500 PetscInt *closure = NULL; 2501 PetscInt closureSize, cl; 2502 PetscInt cellLoc, numCorners = 0, faceSize = 0; 2503 2504 PetscCall(DMLabelGetValue(subpointMap, cell, &cellLoc)); 2505 if (cellLoc == 2) continue; 2506 PetscCheck(cellLoc < 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " has dimension %" PetscInt_FMT " in the surface label", cell, cellLoc); 2507 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 2508 for (cl = 0; cl < closureSize * 2; cl += 2) { 2509 const PetscInt point = closure[cl]; 2510 PetscInt vertexLoc; 2511 2512 if (point >= pStart[0] && point < pEnd[0]) { 2513 ++numCorners; 2514 PetscCall(DMLabelGetValue(vertexLabel, point, &vertexLoc)); 2515 if (vertexLoc == value) closure[faceSize++] = point; 2516 } 2517 } 2518 if (!*nFV) PetscCall(DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV)); 2519 PetscCheck(faceSize <= *nFV, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %" PetscInt_FMT " of an element on the surface", faceSize); 2520 if (faceSize == *nFV) { 2521 const PetscInt *cells = NULL; 2522 PetscInt numCells, nc; 2523 2524 ++(*numFaces); 2525 for (cl = 0; cl < faceSize; ++cl) PetscCall(DMLabelSetValue(subpointMap, closure[cl], 0)); 2526 PetscCall(DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells)); 2527 for (nc = 0; nc < numCells; ++nc) PetscCall(DMLabelSetValue(subpointMap, cells[nc], 2)); 2528 PetscCall(DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells)); 2529 } 2530 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 2531 } 2532 PetscCall(DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star)); 2533 } 2534 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subvertices)); 2535 PetscCall(ISDestroy(&subvertexIS)); 2536 PetscCall(PetscFree2(pStart, pEnd)); 2537 PetscFunctionReturn(PETSC_SUCCESS); 2538 } 2539 2540 PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, PetscBool addCells, DMLabel subpointMap, DM subdm) 2541 { 2542 IS subvertexIS = NULL; 2543 const PetscInt *subvertices; 2544 PetscInt *pStart, *pEnd; 2545 PetscInt dim, d, numSubVerticesInitial = 0, v; 2546 2547 PetscFunctionBegin; 2548 PetscCall(DMGetDimension(dm, &dim)); 2549 PetscCall(PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd)); 2550 for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetSimplexOrBoxCells(dm, dim - d, &pStart[d], &pEnd[d])); 2551 /* Loop over initial vertices and mark all faces in the collective star() */ 2552 if (vertexLabel) { 2553 PetscCall(DMLabelGetStratumIS(vertexLabel, value, &subvertexIS)); 2554 if (subvertexIS) { 2555 PetscCall(ISGetSize(subvertexIS, &numSubVerticesInitial)); 2556 PetscCall(ISGetIndices(subvertexIS, &subvertices)); 2557 } 2558 } 2559 for (v = 0; v < numSubVerticesInitial; ++v) { 2560 const PetscInt vertex = subvertices[v]; 2561 PetscInt *star = NULL; 2562 PetscInt starSize, s, numFaces = 0, f; 2563 2564 PetscCall(DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star)); 2565 for (s = 0; s < starSize * 2; s += 2) { 2566 const PetscInt point = star[s]; 2567 PetscInt faceLoc; 2568 2569 if (point >= pStart[dim - 1] && point < pEnd[dim - 1]) { 2570 if (markedFaces) { 2571 PetscCall(DMLabelGetValue(vertexLabel, point, &faceLoc)); 2572 if (faceLoc < 0) continue; 2573 } 2574 star[numFaces++] = point; 2575 } 2576 } 2577 for (f = 0; f < numFaces; ++f) { 2578 const PetscInt face = star[f]; 2579 PetscInt *closure = NULL; 2580 PetscInt closureSize, c; 2581 PetscInt faceLoc; 2582 2583 PetscCall(DMLabelGetValue(subpointMap, face, &faceLoc)); 2584 if (faceLoc == dim - 1) continue; 2585 PetscCheck(faceLoc < 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %" PetscInt_FMT " has dimension %" PetscInt_FMT " in the surface label", face, faceLoc); 2586 PetscCall(DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure)); 2587 for (c = 0; c < closureSize * 2; c += 2) { 2588 const PetscInt point = closure[c]; 2589 PetscInt vertexLoc; 2590 2591 if (point >= pStart[0] && point < pEnd[0]) { 2592 PetscCall(DMLabelGetValue(vertexLabel, point, &vertexLoc)); 2593 if (vertexLoc != value) break; 2594 } 2595 } 2596 if (c == closureSize * 2) { 2597 const PetscInt *support; 2598 PetscInt supportSize, s; 2599 2600 for (c = 0; c < closureSize * 2; c += 2) { 2601 const PetscInt point = closure[c]; 2602 2603 for (d = 0; d < dim; ++d) { 2604 if (point >= pStart[d] && point < pEnd[d]) { 2605 PetscCall(DMLabelSetValue(subpointMap, point, d)); 2606 break; 2607 } 2608 } 2609 } 2610 PetscCall(DMPlexGetSupportSize(dm, face, &supportSize)); 2611 PetscCall(DMPlexGetSupport(dm, face, &support)); 2612 if (addCells) 2613 for (s = 0; s < supportSize; ++s) PetscCall(DMLabelSetValue(subpointMap, support[s], dim)); 2614 } 2615 PetscCall(DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure)); 2616 } 2617 PetscCall(DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star)); 2618 } 2619 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subvertices)); 2620 PetscCall(ISDestroy(&subvertexIS)); 2621 PetscCall(PetscFree2(pStart, pEnd)); 2622 PetscFunctionReturn(PETSC_SUCCESS); 2623 } 2624 2625 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 2626 { 2627 DMLabel label = NULL; 2628 const PetscInt *cone; 2629 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1; 2630 2631 PetscFunctionBegin; 2632 *numFaces = 0; 2633 *nFV = 0; 2634 if (labelname) PetscCall(DMGetLabel(dm, labelname, &label)); 2635 *subCells = NULL; 2636 PetscCall(DMGetDimension(dm, &dim)); 2637 PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd)); 2638 if (cMax < 0) PetscFunctionReturn(PETSC_SUCCESS); 2639 if (label) { 2640 for (c = cMax; c < cEnd; ++c) { 2641 PetscInt val; 2642 2643 PetscCall(DMLabelGetValue(label, c, &val)); 2644 if (val == value) { 2645 ++(*numFaces); 2646 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 2647 } 2648 } 2649 } else { 2650 *numFaces = cEnd - cMax; 2651 PetscCall(DMPlexGetConeSize(dm, cMax, &coneSize)); 2652 } 2653 PetscCall(PetscMalloc1(*numFaces * 2, subCells)); 2654 if (!*numFaces) PetscFunctionReturn(PETSC_SUCCESS); 2655 *nFV = hasLagrange ? coneSize / 3 : coneSize / 2; 2656 for (c = cMax; c < cEnd; ++c) { 2657 const PetscInt *cells; 2658 PetscInt numCells; 2659 2660 if (label) { 2661 PetscInt val; 2662 2663 PetscCall(DMLabelGetValue(label, c, &val)); 2664 if (val != value) continue; 2665 } 2666 PetscCall(DMPlexGetCone(dm, c, &cone)); 2667 for (p = 0; p < *nFV; ++p) PetscCall(DMLabelSetValue(subpointMap, cone[p], 0)); 2668 /* Negative face */ 2669 PetscCall(DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells)); 2670 /* Not true in parallel 2671 PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2672 for (p = 0; p < numCells; ++p) { 2673 PetscCall(DMLabelSetValue(subpointMap, cells[p], 2)); 2674 (*subCells)[subc++] = cells[p]; 2675 } 2676 PetscCall(DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells)); 2677 /* Positive face is not included */ 2678 } 2679 PetscFunctionReturn(PETSC_SUCCESS); 2680 } 2681 2682 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 2683 { 2684 PetscInt *pStart, *pEnd; 2685 PetscInt dim, cMax, cEnd, c, d; 2686 2687 PetscFunctionBegin; 2688 PetscCall(DMGetDimension(dm, &dim)); 2689 PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd)); 2690 if (cMax < 0) PetscFunctionReturn(PETSC_SUCCESS); 2691 PetscCall(PetscMalloc2(dim + 1, &pStart, dim + 1, &pEnd)); 2692 for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d])); 2693 for (c = cMax; c < cEnd; ++c) { 2694 const PetscInt *cone; 2695 PetscInt *closure = NULL; 2696 PetscInt fconeSize, coneSize, closureSize, cl, val; 2697 2698 if (label) { 2699 PetscCall(DMLabelGetValue(label, c, &val)); 2700 if (val != value) continue; 2701 } 2702 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 2703 PetscCall(DMPlexGetCone(dm, c, &cone)); 2704 PetscCall(DMPlexGetConeSize(dm, cone[0], &fconeSize)); 2705 PetscCheck(coneSize == (fconeSize ? fconeSize : 1) + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2706 /* Negative face */ 2707 PetscCall(DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure)); 2708 for (cl = 0; cl < closureSize * 2; cl += 2) { 2709 const PetscInt point = closure[cl]; 2710 2711 for (d = 0; d <= dim; ++d) { 2712 if (point >= pStart[d] && point < pEnd[d]) { 2713 PetscCall(DMLabelSetValue(subpointMap, point, d)); 2714 break; 2715 } 2716 } 2717 } 2718 PetscCall(DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure)); 2719 /* Cells -- positive face is not included */ 2720 for (cl = 0; cl < 1; ++cl) { 2721 const PetscInt *support; 2722 PetscInt supportSize, s; 2723 2724 PetscCall(DMPlexGetSupportSize(dm, cone[cl], &supportSize)); 2725 /* PetscCheck(supportSize == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2726 PetscCall(DMPlexGetSupport(dm, cone[cl], &support)); 2727 for (s = 0; s < supportSize; ++s) PetscCall(DMLabelSetValue(subpointMap, support[s], dim)); 2728 } 2729 } 2730 PetscCall(PetscFree2(pStart, pEnd)); 2731 PetscFunctionReturn(PETSC_SUCCESS); 2732 } 2733 2734 static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2735 { 2736 MPI_Comm comm; 2737 PetscBool posOrient = PETSC_FALSE; 2738 const PetscInt debug = 0; 2739 PetscInt cellDim, faceSize, f; 2740 2741 PetscFunctionBegin; 2742 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2743 PetscCall(DMGetDimension(dm, &cellDim)); 2744 if (debug) PetscCall(PetscPrintf(comm, "cellDim: %" PetscInt_FMT " numCorners: %" PetscInt_FMT "\n", cellDim, numCorners)); 2745 2746 if (cellDim == 1 && numCorners == 2) { 2747 /* Triangle */ 2748 faceSize = numCorners - 1; 2749 posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE; 2750 } else if (cellDim == 2 && numCorners == 3) { 2751 /* Triangle */ 2752 faceSize = numCorners - 1; 2753 posOrient = !(oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE; 2754 } else if (cellDim == 3 && numCorners == 4) { 2755 /* Tetrahedron */ 2756 faceSize = numCorners - 1; 2757 posOrient = (oppositeVertex % 2) ? PETSC_TRUE : PETSC_FALSE; 2758 } else if (cellDim == 1 && numCorners == 3) { 2759 /* Quadratic line */ 2760 faceSize = 1; 2761 posOrient = PETSC_TRUE; 2762 } else if (cellDim == 2 && numCorners == 4) { 2763 /* Quads */ 2764 faceSize = 2; 2765 if (indices[1] > indices[0] && indices[1] - indices[0] == 1) { 2766 posOrient = PETSC_TRUE; 2767 } else if (indices[0] == 3 && indices[1] == 0) { 2768 posOrient = PETSC_TRUE; 2769 } else { 2770 PetscCheck((indices[0] > indices[1] && indices[0] - indices[1] == 1) || (indices[0] == 0 && indices[1] == 3), comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2771 posOrient = PETSC_FALSE; 2772 } 2773 } else if (cellDim == 2 && numCorners == 6) { 2774 /* Quadratic triangle (I hate this) */ 2775 /* Edges are determined by the first 2 vertices (corners of edges) */ 2776 const PetscInt faceSizeTri = 3; 2777 PetscInt sortedIndices[3], i, iFace; 2778 PetscBool found = PETSC_FALSE; 2779 PetscInt faceVerticesTriSorted[9] = { 2780 0, 3, 4, /* bottom */ 2781 1, 4, 5, /* right */ 2782 2, 3, 5, /* left */ 2783 }; 2784 PetscInt faceVerticesTri[9] = { 2785 0, 3, 4, /* bottom */ 2786 1, 4, 5, /* right */ 2787 2, 5, 3, /* left */ 2788 }; 2789 2790 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2791 PetscCall(PetscSortInt(faceSizeTri, sortedIndices)); 2792 for (iFace = 0; iFace < 3; ++iFace) { 2793 const PetscInt ii = iFace * faceSizeTri; 2794 PetscInt fVertex, cVertex; 2795 2796 if (sortedIndices[0] == faceVerticesTriSorted[ii + 0] && sortedIndices[1] == faceVerticesTriSorted[ii + 1]) { 2797 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2798 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2799 if (indices[cVertex] == faceVerticesTri[ii + fVertex]) { 2800 faceVertices[fVertex] = origVertices[cVertex]; 2801 break; 2802 } 2803 } 2804 } 2805 found = PETSC_TRUE; 2806 break; 2807 } 2808 } 2809 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2810 if (posOriented) *posOriented = PETSC_TRUE; 2811 PetscFunctionReturn(PETSC_SUCCESS); 2812 } else if (cellDim == 2 && numCorners == 9) { 2813 /* Quadratic quad (I hate this) */ 2814 /* Edges are determined by the first 2 vertices (corners of edges) */ 2815 const PetscInt faceSizeQuad = 3; 2816 PetscInt sortedIndices[3], i, iFace; 2817 PetscBool found = PETSC_FALSE; 2818 PetscInt faceVerticesQuadSorted[12] = { 2819 0, 1, 4, /* bottom */ 2820 1, 2, 5, /* right */ 2821 2, 3, 6, /* top */ 2822 0, 3, 7, /* left */ 2823 }; 2824 PetscInt faceVerticesQuad[12] = { 2825 0, 1, 4, /* bottom */ 2826 1, 2, 5, /* right */ 2827 2, 3, 6, /* top */ 2828 3, 0, 7, /* left */ 2829 }; 2830 2831 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2832 PetscCall(PetscSortInt(faceSizeQuad, sortedIndices)); 2833 for (iFace = 0; iFace < 4; ++iFace) { 2834 const PetscInt ii = iFace * faceSizeQuad; 2835 PetscInt fVertex, cVertex; 2836 2837 if (sortedIndices[0] == faceVerticesQuadSorted[ii + 0] && sortedIndices[1] == faceVerticesQuadSorted[ii + 1]) { 2838 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2839 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2840 if (indices[cVertex] == faceVerticesQuad[ii + fVertex]) { 2841 faceVertices[fVertex] = origVertices[cVertex]; 2842 break; 2843 } 2844 } 2845 } 2846 found = PETSC_TRUE; 2847 break; 2848 } 2849 } 2850 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2851 if (posOriented) *posOriented = PETSC_TRUE; 2852 PetscFunctionReturn(PETSC_SUCCESS); 2853 } else if (cellDim == 3 && numCorners == 8) { 2854 /* Hexes 2855 A hex is two oriented quads with the normal of the first 2856 pointing up at the second. 2857 2858 7---6 2859 /| /| 2860 4---5 | 2861 | 1-|-2 2862 |/ |/ 2863 0---3 2864 2865 Faces are determined by the first 4 vertices (corners of faces) */ 2866 const PetscInt faceSizeHex = 4; 2867 PetscInt sortedIndices[4], i, iFace; 2868 PetscBool found = PETSC_FALSE; 2869 PetscInt faceVerticesHexSorted[24] = { 2870 0, 1, 2, 3, /* bottom */ 2871 4, 5, 6, 7, /* top */ 2872 0, 3, 4, 5, /* front */ 2873 2, 3, 5, 6, /* right */ 2874 1, 2, 6, 7, /* back */ 2875 0, 1, 4, 7, /* left */ 2876 }; 2877 PetscInt faceVerticesHex[24] = { 2878 1, 2, 3, 0, /* bottom */ 2879 4, 5, 6, 7, /* top */ 2880 0, 3, 5, 4, /* front */ 2881 3, 2, 6, 5, /* right */ 2882 2, 1, 7, 6, /* back */ 2883 1, 0, 4, 7, /* left */ 2884 }; 2885 2886 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2887 PetscCall(PetscSortInt(faceSizeHex, sortedIndices)); 2888 for (iFace = 0; iFace < 6; ++iFace) { 2889 const PetscInt ii = iFace * faceSizeHex; 2890 PetscInt fVertex, cVertex; 2891 2892 if (sortedIndices[0] == faceVerticesHexSorted[ii + 0] && sortedIndices[1] == faceVerticesHexSorted[ii + 1] && sortedIndices[2] == faceVerticesHexSorted[ii + 2] && sortedIndices[3] == faceVerticesHexSorted[ii + 3]) { 2893 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2894 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2895 if (indices[cVertex] == faceVerticesHex[ii + fVertex]) { 2896 faceVertices[fVertex] = origVertices[cVertex]; 2897 break; 2898 } 2899 } 2900 } 2901 found = PETSC_TRUE; 2902 break; 2903 } 2904 } 2905 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2906 if (posOriented) *posOriented = PETSC_TRUE; 2907 PetscFunctionReturn(PETSC_SUCCESS); 2908 } else if (cellDim == 3 && numCorners == 10) { 2909 /* Quadratic tet */ 2910 /* Faces are determined by the first 3 vertices (corners of faces) */ 2911 const PetscInt faceSizeTet = 6; 2912 PetscInt sortedIndices[6], i, iFace; 2913 PetscBool found = PETSC_FALSE; 2914 PetscInt faceVerticesTetSorted[24] = { 2915 0, 1, 2, 6, 7, 8, /* bottom */ 2916 0, 3, 4, 6, 7, 9, /* front */ 2917 1, 4, 5, 7, 8, 9, /* right */ 2918 2, 3, 5, 6, 8, 9, /* left */ 2919 }; 2920 PetscInt faceVerticesTet[24] = { 2921 0, 1, 2, 6, 7, 8, /* bottom */ 2922 0, 4, 3, 6, 7, 9, /* front */ 2923 1, 5, 4, 7, 8, 9, /* right */ 2924 2, 3, 5, 8, 6, 9, /* left */ 2925 }; 2926 2927 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2928 PetscCall(PetscSortInt(faceSizeTet, sortedIndices)); 2929 for (iFace = 0; iFace < 4; ++iFace) { 2930 const PetscInt ii = iFace * faceSizeTet; 2931 PetscInt fVertex, cVertex; 2932 2933 if (sortedIndices[0] == faceVerticesTetSorted[ii + 0] && sortedIndices[1] == faceVerticesTetSorted[ii + 1] && sortedIndices[2] == faceVerticesTetSorted[ii + 2] && sortedIndices[3] == faceVerticesTetSorted[ii + 3]) { 2934 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2935 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2936 if (indices[cVertex] == faceVerticesTet[ii + fVertex]) { 2937 faceVertices[fVertex] = origVertices[cVertex]; 2938 break; 2939 } 2940 } 2941 } 2942 found = PETSC_TRUE; 2943 break; 2944 } 2945 } 2946 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2947 if (posOriented) *posOriented = PETSC_TRUE; 2948 PetscFunctionReturn(PETSC_SUCCESS); 2949 } else if (cellDim == 3 && numCorners == 27) { 2950 /* Quadratic hexes (I hate this) 2951 A hex is two oriented quads with the normal of the first 2952 pointing up at the second. 2953 2954 7---6 2955 /| /| 2956 4---5 | 2957 | 3-|-2 2958 |/ |/ 2959 0---1 2960 2961 Faces are determined by the first 4 vertices (corners of faces) */ 2962 const PetscInt faceSizeQuadHex = 9; 2963 PetscInt sortedIndices[9], i, iFace; 2964 PetscBool found = PETSC_FALSE; 2965 PetscInt faceVerticesQuadHexSorted[54] = { 2966 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2967 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2968 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2969 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2970 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2971 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2972 }; 2973 PetscInt faceVerticesQuadHex[54] = { 2974 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2975 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2976 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2977 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2978 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2979 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2980 }; 2981 2982 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2983 PetscCall(PetscSortInt(faceSizeQuadHex, sortedIndices)); 2984 for (iFace = 0; iFace < 6; ++iFace) { 2985 const PetscInt ii = iFace * faceSizeQuadHex; 2986 PetscInt fVertex, cVertex; 2987 2988 if (sortedIndices[0] == faceVerticesQuadHexSorted[ii + 0] && sortedIndices[1] == faceVerticesQuadHexSorted[ii + 1] && sortedIndices[2] == faceVerticesQuadHexSorted[ii + 2] && sortedIndices[3] == faceVerticesQuadHexSorted[ii + 3]) { 2989 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2990 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2991 if (indices[cVertex] == faceVerticesQuadHex[ii + fVertex]) { 2992 faceVertices[fVertex] = origVertices[cVertex]; 2993 break; 2994 } 2995 } 2996 } 2997 found = PETSC_TRUE; 2998 break; 2999 } 3000 } 3001 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 3002 if (posOriented) *posOriented = PETSC_TRUE; 3003 PetscFunctionReturn(PETSC_SUCCESS); 3004 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 3005 if (!posOrient) { 3006 if (debug) PetscCall(PetscPrintf(comm, " Reversing initial face orientation\n")); 3007 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize - 1 - f]; 3008 } else { 3009 if (debug) PetscCall(PetscPrintf(comm, " Keeping initial face orientation\n")); 3010 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 3011 } 3012 if (posOriented) *posOriented = posOrient; 3013 PetscFunctionReturn(PETSC_SUCCESS); 3014 } 3015 3016 /*@ 3017 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 3018 in faceVertices. The orientation is such that the face normal points out of the cell 3019 3020 Not Collective 3021 3022 Input Parameters: 3023 + dm - The original mesh 3024 . cell - The cell mesh point 3025 . faceSize - The number of vertices on the face 3026 . face - The face vertices 3027 . numCorners - The number of vertices on the cell 3028 . indices - Local numbering of face vertices in cell cone 3029 - origVertices - Original face vertices 3030 3031 Output Parameters: 3032 + faceVertices - The face vertices properly oriented 3033 - posOriented - `PETSC_TRUE` if the face was oriented with outward normal 3034 3035 Level: developer 3036 3037 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetCone()` 3038 @*/ 3039 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 3040 { 3041 const PetscInt *cone = NULL; 3042 PetscInt coneSize, v, f, v2; 3043 PetscInt oppositeVertex = -1; 3044 3045 PetscFunctionBegin; 3046 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3047 PetscCall(DMPlexGetCone(dm, cell, &cone)); 3048 for (v = 0, v2 = 0; v < coneSize; ++v) { 3049 PetscBool found = PETSC_FALSE; 3050 3051 for (f = 0; f < faceSize; ++f) { 3052 if (face[f] == cone[v]) { 3053 found = PETSC_TRUE; 3054 break; 3055 } 3056 } 3057 if (found) { 3058 indices[v2] = v; 3059 origVertices[v2] = cone[v]; 3060 ++v2; 3061 } else { 3062 oppositeVertex = v; 3063 } 3064 } 3065 PetscCall(DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented)); 3066 PetscFunctionReturn(PETSC_SUCCESS); 3067 } 3068 3069 /* 3070 DMPlexInsertFace_Internal - Puts a face into the mesh 3071 3072 Not Collective 3073 3074 Input Parameters: 3075 + dm - The `DMPLEX` 3076 . numFaceVertex - The number of vertices in the face 3077 . faceVertices - The vertices in the face for `dm` 3078 . subfaceVertices - The vertices in the face for subdm 3079 . numCorners - The number of vertices in the `cell` 3080 . cell - A cell in `dm` containing the face 3081 . subcell - A cell in subdm containing the face 3082 . firstFace - First face in the mesh 3083 - newFacePoint - Next face in the mesh 3084 3085 Output Parameter: 3086 . newFacePoint - Contains next face point number on input, updated on output 3087 3088 Level: developer 3089 */ 3090 static PetscErrorCode DMPlexInsertFace_Internal(DM dm, DM subdm, PetscInt numFaceVertices, const PetscInt faceVertices[], const PetscInt subfaceVertices[], PetscInt numCorners, PetscInt cell, PetscInt subcell, PetscInt firstFace, PetscInt *newFacePoint) 3091 { 3092 MPI_Comm comm; 3093 DM_Plex *submesh = (DM_Plex *)subdm->data; 3094 const PetscInt *faces; 3095 PetscInt numFaces, coneSize; 3096 3097 PetscFunctionBegin; 3098 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3099 PetscCall(DMPlexGetConeSize(subdm, subcell, &coneSize)); 3100 PetscCheck(coneSize == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %" PetscInt_FMT " is %" PetscInt_FMT " != 1", cell, coneSize); 3101 #if 0 3102 /* Cannot use this because support() has not been constructed yet */ 3103 PetscCall(DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces)); 3104 #else 3105 { 3106 PetscInt f; 3107 3108 numFaces = 0; 3109 PetscCall(DMGetWorkArray(subdm, 1, MPIU_INT, (void **)&faces)); 3110 for (f = firstFace; f < *newFacePoint; ++f) { 3111 PetscInt dof, off, d; 3112 3113 PetscCall(PetscSectionGetDof(submesh->coneSection, f, &dof)); 3114 PetscCall(PetscSectionGetOffset(submesh->coneSection, f, &off)); 3115 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 3116 for (d = 0; d < dof; ++d) { 3117 const PetscInt p = submesh->cones[off + d]; 3118 PetscInt v; 3119 3120 for (v = 0; v < numFaceVertices; ++v) { 3121 if (subfaceVertices[v] == p) break; 3122 } 3123 if (v == numFaceVertices) break; 3124 } 3125 if (d == dof) { 3126 numFaces = 1; 3127 ((PetscInt *)faces)[0] = f; 3128 } 3129 } 3130 } 3131 #endif 3132 PetscCheck(numFaces <= 1, comm, PETSC_ERR_ARG_WRONG, "Vertex set had %" PetscInt_FMT " faces, not one", numFaces); 3133 if (numFaces == 1) { 3134 /* Add the other cell neighbor for this face */ 3135 PetscCall(DMPlexSetCone(subdm, subcell, faces)); 3136 } else { 3137 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 3138 PetscBool posOriented; 3139 3140 PetscCall(DMGetWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices)); 3141 origVertices = &orientedVertices[numFaceVertices]; 3142 indices = &orientedVertices[numFaceVertices * 2]; 3143 orientedSubVertices = &orientedVertices[numFaceVertices * 3]; 3144 PetscCall(DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented)); 3145 /* TODO: I know that routine should return a permutation, not the indices */ 3146 for (v = 0; v < numFaceVertices; ++v) { 3147 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 3148 for (ov = 0; ov < numFaceVertices; ++ov) { 3149 if (orientedVertices[ov] == vertex) { 3150 orientedSubVertices[ov] = subvertex; 3151 break; 3152 } 3153 } 3154 PetscCheck(ov != numFaceVertices, comm, PETSC_ERR_PLIB, "Could not find face vertex %" PetscInt_FMT " in orientated set", vertex); 3155 } 3156 PetscCall(DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices)); 3157 PetscCall(DMPlexSetCone(subdm, subcell, newFacePoint)); 3158 PetscCall(DMRestoreWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices)); 3159 ++(*newFacePoint); 3160 } 3161 #if 0 3162 PetscCall(DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces)); 3163 #else 3164 PetscCall(DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **)&faces)); 3165 #endif 3166 PetscFunctionReturn(PETSC_SUCCESS); 3167 } 3168 3169 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 3170 { 3171 MPI_Comm comm; 3172 DMLabel subpointMap; 3173 IS subvertexIS, subcellIS; 3174 const PetscInt *subVertices, *subCells; 3175 PetscInt numSubVertices, firstSubVertex, numSubCells; 3176 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 3177 PetscInt vStart, vEnd, c, f; 3178 3179 PetscFunctionBegin; 3180 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3181 /* Create subpointMap which marks the submesh */ 3182 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3183 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3184 if (vertexLabel) PetscCall(DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm)); 3185 /* Setup chart */ 3186 PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3187 PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3188 PetscCall(DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices)); 3189 PetscCall(DMPlexSetVTKCellHeight(subdm, 1)); 3190 /* Set cone sizes */ 3191 firstSubVertex = numSubCells; 3192 firstSubFace = numSubCells + numSubVertices; 3193 newFacePoint = firstSubFace; 3194 PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3195 if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices)); 3196 PetscCall(DMLabelGetStratumIS(subpointMap, 2, &subcellIS)); 3197 if (subcellIS) PetscCall(ISGetIndices(subcellIS, &subCells)); 3198 for (c = 0; c < numSubCells; ++c) PetscCall(DMPlexSetConeSize(subdm, c, 1)); 3199 for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) PetscCall(DMPlexSetConeSize(subdm, f, nFV)); 3200 PetscCall(DMSetUp(subdm)); 3201 PetscCall(DMLabelDestroy(&subpointMap)); 3202 /* Create face cones */ 3203 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3204 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3205 PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 3206 for (c = 0; c < numSubCells; ++c) { 3207 const PetscInt cell = subCells[c]; 3208 const PetscInt subcell = c; 3209 PetscInt *closure = NULL; 3210 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 3211 3212 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 3213 for (cl = 0; cl < closureSize * 2; cl += 2) { 3214 const PetscInt point = closure[cl]; 3215 PetscInt subVertex; 3216 3217 if (point >= vStart && point < vEnd) { 3218 ++numCorners; 3219 PetscCall(PetscFindInt(point, numSubVertices, subVertices, &subVertex)); 3220 if (subVertex >= 0) { 3221 closure[faceSize] = point; 3222 subface[faceSize] = firstSubVertex + subVertex; 3223 ++faceSize; 3224 } 3225 } 3226 } 3227 PetscCheck(faceSize <= nFV, comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %" PetscInt_FMT " of an element on the surface", faceSize); 3228 if (faceSize == nFV) PetscCall(DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint)); 3229 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 3230 } 3231 PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 3232 PetscCall(DMPlexSymmetrize(subdm)); 3233 PetscCall(DMPlexStratify(subdm)); 3234 /* Build coordinates */ 3235 { 3236 PetscSection coordSection, subCoordSection; 3237 Vec coordinates, subCoordinates; 3238 PetscScalar *coords, *subCoords; 3239 PetscInt numComp, coordSize, v; 3240 const char *name; 3241 3242 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3243 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3244 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3245 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3246 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3247 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3248 PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices)); 3249 for (v = 0; v < numSubVertices; ++v) { 3250 const PetscInt vertex = subVertices[v]; 3251 const PetscInt subvertex = firstSubVertex + v; 3252 PetscInt dof; 3253 3254 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3255 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3256 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3257 } 3258 PetscCall(PetscSectionSetUp(subCoordSection)); 3259 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3260 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3261 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 3262 PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name)); 3263 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3264 PetscCall(VecSetType(subCoordinates, VECSTANDARD)); 3265 if (coordSize) { 3266 PetscCall(VecGetArray(coordinates, &coords)); 3267 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3268 for (v = 0; v < numSubVertices; ++v) { 3269 const PetscInt vertex = subVertices[v]; 3270 const PetscInt subvertex = firstSubVertex + v; 3271 PetscInt dof, off, sdof, soff, d; 3272 3273 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3274 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 3275 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3276 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3277 PetscCheck(dof == sdof, comm, PETSC_ERR_PLIB, "Coordinate dimension %" PetscInt_FMT " on subvertex %" PetscInt_FMT ", vertex %" PetscInt_FMT " should be %" PetscInt_FMT, sdof, subvertex, vertex, dof); 3278 for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d]; 3279 } 3280 PetscCall(VecRestoreArray(coordinates, &coords)); 3281 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3282 } 3283 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3284 PetscCall(VecDestroy(&subCoordinates)); 3285 } 3286 /* Cleanup */ 3287 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices)); 3288 PetscCall(ISDestroy(&subvertexIS)); 3289 if (subcellIS) PetscCall(ISRestoreIndices(subcellIS, &subCells)); 3290 PetscCall(ISDestroy(&subcellIS)); 3291 PetscFunctionReturn(PETSC_SUCCESS); 3292 } 3293 3294 /* TODO: Fix this to properly propagate up error conditions it may find */ 3295 static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 3296 { 3297 PetscInt subPoint; 3298 PetscErrorCode ierr; 3299 3300 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); 3301 if (ierr) return -1; 3302 return subPoint < 0 ? subPoint : firstSubPoint + subPoint; 3303 } 3304 3305 /* TODO: Fix this to properly propagate up error conditions it may find */ 3306 static inline PetscInt DMPlexFilterPointPerm_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[], const PetscInt subIndices[]) 3307 { 3308 PetscInt subPoint; 3309 PetscErrorCode ierr; 3310 3311 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); 3312 if (ierr) return -1; 3313 return subPoint < 0 ? subPoint : firstSubPoint + (subIndices ? subIndices[subPoint] : subPoint); 3314 } 3315 3316 static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm) 3317 { 3318 DMLabel depthLabel; 3319 PetscInt Nl, l, d; 3320 3321 PetscFunctionBegin; 3322 // Reset depth label for fast lookup 3323 PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); 3324 PetscCall(DMLabelMakeAllInvalid_Internal(depthLabel)); 3325 PetscCall(DMGetNumLabels(dm, &Nl)); 3326 for (l = 0; l < Nl; ++l) { 3327 DMLabel label, newlabel; 3328 const char *lname; 3329 PetscBool isDepth, isDim, isCelltype, isVTK; 3330 IS valueIS; 3331 const PetscInt *values; 3332 PetscInt Nv, v; 3333 3334 PetscCall(DMGetLabelName(dm, l, &lname)); 3335 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 3336 PetscCall(PetscStrcmp(lname, "dim", &isDim)); 3337 PetscCall(PetscStrcmp(lname, "celltype", &isCelltype)); 3338 PetscCall(PetscStrcmp(lname, "vtk", &isVTK)); 3339 if (isDepth || isDim || isCelltype || isVTK) continue; 3340 PetscCall(DMCreateLabel(subdm, lname)); 3341 PetscCall(DMGetLabel(dm, lname, &label)); 3342 PetscCall(DMGetLabel(subdm, lname, &newlabel)); 3343 PetscCall(DMLabelGetDefaultValue(label, &v)); 3344 PetscCall(DMLabelSetDefaultValue(newlabel, v)); 3345 PetscCall(DMLabelGetValueIS(label, &valueIS)); 3346 PetscCall(ISGetLocalSize(valueIS, &Nv)); 3347 PetscCall(ISGetIndices(valueIS, &values)); 3348 for (v = 0; v < Nv; ++v) { 3349 IS pointIS; 3350 const PetscInt *points; 3351 PetscInt Np, p; 3352 3353 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 3354 PetscCall(ISGetLocalSize(pointIS, &Np)); 3355 PetscCall(ISGetIndices(pointIS, &points)); 3356 for (p = 0; p < Np; ++p) { 3357 const PetscInt point = points[p]; 3358 PetscInt subp, subdepth; 3359 3360 PetscCall(DMPlexGetPointDepth(dm, point, &d)); 3361 PetscCall(DMPlexGetDepth(subdm, &subdepth)); 3362 if (d > subdepth) continue; 3363 subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]); 3364 if (subp >= 0) PetscCall(DMLabelSetValue(newlabel, subp, values[v])); 3365 } 3366 PetscCall(ISRestoreIndices(pointIS, &points)); 3367 PetscCall(ISDestroy(&pointIS)); 3368 } 3369 PetscCall(ISRestoreIndices(valueIS, &values)); 3370 PetscCall(ISDestroy(&valueIS)); 3371 } 3372 PetscFunctionReturn(PETSC_SUCCESS); 3373 } 3374 3375 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, PetscBool ignoreLabelHalo, PetscBool sanitizeSubmesh, PetscSF *ownershipTransferSF, DM subdm) 3376 { 3377 MPI_Comm comm; 3378 DMLabel subpointMap; 3379 IS *subpointIS; 3380 const PetscInt **subpoints; 3381 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 3382 PetscInt totSubPoints = 0, maxConeSize, dim, sdim, cdim, p, d, coordinate_type; 3383 PetscMPIInt rank; 3384 3385 PetscFunctionBegin; 3386 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3387 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 3388 /* Create subpointMap which marks the submesh */ 3389 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3390 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3391 if (cellHeight) { 3392 if (isCohesive) PetscCall(DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm)); 3393 else PetscCall(DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, PETSC_TRUE, subpointMap, subdm)); 3394 } else { 3395 DMLabel depth; 3396 IS pointIS; 3397 const PetscInt *points; 3398 PetscInt numPoints = 0, pStart = 0; 3399 PetscBool *isGhost = NULL; 3400 3401 PetscCall(DMPlexGetDepthLabel(dm, &depth)); 3402 PetscCall(DMLabelGetStratumIS(label, value, &pointIS)); 3403 if (pointIS) { 3404 PetscCall(ISGetIndices(pointIS, &points)); 3405 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 3406 } 3407 if (ignoreLabelHalo) { 3408 PetscSF pointSF; 3409 PetscInt nleaves, pEnd; 3410 const PetscInt *ilocal = NULL; 3411 3412 PetscCall(DMGetPointSF(dm, &pointSF)); 3413 PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, NULL)); 3414 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3415 PetscCall(PetscMalloc1(pEnd - pStart, &isGhost)); 3416 for (p = 0; p < pEnd - pStart; ++p) isGhost[p] = PETSC_FALSE; 3417 for (p = 0; p < nleaves; ++p) isGhost[(ilocal ? ilocal[p] : p) - pStart] = PETSC_TRUE; 3418 } 3419 for (p = 0; p < numPoints; ++p) { 3420 PetscInt *closure = NULL; 3421 PetscInt closureSize, c, pdim; 3422 3423 if (ignoreLabelHalo && isGhost[points[p] - pStart]) continue; 3424 PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3425 for (c = 0; c < closureSize * 2; c += 2) { 3426 PetscCall(DMLabelGetValue(depth, closure[c], &pdim)); 3427 PetscCall(DMLabelSetValue(subpointMap, closure[c], pdim)); 3428 } 3429 PetscCall(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3430 } 3431 PetscCall(PetscFree(isGhost)); 3432 if (pointIS) PetscCall(ISRestoreIndices(pointIS, &points)); 3433 PetscCall(ISDestroy(&pointIS)); 3434 } 3435 /* Setup chart */ 3436 PetscCall(DMGetDimension(dm, &dim)); 3437 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3438 PetscCall(PetscMalloc4(dim + 1, &numSubPoints, dim + 1, &firstSubPoint, dim + 1, &subpointIS, dim + 1, &subpoints)); 3439 for (d = 0; d <= dim; ++d) { 3440 PetscCall(DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d])); 3441 totSubPoints += numSubPoints[d]; 3442 } 3443 // Determine submesh dimension 3444 PetscCall(DMGetDimension(subdm, &sdim)); 3445 if (sdim > 0) { 3446 // Calling function knows what dimension to use, and we include neighboring cells as well 3447 sdim = dim; 3448 } else { 3449 // We reset the subdimension based on what is being selected 3450 PetscInt lsdim; 3451 for (lsdim = dim; lsdim >= 0; --lsdim) 3452 if (numSubPoints[lsdim]) break; 3453 PetscCallMPI(MPIU_Allreduce(&lsdim, &sdim, 1, MPIU_INT, MPI_MAX, comm)); 3454 PetscCall(DMSetDimension(subdm, sdim)); 3455 PetscCall(DMSetCoordinateDim(subdm, cdim)); 3456 } 3457 PetscCall(DMPlexSetChart(subdm, 0, totSubPoints)); 3458 PetscCall(DMPlexSetVTKCellHeight(subdm, cellHeight)); 3459 /* Set cone sizes */ 3460 firstSubPoint[sdim] = 0; 3461 firstSubPoint[0] = firstSubPoint[sdim] + numSubPoints[sdim]; 3462 if (sdim > 1) firstSubPoint[sdim - 1] = firstSubPoint[0] + numSubPoints[0]; 3463 if (sdim > 2) firstSubPoint[sdim - 2] = firstSubPoint[sdim - 1] + numSubPoints[sdim - 1]; 3464 for (d = 0; d <= sdim; ++d) { 3465 PetscCall(DMLabelGetStratumIS(subpointMap, d, &subpointIS[d])); 3466 if (subpointIS[d]) PetscCall(ISGetIndices(subpointIS[d], &subpoints[d])); 3467 } 3468 /* We do not want this label automatically computed, instead we compute it here */ 3469 PetscCall(DMCreateLabel(subdm, "celltype")); 3470 for (d = 0; d <= sdim; ++d) { 3471 for (p = 0; p < numSubPoints[d]; ++p) { 3472 const PetscInt point = subpoints[d][p]; 3473 const PetscInt subpoint = firstSubPoint[d] + p; 3474 const PetscInt *cone; 3475 PetscInt coneSize; 3476 3477 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 3478 if (cellHeight && (d == sdim)) { 3479 PetscInt coneSizeNew, c, val; 3480 3481 PetscCall(DMPlexGetCone(dm, point, &cone)); 3482 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3483 PetscCall(DMLabelGetValue(subpointMap, cone[c], &val)); 3484 if (val >= 0) coneSizeNew++; 3485 } 3486 PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSizeNew)); 3487 PetscCall(DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST)); 3488 } else { 3489 DMPolytopeType ct; 3490 3491 PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSize)); 3492 PetscCall(DMPlexGetCellType(dm, point, &ct)); 3493 PetscCall(DMPlexSetCellType(subdm, subpoint, ct)); 3494 } 3495 } 3496 } 3497 PetscCall(DMLabelDestroy(&subpointMap)); 3498 PetscCall(DMSetUp(subdm)); 3499 /* Set cones */ 3500 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3501 PetscCall(PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew)); 3502 for (d = 0; d <= sdim; ++d) { 3503 for (p = 0; p < numSubPoints[d]; ++p) { 3504 const PetscInt point = subpoints[d][p]; 3505 const PetscInt subpoint = firstSubPoint[d] + p; 3506 const PetscInt *cone, *ornt; 3507 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 3508 3509 if (d == sdim - 1) { 3510 const PetscInt *support, *cone, *ornt; 3511 PetscInt supportSize, coneSize, s, subc; 3512 3513 PetscCall(DMPlexGetSupport(dm, point, &support)); 3514 PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); 3515 for (s = 0; s < supportSize; ++s) { 3516 PetscBool isHybrid = PETSC_FALSE; 3517 3518 PetscCall(DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid)); 3519 if (!isHybrid) continue; 3520 PetscCall(PetscFindInt(support[s], numSubPoints[d + 1], subpoints[d + 1], &subc)); 3521 if (subc >= 0) { 3522 const PetscInt ccell = subpoints[d + 1][subc]; 3523 3524 PetscCall(DMPlexGetCone(dm, ccell, &cone)); 3525 PetscCall(DMPlexGetConeSize(dm, ccell, &coneSize)); 3526 PetscCall(DMPlexGetConeOrientation(dm, ccell, &ornt)); 3527 for (c = 0; c < coneSize; ++c) { 3528 if (cone[c] == point) { 3529 fornt = ornt[c]; 3530 break; 3531 } 3532 } 3533 break; 3534 } 3535 } 3536 } 3537 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 3538 PetscCall(DMPlexGetConeSize(subdm, subpoint, &subconeSize)); 3539 PetscCall(DMPlexGetCone(dm, point, &cone)); 3540 PetscCall(DMPlexGetConeOrientation(dm, point, &ornt)); 3541 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3542 PetscCall(PetscFindInt(cone[c], numSubPoints[d - 1], subpoints[d - 1], &subc)); 3543 if (subc >= 0) { 3544 coneNew[coneSizeNew] = firstSubPoint[d - 1] + subc; 3545 orntNew[coneSizeNew] = ornt[c]; 3546 ++coneSizeNew; 3547 } 3548 } 3549 PetscCheck(coneSizeNew == subconeSize, comm, PETSC_ERR_PLIB, "Number of cone points located %" PetscInt_FMT " does not match subcone size %" PetscInt_FMT, coneSizeNew, subconeSize); 3550 PetscCall(DMPlexSetCone(subdm, subpoint, coneNew)); 3551 PetscCall(DMPlexSetConeOrientation(subdm, subpoint, orntNew)); 3552 if (fornt < 0) PetscCall(DMPlexOrientPoint(subdm, subpoint, fornt)); 3553 } 3554 } 3555 PetscCall(PetscFree2(coneNew, orntNew)); 3556 PetscCall(DMPlexSymmetrize(subdm)); 3557 PetscCall(DMPlexStratify(subdm)); 3558 /* Build coordinates */ 3559 for (coordinate_type = 0; coordinate_type < 2; ++coordinate_type) { 3560 DM coordDM = NULL, subCoordDM = NULL; 3561 PetscSection coordSection = NULL, subCoordSection = NULL; 3562 Vec coordinates = NULL, subCoordinates = NULL; 3563 PetscScalar *coords = NULL, *subCoords = NULL; 3564 PetscInt bs, numComp, coordSize, firstP, lastP, firstSubP = totSubPoints, lastSubP = -1, numFields; 3565 const char *name; 3566 PetscBool localized = (PetscBool)coordinate_type; 3567 3568 if (!dm->coordinates[coordinate_type].dm) continue; 3569 if (!localized) { 3570 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3571 PetscCall(DMGetCoordinateDM(subdm, &subCoordDM)); 3572 } else { 3573 { 3574 PetscInt localizationHeight; 3575 PetscBool sparseLocalize; 3576 3577 PetscCall(DMGetSparseLocalize(dm, &sparseLocalize)); 3578 PetscCall(DMSetSparseLocalize(subdm, sparseLocalize)); 3579 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3580 PetscCall(DMGetCoordinateDM(subdm, &subCoordDM)); 3581 PetscCall(DMPlexGetMaxProjectionHeight(coordDM, &localizationHeight)); 3582 PetscCall(DMPlexSetMaxProjectionHeight(subCoordDM, localizationHeight - (dim - sdim) - cellHeight)); 3583 PetscUseTypeMethod(subdm, createcellcoordinatedm, &subCoordDM); 3584 PetscCall(DMSetCellCoordinateDM(subdm, subCoordDM)); 3585 PetscCall(DMDestroy(&subCoordDM)); 3586 } 3587 PetscCall(DMGetCellCoordinateDM(dm, &coordDM)); 3588 PetscCall(DMGetCellCoordinateDM(subdm, &subCoordDM)); 3589 } 3590 PetscCall(DMGetNumFields(coordDM, &numFields)); 3591 if (numFields > 0) { 3592 PetscFE fe = NULL; 3593 PetscSpace P = NULL; 3594 PetscClassId id; 3595 PetscInt degree; 3596 3597 PetscCall(DMGetField(coordDM, 0, NULL, (PetscObject *)&fe)); 3598 PetscCall(PetscObjectGetClassId((PetscObject)fe, &id)); 3599 if (id == PETSCFE_CLASSID) { 3600 if (sdim == dim && cellHeight == 0) { 3601 /* TODO: Handle Field labels correctly */ 3602 PetscCall(DMSetField(subCoordDM, 0, NULL, (PetscObject)fe)); 3603 PetscCall(DMCreateDS(subCoordDM)); 3604 } else { 3605 /* TODO: Reconstruct the lower-dimensional FE more robustly */ 3606 PetscCall(PetscFEGetBasisSpace(fe, &P)); 3607 PetscCall(PetscSpaceGetDegree(P, °ree, NULL)); 3608 PetscCall(DMPlexCreateCoordinateSpace(subdm, degree, localized, PETSC_FALSE)); 3609 } 3610 } 3611 } 3612 if (!localized) { 3613 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3614 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3615 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3616 } else { 3617 PetscCall(DMGetCellCoordinateSection(dm, &coordSection)); 3618 PetscCall(DMGetCellCoordinateSection(subdm, &subCoordSection)); 3619 PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates)); 3620 } 3621 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3622 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3623 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3624 PetscCall(PetscSectionGetChart(coordSection, &firstP, &lastP)); 3625 for (d = 0; d <= sdim; ++d) { 3626 for (p = 0; p < numSubPoints[d]; ++p) { 3627 const PetscInt point = subpoints[d][p]; 3628 const PetscInt subpoint = firstSubPoint[d] + p; 3629 3630 if (point >= firstP && point < lastP) { 3631 firstSubP = PetscMin(firstSubP, subpoint); 3632 lastSubP = PetscMax(lastSubP, subpoint); 3633 } 3634 } 3635 } 3636 lastSubP += 1; 3637 if (firstSubP == totSubPoints) { 3638 /* Zero if there is no coordinate point. */ 3639 firstSubP = 0; 3640 lastSubP = 0; 3641 } 3642 PetscCall(PetscSectionSetChart(subCoordSection, firstSubP, lastSubP)); 3643 for (d = 0; d <= sdim; ++d) { 3644 for (p = 0; p < numSubPoints[d]; ++p) { 3645 const PetscInt point = subpoints[d][p]; 3646 const PetscInt subpoint = firstSubPoint[d] + p; 3647 PetscInt dof; 3648 3649 if (point >= firstP && point < lastP) { 3650 PetscCall(PetscSectionGetDof(coordSection, point, &dof)); 3651 PetscCall(PetscSectionSetDof(subCoordSection, subpoint, dof)); 3652 PetscCall(PetscSectionSetFieldDof(subCoordSection, subpoint, 0, dof)); 3653 } 3654 } 3655 } 3656 PetscCall(PetscSectionSetUp(subCoordSection)); 3657 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3658 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3659 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 3660 PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name)); 3661 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3662 PetscCall(VecGetBlockSize(coordinates, &bs)); 3663 PetscCall(VecSetBlockSize(subCoordinates, bs)); 3664 PetscCall(VecSetType(subCoordinates, VECSTANDARD)); 3665 PetscCall(VecGetArray(coordinates, &coords)); 3666 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3667 for (d = 0; d <= sdim; ++d) { 3668 for (p = 0; p < numSubPoints[d]; ++p) { 3669 const PetscInt point = subpoints[d][p]; 3670 const PetscInt subpoint = firstSubPoint[d] + p; 3671 PetscInt dof, off, sdof, soff, d; 3672 3673 if (point >= firstP && point < lastP) { 3674 PetscCall(PetscSectionGetDof(coordSection, point, &dof)); 3675 PetscCall(PetscSectionGetOffset(coordSection, point, &off)); 3676 PetscCall(PetscSectionGetDof(subCoordSection, subpoint, &sdof)); 3677 PetscCall(PetscSectionGetOffset(subCoordSection, subpoint, &soff)); 3678 PetscCheck(dof == sdof, comm, PETSC_ERR_PLIB, "Coordinate dimension %" PetscInt_FMT " on subpoint %" PetscInt_FMT ", point %" PetscInt_FMT " should be %" PetscInt_FMT, sdof, subpoint, point, dof); 3679 for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d]; 3680 } 3681 } 3682 } 3683 PetscCall(VecRestoreArray(coordinates, &coords)); 3684 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3685 switch (coordinate_type) { 3686 case 0: 3687 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3688 break; 3689 case 1: 3690 PetscCall(DMSetCellCoordinatesLocal(subdm, subCoordinates)); 3691 break; 3692 default: 3693 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "coordinate_type must be {0, 1}"); 3694 } 3695 PetscCall(VecDestroy(&subCoordinates)); 3696 } 3697 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3698 { 3699 PetscSF sfPoint, sfPointSub; 3700 IS subpIS; 3701 const PetscSFNode *remotePoints; 3702 PetscSFNode *sremotePoints = NULL, *newLocalPoints = NULL, *newOwners = NULL; 3703 const PetscInt *localPoints, *subpoints, *rootdegree; 3704 PetscInt *slocalPoints = NULL, *sortedPoints = NULL, *sortedIndices = NULL; 3705 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl = 0, ll = 0, pStart, pEnd, p; 3706 PetscMPIInt rank, size; 3707 3708 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 3709 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 3710 PetscCall(DMGetPointSF(dm, &sfPoint)); 3711 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 3712 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3713 PetscCall(DMPlexGetChart(subdm, NULL, &numSubroots)); 3714 PetscCall(DMPlexGetSubpointIS(subdm, &subpIS)); 3715 if (subpIS) { 3716 PetscBool sorted = PETSC_TRUE; 3717 3718 PetscCall(ISGetIndices(subpIS, &subpoints)); 3719 PetscCall(ISGetLocalSize(subpIS, &numSubpoints)); 3720 for (p = 1; p < numSubpoints; ++p) sorted = sorted && (subpoints[p] >= subpoints[p - 1]) ? PETSC_TRUE : PETSC_FALSE; 3721 if (!sorted) { 3722 PetscCall(PetscMalloc2(numSubpoints, &sortedPoints, numSubpoints, &sortedIndices)); 3723 for (p = 0; p < numSubpoints; ++p) sortedIndices[p] = p; 3724 PetscCall(PetscArraycpy(sortedPoints, subpoints, numSubpoints)); 3725 PetscCall(PetscSortIntWithArray(numSubpoints, sortedPoints, sortedIndices)); 3726 } 3727 } 3728 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3729 if (numRoots >= 0) { 3730 PetscCall(PetscSFComputeDegreeBegin(sfPoint, &rootdegree)); 3731 PetscCall(PetscSFComputeDegreeEnd(sfPoint, &rootdegree)); 3732 PetscCall(PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners)); 3733 for (p = 0; p < pEnd - pStart; ++p) { 3734 newLocalPoints[p].rank = -2; 3735 newLocalPoints[p].index = -2; 3736 } 3737 for (p = pStart; p < pEnd; ++p) { 3738 newOwners[p - pStart].rank = -3; 3739 newOwners[p - pStart].index = -3; 3740 } 3741 if (sanitizeSubmesh) { 3742 /* A subpoint is forced to be owned by a rank that owns */ 3743 /* a subcell that contains the subpoint in its closure. */ 3744 PetscInt cStart, cEnd, c, clSize, cl; 3745 PetscInt *ownedCells, *closure = NULL; 3746 3747 /* claim ownership */ 3748 for (p = 0; p < numSubpoints; ++p) { 3749 const PetscInt point = subpoints[p]; 3750 3751 newLocalPoints[point - pStart].rank = rank; 3752 newLocalPoints[point - pStart].index = p; 3753 } 3754 PetscCall(DMGetDimension(subdm, &sdim)); 3755 PetscCall(DMPlexGetDepthStratum(dm, sdim, &cStart, &cEnd)); 3756 PetscCall(PetscMalloc1(cEnd - cStart, &ownedCells)); 3757 for (c = cStart; c < cEnd; ++c) ownedCells[c - cStart] = 0; 3758 for (p = 0; p < numSubpoints; ++p) { 3759 c = subpoints[p]; 3760 if (c >= cStart && c < cEnd) ownedCells[c - cStart] = 1; 3761 } 3762 for (l = 0; l < numLeaves; ++l) { 3763 c = localPoints ? localPoints[l] : l; 3764 if (c >= cStart && c < cEnd) ownedCells[c - cStart] = 0; 3765 } 3766 for (c = cStart; c < cEnd; ++c) { 3767 if (ownedCells[c - cStart] == 0) continue; 3768 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 3769 for (cl = 0; cl < clSize * 2; cl += 2) { 3770 p = closure[cl]; 3771 if (newLocalPoints[p - pStart].rank < size) newLocalPoints[p - pStart].rank += size; 3772 } 3773 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 3774 } 3775 PetscCall(PetscFree(ownedCells)); 3776 for (p = 0; p < numRoots; ++p) { 3777 newOwners[p].rank = newLocalPoints[p].rank; 3778 newOwners[p].index = newLocalPoints[p].index; 3779 } 3780 } else { 3781 /* Set subleaves */ 3782 for (l = 0; l < numLeaves; ++l) { 3783 const PetscInt point = localPoints[l]; 3784 const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices); 3785 3786 if (subpoint < 0) continue; 3787 newLocalPoints[point - pStart].rank = rank; 3788 newLocalPoints[point - pStart].index = subpoint; 3789 ++numSubleaves; 3790 } 3791 /* Must put in owned subpoints */ 3792 for (p = 0; p < numSubpoints; ++p) { 3793 /* Hold on to currently owned points */ 3794 if (rootdegree[subpoints[p] - pStart]) newOwners[subpoints[p] - pStart].rank = rank + size; 3795 else newOwners[subpoints[p] - pStart].rank = rank; 3796 newOwners[subpoints[p] - pStart].index = p; 3797 } 3798 } 3799 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 3800 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 3801 for (p = pStart; p < pEnd; ++p) 3802 if (newOwners[p - pStart].rank >= size) newOwners[p - pStart].rank -= size; 3803 if (ownershipTransferSF) { 3804 PetscSFNode *iremote1 = NULL, *newOwners1 = NULL; 3805 PetscInt *ilocal1 = NULL; 3806 PetscInt nleaves1, point; 3807 3808 for (p = 0; p < numSubpoints; ++p) { 3809 point = subpoints[p]; 3810 newLocalPoints[point - pStart].index = point; 3811 } 3812 PetscCall(PetscMalloc1(numRoots, &newOwners1)); 3813 for (p = 0; p < numRoots; ++p) { 3814 newOwners1[p].rank = -1; 3815 newOwners1[p].index = -1; 3816 } 3817 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners1, MPI_MAXLOC)); 3818 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners1, MPI_MAXLOC)); 3819 for (p = 0, nleaves1 = 0; p < numRoots; ++p) { 3820 if (newOwners[p].rank >= 0 && newOwners[p].rank != rank) { ++nleaves1; } 3821 } 3822 PetscCall(PetscMalloc1(nleaves1, &ilocal1)); 3823 PetscCall(PetscMalloc1(nleaves1, &iremote1)); 3824 for (p = 0, nleaves1 = 0; p < numRoots; ++p) { 3825 if (newOwners[p].rank >= 0 && newOwners[p].rank != rank) { 3826 ilocal1[nleaves1] = pStart + p; 3827 iremote1[nleaves1].rank = newOwners[p].rank; 3828 iremote1[nleaves1].index = newOwners1[p].index; 3829 ++nleaves1; 3830 } 3831 } 3832 PetscCall(PetscFree(newOwners1)); 3833 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), ownershipTransferSF)); 3834 PetscCall(PetscSFSetFromOptions(*ownershipTransferSF)); 3835 PetscCall(PetscSFSetGraph(*ownershipTransferSF, pEnd - pStart, nleaves1, ilocal1, PETSC_OWN_POINTER, iremote1, PETSC_OWN_POINTER)); 3836 } 3837 if (sanitizeSubmesh) { 3838 for (p = pStart; p < pEnd; ++p) { 3839 newLocalPoints[p - pStart].rank = newOwners[p - pStart].rank; 3840 newLocalPoints[p - pStart].index = newOwners[p - pStart].index; 3841 } 3842 } 3843 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 3844 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 3845 if (sanitizeSubmesh) { 3846 for (p = 0; p < numSubpoints; ++p) { 3847 const PetscInt point = subpoints[p]; 3848 3849 if (newLocalPoints[point - pStart].rank >= 0 && newLocalPoints[point - pStart].rank != rank) ++sl; 3850 } 3851 PetscCall(PetscMalloc1(sl, &slocalPoints)); 3852 PetscCall(PetscMalloc1(sl, &sremotePoints)); 3853 for (p = 0, sl = 0; p < numSubpoints; ++p) { 3854 const PetscInt point = subpoints[p]; 3855 3856 if (newLocalPoints[point - pStart].rank >= 0 && newLocalPoints[point - pStart].rank != rank) { 3857 slocalPoints[sl] = p; 3858 sremotePoints[sl].rank = newLocalPoints[point - pStart].rank; 3859 sremotePoints[sl].index = newLocalPoints[point - pStart].index; 3860 PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank"); 3861 PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint"); 3862 ++sl; 3863 } 3864 } 3865 } else { 3866 PetscCall(PetscMalloc1(numSubleaves, &slocalPoints)); 3867 PetscCall(PetscMalloc1(numSubleaves, &sremotePoints)); 3868 for (l = 0; l < numLeaves; ++l) { 3869 const PetscInt point = localPoints[l]; 3870 const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices); 3871 3872 if (subpoint < 0) continue; 3873 if (newLocalPoints[point].rank == rank) { 3874 ++ll; 3875 continue; 3876 } 3877 slocalPoints[sl] = subpoint; 3878 sremotePoints[sl].rank = newLocalPoints[point].rank; 3879 sremotePoints[sl].index = newLocalPoints[point].index; 3880 PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 3881 PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 3882 ++sl; 3883 } 3884 PetscCheck(sl + ll == numSubleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubleaves); 3885 } 3886 PetscCall(PetscFree2(newLocalPoints, newOwners)); 3887 PetscCall(PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3888 } 3889 if (subpIS) PetscCall(ISRestoreIndices(subpIS, &subpoints)); 3890 PetscCall(PetscFree2(sortedPoints, sortedIndices)); 3891 } 3892 /* Filter labels */ 3893 PetscCall(DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm)); 3894 /* Cleanup */ 3895 for (d = 0; d <= sdim; ++d) { 3896 if (subpointIS[d]) PetscCall(ISRestoreIndices(subpointIS[d], &subpoints[d])); 3897 PetscCall(ISDestroy(&subpointIS[d])); 3898 } 3899 PetscCall(PetscFree4(numSubPoints, firstSubPoint, subpointIS, subpoints)); 3900 PetscFunctionReturn(PETSC_SUCCESS); 3901 } 3902 3903 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm) 3904 { 3905 PetscFunctionBegin; 3906 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm)); 3907 PetscFunctionReturn(PETSC_SUCCESS); 3908 } 3909 3910 /*@ 3911 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3912 3913 Input Parameters: 3914 + dm - The original mesh 3915 . vertexLabel - The `DMLabel` marking points contained in the surface 3916 . value - The label value to use 3917 - markedFaces - `PETSC_TRUE` if surface faces are marked in addition to vertices, `PETSC_FALSE` if only vertices are marked 3918 3919 Output Parameter: 3920 . subdm - The surface mesh 3921 3922 Level: developer 3923 3924 Note: 3925 This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`. 3926 3927 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()` 3928 @*/ 3929 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm) 3930 { 3931 DMPlexInterpolatedFlag interpolated; 3932 PetscInt dim, cdim; 3933 3934 PetscFunctionBegin; 3935 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3936 PetscAssertPointer(subdm, 5); 3937 PetscCall(DMGetDimension(dm, &dim)); 3938 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3939 PetscCall(DMSetType(*subdm, DMPLEX)); 3940 PetscCall(DMSetDimension(*subdm, dim - 1)); 3941 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3942 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 3943 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3944 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 3945 if (interpolated) { 3946 PetscCall(DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm)); 3947 } else { 3948 PetscCall(DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm)); 3949 } 3950 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3951 PetscFunctionReturn(PETSC_SUCCESS); 3952 } 3953 3954 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3955 { 3956 MPI_Comm comm; 3957 DMLabel subpointMap; 3958 IS subvertexIS; 3959 const PetscInt *subVertices; 3960 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3961 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3962 PetscInt c, f; 3963 3964 PetscFunctionBegin; 3965 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3966 /* Create subpointMap which marks the submesh */ 3967 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3968 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3969 PetscCall(DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm)); 3970 /* Setup chart */ 3971 PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3972 PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3973 PetscCall(DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices)); 3974 PetscCall(DMPlexSetVTKCellHeight(subdm, 1)); 3975 /* Set cone sizes */ 3976 firstSubVertex = numSubCells; 3977 firstSubFace = numSubCells + numSubVertices; 3978 newFacePoint = firstSubFace; 3979 PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3980 if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices)); 3981 for (c = 0; c < numSubCells; ++c) PetscCall(DMPlexSetConeSize(subdm, c, 1)); 3982 for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) PetscCall(DMPlexSetConeSize(subdm, f, nFV)); 3983 PetscCall(DMSetUp(subdm)); 3984 PetscCall(DMLabelDestroy(&subpointMap)); 3985 /* Create face cones */ 3986 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3987 PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 3988 for (c = 0; c < numSubCells; ++c) { 3989 const PetscInt cell = subCells[c]; 3990 const PetscInt subcell = c; 3991 const PetscInt *cone, *cells; 3992 PetscBool isHybrid = PETSC_FALSE; 3993 PetscInt numCells, subVertex, p, v; 3994 3995 PetscCall(DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid)); 3996 if (!isHybrid) continue; 3997 PetscCall(DMPlexGetCone(dm, cell, &cone)); 3998 for (v = 0; v < nFV; ++v) { 3999 PetscCall(PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex)); 4000 subface[v] = firstSubVertex + subVertex; 4001 } 4002 PetscCall(DMPlexSetCone(subdm, newFacePoint, subface)); 4003 PetscCall(DMPlexSetCone(subdm, subcell, &newFacePoint)); 4004 PetscCall(DMPlexGetJoin(dm, nFV, cone, &numCells, &cells)); 4005 /* Not true in parallel 4006 PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 4007 for (p = 0; p < numCells; ++p) { 4008 PetscInt negsubcell; 4009 PetscBool isHybrid = PETSC_FALSE; 4010 4011 PetscCall(DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid)); 4012 if (isHybrid) continue; 4013 /* I know this is a crap search */ 4014 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 4015 if (subCells[negsubcell] == cells[p]) break; 4016 } 4017 PetscCheck(negsubcell != numSubCells, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %" PetscInt_FMT, cell); 4018 PetscCall(DMPlexSetCone(subdm, negsubcell, &newFacePoint)); 4019 } 4020 PetscCall(DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells)); 4021 ++newFacePoint; 4022 } 4023 PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 4024 PetscCall(DMPlexSymmetrize(subdm)); 4025 PetscCall(DMPlexStratify(subdm)); 4026 /* Build coordinates */ 4027 { 4028 PetscSection coordSection, subCoordSection; 4029 Vec coordinates, subCoordinates; 4030 PetscScalar *coords, *subCoords; 4031 PetscInt cdim, numComp, coordSize, v; 4032 const char *name; 4033 4034 PetscCall(DMGetCoordinateDim(dm, &cdim)); 4035 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4036 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 4037 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 4038 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 4039 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 4040 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 4041 PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices)); 4042 for (v = 0; v < numSubVertices; ++v) { 4043 const PetscInt vertex = subVertices[v]; 4044 const PetscInt subvertex = firstSubVertex + v; 4045 PetscInt dof; 4046 4047 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 4048 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 4049 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 4050 } 4051 PetscCall(PetscSectionSetUp(subCoordSection)); 4052 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 4053 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 4054 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 4055 PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name)); 4056 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 4057 PetscCall(VecSetBlockSize(subCoordinates, cdim)); 4058 PetscCall(VecSetType(subCoordinates, VECSTANDARD)); 4059 PetscCall(VecGetArray(coordinates, &coords)); 4060 PetscCall(VecGetArray(subCoordinates, &subCoords)); 4061 for (v = 0; v < numSubVertices; ++v) { 4062 const PetscInt vertex = subVertices[v]; 4063 const PetscInt subvertex = firstSubVertex + v; 4064 PetscInt dof, off, sdof, soff, d; 4065 4066 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 4067 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 4068 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 4069 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 4070 PetscCheck(dof == sdof, comm, PETSC_ERR_PLIB, "Coordinate dimension %" PetscInt_FMT " on subvertex %" PetscInt_FMT ", vertex %" PetscInt_FMT " should be %" PetscInt_FMT, sdof, subvertex, vertex, dof); 4071 for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d]; 4072 } 4073 PetscCall(VecRestoreArray(coordinates, &coords)); 4074 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 4075 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 4076 PetscCall(VecDestroy(&subCoordinates)); 4077 } 4078 /* Build SF */ 4079 CHKMEMQ; 4080 { 4081 PetscSF sfPoint, sfPointSub; 4082 const PetscSFNode *remotePoints; 4083 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 4084 const PetscInt *localPoints; 4085 PetscInt *slocalPoints; 4086 PetscInt numRoots, numLeaves, numSubRoots = numSubCells + numSubFaces + numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 4087 PetscMPIInt rank; 4088 4089 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 4090 PetscCall(DMGetPointSF(dm, &sfPoint)); 4091 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 4092 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 4093 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4094 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 4095 if (numRoots >= 0) { 4096 /* Only vertices should be shared */ 4097 PetscCall(PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners)); 4098 for (p = 0; p < pEnd - pStart; ++p) { 4099 newLocalPoints[p].rank = -2; 4100 newLocalPoints[p].index = -2; 4101 } 4102 /* Set subleaves */ 4103 for (l = 0; l < numLeaves; ++l) { 4104 const PetscInt point = localPoints[l]; 4105 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 4106 4107 PetscCheck(!(point < vStart) || !(point >= vEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %" PetscInt_FMT, point); 4108 if (subPoint < 0) continue; 4109 newLocalPoints[point - pStart].rank = rank; 4110 newLocalPoints[point - pStart].index = subPoint; 4111 ++numSubLeaves; 4112 } 4113 /* Must put in owned subpoints */ 4114 for (p = pStart; p < pEnd; ++p) { 4115 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 4116 4117 if (subPoint < 0) { 4118 newOwners[p - pStart].rank = -3; 4119 newOwners[p - pStart].index = -3; 4120 } else { 4121 newOwners[p - pStart].rank = rank; 4122 newOwners[p - pStart].index = subPoint; 4123 } 4124 } 4125 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 4126 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 4127 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 4128 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 4129 PetscCall(PetscMalloc1(numSubLeaves, &slocalPoints)); 4130 PetscCall(PetscMalloc1(numSubLeaves, &sremotePoints)); 4131 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 4132 const PetscInt point = localPoints[l]; 4133 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 4134 4135 if (subPoint < 0) continue; 4136 if (newLocalPoints[point].rank == rank) { 4137 ++ll; 4138 continue; 4139 } 4140 slocalPoints[sl] = subPoint; 4141 sremotePoints[sl].rank = newLocalPoints[point].rank; 4142 sremotePoints[sl].index = newLocalPoints[point].index; 4143 PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 4144 PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 4145 ++sl; 4146 } 4147 PetscCall(PetscFree2(newLocalPoints, newOwners)); 4148 PetscCheck(sl + ll == numSubLeaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubLeaves); 4149 PetscCall(PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 4150 } 4151 } 4152 CHKMEMQ; 4153 /* Cleanup */ 4154 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices)); 4155 PetscCall(ISDestroy(&subvertexIS)); 4156 PetscCall(PetscFree(subCells)); 4157 PetscFunctionReturn(PETSC_SUCCESS); 4158 } 4159 4160 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 4161 { 4162 DMLabel label = NULL; 4163 4164 PetscFunctionBegin; 4165 if (labelname) PetscCall(DMGetLabel(dm, labelname, &label)); 4166 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm)); 4167 PetscFunctionReturn(PETSC_SUCCESS); 4168 } 4169 4170 /*@ 4171 DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a label can be given to restrict the cells. 4172 4173 Input Parameters: 4174 + dm - The original mesh 4175 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 4176 . label - A label name, or `NULL` 4177 - value - A label value 4178 4179 Output Parameter: 4180 . subdm - The surface mesh 4181 4182 Level: developer 4183 4184 Note: 4185 This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`. 4186 4187 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()` 4188 @*/ 4189 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 4190 { 4191 PetscInt dim, cdim, depth; 4192 4193 PetscFunctionBegin; 4194 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4195 PetscAssertPointer(subdm, 5); 4196 PetscCall(DMGetDimension(dm, &dim)); 4197 PetscCall(DMPlexGetDepth(dm, &depth)); 4198 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 4199 PetscCall(DMSetType(*subdm, DMPLEX)); 4200 PetscCall(DMSetDimension(*subdm, dim - 1)); 4201 PetscCall(DMGetCoordinateDim(dm, &cdim)); 4202 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 4203 if (depth == dim) { 4204 PetscCall(DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm)); 4205 } else { 4206 PetscCall(DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm)); 4207 } 4208 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 4209 PetscFunctionReturn(PETSC_SUCCESS); 4210 } 4211 4212 /*@ 4213 DMPlexReorderCohesiveSupports - Ensure that face supports for cohesive end caps are ordered 4214 4215 Not Collective 4216 4217 Input Parameter: 4218 . dm - The `DM` containing cohesive cells 4219 4220 Level: developer 4221 4222 Note: 4223 For the negative size (first) face, the cohesive cell should be first in the support, and for 4224 the positive side (second) face, the cohesive cell should be second in the support. 4225 4226 .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexCreateCohesiveSubmesh()` 4227 @*/ 4228 PetscErrorCode DMPlexReorderCohesiveSupports(DM dm) 4229 { 4230 PetscInt dim, cStart, cEnd; 4231 4232 PetscFunctionBegin; 4233 PetscCall(DMGetDimension(dm, &dim)); 4234 PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cStart, &cEnd)); 4235 for (PetscInt c = cStart; c < cEnd; ++c) { 4236 const PetscInt *cone; 4237 PetscInt coneSize; 4238 4239 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 4240 PetscCall(DMPlexGetCone(dm, c, &cone)); 4241 PetscCheck(coneSize >= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hybrid cell %" PetscInt_FMT " cone size %" PetscInt_FMT " must be >= 2", c, coneSize); 4242 for (PetscInt s = 0; s < 2; ++s) { 4243 const PetscInt *supp; 4244 PetscInt suppSize, newsupp[2]; 4245 4246 PetscCall(DMPlexGetSupportSize(dm, cone[s], &suppSize)); 4247 PetscCall(DMPlexGetSupport(dm, cone[s], &supp)); 4248 if (suppSize == 2) { 4249 /* Reorder hybrid end cap support */ 4250 if (supp[s] == c) { 4251 newsupp[0] = supp[1]; 4252 newsupp[1] = supp[0]; 4253 } else { 4254 newsupp[0] = supp[0]; 4255 newsupp[1] = supp[1]; 4256 } 4257 PetscCheck(newsupp[1 - s] == c, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hybrid end cap %" PetscInt_FMT " support entry %" PetscInt_FMT " != %" PetscInt_FMT " hybrid cell", cone[s], newsupp[1 - s], c); 4258 PetscCall(DMPlexSetSupport(dm, cone[s], newsupp)); 4259 } 4260 } 4261 } 4262 PetscFunctionReturn(PETSC_SUCCESS); 4263 } 4264 4265 /*@ 4266 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 4267 4268 Input Parameters: 4269 + dm - The original mesh 4270 . cellLabel - The `DMLabel` marking cells contained in the new mesh 4271 . value - The label value to use 4272 . ignoreLabelHalo - The flag indicating if labeled points that are in the halo are ignored 4273 - sanitizeSubmesh - The flag indicating if a subpoint is forced to be owned by a rank that owns a subcell that contains that point in its closure 4274 4275 Output Parameters: 4276 + ownershipTransferSF - The `PetscSF` representing the ownership transfers between parent local meshes due to submeshing. 4277 - subdm - The new mesh 4278 4279 Level: developer 4280 4281 Note: 4282 This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`. 4283 4284 On a given rank, the leaves of the ownershipTransferSF are the points in the local mesh that this rank gives up ownership of (in the submesh), and 4285 the remote locations for these leaves are tuples of rank and point that represent the new owners. 4286 4287 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`, `DMPlexCreateSubmesh()` 4288 @*/ 4289 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, PetscBool ignoreLabelHalo, PetscBool sanitizeSubmesh, PetscSF *ownershipTransferSF, DM *subdm) 4290 { 4291 PetscInt dim, overlap; 4292 4293 PetscFunctionBegin; 4294 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4295 if (ownershipTransferSF) PetscAssertPointer(ownershipTransferSF, 6); 4296 PetscAssertPointer(subdm, 7); 4297 PetscCall(DMGetDimension(dm, &dim)); 4298 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 4299 PetscCall(DMSetType(*subdm, DMPLEX)); 4300 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 4301 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, ignoreLabelHalo, sanitizeSubmesh, ownershipTransferSF, *subdm)); 4302 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 4303 // It is possible to obtain a surface mesh where some faces are in SF 4304 // We should either mark the mesh as having an overlap, or delete these from the SF 4305 PetscCall(DMPlexGetOverlap(dm, &overlap)); 4306 if (!overlap) { 4307 PetscSF sf; 4308 const PetscInt *leaves; 4309 PetscInt cStart, cEnd, Nl; 4310 PetscBool hasSubcell = PETSC_FALSE, ghasSubcell; 4311 4312 PetscCall(DMPlexGetHeightStratum(*subdm, 0, &cStart, &cEnd)); 4313 PetscCall(DMGetPointSF(*subdm, &sf)); 4314 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL)); 4315 for (PetscInt l = 0; l < Nl; ++l) { 4316 const PetscInt point = leaves ? leaves[l] : l; 4317 4318 if (point >= cStart && point < cEnd) { 4319 hasSubcell = PETSC_TRUE; 4320 break; 4321 } 4322 } 4323 PetscCallMPI(MPIU_Allreduce(&hasSubcell, &ghasSubcell, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm))); 4324 if (ghasSubcell) PetscCall(DMPlexSetOverlap(*subdm, NULL, 1)); 4325 } 4326 PetscFunctionReturn(PETSC_SUCCESS); 4327 } 4328 4329 /*@ 4330 DMPlexGetSubpointMap - Returns a `DMLabel` with point dimension as values 4331 4332 Input Parameter: 4333 . dm - The submesh `DM` 4334 4335 Output Parameter: 4336 . subpointMap - The `DMLabel` of all the points from the original mesh in this submesh, or `NULL` if this is not a submesh 4337 4338 Level: developer 4339 4340 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 4341 @*/ 4342 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 4343 { 4344 PetscFunctionBegin; 4345 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4346 PetscAssertPointer(subpointMap, 2); 4347 *subpointMap = ((DM_Plex *)dm->data)->subpointMap; 4348 PetscFunctionReturn(PETSC_SUCCESS); 4349 } 4350 4351 /*@ 4352 DMPlexSetSubpointMap - Sets the `DMLabel` with point dimension as values 4353 4354 Input Parameters: 4355 + dm - The submesh `DM` 4356 - subpointMap - The `DMLabel` of all the points from the original mesh in this submesh 4357 4358 Level: developer 4359 4360 Note: 4361 Should normally not be called by the user, since it is set in `DMPlexCreateSubmesh()` 4362 4363 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 4364 @*/ 4365 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 4366 { 4367 DM_Plex *mesh = (DM_Plex *)dm->data; 4368 DMLabel tmp; 4369 4370 PetscFunctionBegin; 4371 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4372 tmp = mesh->subpointMap; 4373 mesh->subpointMap = subpointMap; 4374 PetscCall(PetscObjectReference((PetscObject)mesh->subpointMap)); 4375 PetscCall(DMLabelDestroy(&tmp)); 4376 PetscFunctionReturn(PETSC_SUCCESS); 4377 } 4378 4379 static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS) 4380 { 4381 DMLabel spmap; 4382 PetscInt depth, d; 4383 4384 PetscFunctionBegin; 4385 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 4386 PetscCall(DMPlexGetDepth(dm, &depth)); 4387 if (spmap && depth >= 0) { 4388 DM_Plex *mesh = (DM_Plex *)dm->data; 4389 PetscInt *points, *depths; 4390 PetscInt pStart, pEnd, p, off; 4391 4392 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 4393 PetscCheck(!pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %" PetscInt_FMT, pStart); 4394 PetscCall(PetscMalloc1(pEnd, &points)); 4395 PetscCall(DMGetWorkArray(dm, depth + 1, MPIU_INT, &depths)); 4396 depths[0] = depth; 4397 depths[1] = 0; 4398 for (d = 2; d <= depth; ++d) depths[d] = depth + 1 - d; 4399 for (d = 0, off = 0; d <= depth; ++d) { 4400 const PetscInt dep = depths[d]; 4401 PetscInt depStart, depEnd, n; 4402 4403 PetscCall(DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd)); 4404 PetscCall(DMLabelGetStratumSize(spmap, dep, &n)); 4405 if ((d < 2 && depth > 1) || d == 1) { /* Only check vertices and cells for now since the map is broken for others */ 4406 PetscCheck(n == depEnd - depStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %" PetscInt_FMT " at depth %" PetscInt_FMT " should be %" PetscInt_FMT, n, dep, depEnd - depStart); 4407 } else { 4408 if (!n) { 4409 if (d == 0) { 4410 /* Missing cells */ 4411 for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = -1; 4412 } else { 4413 /* Missing faces */ 4414 for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = PETSC_INT_MAX; 4415 } 4416 } 4417 } 4418 if (n) { 4419 IS is; 4420 const PetscInt *opoints; 4421 4422 PetscCall(DMLabelGetStratumIS(spmap, dep, &is)); 4423 PetscCall(ISGetIndices(is, &opoints)); 4424 for (p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 4425 PetscCall(ISRestoreIndices(is, &opoints)); 4426 PetscCall(ISDestroy(&is)); 4427 } 4428 } 4429 PetscCall(DMRestoreWorkArray(dm, depth + 1, MPIU_INT, &depths)); 4430 PetscCheck(off == pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %" PetscInt_FMT " should be %" PetscInt_FMT, off, pEnd); 4431 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS)); 4432 PetscCall(PetscObjectStateGet((PetscObject)spmap, &mesh->subpointState)); 4433 } 4434 PetscFunctionReturn(PETSC_SUCCESS); 4435 } 4436 4437 /*@ 4438 DMPlexGetSubpointIS - Returns an `IS` covering the entire subdm chart with the original points as data 4439 4440 Input Parameter: 4441 . dm - The submesh `DM` 4442 4443 Output Parameter: 4444 . subpointIS - The `IS` of all the points from the original mesh in this submesh, or `NULL` if this is not a submesh 4445 4446 Level: developer 4447 4448 Note: 4449 This `IS` is guaranteed to be sorted by the construction of the submesh. However, if the filtering operation removes an entire stratum, then the strata in the submesh can be in a different order, and the `subpointIS` will only be sorted within each stratum. 4450 4451 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()` 4452 @*/ 4453 PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS) 4454 { 4455 DM_Plex *mesh = (DM_Plex *)dm->data; 4456 DMLabel spmap; 4457 PetscObjectState state; 4458 4459 PetscFunctionBegin; 4460 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4461 PetscAssertPointer(subpointIS, 2); 4462 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 4463 PetscCall(PetscObjectStateGet((PetscObject)spmap, &state)); 4464 if (state != mesh->subpointState || !mesh->subpointIS) PetscCall(DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS)); 4465 *subpointIS = mesh->subpointIS; 4466 PetscFunctionReturn(PETSC_SUCCESS); 4467 } 4468 4469 /*@ 4470 DMGetEnclosureRelation - Get the relationship between `dmA` and `dmB` 4471 4472 Input Parameters: 4473 + dmA - The first `DM` 4474 - dmB - The second `DM` 4475 4476 Output Parameter: 4477 . rel - The relation of `dmA` to `dmB` 4478 4479 Level: intermediate 4480 4481 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMGetEnclosurePoint()` 4482 @*/ 4483 PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel) 4484 { 4485 DM plexA, plexB, sdm; 4486 DMLabel spmap; 4487 PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB; 4488 4489 PetscFunctionBegin; 4490 PetscAssertPointer(rel, 3); 4491 *rel = DM_ENC_NONE; 4492 if (!dmA || !dmB) PetscFunctionReturn(PETSC_SUCCESS); 4493 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 4494 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 4495 if (dmA == dmB) { 4496 *rel = DM_ENC_EQUALITY; 4497 PetscFunctionReturn(PETSC_SUCCESS); 4498 } 4499 PetscCall(DMConvert(dmA, DMPLEX, &plexA)); 4500 PetscCall(DMConvert(dmB, DMPLEX, &plexB)); 4501 PetscCall(DMPlexGetChart(plexA, &pStartA, &pEndA)); 4502 PetscCall(DMPlexGetChart(plexB, &pStartB, &pEndB)); 4503 /* Assumption 1: subDMs have smaller charts than the DMs that they originate from 4504 - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */ 4505 if (pStartA == pStartB && pEndA == pEndB) { 4506 *rel = DM_ENC_EQUALITY; 4507 goto end; 4508 } 4509 NpA = pEndA - pStartA; 4510 NpB = pEndB - pStartB; 4511 if (NpA == NpB) goto end; 4512 sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */ 4513 PetscCall(DMPlexGetSubpointMap(sdm, &spmap)); 4514 if (!spmap) goto end; 4515 /* TODO Check the space mapped to by subpointMap is same size as dm */ 4516 if (NpA > NpB) { 4517 *rel = DM_ENC_SUPERMESH; 4518 } else { 4519 *rel = DM_ENC_SUBMESH; 4520 } 4521 end: 4522 PetscCall(DMDestroy(&plexA)); 4523 PetscCall(DMDestroy(&plexB)); 4524 PetscFunctionReturn(PETSC_SUCCESS); 4525 } 4526 4527 /*@ 4528 DMGetEnclosurePoint - Get the point `pA` in `dmA` which corresponds to the point `pB` in `dmB` 4529 4530 Input Parameters: 4531 + dmA - The first `DM` 4532 . dmB - The second `DM` 4533 . etype - The type of enclosure relation that `dmA` has to `dmB` 4534 - pB - A point of `dmB` 4535 4536 Output Parameter: 4537 . pA - The corresponding point of `dmA` 4538 4539 Level: intermediate 4540 4541 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMGetEnclosureRelation()` 4542 @*/ 4543 PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA) 4544 { 4545 DM sdm; 4546 IS subpointIS; 4547 const PetscInt *subpoints; 4548 PetscInt numSubpoints; 4549 4550 PetscFunctionBegin; 4551 /* TODO Cache the IS, making it look like an index */ 4552 switch (etype) { 4553 case DM_ENC_SUPERMESH: 4554 sdm = dmB; 4555 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4556 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4557 *pA = subpoints[pB]; 4558 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4559 break; 4560 case DM_ENC_SUBMESH: 4561 sdm = dmA; 4562 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4563 PetscCall(ISGetLocalSize(subpointIS, &numSubpoints)); 4564 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4565 PetscCall(PetscFindInt(pB, numSubpoints, subpoints, pA)); 4566 if (*pA < 0) { 4567 PetscCall(DMViewFromOptions(dmA, NULL, "-dm_enc_A_view")); 4568 PetscCall(DMViewFromOptions(dmB, NULL, "-dm_enc_B_view")); 4569 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB); 4570 } 4571 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4572 break; 4573 case DM_ENC_EQUALITY: 4574 case DM_ENC_NONE: 4575 *pA = pB; 4576 break; 4577 case DM_ENC_UNKNOWN: { 4578 DMEnclosureType enc; 4579 4580 PetscCall(DMGetEnclosureRelation(dmA, dmB, &enc)); 4581 PetscCall(DMGetEnclosurePoint(dmA, dmB, enc, pB, pA)); 4582 } break; 4583 default: 4584 SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int)etype); 4585 } 4586 PetscFunctionReturn(PETSC_SUCCESS); 4587 } 4588