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