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 PetscInt Nl, l, d; 3184 3185 PetscFunctionBegin; 3186 PetscCall(DMGetNumLabels(dm, &Nl)); 3187 for (l = 0; l < Nl; ++l) { 3188 DMLabel label, newlabel; 3189 const char *lname; 3190 PetscBool isDepth, isDim, isCelltype, isVTK; 3191 IS valueIS; 3192 const PetscInt *values; 3193 PetscInt Nv, v; 3194 3195 PetscCall(DMGetLabelName(dm, l, &lname)); 3196 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 3197 PetscCall(PetscStrcmp(lname, "dim", &isDim)); 3198 PetscCall(PetscStrcmp(lname, "celltype", &isCelltype)); 3199 PetscCall(PetscStrcmp(lname, "vtk", &isVTK)); 3200 if (isDepth || isDim || isCelltype || isVTK) continue; 3201 PetscCall(DMCreateLabel(subdm, lname)); 3202 PetscCall(DMGetLabel(dm, lname, &label)); 3203 PetscCall(DMGetLabel(subdm, lname, &newlabel)); 3204 PetscCall(DMLabelGetDefaultValue(label, &v)); 3205 PetscCall(DMLabelSetDefaultValue(newlabel, v)); 3206 PetscCall(DMLabelGetValueIS(label, &valueIS)); 3207 PetscCall(ISGetLocalSize(valueIS, &Nv)); 3208 PetscCall(ISGetIndices(valueIS, &values)); 3209 for (v = 0; v < Nv; ++v) { 3210 IS pointIS; 3211 const PetscInt *points; 3212 PetscInt Np, p; 3213 3214 PetscCall(DMLabelGetStratumIS(label, values[v], &pointIS)); 3215 PetscCall(ISGetLocalSize(pointIS, &Np)); 3216 PetscCall(ISGetIndices(pointIS, &points)); 3217 for (p = 0; p < Np; ++p) { 3218 const PetscInt point = points[p]; 3219 PetscInt subp; 3220 3221 PetscCall(DMPlexGetPointDepth(dm, point, &d)); 3222 subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]); 3223 if (subp >= 0) PetscCall(DMLabelSetValue(newlabel, subp, values[v])); 3224 } 3225 PetscCall(ISRestoreIndices(pointIS, &points)); 3226 PetscCall(ISDestroy(&pointIS)); 3227 } 3228 PetscCall(ISRestoreIndices(valueIS, &values)); 3229 PetscCall(ISDestroy(&valueIS)); 3230 } 3231 PetscFunctionReturn(0); 3232 } 3233 3234 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 3235 { 3236 MPI_Comm comm; 3237 DMLabel subpointMap; 3238 IS *subpointIS; 3239 const PetscInt **subpoints; 3240 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 3241 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 3242 PetscMPIInt rank; 3243 3244 PetscFunctionBegin; 3245 PetscCall(PetscObjectGetComm((PetscObject)dm,&comm)); 3246 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 3247 /* Create subpointMap which marks the submesh */ 3248 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3249 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3250 if (cellHeight) { 3251 if (isCohesive) PetscCall(DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm)); 3252 else PetscCall(DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm)); 3253 } else { 3254 DMLabel depth; 3255 IS pointIS; 3256 const PetscInt *points; 3257 PetscInt numPoints=0; 3258 3259 PetscCall(DMPlexGetDepthLabel(dm, &depth)); 3260 PetscCall(DMLabelGetStratumIS(label, value, &pointIS)); 3261 if (pointIS) { 3262 PetscCall(ISGetIndices(pointIS, &points)); 3263 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 3264 } 3265 for (p = 0; p < numPoints; ++p) { 3266 PetscInt *closure = NULL; 3267 PetscInt closureSize, c, pdim; 3268 3269 PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3270 for (c = 0; c < closureSize*2; c += 2) { 3271 PetscCall(DMLabelGetValue(depth, closure[c], &pdim)); 3272 PetscCall(DMLabelSetValue(subpointMap, closure[c], pdim)); 3273 } 3274 PetscCall(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3275 } 3276 if (pointIS) PetscCall(ISRestoreIndices(pointIS, &points)); 3277 PetscCall(ISDestroy(&pointIS)); 3278 } 3279 /* Setup chart */ 3280 PetscCall(DMGetDimension(dm, &dim)); 3281 PetscCall(PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints)); 3282 for (d = 0; d <= dim; ++d) { 3283 PetscCall(DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d])); 3284 totSubPoints += numSubPoints[d]; 3285 } 3286 PetscCall(DMPlexSetChart(subdm, 0, totSubPoints)); 3287 PetscCall(DMPlexSetVTKCellHeight(subdm, cellHeight)); 3288 /* Set cone sizes */ 3289 firstSubPoint[dim] = 0; 3290 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 3291 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 3292 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 3293 for (d = 0; d <= dim; ++d) { 3294 PetscCall(DMLabelGetStratumIS(subpointMap, d, &subpointIS[d])); 3295 if (subpointIS[d]) PetscCall(ISGetIndices(subpointIS[d], &subpoints[d])); 3296 } 3297 /* We do not want this label automatically computed, instead we compute it here */ 3298 PetscCall(DMCreateLabel(subdm, "celltype")); 3299 for (d = 0; d <= dim; ++d) { 3300 for (p = 0; p < numSubPoints[d]; ++p) { 3301 const PetscInt point = subpoints[d][p]; 3302 const PetscInt subpoint = firstSubPoint[d] + p; 3303 const PetscInt *cone; 3304 PetscInt coneSize; 3305 3306 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 3307 if (cellHeight && (d == dim)) { 3308 PetscInt coneSizeNew, c, val; 3309 3310 PetscCall(DMPlexGetCone(dm, point, &cone)); 3311 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3312 PetscCall(DMLabelGetValue(subpointMap, cone[c], &val)); 3313 if (val >= 0) coneSizeNew++; 3314 } 3315 PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSizeNew)); 3316 PetscCall(DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST)); 3317 } else { 3318 DMPolytopeType ct; 3319 3320 PetscCall(DMPlexSetConeSize(subdm, subpoint, coneSize)); 3321 PetscCall(DMPlexGetCellType(dm, point, &ct)); 3322 PetscCall(DMPlexSetCellType(subdm, subpoint, ct)); 3323 } 3324 } 3325 } 3326 PetscCall(DMLabelDestroy(&subpointMap)); 3327 PetscCall(DMSetUp(subdm)); 3328 /* Set cones */ 3329 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3330 PetscCall(PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew)); 3331 for (d = 0; d <= dim; ++d) { 3332 for (p = 0; p < numSubPoints[d]; ++p) { 3333 const PetscInt point = subpoints[d][p]; 3334 const PetscInt subpoint = firstSubPoint[d] + p; 3335 const PetscInt *cone, *ornt; 3336 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 3337 3338 if (d == dim-1) { 3339 const PetscInt *support, *cone, *ornt; 3340 PetscInt supportSize, coneSize, s, subc; 3341 3342 PetscCall(DMPlexGetSupport(dm, point, &support)); 3343 PetscCall(DMPlexGetSupportSize(dm, point, &supportSize)); 3344 for (s = 0; s < supportSize; ++s) { 3345 PetscBool isHybrid; 3346 3347 PetscCall(DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid)); 3348 if (!isHybrid) continue; 3349 PetscCall(PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc)); 3350 if (subc >= 0) { 3351 const PetscInt ccell = subpoints[d+1][subc]; 3352 3353 PetscCall(DMPlexGetCone(dm, ccell, &cone)); 3354 PetscCall(DMPlexGetConeSize(dm, ccell, &coneSize)); 3355 PetscCall(DMPlexGetConeOrientation(dm, ccell, &ornt)); 3356 for (c = 0; c < coneSize; ++c) { 3357 if (cone[c] == point) { 3358 fornt = ornt[c]; 3359 break; 3360 } 3361 } 3362 break; 3363 } 3364 } 3365 } 3366 PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 3367 PetscCall(DMPlexGetConeSize(subdm, subpoint, &subconeSize)); 3368 PetscCall(DMPlexGetCone(dm, point, &cone)); 3369 PetscCall(DMPlexGetConeOrientation(dm, point, &ornt)); 3370 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3371 PetscCall(PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc)); 3372 if (subc >= 0) { 3373 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 3374 orntNew[coneSizeNew] = ornt[c]; 3375 ++coneSizeNew; 3376 } 3377 } 3378 PetscCheck(coneSizeNew == subconeSize,comm, PETSC_ERR_PLIB, "Number of cone points located %" PetscInt_FMT " does not match subcone size %" PetscInt_FMT, coneSizeNew, subconeSize); 3379 PetscCall(DMPlexSetCone(subdm, subpoint, coneNew)); 3380 PetscCall(DMPlexSetConeOrientation(subdm, subpoint, orntNew)); 3381 if (fornt < 0) PetscCall(DMPlexOrientPoint(subdm, subpoint, fornt)); 3382 } 3383 } 3384 PetscCall(PetscFree2(coneNew,orntNew)); 3385 PetscCall(DMPlexSymmetrize(subdm)); 3386 PetscCall(DMPlexStratify(subdm)); 3387 /* Build coordinates */ 3388 { 3389 PetscSection coordSection, subCoordSection; 3390 Vec coordinates, subCoordinates; 3391 PetscScalar *coords, *subCoords; 3392 PetscInt cdim, numComp, coordSize; 3393 const char *name; 3394 3395 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3396 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3397 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3398 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3399 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3400 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3401 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3402 PetscCall(PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0])); 3403 for (v = 0; v < numSubPoints[0]; ++v) { 3404 const PetscInt vertex = subpoints[0][v]; 3405 const PetscInt subvertex = firstSubPoint[0]+v; 3406 PetscInt dof; 3407 3408 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3409 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3410 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3411 } 3412 PetscCall(PetscSectionSetUp(subCoordSection)); 3413 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3414 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3415 PetscCall(PetscObjectGetName((PetscObject)coordinates,&name)); 3416 PetscCall(PetscObjectSetName((PetscObject)subCoordinates,name)); 3417 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3418 PetscCall(VecSetBlockSize(subCoordinates, cdim)); 3419 PetscCall(VecSetType(subCoordinates,VECSTANDARD)); 3420 PetscCall(VecGetArray(coordinates, &coords)); 3421 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3422 for (v = 0; v < numSubPoints[0]; ++v) { 3423 const PetscInt vertex = subpoints[0][v]; 3424 const PetscInt subvertex = firstSubPoint[0]+v; 3425 PetscInt dof, off, sdof, soff, d; 3426 3427 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3428 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 3429 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3430 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3431 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); 3432 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3433 } 3434 PetscCall(VecRestoreArray(coordinates, &coords)); 3435 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3436 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3437 PetscCall(VecDestroy(&subCoordinates)); 3438 } 3439 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3440 { 3441 PetscSF sfPoint, sfPointSub; 3442 IS subpIS; 3443 const PetscSFNode *remotePoints; 3444 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3445 const PetscInt *localPoints, *subpoints; 3446 PetscInt *slocalPoints; 3447 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3448 PetscMPIInt rank; 3449 3450 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 3451 PetscCall(DMGetPointSF(dm, &sfPoint)); 3452 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 3453 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3454 PetscCall(DMPlexGetChart(subdm, NULL, &numSubroots)); 3455 PetscCall(DMPlexGetSubpointIS(subdm, &subpIS)); 3456 if (subpIS) { 3457 PetscCall(ISGetIndices(subpIS, &subpoints)); 3458 PetscCall(ISGetLocalSize(subpIS, &numSubpoints)); 3459 } 3460 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3461 if (numRoots >= 0) { 3462 PetscCall(PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners)); 3463 for (p = 0; p < pEnd-pStart; ++p) { 3464 newLocalPoints[p].rank = -2; 3465 newLocalPoints[p].index = -2; 3466 } 3467 /* Set subleaves */ 3468 for (l = 0; l < numLeaves; ++l) { 3469 const PetscInt point = localPoints[l]; 3470 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3471 3472 if (subpoint < 0) continue; 3473 newLocalPoints[point-pStart].rank = rank; 3474 newLocalPoints[point-pStart].index = subpoint; 3475 ++numSubleaves; 3476 } 3477 /* Must put in owned subpoints */ 3478 for (p = pStart; p < pEnd; ++p) { 3479 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3480 3481 if (subpoint < 0) { 3482 newOwners[p-pStart].rank = -3; 3483 newOwners[p-pStart].index = -3; 3484 } else { 3485 newOwners[p-pStart].rank = rank; 3486 newOwners[p-pStart].index = subpoint; 3487 } 3488 } 3489 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3490 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3491 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3492 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3493 PetscCall(PetscMalloc1(numSubleaves, &slocalPoints)); 3494 PetscCall(PetscMalloc1(numSubleaves, &sremotePoints)); 3495 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3496 const PetscInt point = localPoints[l]; 3497 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3498 3499 if (subpoint < 0) continue; 3500 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3501 slocalPoints[sl] = subpoint; 3502 sremotePoints[sl].rank = newLocalPoints[point].rank; 3503 sremotePoints[sl].index = newLocalPoints[point].index; 3504 PetscCheck(sremotePoints[sl].rank >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 3505 PetscCheck(sremotePoints[sl].index >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 3506 ++sl; 3507 } 3508 PetscCheck(sl + ll == numSubleaves,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubleaves); 3509 PetscCall(PetscFree2(newLocalPoints,newOwners)); 3510 PetscCall(PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3511 } 3512 if (subpIS) { 3513 PetscCall(ISRestoreIndices(subpIS, &subpoints)); 3514 } 3515 } 3516 /* Filter labels */ 3517 PetscCall(DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm)); 3518 /* Cleanup */ 3519 for (d = 0; d <= dim; ++d) { 3520 if (subpointIS[d]) PetscCall(ISRestoreIndices(subpointIS[d], &subpoints[d])); 3521 PetscCall(ISDestroy(&subpointIS[d])); 3522 } 3523 PetscCall(PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints)); 3524 PetscFunctionReturn(0); 3525 } 3526 3527 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm) 3528 { 3529 PetscFunctionBegin; 3530 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm)); 3531 PetscFunctionReturn(0); 3532 } 3533 3534 /*@ 3535 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3536 3537 Input Parameters: 3538 + dm - The original mesh 3539 . vertexLabel - The DMLabel marking points contained in the surface 3540 . value - The label value to use 3541 - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked 3542 3543 Output Parameter: 3544 . subdm - The surface mesh 3545 3546 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3547 3548 Level: developer 3549 3550 .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()` 3551 @*/ 3552 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm) 3553 { 3554 DMPlexInterpolatedFlag interpolated; 3555 PetscInt dim, cdim; 3556 3557 PetscFunctionBegin; 3558 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3559 PetscValidPointer(subdm, 5); 3560 PetscCall(DMGetDimension(dm, &dim)); 3561 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3562 PetscCall(DMSetType(*subdm, DMPLEX)); 3563 PetscCall(DMSetDimension(*subdm, dim-1)); 3564 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3565 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 3566 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 3567 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 3568 if (interpolated) { 3569 PetscCall(DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm)); 3570 } else { 3571 PetscCall(DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm)); 3572 } 3573 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3574 PetscFunctionReturn(0); 3575 } 3576 3577 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3578 { 3579 MPI_Comm comm; 3580 DMLabel subpointMap; 3581 IS subvertexIS; 3582 const PetscInt *subVertices; 3583 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3584 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3585 PetscInt c, f; 3586 3587 PetscFunctionBegin; 3588 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3589 /* Create subpointMap which marks the submesh */ 3590 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3591 PetscCall(DMPlexSetSubpointMap(subdm, subpointMap)); 3592 PetscCall(DMLabelDestroy(&subpointMap)); 3593 PetscCall(DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm)); 3594 /* Setup chart */ 3595 PetscCall(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3596 PetscCall(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3597 PetscCall(DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices)); 3598 PetscCall(DMPlexSetVTKCellHeight(subdm, 1)); 3599 /* Set cone sizes */ 3600 firstSubVertex = numSubCells; 3601 firstSubFace = numSubCells+numSubVertices; 3602 newFacePoint = firstSubFace; 3603 PetscCall(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3604 if (subvertexIS) PetscCall(ISGetIndices(subvertexIS, &subVertices)); 3605 for (c = 0; c < numSubCells; ++c) { 3606 PetscCall(DMPlexSetConeSize(subdm, c, 1)); 3607 } 3608 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3609 PetscCall(DMPlexSetConeSize(subdm, f, nFV)); 3610 } 3611 PetscCall(DMSetUp(subdm)); 3612 /* Create face cones */ 3613 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3614 PetscCall(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 3615 for (c = 0; c < numSubCells; ++c) { 3616 const PetscInt cell = subCells[c]; 3617 const PetscInt subcell = c; 3618 const PetscInt *cone, *cells; 3619 PetscBool isHybrid; 3620 PetscInt numCells, subVertex, p, v; 3621 3622 PetscCall(DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid)); 3623 if (!isHybrid) continue; 3624 PetscCall(DMPlexGetCone(dm, cell, &cone)); 3625 for (v = 0; v < nFV; ++v) { 3626 PetscCall(PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex)); 3627 subface[v] = firstSubVertex+subVertex; 3628 } 3629 PetscCall(DMPlexSetCone(subdm, newFacePoint, subface)); 3630 PetscCall(DMPlexSetCone(subdm, subcell, &newFacePoint)); 3631 PetscCall(DMPlexGetJoin(dm, nFV, cone, &numCells, &cells)); 3632 /* Not true in parallel 3633 PetscCheck(numCells == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3634 for (p = 0; p < numCells; ++p) { 3635 PetscInt negsubcell; 3636 PetscBool isHybrid; 3637 3638 PetscCall(DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid)); 3639 if (isHybrid) continue; 3640 /* I know this is a crap search */ 3641 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3642 if (subCells[negsubcell] == cells[p]) break; 3643 } 3644 PetscCheck(negsubcell != numSubCells,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %" PetscInt_FMT, cell); 3645 PetscCall(DMPlexSetCone(subdm, negsubcell, &newFacePoint)); 3646 } 3647 PetscCall(DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells)); 3648 ++newFacePoint; 3649 } 3650 PetscCall(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 3651 PetscCall(DMPlexSymmetrize(subdm)); 3652 PetscCall(DMPlexStratify(subdm)); 3653 /* Build coordinates */ 3654 { 3655 PetscSection coordSection, subCoordSection; 3656 Vec coordinates, subCoordinates; 3657 PetscScalar *coords, *subCoords; 3658 PetscInt cdim, numComp, coordSize, v; 3659 const char *name; 3660 3661 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3662 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 3663 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3664 PetscCall(DMGetCoordinateSection(subdm, &subCoordSection)); 3665 PetscCall(PetscSectionSetNumFields(subCoordSection, 1)); 3666 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3667 PetscCall(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3668 PetscCall(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices)); 3669 for (v = 0; v < numSubVertices; ++v) { 3670 const PetscInt vertex = subVertices[v]; 3671 const PetscInt subvertex = firstSubVertex+v; 3672 PetscInt dof; 3673 3674 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3675 PetscCall(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3676 PetscCall(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3677 } 3678 PetscCall(PetscSectionSetUp(subCoordSection)); 3679 PetscCall(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3680 PetscCall(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3681 PetscCall(PetscObjectGetName((PetscObject)coordinates,&name)); 3682 PetscCall(PetscObjectSetName((PetscObject)subCoordinates,name)); 3683 PetscCall(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3684 PetscCall(VecSetBlockSize(subCoordinates, cdim)); 3685 PetscCall(VecSetType(subCoordinates,VECSTANDARD)); 3686 PetscCall(VecGetArray(coordinates, &coords)); 3687 PetscCall(VecGetArray(subCoordinates, &subCoords)); 3688 for (v = 0; v < numSubVertices; ++v) { 3689 const PetscInt vertex = subVertices[v]; 3690 const PetscInt subvertex = firstSubVertex+v; 3691 PetscInt dof, off, sdof, soff, d; 3692 3693 PetscCall(PetscSectionGetDof(coordSection, vertex, &dof)); 3694 PetscCall(PetscSectionGetOffset(coordSection, vertex, &off)); 3695 PetscCall(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3696 PetscCall(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3697 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); 3698 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3699 } 3700 PetscCall(VecRestoreArray(coordinates, &coords)); 3701 PetscCall(VecRestoreArray(subCoordinates, &subCoords)); 3702 PetscCall(DMSetCoordinatesLocal(subdm, subCoordinates)); 3703 PetscCall(VecDestroy(&subCoordinates)); 3704 } 3705 /* Build SF */ 3706 CHKMEMQ; 3707 { 3708 PetscSF sfPoint, sfPointSub; 3709 const PetscSFNode *remotePoints; 3710 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3711 const PetscInt *localPoints; 3712 PetscInt *slocalPoints; 3713 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3714 PetscMPIInt rank; 3715 3716 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 3717 PetscCall(DMGetPointSF(dm, &sfPoint)); 3718 PetscCall(DMGetPointSF(subdm, &sfPointSub)); 3719 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3720 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3721 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3722 if (numRoots >= 0) { 3723 /* Only vertices should be shared */ 3724 PetscCall(PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners)); 3725 for (p = 0; p < pEnd-pStart; ++p) { 3726 newLocalPoints[p].rank = -2; 3727 newLocalPoints[p].index = -2; 3728 } 3729 /* Set subleaves */ 3730 for (l = 0; l < numLeaves; ++l) { 3731 const PetscInt point = localPoints[l]; 3732 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3733 3734 PetscCheck(!(point < vStart) || !(point >= vEnd),PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %" PetscInt_FMT, point); 3735 if (subPoint < 0) continue; 3736 newLocalPoints[point-pStart].rank = rank; 3737 newLocalPoints[point-pStart].index = subPoint; 3738 ++numSubLeaves; 3739 } 3740 /* Must put in owned subpoints */ 3741 for (p = pStart; p < pEnd; ++p) { 3742 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3743 3744 if (subPoint < 0) { 3745 newOwners[p-pStart].rank = -3; 3746 newOwners[p-pStart].index = -3; 3747 } else { 3748 newOwners[p-pStart].rank = rank; 3749 newOwners[p-pStart].index = subPoint; 3750 } 3751 } 3752 PetscCall(PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3753 PetscCall(PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3754 PetscCall(PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3755 PetscCall(PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3756 PetscCall(PetscMalloc1(numSubLeaves, &slocalPoints)); 3757 PetscCall(PetscMalloc1(numSubLeaves, &sremotePoints)); 3758 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3759 const PetscInt point = localPoints[l]; 3760 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3761 3762 if (subPoint < 0) continue; 3763 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3764 slocalPoints[sl] = subPoint; 3765 sremotePoints[sl].rank = newLocalPoints[point].rank; 3766 sremotePoints[sl].index = newLocalPoints[point].index; 3767 PetscCheck(sremotePoints[sl].rank >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %" PetscInt_FMT, point); 3768 PetscCheck(sremotePoints[sl].index >= 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %" PetscInt_FMT, point); 3769 ++sl; 3770 } 3771 PetscCall(PetscFree2(newLocalPoints,newOwners)); 3772 PetscCheck(sl + ll == numSubLeaves,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %" PetscInt_FMT " + %" PetscInt_FMT " != %" PetscInt_FMT, sl, ll, numSubLeaves); 3773 PetscCall(PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3774 } 3775 } 3776 CHKMEMQ; 3777 /* Cleanup */ 3778 if (subvertexIS) PetscCall(ISRestoreIndices(subvertexIS, &subVertices)); 3779 PetscCall(ISDestroy(&subvertexIS)); 3780 PetscCall(PetscFree(subCells)); 3781 PetscFunctionReturn(0); 3782 } 3783 3784 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3785 { 3786 DMLabel label = NULL; 3787 3788 PetscFunctionBegin; 3789 if (labelname) PetscCall(DMGetLabel(dm, labelname, &label)); 3790 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm)); 3791 PetscFunctionReturn(0); 3792 } 3793 3794 /*@C 3795 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. 3796 3797 Input Parameters: 3798 + dm - The original mesh 3799 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3800 . label - A label name, or NULL 3801 - value - A label value 3802 3803 Output Parameter: 3804 . subdm - The surface mesh 3805 3806 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3807 3808 Level: developer 3809 3810 .seealso: `DMPlexGetSubpointMap()`, `DMPlexCreateSubmesh()` 3811 @*/ 3812 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3813 { 3814 PetscInt dim, cdim, depth; 3815 3816 PetscFunctionBegin; 3817 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3818 PetscValidPointer(subdm, 5); 3819 PetscCall(DMGetDimension(dm, &dim)); 3820 PetscCall(DMPlexGetDepth(dm, &depth)); 3821 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3822 PetscCall(DMSetType(*subdm, DMPLEX)); 3823 PetscCall(DMSetDimension(*subdm, dim-1)); 3824 PetscCall(DMGetCoordinateDim(dm, &cdim)); 3825 PetscCall(DMSetCoordinateDim(*subdm, cdim)); 3826 if (depth == dim) { 3827 PetscCall(DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm)); 3828 } else { 3829 PetscCall(DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm)); 3830 } 3831 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3832 PetscFunctionReturn(0); 3833 } 3834 3835 /*@ 3836 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3837 3838 Input Parameters: 3839 + dm - The original mesh 3840 . cellLabel - The DMLabel marking cells contained in the new mesh 3841 - value - The label value to use 3842 3843 Output Parameter: 3844 . subdm - The new mesh 3845 3846 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3847 3848 Level: developer 3849 3850 .seealso: `DMPlexGetSubpointMap()`, `DMGetLabel()`, `DMLabelSetValue()` 3851 @*/ 3852 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3853 { 3854 PetscInt dim; 3855 3856 PetscFunctionBegin; 3857 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3858 PetscValidPointer(subdm, 4); 3859 PetscCall(DMGetDimension(dm, &dim)); 3860 PetscCall(DMCreate(PetscObjectComm((PetscObject) dm), subdm)); 3861 PetscCall(DMSetType(*subdm, DMPLEX)); 3862 PetscCall(DMSetDimension(*subdm, dim)); 3863 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3864 PetscCall(DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm)); 3865 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *subdm)); 3866 PetscFunctionReturn(0); 3867 } 3868 3869 /*@ 3870 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3871 3872 Input Parameter: 3873 . dm - The submesh DM 3874 3875 Output Parameter: 3876 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3877 3878 Level: developer 3879 3880 .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 3881 @*/ 3882 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3883 { 3884 PetscFunctionBegin; 3885 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3886 PetscValidPointer(subpointMap, 2); 3887 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3888 PetscFunctionReturn(0); 3889 } 3890 3891 /*@ 3892 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3893 3894 Input Parameters: 3895 + dm - The submesh DM 3896 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3897 3898 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3899 3900 Level: developer 3901 3902 .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointIS()` 3903 @*/ 3904 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3905 { 3906 DM_Plex *mesh = (DM_Plex *) dm->data; 3907 DMLabel tmp; 3908 3909 PetscFunctionBegin; 3910 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3911 tmp = mesh->subpointMap; 3912 mesh->subpointMap = subpointMap; 3913 PetscCall(PetscObjectReference((PetscObject) mesh->subpointMap)); 3914 PetscCall(DMLabelDestroy(&tmp)); 3915 PetscFunctionReturn(0); 3916 } 3917 3918 static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS) 3919 { 3920 DMLabel spmap; 3921 PetscInt depth, d; 3922 3923 PetscFunctionBegin; 3924 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 3925 PetscCall(DMPlexGetDepth(dm, &depth)); 3926 if (spmap && depth >= 0) { 3927 DM_Plex *mesh = (DM_Plex *) dm->data; 3928 PetscInt *points, *depths; 3929 PetscInt pStart, pEnd, p, off; 3930 3931 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 3932 PetscCheck(!pStart,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %" PetscInt_FMT, pStart); 3933 PetscCall(PetscMalloc1(pEnd, &points)); 3934 PetscCall(DMGetWorkArray(dm, depth+1, MPIU_INT, &depths)); 3935 depths[0] = depth; 3936 depths[1] = 0; 3937 for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3938 for (d = 0, off = 0; d <= depth; ++d) { 3939 const PetscInt dep = depths[d]; 3940 PetscInt depStart, depEnd, n; 3941 3942 PetscCall(DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd)); 3943 PetscCall(DMLabelGetStratumSize(spmap, dep, &n)); 3944 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3945 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); 3946 } else { 3947 if (!n) { 3948 if (d == 0) { 3949 /* Missing cells */ 3950 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3951 } else { 3952 /* Missing faces */ 3953 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3954 } 3955 } 3956 } 3957 if (n) { 3958 IS is; 3959 const PetscInt *opoints; 3960 3961 PetscCall(DMLabelGetStratumIS(spmap, dep, &is)); 3962 PetscCall(ISGetIndices(is, &opoints)); 3963 for (p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3964 PetscCall(ISRestoreIndices(is, &opoints)); 3965 PetscCall(ISDestroy(&is)); 3966 } 3967 } 3968 PetscCall(DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths)); 3969 PetscCheck(off == pEnd,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %" PetscInt_FMT " should be %" PetscInt_FMT, off, pEnd); 3970 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS)); 3971 PetscCall(PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState)); 3972 } 3973 PetscFunctionReturn(0); 3974 } 3975 3976 /*@ 3977 DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data 3978 3979 Input Parameter: 3980 . dm - The submesh DM 3981 3982 Output Parameter: 3983 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3984 3985 Note: This IS is guaranteed to be sorted by the construction of the submesh 3986 3987 Level: developer 3988 3989 .seealso: `DMPlexCreateSubmesh()`, `DMPlexGetSubpointMap()` 3990 @*/ 3991 PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS) 3992 { 3993 DM_Plex *mesh = (DM_Plex *) dm->data; 3994 DMLabel spmap; 3995 PetscObjectState state; 3996 3997 PetscFunctionBegin; 3998 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3999 PetscValidPointer(subpointIS, 2); 4000 PetscCall(DMPlexGetSubpointMap(dm, &spmap)); 4001 PetscCall(PetscObjectStateGet((PetscObject) spmap, &state)); 4002 if (state != mesh->subpointState || !mesh->subpointIS) PetscCall(DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS)); 4003 *subpointIS = mesh->subpointIS; 4004 PetscFunctionReturn(0); 4005 } 4006 4007 /*@ 4008 DMGetEnclosureRelation - Get the relationship between dmA and dmB 4009 4010 Input Parameters: 4011 + dmA - The first DM 4012 - dmB - The second DM 4013 4014 Output Parameter: 4015 . rel - The relation of dmA to dmB 4016 4017 Level: intermediate 4018 4019 .seealso: `DMGetEnclosurePoint()` 4020 @*/ 4021 PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel) 4022 { 4023 DM plexA, plexB, sdm; 4024 DMLabel spmap; 4025 PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB; 4026 4027 PetscFunctionBegin; 4028 PetscValidPointer(rel, 3); 4029 *rel = DM_ENC_NONE; 4030 if (!dmA || !dmB) PetscFunctionReturn(0); 4031 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 4032 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 4033 if (dmA == dmB) {*rel = DM_ENC_EQUALITY; PetscFunctionReturn(0);} 4034 PetscCall(DMConvert(dmA, DMPLEX, &plexA)); 4035 PetscCall(DMConvert(dmB, DMPLEX, &plexB)); 4036 PetscCall(DMPlexGetChart(plexA, &pStartA, &pEndA)); 4037 PetscCall(DMPlexGetChart(plexB, &pStartB, &pEndB)); 4038 /* Assumption 1: subDMs have smaller charts than the DMs that they originate from 4039 - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */ 4040 if ((pStartA == pStartB) && (pEndA == pEndB)) { 4041 *rel = DM_ENC_EQUALITY; 4042 goto end; 4043 } 4044 NpA = pEndA - pStartA; 4045 NpB = pEndB - pStartB; 4046 if (NpA == NpB) goto end; 4047 sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */ 4048 PetscCall(DMPlexGetSubpointMap(sdm, &spmap)); 4049 if (!spmap) goto end; 4050 /* TODO Check the space mapped to by subpointMap is same size as dm */ 4051 if (NpA > NpB) { 4052 *rel = DM_ENC_SUPERMESH; 4053 } else { 4054 *rel = DM_ENC_SUBMESH; 4055 } 4056 end: 4057 PetscCall(DMDestroy(&plexA)); 4058 PetscCall(DMDestroy(&plexB)); 4059 PetscFunctionReturn(0); 4060 } 4061 4062 /*@ 4063 DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB 4064 4065 Input Parameters: 4066 + dmA - The first DM 4067 . dmB - The second DM 4068 . etype - The type of enclosure relation that dmA has to dmB 4069 - pB - A point of dmB 4070 4071 Output Parameter: 4072 . pA - The corresponding point of dmA 4073 4074 Level: intermediate 4075 4076 .seealso: `DMGetEnclosureRelation()` 4077 @*/ 4078 PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA) 4079 { 4080 DM sdm; 4081 IS subpointIS; 4082 const PetscInt *subpoints; 4083 PetscInt numSubpoints; 4084 4085 PetscFunctionBegin; 4086 /* TODO Cache the IS, making it look like an index */ 4087 switch (etype) { 4088 case DM_ENC_SUPERMESH: 4089 sdm = dmB; 4090 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4091 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4092 *pA = subpoints[pB]; 4093 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4094 break; 4095 case DM_ENC_SUBMESH: 4096 sdm = dmA; 4097 PetscCall(DMPlexGetSubpointIS(sdm, &subpointIS)); 4098 PetscCall(ISGetLocalSize(subpointIS, &numSubpoints)); 4099 PetscCall(ISGetIndices(subpointIS, &subpoints)); 4100 PetscCall(PetscFindInt(pB, numSubpoints, subpoints, pA)); 4101 if (*pA < 0) { 4102 PetscCall(DMViewFromOptions(dmA, NULL, "-dm_enc_A_view")); 4103 PetscCall(DMViewFromOptions(dmB, NULL, "-dm_enc_B_view")); 4104 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not found in submesh", pB); 4105 } 4106 PetscCall(ISRestoreIndices(subpointIS, &subpoints)); 4107 break; 4108 case DM_ENC_EQUALITY: 4109 case DM_ENC_NONE: 4110 *pA = pB;break; 4111 case DM_ENC_UNKNOWN: 4112 { 4113 DMEnclosureType enc; 4114 4115 PetscCall(DMGetEnclosureRelation(dmA, dmB, &enc)); 4116 PetscCall(DMGetEnclosurePoint(dmA, dmB, enc, pB, pA)); 4117 } 4118 break; 4119 default: SETERRQ(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype); 4120 } 4121 PetscFunctionReturn(0); 4122 } 4123