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 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2771 posOrient = PETSC_FALSE; 2772 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2773 } 2774 } else if (cellDim == 2 && numCorners == 6) { 2775 /* Quadratic triangle (I hate this) */ 2776 /* Edges are determined by the first 2 vertices (corners of edges) */ 2777 const PetscInt faceSizeTri = 3; 2778 PetscInt sortedIndices[3], i, iFace; 2779 PetscBool found = PETSC_FALSE; 2780 PetscInt faceVerticesTriSorted[9] = { 2781 0, 3, 4, /* bottom */ 2782 1, 4, 5, /* right */ 2783 2, 3, 5, /* left */ 2784 }; 2785 PetscInt faceVerticesTri[9] = { 2786 0, 3, 4, /* bottom */ 2787 1, 4, 5, /* right */ 2788 2, 5, 3, /* left */ 2789 }; 2790 2791 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2792 PetscCall(PetscSortInt(faceSizeTri, sortedIndices)); 2793 for (iFace = 0; iFace < 3; ++iFace) { 2794 const PetscInt ii = iFace * faceSizeTri; 2795 PetscInt fVertex, cVertex; 2796 2797 if ((sortedIndices[0] == faceVerticesTriSorted[ii + 0]) && (sortedIndices[1] == faceVerticesTriSorted[ii + 1])) { 2798 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2799 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2800 if (indices[cVertex] == faceVerticesTri[ii + fVertex]) { 2801 faceVertices[fVertex] = origVertices[cVertex]; 2802 break; 2803 } 2804 } 2805 } 2806 found = PETSC_TRUE; 2807 break; 2808 } 2809 } 2810 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2811 if (posOriented) *posOriented = PETSC_TRUE; 2812 PetscFunctionReturn(PETSC_SUCCESS); 2813 } else if (cellDim == 2 && numCorners == 9) { 2814 /* Quadratic quad (I hate this) */ 2815 /* Edges are determined by the first 2 vertices (corners of edges) */ 2816 const PetscInt faceSizeQuad = 3; 2817 PetscInt sortedIndices[3], i, iFace; 2818 PetscBool found = PETSC_FALSE; 2819 PetscInt faceVerticesQuadSorted[12] = { 2820 0, 1, 4, /* bottom */ 2821 1, 2, 5, /* right */ 2822 2, 3, 6, /* top */ 2823 0, 3, 7, /* left */ 2824 }; 2825 PetscInt faceVerticesQuad[12] = { 2826 0, 1, 4, /* bottom */ 2827 1, 2, 5, /* right */ 2828 2, 3, 6, /* top */ 2829 3, 0, 7, /* left */ 2830 }; 2831 2832 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2833 PetscCall(PetscSortInt(faceSizeQuad, sortedIndices)); 2834 for (iFace = 0; iFace < 4; ++iFace) { 2835 const PetscInt ii = iFace * faceSizeQuad; 2836 PetscInt fVertex, cVertex; 2837 2838 if ((sortedIndices[0] == faceVerticesQuadSorted[ii + 0]) && (sortedIndices[1] == faceVerticesQuadSorted[ii + 1])) { 2839 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2840 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2841 if (indices[cVertex] == faceVerticesQuad[ii + fVertex]) { 2842 faceVertices[fVertex] = origVertices[cVertex]; 2843 break; 2844 } 2845 } 2846 } 2847 found = PETSC_TRUE; 2848 break; 2849 } 2850 } 2851 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2852 if (posOriented) *posOriented = PETSC_TRUE; 2853 PetscFunctionReturn(PETSC_SUCCESS); 2854 } else if (cellDim == 3 && numCorners == 8) { 2855 /* Hexes 2856 A hex is two oriented quads with the normal of the first 2857 pointing up at the second. 2858 2859 7---6 2860 /| /| 2861 4---5 | 2862 | 1-|-2 2863 |/ |/ 2864 0---3 2865 2866 Faces are determined by the first 4 vertices (corners of faces) */ 2867 const PetscInt faceSizeHex = 4; 2868 PetscInt sortedIndices[4], i, iFace; 2869 PetscBool found = PETSC_FALSE; 2870 PetscInt faceVerticesHexSorted[24] = { 2871 0, 1, 2, 3, /* bottom */ 2872 4, 5, 6, 7, /* top */ 2873 0, 3, 4, 5, /* front */ 2874 2, 3, 5, 6, /* right */ 2875 1, 2, 6, 7, /* back */ 2876 0, 1, 4, 7, /* left */ 2877 }; 2878 PetscInt faceVerticesHex[24] = { 2879 1, 2, 3, 0, /* bottom */ 2880 4, 5, 6, 7, /* top */ 2881 0, 3, 5, 4, /* front */ 2882 3, 2, 6, 5, /* right */ 2883 2, 1, 7, 6, /* back */ 2884 1, 0, 4, 7, /* left */ 2885 }; 2886 2887 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2888 PetscCall(PetscSortInt(faceSizeHex, sortedIndices)); 2889 for (iFace = 0; iFace < 6; ++iFace) { 2890 const PetscInt ii = iFace * faceSizeHex; 2891 PetscInt fVertex, cVertex; 2892 2893 if ((sortedIndices[0] == faceVerticesHexSorted[ii + 0]) && (sortedIndices[1] == faceVerticesHexSorted[ii + 1]) && (sortedIndices[2] == faceVerticesHexSorted[ii + 2]) && (sortedIndices[3] == faceVerticesHexSorted[ii + 3])) { 2894 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2895 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2896 if (indices[cVertex] == faceVerticesHex[ii + fVertex]) { 2897 faceVertices[fVertex] = origVertices[cVertex]; 2898 break; 2899 } 2900 } 2901 } 2902 found = PETSC_TRUE; 2903 break; 2904 } 2905 } 2906 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2907 if (posOriented) *posOriented = PETSC_TRUE; 2908 PetscFunctionReturn(PETSC_SUCCESS); 2909 } else if (cellDim == 3 && numCorners == 10) { 2910 /* Quadratic tet */ 2911 /* Faces are determined by the first 3 vertices (corners of faces) */ 2912 const PetscInt faceSizeTet = 6; 2913 PetscInt sortedIndices[6], i, iFace; 2914 PetscBool found = PETSC_FALSE; 2915 PetscInt faceVerticesTetSorted[24] = { 2916 0, 1, 2, 6, 7, 8, /* bottom */ 2917 0, 3, 4, 6, 7, 9, /* front */ 2918 1, 4, 5, 7, 8, 9, /* right */ 2919 2, 3, 5, 6, 8, 9, /* left */ 2920 }; 2921 PetscInt faceVerticesTet[24] = { 2922 0, 1, 2, 6, 7, 8, /* bottom */ 2923 0, 4, 3, 6, 7, 9, /* front */ 2924 1, 5, 4, 7, 8, 9, /* right */ 2925 2, 3, 5, 8, 6, 9, /* left */ 2926 }; 2927 2928 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2929 PetscCall(PetscSortInt(faceSizeTet, sortedIndices)); 2930 for (iFace = 0; iFace < 4; ++iFace) { 2931 const PetscInt ii = iFace * faceSizeTet; 2932 PetscInt fVertex, cVertex; 2933 2934 if ((sortedIndices[0] == faceVerticesTetSorted[ii + 0]) && (sortedIndices[1] == faceVerticesTetSorted[ii + 1]) && (sortedIndices[2] == faceVerticesTetSorted[ii + 2]) && (sortedIndices[3] == faceVerticesTetSorted[ii + 3])) { 2935 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2936 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2937 if (indices[cVertex] == faceVerticesTet[ii + fVertex]) { 2938 faceVertices[fVertex] = origVertices[cVertex]; 2939 break; 2940 } 2941 } 2942 } 2943 found = PETSC_TRUE; 2944 break; 2945 } 2946 } 2947 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2948 if (posOriented) *posOriented = PETSC_TRUE; 2949 PetscFunctionReturn(PETSC_SUCCESS); 2950 } else if (cellDim == 3 && numCorners == 27) { 2951 /* Quadratic hexes (I hate this) 2952 A hex is two oriented quads with the normal of the first 2953 pointing up at the second. 2954 2955 7---6 2956 /| /| 2957 4---5 | 2958 | 3-|-2 2959 |/ |/ 2960 0---1 2961 2962 Faces are determined by the first 4 vertices (corners of faces) */ 2963 const PetscInt faceSizeQuadHex = 9; 2964 PetscInt sortedIndices[9], i, iFace; 2965 PetscBool found = PETSC_FALSE; 2966 PetscInt faceVerticesQuadHexSorted[54] = { 2967 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2968 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2969 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2970 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2971 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2972 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2973 }; 2974 PetscInt faceVerticesQuadHex[54] = { 2975 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2976 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2977 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2978 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2979 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2980 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2981 }; 2982 2983 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2984 PetscCall(PetscSortInt(faceSizeQuadHex, sortedIndices)); 2985 for (iFace = 0; iFace < 6; ++iFace) { 2986 const PetscInt ii = iFace * faceSizeQuadHex; 2987 PetscInt fVertex, cVertex; 2988 2989 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii + 0]) && (sortedIndices[1] == faceVerticesQuadHexSorted[ii + 1]) && (sortedIndices[2] == faceVerticesQuadHexSorted[ii + 2]) && (sortedIndices[3] == faceVerticesQuadHexSorted[ii + 3])) { 2990 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2991 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2992 if (indices[cVertex] == faceVerticesQuadHex[ii + fVertex]) { 2993 faceVertices[fVertex] = origVertices[cVertex]; 2994 break; 2995 } 2996 } 2997 } 2998 found = PETSC_TRUE; 2999 break; 3000 } 3001 } 3002 PetscCheck(found, comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 3003 if (posOriented) *posOriented = PETSC_TRUE; 3004 PetscFunctionReturn(PETSC_SUCCESS); 3005 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 3006 if (!posOrient) { 3007 if (debug) PetscCall(PetscPrintf(comm, " Reversing initial face orientation\n")); 3008 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize - 1 - f]; 3009 } else { 3010 if (debug) PetscCall(PetscPrintf(comm, " Keeping initial face orientation\n")); 3011 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 3012 } 3013 if (posOriented) *posOriented = posOrient; 3014 PetscFunctionReturn(PETSC_SUCCESS); 3015 } 3016 3017 /*@ 3018 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 3019 in faceVertices. The orientation is such that the face normal points out of the cell 3020 3021 Not Collective 3022 3023 Input Parameters: 3024 + dm - The original mesh 3025 . cell - The cell mesh point 3026 . faceSize - The number of vertices on the face 3027 . face - The face vertices 3028 . numCorners - The number of vertices on the cell 3029 . indices - Local numbering of face vertices in cell cone 3030 - origVertices - Original face vertices 3031 3032 Output Parameters: 3033 + faceVertices - The face vertices properly oriented 3034 - posOriented - `PETSC_TRUE` if the face was oriented with outward normal 3035 3036 Level: developer 3037 3038 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetCone()` 3039 @*/ 3040 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 3041 { 3042 const PetscInt *cone = NULL; 3043 PetscInt coneSize, v, f, v2; 3044 PetscInt oppositeVertex = -1; 3045 3046 PetscFunctionBegin; 3047 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 3048 PetscCall(DMPlexGetCone(dm, cell, &cone)); 3049 for (v = 0, v2 = 0; v < coneSize; ++v) { 3050 PetscBool found = PETSC_FALSE; 3051 3052 for (f = 0; f < faceSize; ++f) { 3053 if (face[f] == cone[v]) { 3054 found = PETSC_TRUE; 3055 break; 3056 } 3057 } 3058 if (found) { 3059 indices[v2] = v; 3060 origVertices[v2] = cone[v]; 3061 ++v2; 3062 } else { 3063 oppositeVertex = v; 3064 } 3065 } 3066 PetscCall(DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented)); 3067 PetscFunctionReturn(PETSC_SUCCESS); 3068 } 3069 3070 /* 3071 DMPlexInsertFace_Internal - Puts a face into the mesh 3072 3073 Not Collective 3074 3075 Input Parameters: 3076 + dm - The `DMPLEX` 3077 . numFaceVertex - The number of vertices in the face 3078 . faceVertices - The vertices in the face for `dm` 3079 . subfaceVertices - The vertices in the face for subdm 3080 . numCorners - The number of vertices in the `cell` 3081 . cell - A cell in `dm` containing the face 3082 . subcell - A cell in subdm containing the face 3083 . firstFace - First face in the mesh 3084 - newFacePoint - Next face in the mesh 3085 3086 Output Parameter: 3087 . newFacePoint - Contains next face point number on input, updated on output 3088 3089 Level: developer 3090 */ 3091 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) 3092 { 3093 MPI_Comm comm; 3094 DM_Plex *submesh = (DM_Plex *)subdm->data; 3095 const PetscInt *faces; 3096 PetscInt numFaces, coneSize; 3097 3098 PetscFunctionBegin; 3099 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3100 PetscCall(DMPlexGetConeSize(subdm, subcell, &coneSize)); 3101 PetscCheck(coneSize == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %" PetscInt_FMT " is %" PetscInt_FMT " != 1", cell, coneSize); 3102 #if 0 3103 /* Cannot use this because support() has not been constructed yet */ 3104 PetscCall(DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces)); 3105 #else 3106 { 3107 PetscInt f; 3108 3109 numFaces = 0; 3110 PetscCall(DMGetWorkArray(subdm, 1, MPIU_INT, (void **)&faces)); 3111 for (f = firstFace; f < *newFacePoint; ++f) { 3112 PetscInt dof, off, d; 3113 3114 PetscCall(PetscSectionGetDof(submesh->coneSection, f, &dof)); 3115 PetscCall(PetscSectionGetOffset(submesh->coneSection, f, &off)); 3116 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 3117 for (d = 0; d < dof; ++d) { 3118 const PetscInt p = submesh->cones[off + d]; 3119 PetscInt v; 3120 3121 for (v = 0; v < numFaceVertices; ++v) { 3122 if (subfaceVertices[v] == p) break; 3123 } 3124 if (v == numFaceVertices) break; 3125 } 3126 if (d == dof) { 3127 numFaces = 1; 3128 ((PetscInt *)faces)[0] = f; 3129 } 3130 } 3131 } 3132 #endif 3133 PetscCheck(numFaces <= 1, comm, PETSC_ERR_ARG_WRONG, "Vertex set had %" PetscInt_FMT " faces, not one", numFaces); 3134 if (numFaces == 1) { 3135 /* Add the other cell neighbor for this face */ 3136 PetscCall(DMPlexSetCone(subdm, subcell, faces)); 3137 } else { 3138 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 3139 PetscBool posOriented; 3140 3141 PetscCall(DMGetWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices)); 3142 origVertices = &orientedVertices[numFaceVertices]; 3143 indices = &orientedVertices[numFaceVertices * 2]; 3144 orientedSubVertices = &orientedVertices[numFaceVertices * 3]; 3145 PetscCall(DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented)); 3146 /* TODO: I know that routine should return a permutation, not the indices */ 3147 for (v = 0; v < numFaceVertices; ++v) { 3148 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 3149 for (ov = 0; ov < numFaceVertices; ++ov) { 3150 if (orientedVertices[ov] == vertex) { 3151 orientedSubVertices[ov] = subvertex; 3152 break; 3153 } 3154 } 3155 PetscCheck(ov != numFaceVertices, comm, PETSC_ERR_PLIB, "Could not find face vertex %" PetscInt_FMT " in orientated set", vertex); 3156 } 3157 PetscCall(DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices)); 3158 PetscCall(DMPlexSetCone(subdm, subcell, newFacePoint)); 3159 PetscCall(DMRestoreWorkArray(subdm, 4 * numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices)); 3160 ++(*newFacePoint); 3161 } 3162 #if 0 3163 PetscCall(DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces)); 3164 #else 3165 PetscCall(DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **)&faces)); 3166 #endif 3167 PetscFunctionReturn(PETSC_SUCCESS); 3168 } 3169 3170 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 3171 { 3172 MPI_Comm comm; 3173 DMLabel subpointMap; 3174 IS subvertexIS, subcellIS; 3175 const PetscInt *subVertices, *subCells; 3176 PetscInt numSubVertices, firstSubVertex, numSubCells; 3177 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 3178 PetscInt vStart, vEnd, c, f; 3179 3180 PetscFunctionBegin; 3181 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3182 /* Create subpointMap which marks the submesh */ 3183 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3184 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3185 if (vertexLabel) PetscCall(DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm)); 3186 /* Setup chart */ 3187 PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3188 PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3189 PetscCall(DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices)); 3190 PetscCall(DMPlexSetVTKCellHeight(subdm, 1)); 3191 /* Set cone sizes */ 3192 firstSubVertex = numSubCells; 3193 firstSubFace = numSubCells + numSubVertices; 3194 newFacePoint = firstSubFace; 3195 PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3196 if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices)); 3197 PetscCall(DMLabelGetStratumIS(subpointMap, 2, &subcellIS)); 3198 if (subcellIS) PetscCall(ISGetIndices(subcellIS, &subCells)); 3199 for (c = 0; c < numSubCells; ++c) PetscCall(DMPlexSetConeSize(subdm, c, 1)); 3200 for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) PetscCall(DMPlexSetConeSize(subdm, f, nFV)); 3201 PetscCall(DMSetUp(subdm)); 3202 PetscCall(DMLabelDestroy(&subpointMap)); 3203 /* Create face cones */ 3204 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3205 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3206 PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 3207 for (c = 0; c < numSubCells; ++c) { 3208 const PetscInt cell = subCells[c]; 3209 const PetscInt subcell = c; 3210 PetscInt *closure = NULL; 3211 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 3212 3213 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 3214 for (cl = 0; cl < closureSize * 2; cl += 2) { 3215 const PetscInt point = closure[cl]; 3216 PetscInt subVertex; 3217 3218 if ((point >= vStart) && (point < vEnd)) { 3219 ++numCorners; 3220 PetscCall(PetscFindInt(point, numSubVertices, subVertices, &subVertex)); 3221 if (subVertex >= 0) { 3222 closure[faceSize] = point; 3223 subface[faceSize] = firstSubVertex + subVertex; 3224 ++faceSize; 3225 } 3226 } 3227 } 3228 PetscCheck(faceSize <= nFV, comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %" PetscInt_FMT " of an element on the surface", faceSize); 3229 if (faceSize == nFV) PetscCall(DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint)); 3230 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 3231 } 3232 PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 3233 PetscCall(DMPlexSymmetrize(subdm)); 3234 PetscCall(DMPlexStratify(subdm)); 3235 /* Build coordinates */ 3236 { 3237 PetscSection coordSection, subCoordSection; 3238 Vec coordinates, subCoordinates; 3239 PetscScalar *coords, *subCoords; 3240 PetscInt numComp, coordSize, v; 3241 const char *name; 3242 3243 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3244 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3245 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3246 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3247 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3248 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3249 PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices)); 3250 for (v = 0; v < numSubVertices; ++v) { 3251 const PetscInt vertex = subVertices[v]; 3252 const PetscInt subvertex = firstSubVertex + v; 3253 PetscInt dof; 3254 3255 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3256 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3257 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3258 } 3259 PetscCall(PetscSectionSetUp(subCoordSection)); 3260 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3261 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3262 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 3263 PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name)); 3264 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3265 PetscCall(VecSetType(subCoordinates, VECSTANDARD)); 3266 if (coordSize) { 3267 PetscCall(VecGetArray(coordinates, &coords)); 3268 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3269 for (v = 0; v < numSubVertices; ++v) { 3270 const PetscInt vertex = subVertices[v]; 3271 const PetscInt subvertex = firstSubVertex + v; 3272 PetscInt dof, off, sdof, soff, d; 3273 3274 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3275 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 3276 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3277 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3278 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); 3279 for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d]; 3280 } 3281 PetscCall(VecRestoreArray(coordinates, &coords)); 3282 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3283 } 3284 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3285 PetscCall(VecDestroy(&subCoordinates)); 3286 } 3287 /* Cleanup */ 3288 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices)); 3289 PetscCall(ISDestroy(&subvertexIS)); 3290 if (subcellIS) PetscCall(ISRestoreIndices(subcellIS, &subCells)); 3291 PetscCall(ISDestroy(&subcellIS)); 3292 PetscFunctionReturn(PETSC_SUCCESS); 3293 } 3294 3295 /* TODO: Fix this to properly propagate up error conditions it may find */ 3296 static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 3297 { 3298 PetscInt subPoint; 3299 PetscErrorCode ierr; 3300 3301 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); 3302 if (ierr) return -1; 3303 return subPoint < 0 ? subPoint : firstSubPoint + subPoint; 3304 } 3305 3306 /* TODO: Fix this to properly propagate up error conditions it may find */ 3307 static inline PetscInt DMPlexFilterPointPerm_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[], const PetscInt subIndices[]) 3308 { 3309 PetscInt subPoint; 3310 PetscErrorCode ierr; 3311 3312 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); 3313 if (ierr) return -1; 3314 return subPoint < 0 ? subPoint : firstSubPoint + (subIndices ? subIndices[subPoint] : subPoint); 3315 } 3316 3317 static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm) 3318 { 3319 DMLabel depthLabel; 3320 PetscInt Nl, l, d; 3321 3322 PetscFunctionBegin; 3323 // Reset depth label for fast lookup 3324 PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); 3325 PetscCall(DMLabelMakeAllInvalid_Internal(depthLabel)); 3326 PetscCall(DMGetNumLabels(dm, &Nl)); 3327 for (l = 0; l < Nl; ++l) { 3328 DMLabel label, newlabel; 3329 const char *lname; 3330 PetscBool isDepth, isDim, isCelltype, isVTK; 3331 IS valueIS; 3332 const PetscInt *values; 3333 PetscInt Nv, v; 3334 3335 PetscCall(DMGetLabelName(dm, l, &lname)); 3336 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 3337 PetscCall(PetscStrcmp(lname, "dim", &isDim)); 3338 PetscCall(PetscStrcmp(lname, "celltype", &isCelltype)); 3339 PetscCall(PetscStrcmp(lname, "vtk", &isVTK)); 3340 if (isDepth || isDim || isCelltype || isVTK) continue; 3341 PetscCall(DMCreateLabel(subdm, lname)); 3342 PetscCall(DMGetLabel(dm, lname, &label)); 3343 PetscCall(DMGetLabel(subdm, lname, &newlabel)); 3344 PetscCall(DMLabelGetDefaultValue(label, &v)); 3345 PetscCall(DMLabelSetDefaultValue(newlabel, v)); 3346 PetscCall(DMLabelGetValueIS(label, &valueIS)); 3347 PetscCall(ISGetLocalSize(valueIS, &Nv)); 3348 PetscCall(ISGetIndices(valueIS, &values)); 3349 for (v = 0; v < Nv; ++v) { 3350 IS pointIS; 3351 const PetscInt *points; 3352 PetscInt Np, p; 3353 3354 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 3355 PetscCall(ISGetLocalSize(pointIS, &Np)); 3356 PetscCall(ISGetIndices(pointIS, &points)); 3357 for (p = 0; p < Np; ++p) { 3358 const PetscInt point = points[p]; 3359 PetscInt subp, subdepth; 3360 3361 PetscCall(DMPlexGetPointDepth(dm, point, &d)); 3362 PetscCall(DMPlexGetDepth(subdm, &subdepth)); 3363 if (d > subdepth) continue; 3364 subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]); 3365 if (subp >= 0) PetscCall(DMLabelSetValue(newlabel, subp, values[v])); 3366 } 3367 PetscCall(ISRestoreIndices(pointIS, &points)); 3368 PetscCall(ISDestroy(&pointIS)); 3369 } 3370 PetscCall(ISRestoreIndices(valueIS, &values)); 3371 PetscCall(ISDestroy(&valueIS)); 3372 } 3373 PetscFunctionReturn(PETSC_SUCCESS); 3374 } 3375 3376 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, PetscBool ignoreLabelHalo, PetscBool sanitizeSubmesh, PetscSF *ownershipTransferSF, DM subdm) 3377 { 3378 MPI_Comm comm; 3379 DMLabel subpointMap; 3380 IS *subpointIS; 3381 const PetscInt **subpoints; 3382 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 3383 PetscInt totSubPoints = 0, maxConeSize, dim, sdim, cdim, p, d, coordinate_type; 3384 PetscMPIInt rank; 3385 3386 PetscFunctionBegin; 3387 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3388 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 3389 /* Create subpointMap which marks the submesh */ 3390 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3391 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3392 if (cellHeight) { 3393 if (isCohesive) PetscCall(DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm)); 3394 else PetscCall(DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, PETSC_TRUE, subpointMap, subdm)); 3395 } else { 3396 DMLabel depth; 3397 IS pointIS; 3398 const PetscInt *points; 3399 PetscInt numPoints = 0, pStart = 0; 3400 PetscBool *isGhost = NULL; 3401 3402 PetscCall(DMPlexGetDepthLabel(dm, &depth)); 3403 PetscCall(DMLabelGetStratumIS(label, value, &pointIS)); 3404 if (pointIS) { 3405 PetscCall(ISGetIndices(pointIS, &points)); 3406 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 3407 } 3408 if (ignoreLabelHalo) { 3409 PetscSF pointSF; 3410 PetscInt nleaves, pEnd; 3411 const PetscInt *ilocal = NULL; 3412 3413 PetscCall(DMGetPointSF(dm, &pointSF)); 3414 PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, NULL)); 3415 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3416 PetscCall(PetscMalloc1(pEnd - pStart, &isGhost)); 3417 for (p = 0; p < pEnd - pStart; ++p) isGhost[p] = PETSC_FALSE; 3418 for (p = 0; p < nleaves; ++p) isGhost[(ilocal ? ilocal[p] : p) - pStart] = PETSC_TRUE; 3419 } 3420 for (p = 0; p < numPoints; ++p) { 3421 PetscInt *closure = NULL; 3422 PetscInt closureSize, c, pdim; 3423 3424 if (ignoreLabelHalo && isGhost[points[p] - pStart]) continue; 3425 PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3426 for (c = 0; c < closureSize * 2; c += 2) { 3427 PetscCall(DMLabelGetValue(depth, closure[c], &pdim)); 3428 PetscCall(DMLabelSetValue(subpointMap, closure[c], pdim)); 3429 } 3430 PetscCall(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3431 } 3432 PetscCall(PetscFree(isGhost)); 3433 if (pointIS) PetscCall(ISRestoreIndices(pointIS, &points)); 3434 PetscCall(ISDestroy(&pointIS)); 3435 } 3436 /* Setup chart */ 3437 PetscCall(DMGetDimension(dm, &dim)); 3438 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3439 PetscCall(PetscMalloc4(dim + 1, &numSubPoints, dim + 1, &firstSubPoint, dim + 1, &subpointIS, dim + 1, &subpoints)); 3440 for (d = 0; d <= dim; ++d) { 3441 PetscCall(DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d])); 3442 totSubPoints += numSubPoints[d]; 3443 } 3444 // Determine submesh dimension 3445 PetscCall(DMGetDimension(subdm, &sdim)); 3446 if (sdim > 0) { 3447 // Calling function knows what dimension to use, and we include neighboring cells as well 3448 sdim = dim; 3449 } else { 3450 // We reset the subdimension based on what is being selected 3451 PetscInt lsdim; 3452 for (lsdim = dim; lsdim >= 0; --lsdim) 3453 if (numSubPoints[lsdim]) break; 3454 PetscCallMPI(MPIU_Allreduce(&lsdim, &sdim, 1, MPIU_INT, MPI_MAX, comm)); 3455 PetscCall(DMSetDimension(subdm, sdim)); 3456 PetscCall(DMSetCoordinateDim(subdm, cdim)); 3457 } 3458 PetscCall(DMPlexSetChart(subdm, 0, totSubPoints)); 3459 PetscCall(DMPlexSetVTKCellHeight(subdm, cellHeight)); 3460 /* Set cone sizes */ 3461 firstSubPoint[sdim] = 0; 3462 firstSubPoint[0] = firstSubPoint[sdim] + numSubPoints[sdim]; 3463 if (sdim > 1) firstSubPoint[sdim - 1] = firstSubPoint[0] + numSubPoints[0]; 3464 if (sdim > 2) firstSubPoint[sdim - 2] = firstSubPoint[sdim - 1] + numSubPoints[sdim - 1]; 3465 for (d = 0; d <= sdim; ++d) { 3466 PetscCall(DMLabelGetStratumIS(subpointMap, d, &subpointIS[d])); 3467 if (subpointIS[d]) PetscCall(ISGetIndices(subpointIS[d], &subpoints[d])); 3468 } 3469 /* We do not want this label automatically computed, instead we compute it here */ 3470 PetscCall(DMCreateLabel(subdm, "celltype")); 3471 for (d = 0; d <= sdim; ++d) { 3472 for (p = 0; p < numSubPoints[d]; ++p) { 3473 const PetscInt point = subpoints[d][p]; 3474 const PetscInt subpoint = firstSubPoint[d] + p; 3475 const PetscInt *cone; 3476 PetscInt coneSize; 3477 3478 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 3479 if (cellHeight && (d == sdim)) { 3480 PetscInt coneSizeNew, c, val; 3481 3482 PetscCall(DMPlexGetCone(dm, point, &cone)); 3483 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3484 PetscCall(DMLabelGetValue(subpointMap, cone[c], &val)); 3485 if (val >= 0) coneSizeNew++; 3486 } 3487 PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSizeNew)); 3488 PetscCall(DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST)); 3489 } else { 3490 DMPolytopeType ct; 3491 3492 PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSize)); 3493 PetscCall(DMPlexGetCellType(dm, point, &ct)); 3494 PetscCall(DMPlexSetCellType(subdm, subpoint, ct)); 3495 } 3496 } 3497 } 3498 PetscCall(DMLabelDestroy(&subpointMap)); 3499 PetscCall(DMSetUp(subdm)); 3500 /* Set cones */ 3501 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3502 PetscCall(PetscMalloc2(maxConeSize, &coneNew, maxConeSize, &orntNew)); 3503 for (d = 0; d <= sdim; ++d) { 3504 for (p = 0; p < numSubPoints[d]; ++p) { 3505 const PetscInt point = subpoints[d][p]; 3506 const PetscInt subpoint = firstSubPoint[d] + p; 3507 const PetscInt *cone, *ornt; 3508 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 3509 3510 if (d == sdim - 1) { 3511 const PetscInt *support, *cone, *ornt; 3512 PetscInt supportSize, coneSize, s, subc; 3513 3514 PetscCall(DMPlexGetSupport(dm, point, &support)); 3515 PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); 3516 for (s = 0; s < supportSize; ++s) { 3517 PetscBool isHybrid = PETSC_FALSE; 3518 3519 PetscCall(DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid)); 3520 if (!isHybrid) continue; 3521 PetscCall(PetscFindInt(support[s], numSubPoints[d + 1], subpoints[d + 1], &subc)); 3522 if (subc >= 0) { 3523 const PetscInt ccell = subpoints[d + 1][subc]; 3524 3525 PetscCall(DMPlexGetCone(dm, ccell, &cone)); 3526 PetscCall(DMPlexGetConeSize(dm, ccell, &coneSize)); 3527 PetscCall(DMPlexGetConeOrientation(dm, ccell, &ornt)); 3528 for (c = 0; c < coneSize; ++c) { 3529 if (cone[c] == point) { 3530 fornt = ornt[c]; 3531 break; 3532 } 3533 } 3534 break; 3535 } 3536 } 3537 } 3538 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 3539 PetscCall(DMPlexGetConeSize(subdm, subpoint, &subconeSize)); 3540 PetscCall(DMPlexGetCone(dm, point, &cone)); 3541 PetscCall(DMPlexGetConeOrientation(dm, point, &ornt)); 3542 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3543 PetscCall(PetscFindInt(cone[c], numSubPoints[d - 1], subpoints[d - 1], &subc)); 3544 if (subc >= 0) { 3545 coneNew[coneSizeNew] = firstSubPoint[d - 1] + subc; 3546 orntNew[coneSizeNew] = ornt[c]; 3547 ++coneSizeNew; 3548 } 3549 } 3550 PetscCheck(coneSizeNew == subconeSize, comm, PETSC_ERR_PLIB, "Number of cone points located %" PetscInt_FMT " does not match subcone size %" PetscInt_FMT, coneSizeNew, subconeSize); 3551 PetscCall(DMPlexSetCone(subdm, subpoint, coneNew)); 3552 PetscCall(DMPlexSetConeOrientation(subdm, subpoint, orntNew)); 3553 if (fornt < 0) PetscCall(DMPlexOrientPoint(subdm, subpoint, fornt)); 3554 } 3555 } 3556 PetscCall(PetscFree2(coneNew, orntNew)); 3557 PetscCall(DMPlexSymmetrize(subdm)); 3558 PetscCall(DMPlexStratify(subdm)); 3559 /* Build coordinates */ 3560 for (coordinate_type = 0; coordinate_type < 2; ++coordinate_type) { 3561 DM coordDM = NULL, subCoordDM = NULL; 3562 PetscSection coordSection = NULL, subCoordSection = NULL; 3563 Vec coordinates = NULL, subCoordinates = NULL; 3564 PetscScalar *coords = NULL, *subCoords = NULL; 3565 PetscInt cdim, numComp, coordSize, firstP, lastP, firstSubP = totSubPoints, lastSubP = -1, numFields; 3566 const char *name; 3567 PetscBool localized = (PetscBool)coordinate_type; 3568 3569 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3570 if (!dm->coordinates[coordinate_type].dm) continue; 3571 if (!localized) { 3572 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3573 PetscCall(DMGetCoordinateDM(subdm, &subCoordDM)); 3574 } else { 3575 { 3576 PetscInt localizationHeight; 3577 PetscBool sparseLocalize; 3578 3579 PetscCall(DMGetSparseLocalize(dm, &sparseLocalize)); 3580 PetscCall(DMSetSparseLocalize(subdm, sparseLocalize)); 3581 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 3582 PetscCall(DMGetCoordinateDM(subdm, &subCoordDM)); 3583 PetscCall(DMPlexGetMaxProjectionHeight(coordDM, &localizationHeight)); 3584 PetscCall(DMPlexSetMaxProjectionHeight(subCoordDM, localizationHeight - (dim - sdim) - cellHeight)); 3585 PetscUseTypeMethod(subdm, createcellcoordinatedm, &subCoordDM); 3586 PetscCall(DMSetCellCoordinateDM(subdm, subCoordDM)); 3587 PetscCall(DMDestroy(&subCoordDM)); 3588 } 3589 PetscCall(DMGetCellCoordinateDM(dm, &coordDM)); 3590 PetscCall(DMGetCellCoordinateDM(subdm, &subCoordDM)); 3591 } 3592 PetscCall(DMGetNumFields(coordDM, &numFields)); 3593 if (numFields > 0) { 3594 PetscFE fe = NULL; 3595 PetscSpace P = NULL; 3596 PetscClassId id; 3597 PetscInt degree; 3598 3599 PetscCall(DMGetField(coordDM, 0, NULL, (PetscObject *)&fe)); 3600 PetscCall(PetscObjectGetClassId((PetscObject)fe, &id)); 3601 if (id == PETSCFE_CLASSID) { 3602 if (sdim == dim && cellHeight == 0) { 3603 /* TODO: Handle Field labels correctly */ 3604 PetscCall(DMSetField(subCoordDM, 0, NULL, (PetscObject)fe)); 3605 PetscCall(DMCreateDS(subCoordDM)); 3606 } else { 3607 /* TODO: Reconstruct the lower-dimensional FE more robustly */ 3608 PetscCall(PetscFEGetBasisSpace(fe, &P)); 3609 PetscCall(PetscSpaceGetDegree(P, °ree, NULL)); 3610 PetscCall(DMPlexCreateCoordinateSpace(subdm, degree, localized, PETSC_FALSE)); 3611 } 3612 } 3613 } 3614 if (!localized) { 3615 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3616 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3617 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3618 } else { 3619 PetscCall(DMGetCellCoordinateSection(dm, &coordSection)); 3620 PetscCall(DMGetCellCoordinateSection(subdm, &subCoordSection)); 3621 PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates)); 3622 } 3623 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3624 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3625 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3626 PetscCall(PetscSectionGetChart(coordSection, &firstP, &lastP)); 3627 for (d = 0; d <= sdim; ++d) { 3628 for (p = 0; p < numSubPoints[d]; ++p) { 3629 const PetscInt point = subpoints[d][p]; 3630 const PetscInt subpoint = firstSubPoint[d] + p; 3631 3632 if (point >= firstP && point < lastP) { 3633 firstSubP = PetscMin(firstSubP, subpoint); 3634 lastSubP = PetscMax(lastSubP, subpoint); 3635 } 3636 } 3637 } 3638 lastSubP += 1; 3639 if (firstSubP == totSubPoints) { 3640 /* Zero if there is no coordinate point. */ 3641 firstSubP = 0; 3642 lastSubP = 0; 3643 } 3644 PetscCall(PetscSectionSetChart(subCoordSection, firstSubP, lastSubP)); 3645 for (d = 0; d <= sdim; ++d) { 3646 for (p = 0; p < numSubPoints[d]; ++p) { 3647 const PetscInt point = subpoints[d][p]; 3648 const PetscInt subpoint = firstSubPoint[d] + p; 3649 PetscInt dof; 3650 3651 if (point >= firstP && point < lastP) { 3652 PetscCall(PetscSectionGetDof(coordSection, point, &dof)); 3653 PetscCall(PetscSectionSetDof(subCoordSection, subpoint, dof)); 3654 PetscCall(PetscSectionSetFieldDof(subCoordSection, subpoint, 0, dof)); 3655 } 3656 } 3657 } 3658 PetscCall(PetscSectionSetUp(subCoordSection)); 3659 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3660 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3661 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 3662 PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name)); 3663 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3664 PetscCall(VecSetBlockSize(subCoordinates, cdim)); 3665 PetscCall(VecSetType(subCoordinates, VECSTANDARD)); 3666 PetscCall(VecGetArray(coordinates, &coords)); 3667 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3668 for (d = 0; d <= sdim; ++d) { 3669 for (p = 0; p < numSubPoints[d]; ++p) { 3670 const PetscInt point = subpoints[d][p]; 3671 const PetscInt subpoint = firstSubPoint[d] + p; 3672 PetscInt dof, off, sdof, soff, d; 3673 3674 if (point >= firstP && point < lastP) { 3675 PetscCall(PetscSectionGetDof(coordSection, point, &dof)); 3676 PetscCall(PetscSectionGetOffset(coordSection, point, &off)); 3677 PetscCall(PetscSectionGetDof(subCoordSection, subpoint, &sdof)); 3678 PetscCall(PetscSectionGetOffset(subCoordSection, subpoint, &soff)); 3679 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); 3680 for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d]; 3681 } 3682 } 3683 } 3684 PetscCall(VecRestoreArray(coordinates, &coords)); 3685 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3686 switch (coordinate_type) { 3687 case 0: 3688 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3689 break; 3690 case 1: 3691 PetscCall(DMSetCellCoordinatesLocal(subdm, subCoordinates)); 3692 break; 3693 default: 3694 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "coordinate_type must be {0, 1}"); 3695 } 3696 PetscCall(VecDestroy(&subCoordinates)); 3697 } 3698 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3699 { 3700 PetscSF sfPoint, sfPointSub; 3701 IS subpIS; 3702 const PetscSFNode *remotePoints; 3703 PetscSFNode *sremotePoints = NULL, *newLocalPoints = NULL, *newOwners = NULL; 3704 const PetscInt *localPoints, *subpoints, *rootdegree; 3705 PetscInt *slocalPoints = NULL, *sortedPoints = NULL, *sortedIndices = NULL; 3706 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl = 0, ll = 0, pStart, pEnd, p; 3707 PetscMPIInt rank, size; 3708 3709 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 3710 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 3711 PetscCall(DMGetPointSF(dm, &sfPoint)); 3712 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 3713 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3714 PetscCall(DMPlexGetChart(subdm, NULL, &numSubroots)); 3715 PetscCall(DMPlexGetSubpointIS(subdm, &subpIS)); 3716 if (subpIS) { 3717 PetscBool sorted = PETSC_TRUE; 3718 3719 PetscCall(ISGetIndices(subpIS, &subpoints)); 3720 PetscCall(ISGetLocalSize(subpIS, &numSubpoints)); 3721 for (p = 1; p < numSubpoints; ++p) sorted = sorted && (subpoints[p] >= subpoints[p - 1]) ? PETSC_TRUE : PETSC_FALSE; 3722 if (!sorted) { 3723 PetscCall(PetscMalloc2(numSubpoints, &sortedPoints, numSubpoints, &sortedIndices)); 3724 for (p = 0; p < numSubpoints; ++p) sortedIndices[p] = p; 3725 PetscCall(PetscArraycpy(sortedPoints, subpoints, numSubpoints)); 3726 PetscCall(PetscSortIntWithArray(numSubpoints, sortedPoints, sortedIndices)); 3727 } 3728 } 3729 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3730 if (numRoots >= 0) { 3731 PetscCall(PetscSFComputeDegreeBegin(sfPoint, &rootdegree)); 3732 PetscCall(PetscSFComputeDegreeEnd(sfPoint, &rootdegree)); 3733 PetscCall(PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners)); 3734 for (p = 0; p < pEnd - pStart; ++p) { 3735 newLocalPoints[p].rank = -2; 3736 newLocalPoints[p].index = -2; 3737 } 3738 for (p = pStart; p < pEnd; ++p) { 3739 newOwners[p - pStart].rank = -3; 3740 newOwners[p - pStart].index = -3; 3741 } 3742 if (sanitizeSubmesh) { 3743 /* A subpoint is forced to be owned by a rank that owns */ 3744 /* a subcell that contains the subpoint in its closure. */ 3745 PetscInt cStart, cEnd, c, clSize, cl; 3746 PetscInt *ownedCells, *closure = NULL; 3747 3748 /* claim ownership */ 3749 for (p = 0; p < numSubpoints; ++p) { 3750 const PetscInt point = subpoints[p]; 3751 3752 newLocalPoints[point - pStart].rank = rank; 3753 newLocalPoints[point - pStart].index = p; 3754 } 3755 PetscCall(DMGetDimension(subdm, &sdim)); 3756 PetscCall(DMPlexGetDepthStratum(dm, sdim, &cStart, &cEnd)); 3757 PetscCall(PetscMalloc1(cEnd - cStart, &ownedCells)); 3758 for (c = cStart; c < cEnd; ++c) ownedCells[c - cStart] = 0; 3759 for (p = 0; p < numSubpoints; ++p) { 3760 c = subpoints[p]; 3761 if (c >= cStart && c < cEnd) ownedCells[c - cStart] = 1; 3762 } 3763 for (l = 0; l < numLeaves; ++l) { 3764 c = localPoints ? localPoints[l] : l; 3765 if (c >= cStart && c < cEnd) ownedCells[c - cStart] = 0; 3766 } 3767 for (c = cStart; c < cEnd; ++c) { 3768 if (ownedCells[c - cStart] == 0) continue; 3769 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 3770 for (cl = 0; cl < clSize * 2; cl += 2) { 3771 p = closure[cl]; 3772 if (newLocalPoints[p - pStart].rank < size) newLocalPoints[p - pStart].rank += size; 3773 } 3774 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 3775 } 3776 PetscCall(PetscFree(ownedCells)); 3777 for (p = 0; p < numRoots; ++p) { 3778 newOwners[p].rank = newLocalPoints[p].rank; 3779 newOwners[p].index = newLocalPoints[p].index; 3780 } 3781 } else { 3782 /* Set subleaves */ 3783 for (l = 0; l < numLeaves; ++l) { 3784 const PetscInt point = localPoints[l]; 3785 const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices); 3786 3787 if (subpoint < 0) continue; 3788 newLocalPoints[point - pStart].rank = rank; 3789 newLocalPoints[point - pStart].index = subpoint; 3790 ++numSubleaves; 3791 } 3792 /* Must put in owned subpoints */ 3793 for (p = 0; p < numSubpoints; ++p) { 3794 /* Hold on to currently owned points */ 3795 if (rootdegree[subpoints[p] - pStart]) newOwners[subpoints[p] - pStart].rank = rank + size; 3796 else newOwners[subpoints[p] - pStart].rank = rank; 3797 newOwners[subpoints[p] - pStart].index = p; 3798 } 3799 } 3800 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 3801 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 3802 for (p = pStart; p < pEnd; ++p) 3803 if (newOwners[p - pStart].rank >= size) newOwners[p - pStart].rank -= size; 3804 if (ownershipTransferSF) { 3805 PetscSFNode *iremote1 = NULL, *newOwners1 = NULL; 3806 PetscInt *ilocal1 = NULL; 3807 PetscInt nleaves1, point; 3808 3809 for (p = 0; p < numSubpoints; ++p) { 3810 point = subpoints[p]; 3811 newLocalPoints[point - pStart].index = point; 3812 } 3813 PetscCall(PetscMalloc1(numRoots, &newOwners1)); 3814 for (p = 0; p < numRoots; ++p) { 3815 newOwners1[p].rank = -1; 3816 newOwners1[p].index = -1; 3817 } 3818 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners1, MPI_MAXLOC)); 3819 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners1, MPI_MAXLOC)); 3820 for (p = 0, nleaves1 = 0; p < numRoots; ++p) { 3821 if (newOwners[p].rank >= 0 && newOwners[p].rank != rank) { ++nleaves1; } 3822 } 3823 PetscCall(PetscMalloc1(nleaves1, &ilocal1)); 3824 PetscCall(PetscMalloc1(nleaves1, &iremote1)); 3825 for (p = 0, nleaves1 = 0; p < numRoots; ++p) { 3826 if (newOwners[p].rank >= 0 && newOwners[p].rank != rank) { 3827 ilocal1[nleaves1] = pStart + p; 3828 iremote1[nleaves1].rank = newOwners[p].rank; 3829 iremote1[nleaves1].index = newOwners1[p].index; 3830 ++nleaves1; 3831 } 3832 } 3833 PetscCall(PetscFree(newOwners1)); 3834 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), ownershipTransferSF)); 3835 PetscCall(PetscSFSetFromOptions(*ownershipTransferSF)); 3836 PetscCall(PetscSFSetGraph(*ownershipTransferSF, pEnd - pStart, nleaves1, ilocal1, PETSC_OWN_POINTER, iremote1, PETSC_OWN_POINTER)); 3837 } 3838 if (sanitizeSubmesh) { 3839 for (p = pStart; p < pEnd; ++p) { 3840 newLocalPoints[p - pStart].rank = newOwners[p - pStart].rank; 3841 newLocalPoints[p - pStart].index = newOwners[p - pStart].index; 3842 } 3843 } 3844 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 3845 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 3846 if (sanitizeSubmesh) { 3847 for (p = 0; p < numSubpoints; ++p) { 3848 const PetscInt point = subpoints[p]; 3849 3850 if (newLocalPoints[point - pStart].rank >= 0 && newLocalPoints[point - pStart].rank != rank) ++sl; 3851 } 3852 PetscCall(PetscMalloc1(sl, &slocalPoints)); 3853 PetscCall(PetscMalloc1(sl, &sremotePoints)); 3854 for (p = 0, sl = 0; p < numSubpoints; ++p) { 3855 const PetscInt point = subpoints[p]; 3856 3857 if (newLocalPoints[point - pStart].rank >= 0 && newLocalPoints[point - pStart].rank != rank) { 3858 slocalPoints[sl] = p; 3859 sremotePoints[sl].rank = newLocalPoints[point - pStart].rank; 3860 sremotePoints[sl].index = newLocalPoints[point - pStart].index; 3861 PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank"); 3862 PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint"); 3863 ++sl; 3864 } 3865 } 3866 } else { 3867 PetscCall(PetscMalloc1(numSubleaves, &slocalPoints)); 3868 PetscCall(PetscMalloc1(numSubleaves, &sremotePoints)); 3869 for (l = 0; l < numLeaves; ++l) { 3870 const PetscInt point = localPoints[l]; 3871 const PetscInt subpoint = DMPlexFilterPointPerm_Internal(point, 0, numSubpoints, sortedPoints ? sortedPoints : subpoints, sortedIndices); 3872 3873 if (subpoint < 0) continue; 3874 if (newLocalPoints[point].rank == rank) { 3875 ++ll; 3876 continue; 3877 } 3878 slocalPoints[sl] = subpoint; 3879 sremotePoints[sl].rank = newLocalPoints[point].rank; 3880 sremotePoints[sl].index = newLocalPoints[point].index; 3881 PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 3882 PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 3883 ++sl; 3884 } 3885 PetscCheck(sl + ll == numSubleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubleaves); 3886 } 3887 PetscCall(PetscFree2(newLocalPoints, newOwners)); 3888 PetscCall(PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3889 } 3890 if (subpIS) PetscCall(ISRestoreIndices(subpIS, &subpoints)); 3891 PetscCall(PetscFree2(sortedPoints, sortedIndices)); 3892 } 3893 /* Filter labels */ 3894 PetscCall(DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm)); 3895 /* Cleanup */ 3896 for (d = 0; d <= sdim; ++d) { 3897 if (subpointIS[d]) PetscCall(ISRestoreIndices(subpointIS[d], &subpoints[d])); 3898 PetscCall(ISDestroy(&subpointIS[d])); 3899 } 3900 PetscCall(PetscFree4(numSubPoints, firstSubPoint, subpointIS, subpoints)); 3901 PetscFunctionReturn(PETSC_SUCCESS); 3902 } 3903 3904 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm) 3905 { 3906 PetscFunctionBegin; 3907 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm)); 3908 PetscFunctionReturn(PETSC_SUCCESS); 3909 } 3910 3911 /*@ 3912 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3913 3914 Input Parameters: 3915 + dm - The original mesh 3916 . vertexLabel - The `DMLabel` marking points contained in the surface 3917 . value - The label value to use 3918 - markedFaces - `PETSC_TRUE` if surface faces are marked in addition to vertices, `PETSC_FALSE` if only vertices are marked 3919 3920 Output Parameter: 3921 . subdm - The surface mesh 3922 3923 Level: developer 3924 3925 Note: 3926 This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`. 3927 3928 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()` 3929 @*/ 3930 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm) 3931 { 3932 DMPlexInterpolatedFlag interpolated; 3933 PetscInt dim, cdim; 3934 3935 PetscFunctionBegin; 3936 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3937 PetscAssertPointer(subdm, 5); 3938 PetscCall(DMGetDimension(dm, &dim)); 3939 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3940 PetscCall(DMSetType(*subdm, DMPLEX)); 3941 PetscCall(DMSetDimension(*subdm, dim - 1)); 3942 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3943 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 3944 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3945 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 3946 if (interpolated) { 3947 PetscCall(DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm)); 3948 } else { 3949 PetscCall(DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm)); 3950 } 3951 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3952 PetscFunctionReturn(PETSC_SUCCESS); 3953 } 3954 3955 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3956 { 3957 MPI_Comm comm; 3958 DMLabel subpointMap; 3959 IS subvertexIS; 3960 const PetscInt *subVertices; 3961 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3962 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3963 PetscInt c, f; 3964 3965 PetscFunctionBegin; 3966 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3967 /* Create subpointMap which marks the submesh */ 3968 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3969 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3970 PetscCall(DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm)); 3971 /* Setup chart */ 3972 PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3973 PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3974 PetscCall(DMPlexSetChart(subdm, 0, numSubCells + numSubFaces + numSubVertices)); 3975 PetscCall(DMPlexSetVTKCellHeight(subdm, 1)); 3976 /* Set cone sizes */ 3977 firstSubVertex = numSubCells; 3978 firstSubFace = numSubCells + numSubVertices; 3979 newFacePoint = firstSubFace; 3980 PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3981 if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices)); 3982 for (c = 0; c < numSubCells; ++c) PetscCall(DMPlexSetConeSize(subdm, c, 1)); 3983 for (f = firstSubFace; f < firstSubFace + numSubFaces; ++f) PetscCall(DMPlexSetConeSize(subdm, f, nFV)); 3984 PetscCall(DMSetUp(subdm)); 3985 PetscCall(DMLabelDestroy(&subpointMap)); 3986 /* Create face cones */ 3987 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3988 PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 3989 for (c = 0; c < numSubCells; ++c) { 3990 const PetscInt cell = subCells[c]; 3991 const PetscInt subcell = c; 3992 const PetscInt *cone, *cells; 3993 PetscBool isHybrid = PETSC_FALSE; 3994 PetscInt numCells, subVertex, p, v; 3995 3996 PetscCall(DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid)); 3997 if (!isHybrid) continue; 3998 PetscCall(DMPlexGetCone(dm, cell, &cone)); 3999 for (v = 0; v < nFV; ++v) { 4000 PetscCall(PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex)); 4001 subface[v] = firstSubVertex + subVertex; 4002 } 4003 PetscCall(DMPlexSetCone(subdm, newFacePoint, subface)); 4004 PetscCall(DMPlexSetCone(subdm, subcell, &newFacePoint)); 4005 PetscCall(DMPlexGetJoin(dm, nFV, cone, &numCells, &cells)); 4006 /* Not true in parallel 4007 PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 4008 for (p = 0; p < numCells; ++p) { 4009 PetscInt negsubcell; 4010 PetscBool isHybrid = PETSC_FALSE; 4011 4012 PetscCall(DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid)); 4013 if (isHybrid) continue; 4014 /* I know this is a crap search */ 4015 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 4016 if (subCells[negsubcell] == cells[p]) break; 4017 } 4018 PetscCheck(negsubcell != numSubCells, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %" PetscInt_FMT, cell); 4019 PetscCall(DMPlexSetCone(subdm, negsubcell, &newFacePoint)); 4020 } 4021 PetscCall(DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells)); 4022 ++newFacePoint; 4023 } 4024 PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void **)&subface)); 4025 PetscCall(DMPlexSymmetrize(subdm)); 4026 PetscCall(DMPlexStratify(subdm)); 4027 /* Build coordinates */ 4028 { 4029 PetscSection coordSection, subCoordSection; 4030 Vec coordinates, subCoordinates; 4031 PetscScalar *coords, *subCoords; 4032 PetscInt cdim, numComp, coordSize, v; 4033 const char *name; 4034 4035 PetscCall(DMGetCoordinateDim(dm, &cdim)); 4036 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 4037 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 4038 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 4039 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 4040 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 4041 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 4042 PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex + numSubVertices)); 4043 for (v = 0; v < numSubVertices; ++v) { 4044 const PetscInt vertex = subVertices[v]; 4045 const PetscInt subvertex = firstSubVertex + v; 4046 PetscInt dof; 4047 4048 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 4049 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 4050 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 4051 } 4052 PetscCall(PetscSectionSetUp(subCoordSection)); 4053 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 4054 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 4055 PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 4056 PetscCall(PetscObjectSetName((PetscObject)subCoordinates, name)); 4057 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 4058 PetscCall(VecSetBlockSize(subCoordinates, cdim)); 4059 PetscCall(VecSetType(subCoordinates, VECSTANDARD)); 4060 PetscCall(VecGetArray(coordinates, &coords)); 4061 PetscCall(VecGetArray(subCoordinates, &subCoords)); 4062 for (v = 0; v < numSubVertices; ++v) { 4063 const PetscInt vertex = subVertices[v]; 4064 const PetscInt subvertex = firstSubVertex + v; 4065 PetscInt dof, off, sdof, soff, d; 4066 4067 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 4068 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 4069 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 4070 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 4071 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); 4072 for (d = 0; d < dof; ++d) subCoords[soff + d] = coords[off + d]; 4073 } 4074 PetscCall(VecRestoreArray(coordinates, &coords)); 4075 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 4076 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 4077 PetscCall(VecDestroy(&subCoordinates)); 4078 } 4079 /* Build SF */ 4080 CHKMEMQ; 4081 { 4082 PetscSF sfPoint, sfPointSub; 4083 const PetscSFNode *remotePoints; 4084 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 4085 const PetscInt *localPoints; 4086 PetscInt *slocalPoints; 4087 PetscInt numRoots, numLeaves, numSubRoots = numSubCells + numSubFaces + numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 4088 PetscMPIInt rank; 4089 4090 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 4091 PetscCall(DMGetPointSF(dm, &sfPoint)); 4092 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 4093 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 4094 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 4095 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 4096 if (numRoots >= 0) { 4097 /* Only vertices should be shared */ 4098 PetscCall(PetscMalloc2(pEnd - pStart, &newLocalPoints, numRoots, &newOwners)); 4099 for (p = 0; p < pEnd - pStart; ++p) { 4100 newLocalPoints[p].rank = -2; 4101 newLocalPoints[p].index = -2; 4102 } 4103 /* Set subleaves */ 4104 for (l = 0; l < numLeaves; ++l) { 4105 const PetscInt point = localPoints[l]; 4106 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 4107 4108 PetscCheck(!(point < vStart) || !(point >= vEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %" PetscInt_FMT, point); 4109 if (subPoint < 0) continue; 4110 newLocalPoints[point - pStart].rank = rank; 4111 newLocalPoints[point - pStart].index = subPoint; 4112 ++numSubLeaves; 4113 } 4114 /* Must put in owned subpoints */ 4115 for (p = pStart; p < pEnd; ++p) { 4116 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 4117 4118 if (subPoint < 0) { 4119 newOwners[p - pStart].rank = -3; 4120 newOwners[p - pStart].index = -3; 4121 } else { 4122 newOwners[p - pStart].rank = rank; 4123 newOwners[p - pStart].index = subPoint; 4124 } 4125 } 4126 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 4127 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_SF_NODE, newLocalPoints, newOwners, MPI_MAXLOC)); 4128 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 4129 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_SF_NODE, newOwners, newLocalPoints, MPI_REPLACE)); 4130 PetscCall(PetscMalloc1(numSubLeaves, &slocalPoints)); 4131 PetscCall(PetscMalloc1(numSubLeaves, &sremotePoints)); 4132 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 4133 const PetscInt point = localPoints[l]; 4134 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 4135 4136 if (subPoint < 0) continue; 4137 if (newLocalPoints[point].rank == rank) { 4138 ++ll; 4139 continue; 4140 } 4141 slocalPoints[sl] = subPoint; 4142 sremotePoints[sl].rank = newLocalPoints[point].rank; 4143 sremotePoints[sl].index = newLocalPoints[point].index; 4144 PetscCheck(sremotePoints[sl].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 4145 PetscCheck(sremotePoints[sl].index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 4146 ++sl; 4147 } 4148 PetscCall(PetscFree2(newLocalPoints, newOwners)); 4149 PetscCheck(sl + ll == numSubLeaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubLeaves); 4150 PetscCall(PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 4151 } 4152 } 4153 CHKMEMQ; 4154 /* Cleanup */ 4155 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices)); 4156 PetscCall(ISDestroy(&subvertexIS)); 4157 PetscCall(PetscFree(subCells)); 4158 PetscFunctionReturn(PETSC_SUCCESS); 4159 } 4160 4161 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 4162 { 4163 DMLabel label = NULL; 4164 4165 PetscFunctionBegin; 4166 if (labelname) PetscCall(DMGetLabel(dm, labelname, &label)); 4167 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, PETSC_FALSE, PETSC_FALSE, NULL, subdm)); 4168 PetscFunctionReturn(PETSC_SUCCESS); 4169 } 4170 4171 /*@ 4172 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. 4173 4174 Input Parameters: 4175 + dm - The original mesh 4176 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 4177 . label - A label name, or `NULL` 4178 - value - A label value 4179 4180 Output Parameter: 4181 . subdm - The surface mesh 4182 4183 Level: developer 4184 4185 Note: 4186 This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`. 4187 4188 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()` 4189 @*/ 4190 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 4191 { 4192 PetscInt dim, cdim, depth; 4193 4194 PetscFunctionBegin; 4195 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4196 PetscAssertPointer(subdm, 5); 4197 PetscCall(DMGetDimension(dm, &dim)); 4198 PetscCall(DMPlexGetDepth(dm, &depth)); 4199 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 4200 PetscCall(DMSetType(*subdm, DMPLEX)); 4201 PetscCall(DMSetDimension(*subdm, dim - 1)); 4202 PetscCall(DMGetCoordinateDim(dm, &cdim)); 4203 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 4204 if (depth == dim) { 4205 PetscCall(DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm)); 4206 } else { 4207 PetscCall(DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm)); 4208 } 4209 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 4210 PetscFunctionReturn(PETSC_SUCCESS); 4211 } 4212 4213 /*@ 4214 DMPlexReorderCohesiveSupports - Ensure that face supports for cohesive end caps are ordered 4215 4216 Not Collective 4217 4218 Input Parameter: 4219 . dm - The `DM` containing cohesive cells 4220 4221 Level: developer 4222 4223 Note: 4224 For the negative size (first) face, the cohesive cell should be first in the support, and for 4225 the positive side (second) face, the cohesive cell should be second in the support. 4226 4227 .seealso: `DMPlexConstructCohesiveCells()`, `DMPlexCreateCohesiveSubmesh()` 4228 @*/ 4229 PetscErrorCode DMPlexReorderCohesiveSupports(DM dm) 4230 { 4231 PetscInt dim, cStart, cEnd; 4232 4233 PetscFunctionBegin; 4234 PetscCall(DMGetDimension(dm, &dim)); 4235 PetscCall(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cStart, &cEnd)); 4236 for (PetscInt c = cStart; c < cEnd; ++c) { 4237 const PetscInt *cone; 4238 PetscInt coneSize; 4239 4240 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 4241 PetscCall(DMPlexGetCone(dm, c, &cone)); 4242 PetscCheck(coneSize >= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hybrid cell %" PetscInt_FMT " cone size %" PetscInt_FMT " must be >= 2", c, coneSize); 4243 for (PetscInt s = 0; s < 2; ++s) { 4244 const PetscInt *supp; 4245 PetscInt suppSize, newsupp[2]; 4246 4247 PetscCall(DMPlexGetSupportSize(dm, cone[s], &suppSize)); 4248 PetscCall(DMPlexGetSupport(dm, cone[s], &supp)); 4249 if (suppSize == 2) { 4250 /* Reorder hybrid end cap support */ 4251 if (supp[s] == c) { 4252 newsupp[0] = supp[1]; 4253 newsupp[1] = supp[0]; 4254 } else { 4255 newsupp[0] = supp[0]; 4256 newsupp[1] = supp[1]; 4257 } 4258 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); 4259 PetscCall(DMPlexSetSupport(dm, cone[s], newsupp)); 4260 } 4261 } 4262 } 4263 PetscFunctionReturn(PETSC_SUCCESS); 4264 } 4265 4266 /*@ 4267 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 4268 4269 Input Parameters: 4270 + dm - The original mesh 4271 . cellLabel - The `DMLabel` marking cells contained in the new mesh 4272 . value - The label value to use 4273 . ignoreLabelHalo - The flag indicating if labeled points that are in the halo are ignored 4274 - 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 4275 4276 Output Parameters: 4277 + ownershipTransferSF - The `PetscSF` representing the ownership transfers between parent local meshes due to submeshing. 4278 - subdm - The new mesh 4279 4280 Level: developer 4281 4282 Note: 4283 This function produces a `DMLabel` mapping original points in the submesh to their depth. This can be obtained using `DMPlexGetSubpointMap()`. 4284 4285 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 4286 the remote locations for these leaves are tuples of rank and point that represent the new owners. 4287 4288 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()`, `DMPlexCreateSubmesh()` 4289 @*/ 4290 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, PetscBool ignoreLabelHalo, PetscBool sanitizeSubmesh, PetscSF *ownershipTransferSF, DM *subdm) 4291 { 4292 PetscInt dim, overlap; 4293 4294 PetscFunctionBegin; 4295 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4296 if (ownershipTransferSF) PetscAssertPointer(ownershipTransferSF, 6); 4297 PetscAssertPointer(subdm, 7); 4298 PetscCall(DMGetDimension(dm, &dim)); 4299 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 4300 PetscCall(DMSetType(*subdm, DMPLEX)); 4301 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 4302 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, ignoreLabelHalo, sanitizeSubmesh, ownershipTransferSF, *subdm)); 4303 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 4304 // It is possible to obtain a surface mesh where some faces are in SF 4305 // We should either mark the mesh as having an overlap, or delete these from the SF 4306 PetscCall(DMPlexGetOverlap(dm, &overlap)); 4307 if (!overlap) { 4308 PetscSF sf; 4309 const PetscInt *leaves; 4310 PetscInt cStart, cEnd, Nl; 4311 PetscBool hasSubcell = PETSC_FALSE, ghasSubcell; 4312 4313 PetscCall(DMPlexGetHeightStratum(*subdm, 0, &cStart, &cEnd)); 4314 PetscCall(DMGetPointSF(*subdm, &sf)); 4315 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL)); 4316 for (PetscInt l = 0; l < Nl; ++l) { 4317 const PetscInt point = leaves ? leaves[l] : l; 4318 4319 if (point >= cStart && point < cEnd) { 4320 hasSubcell = PETSC_TRUE; 4321 break; 4322 } 4323 } 4324 PetscCallMPI(MPIU_Allreduce(&hasSubcell, &ghasSubcell, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm))); 4325 if (ghasSubcell) PetscCall(DMPlexSetOverlap(*subdm, NULL, 1)); 4326 } 4327 PetscFunctionReturn(PETSC_SUCCESS); 4328 } 4329 4330 /*@ 4331 DMPlexGetSubpointMap - Returns a `DMLabel` with point dimension as values 4332 4333 Input Parameter: 4334 . dm - The submesh `DM` 4335 4336 Output Parameter: 4337 . subpointMap - The `DMLabel` of all the points from the original mesh in this submesh, or `NULL` if this is not a submesh 4338 4339 Level: developer 4340 4341 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 4342 @*/ 4343 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 4344 { 4345 PetscFunctionBegin; 4346 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4347 PetscAssertPointer(subpointMap, 2); 4348 *subpointMap = ((DM_Plex *)dm->data)->subpointMap; 4349 PetscFunctionReturn(PETSC_SUCCESS); 4350 } 4351 4352 /*@ 4353 DMPlexSetSubpointMap - Sets the `DMLabel` with point dimension as values 4354 4355 Input Parameters: 4356 + dm - The submesh `DM` 4357 - subpointMap - The `DMLabel` of all the points from the original mesh in this submesh 4358 4359 Level: developer 4360 4361 Note: 4362 Should normally not be called by the user, since it is set in `DMPlexCreateSubmesh()` 4363 4364 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 4365 @*/ 4366 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 4367 { 4368 DM_Plex *mesh = (DM_Plex *)dm->data; 4369 DMLabel tmp; 4370 4371 PetscFunctionBegin; 4372 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4373 tmp = mesh->subpointMap; 4374 mesh->subpointMap = subpointMap; 4375 PetscCall(PetscObjectReference((PetscObject)mesh->subpointMap)); 4376 PetscCall(DMLabelDestroy(&tmp)); 4377 PetscFunctionReturn(PETSC_SUCCESS); 4378 } 4379 4380 static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS) 4381 { 4382 DMLabel spmap; 4383 PetscInt depth, d; 4384 4385 PetscFunctionBegin; 4386 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 4387 PetscCall(DMPlexGetDepth(dm, &depth)); 4388 if (spmap && depth >= 0) { 4389 DM_Plex *mesh = (DM_Plex *)dm->data; 4390 PetscInt *points, *depths; 4391 PetscInt pStart, pEnd, p, off; 4392 4393 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 4394 PetscCheck(!pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %" PetscInt_FMT, pStart); 4395 PetscCall(PetscMalloc1(pEnd, &points)); 4396 PetscCall(DMGetWorkArray(dm, depth + 1, MPIU_INT, &depths)); 4397 depths[0] = depth; 4398 depths[1] = 0; 4399 for (d = 2; d <= depth; ++d) depths[d] = depth + 1 - d; 4400 for (d = 0, off = 0; d <= depth; ++d) { 4401 const PetscInt dep = depths[d]; 4402 PetscInt depStart, depEnd, n; 4403 4404 PetscCall(DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd)); 4405 PetscCall(DMLabelGetStratumSize(spmap, dep, &n)); 4406 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 4407 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); 4408 } else { 4409 if (!n) { 4410 if (d == 0) { 4411 /* Missing cells */ 4412 for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = -1; 4413 } else { 4414 /* Missing faces */ 4415 for (p = 0; p < depEnd - depStart; ++p, ++off) points[off] = PETSC_INT_MAX; 4416 } 4417 } 4418 } 4419 if (n) { 4420 IS is; 4421 const PetscInt *opoints; 4422 4423 PetscCall(DMLabelGetStratumIS(spmap, dep, &is)); 4424 PetscCall(ISGetIndices(is, &opoints)); 4425 for (p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 4426 PetscCall(ISRestoreIndices(is, &opoints)); 4427 PetscCall(ISDestroy(&is)); 4428 } 4429 } 4430 PetscCall(DMRestoreWorkArray(dm, depth + 1, MPIU_INT, &depths)); 4431 PetscCheck(off == pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %" PetscInt_FMT " should be %" PetscInt_FMT, off, pEnd); 4432 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS)); 4433 PetscCall(PetscObjectStateGet((PetscObject)spmap, &mesh->subpointState)); 4434 } 4435 PetscFunctionReturn(PETSC_SUCCESS); 4436 } 4437 4438 /*@ 4439 DMPlexGetSubpointIS - Returns an `IS` covering the entire subdm chart with the original points as data 4440 4441 Input Parameter: 4442 . dm - The submesh `DM` 4443 4444 Output Parameter: 4445 . subpointIS - The `IS` of all the points from the original mesh in this submesh, or `NULL` if this is not a submesh 4446 4447 Level: developer 4448 4449 Note: 4450 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. 4451 4452 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()` 4453 @*/ 4454 PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS) 4455 { 4456 DM_Plex *mesh = (DM_Plex *)dm->data; 4457 DMLabel spmap; 4458 PetscObjectState state; 4459 4460 PetscFunctionBegin; 4461 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4462 PetscAssertPointer(subpointIS, 2); 4463 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 4464 PetscCall(PetscObjectStateGet((PetscObject)spmap, &state)); 4465 if (state != mesh->subpointState || !mesh->subpointIS) PetscCall(DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS)); 4466 *subpointIS = mesh->subpointIS; 4467 PetscFunctionReturn(PETSC_SUCCESS); 4468 } 4469 4470 /*@ 4471 DMGetEnclosureRelation - Get the relationship between `dmA` and `dmB` 4472 4473 Input Parameters: 4474 + dmA - The first `DM` 4475 - dmB - The second `DM` 4476 4477 Output Parameter: 4478 . rel - The relation of `dmA` to `dmB` 4479 4480 Level: intermediate 4481 4482 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMGetEnclosurePoint()` 4483 @*/ 4484 PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel) 4485 { 4486 DM plexA, plexB, sdm; 4487 DMLabel spmap; 4488 PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB; 4489 4490 PetscFunctionBegin; 4491 PetscAssertPointer(rel, 3); 4492 *rel = DM_ENC_NONE; 4493 if (!dmA || !dmB) PetscFunctionReturn(PETSC_SUCCESS); 4494 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 4495 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 4496 if (dmA == dmB) { 4497 *rel = DM_ENC_EQUALITY; 4498 PetscFunctionReturn(PETSC_SUCCESS); 4499 } 4500 PetscCall(DMConvert(dmA, DMPLEX, &plexA)); 4501 PetscCall(DMConvert(dmB, DMPLEX, &plexB)); 4502 PetscCall(DMPlexGetChart(plexA, &pStartA, &pEndA)); 4503 PetscCall(DMPlexGetChart(plexB, &pStartB, &pEndB)); 4504 /* Assumption 1: subDMs have smaller charts than the DMs that they originate from 4505 - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */ 4506 if ((pStartA == pStartB) && (pEndA == pEndB)) { 4507 *rel = DM_ENC_EQUALITY; 4508 goto end; 4509 } 4510 NpA = pEndA - pStartA; 4511 NpB = pEndB - pStartB; 4512 if (NpA == NpB) goto end; 4513 sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */ 4514 PetscCall(DMPlexGetSubpointMap(sdm, &spmap)); 4515 if (!spmap) goto end; 4516 /* TODO Check the space mapped to by subpointMap is same size as dm */ 4517 if (NpA > NpB) { 4518 *rel = DM_ENC_SUPERMESH; 4519 } else { 4520 *rel = DM_ENC_SUBMESH; 4521 } 4522 end: 4523 PetscCall(DMDestroy(&plexA)); 4524 PetscCall(DMDestroy(&plexB)); 4525 PetscFunctionReturn(PETSC_SUCCESS); 4526 } 4527 4528 /*@ 4529 DMGetEnclosurePoint - Get the point `pA` in `dmA` which corresponds to the point `pB` in `dmB` 4530 4531 Input Parameters: 4532 + dmA - The first `DM` 4533 . dmB - The second `DM` 4534 . etype - The type of enclosure relation that `dmA` has to `dmB` 4535 - pB - A point of `dmB` 4536 4537 Output Parameter: 4538 . pA - The corresponding point of `dmA` 4539 4540 Level: intermediate 4541 4542 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMGetEnclosureRelation()` 4543 @*/ 4544 PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA) 4545 { 4546 DM sdm; 4547 IS subpointIS; 4548 const PetscInt *subpoints; 4549 PetscInt numSubpoints; 4550 4551 PetscFunctionBegin; 4552 /* TODO Cache the IS, making it look like an index */ 4553 switch (etype) { 4554 case DM_ENC_SUPERMESH: 4555 sdm = dmB; 4556 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4557 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4558 *pA = subpoints[pB]; 4559 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4560 break; 4561 case DM_ENC_SUBMESH: 4562 sdm = dmA; 4563 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4564 PetscCall(ISGetLocalSize(subpointIS, &numSubpoints)); 4565 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4566 PetscCall(PetscFindInt(pB, numSubpoints, subpoints, pA)); 4567 if (*pA < 0) { 4568 PetscCall(DMViewFromOptions(dmA, NULL, "-dm_enc_A_view")); 4569 PetscCall(DMViewFromOptions(dmB, NULL, "-dm_enc_B_view")); 4570 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB); 4571 } 4572 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4573 break; 4574 case DM_ENC_EQUALITY: 4575 case DM_ENC_NONE: 4576 *pA = pB; 4577 break; 4578 case DM_ENC_UNKNOWN: { 4579 DMEnclosureType enc; 4580 4581 PetscCall(DMGetEnclosureRelation(dmA, dmB, &enc)); 4582 PetscCall(DMGetEnclosurePoint(dmA, dmB, enc, pB, pA)); 4583 } break; 4584 default: 4585 SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int)etype); 4586 } 4587 PetscFunctionReturn(PETSC_SUCCESS); 4588 } 4589