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