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