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