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, *newLocalPoints, *newOwners; 3449 const PetscInt *localPoints, *subpoints; 3450 PetscInt *slocalPoints; 3451 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3452 PetscMPIInt rank; 3453 3454 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 3455 PetscCall(DMGetPointSF(dm, &sfPoint)); 3456 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 3457 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3458 PetscCall(DMPlexGetChart(subdm, NULL, &numSubroots)); 3459 PetscCall(DMPlexGetSubpointIS(subdm, &subpIS)); 3460 if (subpIS) { 3461 PetscCall(ISGetIndices(subpIS, &subpoints)); 3462 PetscCall(ISGetLocalSize(subpIS, &numSubpoints)); 3463 } 3464 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3465 if (numRoots >= 0) { 3466 PetscCall(PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners)); 3467 for (p = 0; p < pEnd-pStart; ++p) { 3468 newLocalPoints[p].rank = -2; 3469 newLocalPoints[p].index = -2; 3470 } 3471 /* Set subleaves */ 3472 for (l = 0; l < numLeaves; ++l) { 3473 const PetscInt point = localPoints[l]; 3474 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3475 3476 if (subpoint < 0) continue; 3477 newLocalPoints[point-pStart].rank = rank; 3478 newLocalPoints[point-pStart].index = subpoint; 3479 ++numSubleaves; 3480 } 3481 /* Must put in owned subpoints */ 3482 for (p = pStart; p < pEnd; ++p) { 3483 newOwners[p-pStart].rank = -3; 3484 newOwners[p-pStart].index = -3; 3485 } 3486 for (p = 0; p < numSubpoints; ++p) { 3487 newOwners[subpoints[p]-pStart].rank = rank; 3488 newOwners[subpoints[p]-pStart].index = p; 3489 } 3490 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3491 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3492 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3493 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3494 PetscCall(PetscMalloc1(numSubleaves, &slocalPoints)); 3495 PetscCall(PetscMalloc1(numSubleaves, &sremotePoints)); 3496 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3497 const PetscInt point = localPoints[l]; 3498 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3499 3500 if (subpoint < 0) continue; 3501 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3502 slocalPoints[sl] = subpoint; 3503 sremotePoints[sl].rank = newLocalPoints[point].rank; 3504 sremotePoints[sl].index = newLocalPoints[point].index; 3505 PetscCheck(sremotePoints[sl].rank >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 3506 PetscCheck(sremotePoints[sl].index >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 3507 ++sl; 3508 } 3509 PetscCheck(sl + ll == numSubleaves,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubleaves); 3510 PetscCall(PetscFree2(newLocalPoints,newOwners)); 3511 PetscCall(PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3512 } 3513 if (subpIS) { 3514 PetscCall(ISRestoreIndices(subpIS, &subpoints)); 3515 } 3516 } 3517 /* Filter labels */ 3518 PetscCall(DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm)); 3519 /* Cleanup */ 3520 for (d = 0; d <= dim; ++d) { 3521 if (subpointIS[d]) PetscCall(ISRestoreIndices(subpointIS[d], &subpoints[d])); 3522 PetscCall(ISDestroy(&subpointIS[d])); 3523 } 3524 PetscCall(PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints)); 3525 PetscFunctionReturn(0); 3526 } 3527 3528 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm) 3529 { 3530 PetscFunctionBegin; 3531 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm)); 3532 PetscFunctionReturn(0); 3533 } 3534 3535 /*@ 3536 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3537 3538 Input Parameters: 3539 + dm - The original mesh 3540 . vertexLabel - The DMLabel marking points contained in the surface 3541 . value - The label value to use 3542 - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked 3543 3544 Output Parameter: 3545 . subdm - The surface mesh 3546 3547 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3548 3549 Level: developer 3550 3551 .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()` 3552 @*/ 3553 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm) 3554 { 3555 DMPlexInterpolatedFlag interpolated; 3556 PetscInt dim, cdim; 3557 3558 PetscFunctionBegin; 3559 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3560 PetscValidPointer(subdm, 5); 3561 PetscCall(DMGetDimension(dm, &dim)); 3562 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3563 PetscCall(DMSetType(*subdm, DMPLEX)); 3564 PetscCall(DMSetDimension(*subdm, dim-1)); 3565 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3566 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 3567 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3568 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 3569 if (interpolated) { 3570 PetscCall(DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm)); 3571 } else { 3572 PetscCall(DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm)); 3573 } 3574 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3575 PetscFunctionReturn(0); 3576 } 3577 3578 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3579 { 3580 MPI_Comm comm; 3581 DMLabel subpointMap; 3582 IS subvertexIS; 3583 const PetscInt *subVertices; 3584 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3585 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3586 PetscInt c, f; 3587 3588 PetscFunctionBegin; 3589 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3590 /* Create subpointMap which marks the submesh */ 3591 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3592 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3593 PetscCall(DMLabelDestroy(&subpointMap)); 3594 PetscCall(DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm)); 3595 /* Setup chart */ 3596 PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3597 PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3598 PetscCall(DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices)); 3599 PetscCall(DMPlexSetVTKCellHeight(subdm, 1)); 3600 /* Set cone sizes */ 3601 firstSubVertex = numSubCells; 3602 firstSubFace = numSubCells+numSubVertices; 3603 newFacePoint = firstSubFace; 3604 PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3605 if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices)); 3606 for (c = 0; c < numSubCells; ++c) { 3607 PetscCall(DMPlexSetConeSize(subdm, c, 1)); 3608 } 3609 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3610 PetscCall(DMPlexSetConeSize(subdm, f, nFV)); 3611 } 3612 PetscCall(DMSetUp(subdm)); 3613 /* Create face cones */ 3614 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3615 PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 3616 for (c = 0; c < numSubCells; ++c) { 3617 const PetscInt cell = subCells[c]; 3618 const PetscInt subcell = c; 3619 const PetscInt *cone, *cells; 3620 PetscBool isHybrid; 3621 PetscInt numCells, subVertex, p, v; 3622 3623 PetscCall(DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid)); 3624 if (!isHybrid) continue; 3625 PetscCall(DMPlexGetCone(dm, cell, &cone)); 3626 for (v = 0; v < nFV; ++v) { 3627 PetscCall(PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex)); 3628 subface[v] = firstSubVertex+subVertex; 3629 } 3630 PetscCall(DMPlexSetCone(subdm, newFacePoint, subface)); 3631 PetscCall(DMPlexSetCone(subdm, subcell, &newFacePoint)); 3632 PetscCall(DMPlexGetJoin(dm, nFV, cone, &numCells, &cells)); 3633 /* Not true in parallel 3634 PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3635 for (p = 0; p < numCells; ++p) { 3636 PetscInt negsubcell; 3637 PetscBool isHybrid; 3638 3639 PetscCall(DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid)); 3640 if (isHybrid) continue; 3641 /* I know this is a crap search */ 3642 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3643 if (subCells[negsubcell] == cells[p]) break; 3644 } 3645 PetscCheck(negsubcell != numSubCells,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %" PetscInt_FMT, cell); 3646 PetscCall(DMPlexSetCone(subdm, negsubcell, &newFacePoint)); 3647 } 3648 PetscCall(DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells)); 3649 ++newFacePoint; 3650 } 3651 PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 3652 PetscCall(DMPlexSymmetrize(subdm)); 3653 PetscCall(DMPlexStratify(subdm)); 3654 /* Build coordinates */ 3655 { 3656 PetscSection coordSection, subCoordSection; 3657 Vec coordinates, subCoordinates; 3658 PetscScalar *coords, *subCoords; 3659 PetscInt cdim, numComp, coordSize, v; 3660 const char *name; 3661 3662 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3663 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3664 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3665 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3666 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3667 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3668 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3669 PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices)); 3670 for (v = 0; v < numSubVertices; ++v) { 3671 const PetscInt vertex = subVertices[v]; 3672 const PetscInt subvertex = firstSubVertex+v; 3673 PetscInt dof; 3674 3675 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3676 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3677 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3678 } 3679 PetscCall(PetscSectionSetUp(subCoordSection)); 3680 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3681 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3682 PetscCall(PetscObjectGetName((PetscObject)coordinates,&name)); 3683 PetscCall(PetscObjectSetName((PetscObject)subCoordinates,name)); 3684 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3685 PetscCall(VecSetBlockSize(subCoordinates, cdim)); 3686 PetscCall(VecSetType(subCoordinates,VECSTANDARD)); 3687 PetscCall(VecGetArray(coordinates, &coords)); 3688 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3689 for (v = 0; v < numSubVertices; ++v) { 3690 const PetscInt vertex = subVertices[v]; 3691 const PetscInt subvertex = firstSubVertex+v; 3692 PetscInt dof, off, sdof, soff, d; 3693 3694 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3695 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 3696 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3697 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3698 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); 3699 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3700 } 3701 PetscCall(VecRestoreArray(coordinates, &coords)); 3702 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3703 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3704 PetscCall(VecDestroy(&subCoordinates)); 3705 } 3706 /* Build SF */ 3707 CHKMEMQ; 3708 { 3709 PetscSF sfPoint, sfPointSub; 3710 const PetscSFNode *remotePoints; 3711 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3712 const PetscInt *localPoints; 3713 PetscInt *slocalPoints; 3714 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3715 PetscMPIInt rank; 3716 3717 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 3718 PetscCall(DMGetPointSF(dm, &sfPoint)); 3719 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 3720 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3721 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3722 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3723 if (numRoots >= 0) { 3724 /* Only vertices should be shared */ 3725 PetscCall(PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners)); 3726 for (p = 0; p < pEnd-pStart; ++p) { 3727 newLocalPoints[p].rank = -2; 3728 newLocalPoints[p].index = -2; 3729 } 3730 /* Set subleaves */ 3731 for (l = 0; l < numLeaves; ++l) { 3732 const PetscInt point = localPoints[l]; 3733 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3734 3735 PetscCheck(!(point < vStart) || !(point >= vEnd),PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %" PetscInt_FMT, point); 3736 if (subPoint < 0) continue; 3737 newLocalPoints[point-pStart].rank = rank; 3738 newLocalPoints[point-pStart].index = subPoint; 3739 ++numSubLeaves; 3740 } 3741 /* Must put in owned subpoints */ 3742 for (p = pStart; p < pEnd; ++p) { 3743 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3744 3745 if (subPoint < 0) { 3746 newOwners[p-pStart].rank = -3; 3747 newOwners[p-pStart].index = -3; 3748 } else { 3749 newOwners[p-pStart].rank = rank; 3750 newOwners[p-pStart].index = subPoint; 3751 } 3752 } 3753 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3754 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3755 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3756 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3757 PetscCall(PetscMalloc1(numSubLeaves, &slocalPoints)); 3758 PetscCall(PetscMalloc1(numSubLeaves, &sremotePoints)); 3759 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3760 const PetscInt point = localPoints[l]; 3761 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3762 3763 if (subPoint < 0) continue; 3764 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3765 slocalPoints[sl] = subPoint; 3766 sremotePoints[sl].rank = newLocalPoints[point].rank; 3767 sremotePoints[sl].index = newLocalPoints[point].index; 3768 PetscCheck(sremotePoints[sl].rank >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 3769 PetscCheck(sremotePoints[sl].index >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 3770 ++sl; 3771 } 3772 PetscCall(PetscFree2(newLocalPoints,newOwners)); 3773 PetscCheck(sl + ll == numSubLeaves,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubLeaves); 3774 PetscCall(PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3775 } 3776 } 3777 CHKMEMQ; 3778 /* Cleanup */ 3779 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices)); 3780 PetscCall(ISDestroy(&subvertexIS)); 3781 PetscCall(PetscFree(subCells)); 3782 PetscFunctionReturn(0); 3783 } 3784 3785 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3786 { 3787 DMLabel label = NULL; 3788 3789 PetscFunctionBegin; 3790 if (labelname) PetscCall(DMGetLabel(dm, labelname, &label)); 3791 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm)); 3792 PetscFunctionReturn(0); 3793 } 3794 3795 /*@C 3796 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. 3797 3798 Input Parameters: 3799 + dm - The original mesh 3800 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3801 . label - A label name, or NULL 3802 - value - A label value 3803 3804 Output Parameter: 3805 . subdm - The surface mesh 3806 3807 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3808 3809 Level: developer 3810 3811 .seealso: `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()` 3812 @*/ 3813 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3814 { 3815 PetscInt dim, cdim, depth; 3816 3817 PetscFunctionBegin; 3818 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3819 PetscValidPointer(subdm, 5); 3820 PetscCall(DMGetDimension(dm, &dim)); 3821 PetscCall(DMPlexGetDepth(dm, &depth)); 3822 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3823 PetscCall(DMSetType(*subdm, DMPLEX)); 3824 PetscCall(DMSetDimension(*subdm, dim-1)); 3825 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3826 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 3827 if (depth == dim) { 3828 PetscCall(DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm)); 3829 } else { 3830 PetscCall(DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm)); 3831 } 3832 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3833 PetscFunctionReturn(0); 3834 } 3835 3836 /*@ 3837 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3838 3839 Input Parameters: 3840 + dm - The original mesh 3841 . cellLabel - The DMLabel marking cells contained in the new mesh 3842 - value - The label value to use 3843 3844 Output Parameter: 3845 . subdm - The new mesh 3846 3847 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3848 3849 Level: developer 3850 3851 .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()` 3852 @*/ 3853 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3854 { 3855 PetscInt dim; 3856 3857 PetscFunctionBegin; 3858 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3859 PetscValidPointer(subdm, 4); 3860 PetscCall(DMGetDimension(dm, &dim)); 3861 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), subdm)); 3862 PetscCall(DMSetType(*subdm, DMPLEX)); 3863 PetscCall(DMSetDimension(*subdm, dim)); 3864 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3865 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm)); 3866 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3867 PetscFunctionReturn(0); 3868 } 3869 3870 /*@ 3871 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3872 3873 Input Parameter: 3874 . dm - The submesh DM 3875 3876 Output Parameter: 3877 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3878 3879 Level: developer 3880 3881 .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 3882 @*/ 3883 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3884 { 3885 PetscFunctionBegin; 3886 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3887 PetscValidPointer(subpointMap, 2); 3888 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3889 PetscFunctionReturn(0); 3890 } 3891 3892 /*@ 3893 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3894 3895 Input Parameters: 3896 + dm - The submesh DM 3897 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3898 3899 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3900 3901 Level: developer 3902 3903 .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 3904 @*/ 3905 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3906 { 3907 DM_Plex *mesh = (DM_Plex *) dm->data; 3908 DMLabel tmp; 3909 3910 PetscFunctionBegin; 3911 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3912 tmp = mesh->subpointMap; 3913 mesh->subpointMap = subpointMap; 3914 PetscCall(PetscObjectReference((PetscObject) mesh->subpointMap)); 3915 PetscCall(DMLabelDestroy(&tmp)); 3916 PetscFunctionReturn(0); 3917 } 3918 3919 static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS) 3920 { 3921 DMLabel spmap; 3922 PetscInt depth, d; 3923 3924 PetscFunctionBegin; 3925 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 3926 PetscCall(DMPlexGetDepth(dm, &depth)); 3927 if (spmap && depth >= 0) { 3928 DM_Plex *mesh = (DM_Plex *) dm->data; 3929 PetscInt *points, *depths; 3930 PetscInt pStart, pEnd, p, off; 3931 3932 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3933 PetscCheck(!pStart,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %" PetscInt_FMT, pStart); 3934 PetscCall(PetscMalloc1(pEnd, &points)); 3935 PetscCall(DMGetWorkArray(dm, depth+1, MPIU_INT, &depths)); 3936 depths[0] = depth; 3937 depths[1] = 0; 3938 for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3939 for (d = 0, off = 0; d <= depth; ++d) { 3940 const PetscInt dep = depths[d]; 3941 PetscInt depStart, depEnd, n; 3942 3943 PetscCall(DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd)); 3944 PetscCall(DMLabelGetStratumSize(spmap, dep, &n)); 3945 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3946 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); 3947 } else { 3948 if (!n) { 3949 if (d == 0) { 3950 /* Missing cells */ 3951 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3952 } else { 3953 /* Missing faces */ 3954 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3955 } 3956 } 3957 } 3958 if (n) { 3959 IS is; 3960 const PetscInt *opoints; 3961 3962 PetscCall(DMLabelGetStratumIS(spmap, dep, &is)); 3963 PetscCall(ISGetIndices(is, &opoints)); 3964 for (p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3965 PetscCall(ISRestoreIndices(is, &opoints)); 3966 PetscCall(ISDestroy(&is)); 3967 } 3968 } 3969 PetscCall(DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths)); 3970 PetscCheck(off == pEnd,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %" PetscInt_FMT " should be %" PetscInt_FMT, off, pEnd); 3971 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS)); 3972 PetscCall(PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState)); 3973 } 3974 PetscFunctionReturn(0); 3975 } 3976 3977 /*@ 3978 DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data 3979 3980 Input Parameter: 3981 . dm - The submesh DM 3982 3983 Output Parameter: 3984 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3985 3986 Note: This IS is guaranteed to be sorted by the construction of the submesh 3987 3988 Level: developer 3989 3990 .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()` 3991 @*/ 3992 PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS) 3993 { 3994 DM_Plex *mesh = (DM_Plex *) dm->data; 3995 DMLabel spmap; 3996 PetscObjectState state; 3997 3998 PetscFunctionBegin; 3999 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4000 PetscValidPointer(subpointIS, 2); 4001 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 4002 PetscCall(PetscObjectStateGet((PetscObject) spmap, &state)); 4003 if (state != mesh->subpointState || !mesh->subpointIS) PetscCall(DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS)); 4004 *subpointIS = mesh->subpointIS; 4005 PetscFunctionReturn(0); 4006 } 4007 4008 /*@ 4009 DMGetEnclosureRelation - Get the relationship between dmA and dmB 4010 4011 Input Parameters: 4012 + dmA - The first DM 4013 - dmB - The second DM 4014 4015 Output Parameter: 4016 . rel - The relation of dmA to dmB 4017 4018 Level: intermediate 4019 4020 .seealso: `DMGetEnclosurePoint()` 4021 @*/ 4022 PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel) 4023 { 4024 DM plexA, plexB, sdm; 4025 DMLabel spmap; 4026 PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB; 4027 4028 PetscFunctionBegin; 4029 PetscValidPointer(rel, 3); 4030 *rel = DM_ENC_NONE; 4031 if (!dmA || !dmB) PetscFunctionReturn(0); 4032 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 4033 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 4034 if (dmA == dmB) {*rel = DM_ENC_EQUALITY; PetscFunctionReturn(0);} 4035 PetscCall(DMConvert(dmA, DMPLEX, &plexA)); 4036 PetscCall(DMConvert(dmB, DMPLEX, &plexB)); 4037 PetscCall(DMPlexGetChart(plexA, &pStartA, &pEndA)); 4038 PetscCall(DMPlexGetChart(plexB, &pStartB, &pEndB)); 4039 /* Assumption 1: subDMs have smaller charts than the DMs that they originate from 4040 - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */ 4041 if ((pStartA == pStartB) && (pEndA == pEndB)) { 4042 *rel = DM_ENC_EQUALITY; 4043 goto end; 4044 } 4045 NpA = pEndA - pStartA; 4046 NpB = pEndB - pStartB; 4047 if (NpA == NpB) goto end; 4048 sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */ 4049 PetscCall(DMPlexGetSubpointMap(sdm, &spmap)); 4050 if (!spmap) goto end; 4051 /* TODO Check the space mapped to by subpointMap is same size as dm */ 4052 if (NpA > NpB) { 4053 *rel = DM_ENC_SUPERMESH; 4054 } else { 4055 *rel = DM_ENC_SUBMESH; 4056 } 4057 end: 4058 PetscCall(DMDestroy(&plexA)); 4059 PetscCall(DMDestroy(&plexB)); 4060 PetscFunctionReturn(0); 4061 } 4062 4063 /*@ 4064 DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB 4065 4066 Input Parameters: 4067 + dmA - The first DM 4068 . dmB - The second DM 4069 . etype - The type of enclosure relation that dmA has to dmB 4070 - pB - A point of dmB 4071 4072 Output Parameter: 4073 . pA - The corresponding point of dmA 4074 4075 Level: intermediate 4076 4077 .seealso: `DMGetEnclosureRelation()` 4078 @*/ 4079 PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA) 4080 { 4081 DM sdm; 4082 IS subpointIS; 4083 const PetscInt *subpoints; 4084 PetscInt numSubpoints; 4085 4086 PetscFunctionBegin; 4087 /* TODO Cache the IS, making it look like an index */ 4088 switch (etype) { 4089 case DM_ENC_SUPERMESH: 4090 sdm = dmB; 4091 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4092 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4093 *pA = subpoints[pB]; 4094 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4095 break; 4096 case DM_ENC_SUBMESH: 4097 sdm = dmA; 4098 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4099 PetscCall(ISGetLocalSize(subpointIS, &numSubpoints)); 4100 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4101 PetscCall(PetscFindInt(pB, numSubpoints, subpoints, pA)); 4102 if (*pA < 0) { 4103 PetscCall(DMViewFromOptions(dmA, NULL, "-dm_enc_A_view")); 4104 PetscCall(DMViewFromOptions(dmB, NULL, "-dm_enc_B_view")); 4105 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB); 4106 } 4107 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4108 break; 4109 case DM_ENC_EQUALITY: 4110 case DM_ENC_NONE: 4111 *pA = pB;break; 4112 case DM_ENC_UNKNOWN: 4113 { 4114 DMEnclosureType enc; 4115 4116 PetscCall(DMGetEnclosureRelation(dmA, dmB, &enc)); 4117 PetscCall(DMGetEnclosurePoint(dmA, dmB, enc, pB, pA)); 4118 } 4119 break; 4120 default: SETERRQ(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype); 4121 } 4122 PetscFunctionReturn(0); 4123 } 4124