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