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 CHKERRQ(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 CHKERRQ(DMPlexGetCellTypeLabel(dm, &ctLabel)); 30 switch (dim) { 31 case 1: CHKERRQ(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_POINT_PRISM_TENSOR, cStart, cEnd));break; 32 case 2: CHKERRQ(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_SEG_PRISM_TENSOR, cStart, cEnd));break; 33 case 3: 34 CHKERRQ(DMLabelGetStratumBounds(ctLabel, DM_POLYTOPE_TRI_PRISM_TENSOR, cStart, cEnd)); 35 if (*cStart < 0) CHKERRQ(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 CHKERRQ(DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd)); 48 for (f = fStart; f < fEnd; ++f) { 49 PetscInt supportSize; 50 51 CHKERRQ(DMPlexGetSupportSize(dm, f, &supportSize)); 52 if (supportSize == 1) { 53 if (val < 0) { 54 PetscInt *closure = NULL; 55 PetscInt clSize, cl, cval; 56 57 CHKERRQ(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure)); 58 for (cl = 0; cl < clSize*2; cl += 2) { 59 CHKERRQ(DMLabelGetValue(label, closure[cl], &cval)); 60 if (cval < 0) continue; 61 CHKERRQ(DMLabelSetValue(label, f, cval)); 62 break; 63 } 64 if (cl == clSize*2) CHKERRQ(DMLabelSetValue(label, f, 1)); 65 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &clSize, &closure)); 66 } else { 67 CHKERRQ(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 CHKERRQ(DMPlexIsInterpolated(dm, &flg)); 97 PetscCheckFalse(flg != DMPLEX_INTERPOLATED_FULL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not fully interpolated on this rank"); 98 CHKERRQ(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 CHKERRQ(DMLabelGetNumValues(label, &numValues)); 111 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 112 CHKERRQ(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 113 CHKERRQ(ISGetIndices(valueIS, &values)); 114 for (v = 0; v < numValues; ++v) { 115 IS pointIS; 116 const PetscInt *points; 117 PetscInt numPoints, p; 118 119 CHKERRQ(DMLabelGetStratumSize(label, values[v], &numPoints)); 120 CHKERRQ(DMLabelGetStratumIS(label, values[v], &pointIS)); 121 CHKERRQ(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 CHKERRQ(DMPlexGetTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure)); 131 for (c = 0; c < closureSize*2; c += 2) { 132 CHKERRQ(DMLabelSetValue(label, closure[c], values[v])); 133 } 134 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, q, PETSC_TRUE, &closureSize, &closure)); 135 } 136 CHKERRQ(ISRestoreIndices(pointIS, &points)); 137 CHKERRQ(ISDestroy(&pointIS)); 138 } 139 CHKERRQ(ISRestoreIndices(valueIS, &values)); 140 CHKERRQ(ISDestroy(&valueIS)); 141 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 142 CHKERRQ(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 CHKERRQ(DMLabelGather(label, sfPoint, &lblLeaves)); 151 CHKERRQ(DMLabelGetValueIS(lblLeaves, &valueIS)); 152 CHKERRQ(ISGetLocalSize(valueIS, &numValues)); 153 CHKERRQ(ISGetIndices(valueIS, &values)); 154 for (v = 0; v < numValues; ++v) { 155 const PetscInt value = values[v]; 156 157 CHKERRQ(DMLabelGetStratumIS(lblLeaves, value, &pointIS)); 158 CHKERRQ(DMLabelInsertIS(label, pointIS, value)); 159 CHKERRQ(ISDestroy(&pointIS)); 160 } 161 CHKERRQ(ISRestoreIndices(valueIS, &values)); 162 CHKERRQ(ISDestroy(&valueIS)); 163 CHKERRQ(DMLabelDestroy(&lblLeaves)); 164 /* Push point contributions from roots into remote leaves */ 165 CHKERRQ(DMLabelDistribute(label, sfPoint, &lblRoots)); 166 CHKERRQ(DMLabelGetValueIS(lblRoots, &valueIS)); 167 CHKERRQ(ISGetLocalSize(valueIS, &numValues)); 168 CHKERRQ(ISGetIndices(valueIS, &values)); 169 for (v = 0; v < numValues; ++v) { 170 const PetscInt value = values[v]; 171 172 CHKERRQ(DMLabelGetStratumIS(lblRoots, value, &pointIS)); 173 CHKERRQ(DMLabelInsertIS(label, pointIS, value)); 174 CHKERRQ(ISDestroy(&pointIS)); 175 } 176 CHKERRQ(ISRestoreIndices(valueIS, &values)); 177 CHKERRQ(ISDestroy(&valueIS)); 178 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 228 CHKERRQ(DMLabelGetNumValues(label, &numValues)); 229 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 230 CHKERRQ(ISGetIndices(valueIS, &values)); 231 for (v = 0; v < numValues; ++v) { 232 IS pointIS; 233 const PetscInt *points; 234 PetscInt numPoints, p; 235 236 CHKERRQ(DMLabelGetStratumSize(label, values[v], &numPoints)); 237 CHKERRQ(DMLabelGetStratumIS(label, values[v], &pointIS)); 238 CHKERRQ(ISGetIndices(pointIS, &points)); 239 for (p = 0; p < numPoints; ++p) { 240 PetscInt *closure = NULL; 241 PetscInt closureSize, cl; 242 243 CHKERRQ(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)) {CHKERRQ(DMLabelSetValue(label, cell, values[v])); break;} 247 } 248 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closure)); 249 } 250 CHKERRQ(ISRestoreIndices(pointIS, &points)); 251 CHKERRQ(ISDestroy(&pointIS)); 252 } 253 CHKERRQ(ISRestoreIndices(valueIS, &values)); 254 CHKERRQ(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 CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 282 CHKERRQ(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 283 CHKERRQ(DMLabelGetNumValues(label, &numValues)); 284 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 285 CHKERRQ(ISGetIndices(valueIS, &values)); 286 for (v = 0; v < numValues; ++v) { 287 IS pointIS; 288 const PetscInt *points; 289 PetscInt numPoints, p; 290 291 CHKERRQ(DMLabelGetStratumSize(label, values[v], &numPoints)); 292 CHKERRQ(DMLabelGetStratumIS(label, values[v], &pointIS)); 293 CHKERRQ(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 CHKERRQ(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)) {CHKERRQ(DMLabelSetValue(label, cell, values[v])); break;} 304 } 305 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, face, PETSC_FALSE, &closureSize, &closure)); 306 } 307 CHKERRQ(ISRestoreIndices(pointIS, &points)); 308 CHKERRQ(ISDestroy(&pointIS)); 309 } 310 CHKERRQ(ISRestoreIndices(valueIS, &values)); 311 CHKERRQ(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 CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 339 CHKERRQ(DMLabelGetNumValues(label, &numValues)); 340 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 341 CHKERRQ(ISGetIndices(valueIS, &values)); 342 for (v = 0; v < numValues; ++v) { 343 IS pointIS; 344 const PetscInt *points; 345 PetscInt numPoints, p; 346 347 CHKERRQ(DMLabelGetStratumSize(label, values[v], &numPoints)); 348 CHKERRQ(DMLabelGetStratumIS(label, values[v], &pointIS)); 349 CHKERRQ(ISGetIndices(pointIS, &points)); 350 for (p = 0; p < numPoints; ++p) { 351 PetscInt point = points[p]; 352 353 if (point >= cStart && point < cEnd) { 354 CHKERRQ(DMLabelClearValue(label,point,values[v])); 355 } 356 } 357 CHKERRQ(ISRestoreIndices(pointIS, &points)); 358 CHKERRQ(ISDestroy(&pointIS)); 359 } 360 CHKERRQ(ISRestoreIndices(valueIS, &values)); 361 CHKERRQ(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 437 if (depth < 0) PetscFunctionReturn(0); 438 /* Step 1: Expand chart */ 439 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 440 pEnd = DMPlexShiftPoint_Internal(pEnd,depth,depthShift); 441 CHKERRQ(DMPlexSetChart(dmNew, pStart, pEnd)); 442 CHKERRQ(DMCreateLabel(dmNew,"depth")); 443 CHKERRQ(DMPlexGetDepthLabel(dmNew,&depthLabel)); 444 CHKERRQ(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 CHKERRQ(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 451 pStartNew = DMPlexShiftPoint_Internal(pStart, depth, depthShift); 452 pEndNew = DMPlexShiftPoint_Internal(pEnd, depth, depthShift); 453 CHKERRQ(ISCreateStride(PETSC_COMM_SELF, pEndNew - pStartNew, pStartNew, 1, &pIS)); 454 CHKERRQ(DMLabelSetStratumIS(depthLabel, d, pIS)); 455 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, p, &size)); 462 CHKERRQ(DMPlexSetConeSize(dmNew, newp, size)); 463 CHKERRQ(DMPlexGetSupportSize(dm, p, &size)); 464 CHKERRQ(DMPlexSetSupportSize(dmNew, newp, size)); 465 CHKERRQ(DMPlexGetCellType(dm, p, &ct)); 466 CHKERRQ(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 479 if (depth < 0) PetscFunctionReturn(0); 480 CHKERRQ(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 481 CHKERRQ(DMPlexGetMaxSizes(dmNew, &maxConeSizeNew, &maxSupportSizeNew)); 482 CHKERRQ(PetscMalloc1(PetscMax(PetscMax(maxConeSize, maxSupportSize), PetscMax(maxConeSizeNew, maxSupportSizeNew)),&newpoints)); 483 /* Step 5: Set cones and supports */ 484 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, p, &size)); 490 CHKERRQ(DMPlexGetCone(dm, p, &points)); 491 CHKERRQ(DMPlexGetConeOrientation(dm, p, &orientations)); 492 for (i = 0; i < size; ++i) { 493 newpoints[i] = DMPlexShiftPoint_Internal(points[i], depth, depthShift); 494 } 495 CHKERRQ(DMPlexSetCone(dmNew, newp, newpoints)); 496 CHKERRQ(DMPlexSetConeOrientation(dmNew, newp, orientations)); 497 CHKERRQ(DMPlexGetSupportSize(dm, p, &size)); 498 CHKERRQ(DMPlexGetSupportSize(dmNew, newp, &sizeNew)); 499 CHKERRQ(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 CHKERRQ(DMPlexSetSupport(dmNew, newp, newpoints)); 505 } 506 CHKERRQ(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 CHKERRQ(DMGetCoordinateDim(dm, &dim)); 521 CHKERRQ(DMSetCoordinateDim(dmNew, dim)); 522 CHKERRQ(DMPlexGetDepth(dm, &depth)); 523 /* Step 8: Convert coordinates */ 524 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 525 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 526 CHKERRQ(DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew)); 527 CHKERRQ(DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew)); 528 CHKERRQ(DMGetCoordinateSection(dm, &coordSection)); 529 CHKERRQ(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection)); 530 CHKERRQ(PetscSectionSetNumFields(newCoordSection, 1)); 531 CHKERRQ(PetscSectionSetFieldComponents(newCoordSection, 0, dim)); 532 CHKERRQ(PetscSectionGetChart(coordSection, &sStart, &sEnd)); 533 hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE; 534 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSection, c, &dof)); 540 CHKERRQ(PetscSectionSetDof(newCoordSection, cNew, dof)); 541 CHKERRQ(PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof)); 542 } 543 } 544 for (v = vStartNew; v < vEndNew; ++v) { 545 CHKERRQ(PetscSectionSetDof(newCoordSection, v, dim)); 546 CHKERRQ(PetscSectionSetFieldDof(newCoordSection, v, 0, dim)); 547 } 548 CHKERRQ(PetscSectionSetUp(newCoordSection)); 549 CHKERRQ(DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection)); 550 CHKERRQ(PetscSectionGetStorageSize(newCoordSection, &coordSize)); 551 CHKERRQ(VecCreate(PETSC_COMM_SELF, &newCoordinates)); 552 CHKERRQ(PetscObjectSetName((PetscObject) newCoordinates, "coordinates")); 553 CHKERRQ(VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE)); 554 CHKERRQ(VecSetBlockSize(newCoordinates, dim)); 555 CHKERRQ(VecSetType(newCoordinates,VECSTANDARD)); 556 CHKERRQ(DMSetCoordinatesLocal(dmNew, newCoordinates)); 557 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates)); 558 CHKERRQ(VecGetArray(coordinates, &coords)); 559 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSection, c, &dof)); 565 CHKERRQ(PetscSectionGetOffset(coordSection, c, &off)); 566 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSection, v, &dof)); 574 CHKERRQ(PetscSectionGetOffset(coordSection, v, &off)); 575 CHKERRQ(PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff)); 576 for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d]; 577 } 578 CHKERRQ(VecRestoreArray(coordinates, &coords)); 579 CHKERRQ(VecRestoreArray(newCoordinates, &newCoords)); 580 CHKERRQ(VecDestroy(&newCoordinates)); 581 CHKERRQ(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 595 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 596 CHKERRQ(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 597 totShift = DMPlexShiftPoint_Internal(pEnd, depth, depthShift) - pEnd; 598 if (numRoots >= 0) { 599 CHKERRQ(PetscMalloc2(numRoots, &newLocation, pEnd-pStart, &newRemoteLocation)); 600 for (l = 0; l < numRoots; ++l) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift); 601 CHKERRQ(PetscSFBcastBegin(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE)); 602 CHKERRQ(PetscSFBcastEnd(sf, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE)); 603 CHKERRQ(PetscMalloc1(numLeaves, &glocalPoints)); 604 CHKERRQ(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 CHKERRQ(PetscFree2(newLocation, newRemoteLocation)); 611 CHKERRQ(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 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 624 CHKERRQ(DMGetPointSF(dmNew, &sfPointNew)); 625 CHKERRQ(DMPlexShiftSF_Single(dm, depthShift, sfPoint, sfPointNew)); 626 /* Step 9b: Convert naturalSF */ 627 CHKERRQ(DMGetUseNatural(dm, &useNatural)); 628 if (useNatural) { 629 PetscSF sfNat, sfNatNew; 630 631 CHKERRQ(DMSetUseNatural(dmNew, useNatural)); 632 CHKERRQ(DMGetNaturalSF(dm, &sfNat)); 633 CHKERRQ(DMGetNaturalSF(dmNew, &sfNatNew)); 634 CHKERRQ(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 645 /* Step 10: Convert labels */ 646 CHKERRQ(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 CHKERRQ(DMGetLabelName(dm, l, &lname)); 656 CHKERRQ(PetscStrcmp(lname, "depth", &isDepth)); 657 if (isDepth) continue; 658 CHKERRQ(PetscStrcmp(lname, "dim", &isDim)); 659 if (isDim) continue; 660 CHKERRQ(DMCreateLabel(dmNew, lname)); 661 CHKERRQ(DMGetLabel(dm, lname, &label)); 662 CHKERRQ(DMGetLabel(dmNew, lname, &newlabel)); 663 CHKERRQ(DMLabelGetDefaultValue(label,&val)); 664 CHKERRQ(DMLabelSetDefaultValue(newlabel,val)); 665 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 666 CHKERRQ(ISGetLocalSize(valueIS, &numValues)); 667 CHKERRQ(ISGetIndices(valueIS, &values)); 668 for (val = 0; val < numValues; ++val) { 669 IS pointIS; 670 const PetscInt *points; 671 PetscInt numPoints, p; 672 673 CHKERRQ(DMLabelGetStratumIS(label, values[val], &pointIS)); 674 CHKERRQ(ISGetLocalSize(pointIS, &numPoints)); 675 CHKERRQ(ISGetIndices(pointIS, &points)); 676 for (p = 0; p < numPoints; ++p) { 677 const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift); 678 679 CHKERRQ(DMLabelSetValue(newlabel, newpoint, values[val])); 680 } 681 CHKERRQ(ISRestoreIndices(pointIS, &points)); 682 CHKERRQ(ISDestroy(&pointIS)); 683 } 684 CHKERRQ(ISRestoreIndices(valueIS, &values)); 685 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 702 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 703 CHKERRQ(DMPlexGetVTKCellHeight(dmNew, &cellHeight)); 704 CHKERRQ(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 705 CHKERRQ(PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote)); 706 CHKERRQ(DMCreateLabel(dmNew, "vtk")); 707 CHKERRQ(DMGetLabel(dmNew, "vtk", &vtkLabel)); 708 if (createGhostLabel) { 709 CHKERRQ(DMCreateLabel(dmNew, "ghost")); 710 CHKERRQ(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 CHKERRQ(DMLabelSetValue(vtkLabel, c, 1)); 715 } 716 if (leafLocal[l] >= cEnd) break; 717 if (leafRemote[l].rank == rank) { 718 CHKERRQ(DMLabelSetValue(vtkLabel, c, 1)); 719 } else if (ghostLabel) { 720 CHKERRQ(DMLabelSetValue(ghostLabel, c, 2)); 721 } 722 } 723 for (; c < cEnd; ++c) { 724 CHKERRQ(DMLabelSetValue(vtkLabel, c, 1)); 725 } 726 if (ghostLabel) { 727 CHKERRQ(DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd)); 728 for (f = fStart; f < fEnd; ++f) { 729 PetscInt numCells; 730 731 CHKERRQ(DMPlexGetSupportSize(dmNew, f, &numCells)); 732 if (numCells < 2) { 733 CHKERRQ(DMLabelSetValue(ghostLabel, f, 1)); 734 } else { 735 const PetscInt *cells = NULL; 736 PetscInt vA, vB; 737 738 CHKERRQ(DMPlexGetSupport(dmNew, f, &cells)); 739 CHKERRQ(DMLabelGetValue(vtkLabel, cells[0], &vA)); 740 CHKERRQ(DMLabelGetValue(vtkLabel, cells[1], &vB)); 741 if (vA != 1 && vB != 1) CHKERRQ(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 CHKERRQ(DMPlexGetReferenceTree(dm,&refTree)); 756 CHKERRQ(DMPlexSetReferenceTree(dmNew,refTree)); 757 CHKERRQ(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 CHKERRQ(PetscSectionGetChart(pSec,&pStart,&pEnd)); 764 CHKERRQ(DMPlexGetDepth(dm,&depth)); 765 pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift); 766 pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift); 767 CHKERRQ(PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted)); 768 CHKERRQ(PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted)); 769 CHKERRQ(PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted)); 770 for (p = pStartShifted; p < pEndShifted; p++) { 771 /* start off assuming no children */ 772 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(pSec,p,&dof)); 779 CHKERRQ(PetscSectionSetDof(pSecShifted,pNew,dof)); 780 } 781 CHKERRQ(PetscSectionSetUp(pSecShifted)); 782 for (p = pStart; p < pEnd; p++) { 783 PetscInt dof; 784 PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift); 785 786 CHKERRQ(PetscSectionGetDof(pSec,p,&dof)); 787 if (dof) { 788 PetscInt off, offNew; 789 790 CHKERRQ(PetscSectionGetOffset(pSec,p,&off)); 791 CHKERRQ(PetscSectionGetOffset(pSecShifted,pNew,&offNew)); 792 parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift); 793 childIDsShifted[offNew] = childIDs[off]; 794 } 795 } 796 CHKERRQ(DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted)); 797 CHKERRQ(PetscFree2(parentsShifted,childIDsShifted)); 798 CHKERRQ(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 CHKERRQ(DMGetPointSF(dm, &sf)); 816 CHKERRQ(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL)); 817 nleaves = PetscMax(0, nleaves); 818 CHKERRQ(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 819 /* Count ghost cells */ 820 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 821 CHKERRQ(ISGetLocalSize(valueIS, &numFS)); 822 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, values[fs], &faceIS)); 830 CHKERRQ(ISGetLocalSize(faceIS, &numFaces)); 831 CHKERRQ(ISGetIndices(faceIS, &faces)); 832 for (f = 0; f < numFaces; ++f) { 833 PetscInt numChildren; 834 835 CHKERRQ(PetscFindInt(faces[f], nleaves, leaves, &loc)); 836 CHKERRQ(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 CHKERRQ(ISRestoreIndices(faceIS, &faces)); 843 CHKERRQ(ISDestroy(&faceIS)); 844 } 845 CHKERRQ(DMPlexGetDepth(dm, &depth)); 846 CHKERRQ(PetscMalloc1(2*(depth+1), &depthShift)); 847 for (d = 0; d <= depth; d++) { 848 PetscInt dEnd; 849 850 CHKERRQ(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 CHKERRQ(DMPlexShiftPointSetUp_Internal(depth,depthShift)); 856 CHKERRQ(DMPlexShiftSizes_Internal(dm, depthShift, gdm)); 857 /* Step 3: Set cone/support sizes for new points */ 858 CHKERRQ(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 859 for (c = cEnd; c < cEnd + Ng; ++c) { 860 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, values[fs], &faceIS)); 868 CHKERRQ(ISGetLocalSize(faceIS, &numFaces)); 869 CHKERRQ(ISGetIndices(faceIS, &faces)); 870 for (f = 0; f < numFaces; ++f) { 871 PetscInt size, numChildren; 872 873 CHKERRQ(PetscFindInt(faces[f], nleaves, leaves, &loc)); 874 CHKERRQ(DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL)); 875 if (loc >= 0 || numChildren) continue; 876 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 877 CHKERRQ(DMPlexGetSupportSize(dm, faces[f], &size)); 878 PetscCheckFalse(size != 1,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size); 879 CHKERRQ(DMPlexSetSupportSize(gdm, faces[f] + Ng, 2)); 880 } 881 CHKERRQ(ISRestoreIndices(faceIS, &faces)); 882 CHKERRQ(ISDestroy(&faceIS)); 883 } 884 /* Step 4: Setup ghosted DM */ 885 CHKERRQ(DMSetUp(gdm)); 886 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, values[fs], &faceIS)); 895 CHKERRQ(ISGetLocalSize(faceIS, &numFaces)); 896 CHKERRQ(ISGetIndices(faceIS, &faces)); 897 for (f = 0; f < numFaces; ++f) { 898 PetscInt newFace = faces[f] + Ng, numChildren; 899 900 CHKERRQ(PetscFindInt(faces[f], nleaves, leaves, &loc)); 901 CHKERRQ(DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL)); 902 if (loc >= 0 || numChildren) continue; 903 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 904 CHKERRQ(DMPlexSetCone(gdm, ghostCell, &newFace)); 905 CHKERRQ(DMPlexInsertSupport(gdm, newFace, 1, ghostCell)); 906 ++ghostCell; 907 } 908 CHKERRQ(ISRestoreIndices(faceIS, &faces)); 909 CHKERRQ(ISDestroy(&faceIS)); 910 } 911 CHKERRQ(ISRestoreIndices(valueIS, &values)); 912 CHKERRQ(ISDestroy(&valueIS)); 913 CHKERRQ(DMPlexShiftCoordinates_Internal(dm, depthShift, gdm)); 914 CHKERRQ(DMPlexShiftSF_Internal(dm, depthShift, gdm)); 915 CHKERRQ(DMPlexShiftLabels_Internal(dm, depthShift, gdm)); 916 CHKERRQ(DMPlexCreateVTKLabel_Internal(dm, PETSC_TRUE, gdm)); 917 CHKERRQ(DMPlexShiftTree_Internal(dm, depthShift, gdm)); 918 CHKERRQ(PetscFree(depthShift)); 919 for (c = cEnd; c < cEnd + Ng; ++c) { 920 CHKERRQ(DMPlexSetCellType(gdm, c, DM_POLYTOPE_FV_GHOST)); 921 } 922 /* Step 7: Periodicity */ 923 CHKERRQ(DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd)); 924 CHKERRQ(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 CHKERRQ(DMCreate(PetscObjectComm((PetscObject)dm), &gdm)); 961 CHKERRQ(DMSetType(gdm, DMPLEX)); 962 CHKERRQ(DMGetDimension(dm, &dim)); 963 CHKERRQ(DMSetDimension(gdm, dim)); 964 CHKERRQ(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 965 CHKERRQ(DMSetBasicAdjacency(gdm, useCone, useClosure)); 966 CHKERRQ(DMGetLabel(dm, name, &label)); 967 if (!label) { 968 /* Get label for boundary faces */ 969 CHKERRQ(DMCreateLabel(dm, name)); 970 CHKERRQ(DMGetLabel(dm, name, &label)); 971 CHKERRQ(DMPlexMarkBoundaryFaces(dm, 1, label)); 972 } 973 CHKERRQ(DMPlexConstructGhostCells_Internal(dm, label, &Ng, gdm)); 974 CHKERRQ(DMCopyDisc(dm, gdm)); 975 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 1017 CHKERRQ(DMGetDimension(dm, &dim)); 1018 CHKERRQ(DMPlexGetDepth(dm, &depth)); 1019 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1020 /* We do not want this label automatically computed, instead we compute it here */ 1021 CHKERRQ(DMCreateLabel(sdm, "celltype")); 1022 /* Count split points and add cohesive cells */ 1023 CHKERRQ(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 1024 CHKERRQ(PetscMalloc5(depth+1,&depthMax,depth+1,&depthEnd,2*(depth+1),&depthShift,depth+1,&pMaxNew,depth+1,&numHybridPointsOld)); 1025 CHKERRQ(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 CHKERRQ(DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d])); 1028 CHKERRQ(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 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 1046 CHKERRQ(ISGetLocalSize(valueIS, &numSP)); 1047 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, dep, &splitIS[dep])); 1054 if (splitIS[dep]) { 1055 CHKERRQ(ISGetLocalSize(splitIS[dep], &numSplitPoints[dep])); 1056 CHKERRQ(ISGetIndices(splitIS[dep], &splitPoints[dep])); 1057 } 1058 CHKERRQ(DMLabelGetStratumIS(label, shift2+dep, &unsplitIS[dep])); 1059 if (unsplitIS[dep]) { 1060 CHKERRQ(ISGetLocalSize(unsplitIS[dep], &numUnsplitPoints[dep])); 1061 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, oldp, &coneSize)); 1084 CHKERRQ(DMPlexSetConeSize(sdm, splitp, coneSize)); 1085 CHKERRQ(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1086 CHKERRQ(DMPlexSetSupportSize(sdm, splitp, supportSize)); 1087 CHKERRQ(DMPlexGetCellType(dm, oldp, &ct)); 1088 CHKERRQ(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 CHKERRQ(DMPlexSetConeSize(sdm, hybcell, 2 + coneSize)); 1094 switch (coneSize) { 1095 case 2: CHKERRQ(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_SEG_PRISM_TENSOR));break; 1096 case 3: CHKERRQ(DMPlexSetCellType(sdm, hybcell, DM_POLYTOPE_TRI_PRISM_TENSOR));break; 1097 case 4: CHKERRQ(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 CHKERRQ(DMPlexGetSupport(dm, oldp, &support)); 1103 for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) { 1104 PetscInt val; 1105 1106 CHKERRQ(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 CHKERRQ(DMPlexSetSupportSize(sdm, newp, qn+1)); 1113 /* Split new vertex: Edges into split vertex and new cohesive edge */ 1114 CHKERRQ(DMPlexSetSupportSize(sdm, splitp, qp+1)); 1115 /* Add hybrid edge */ 1116 CHKERRQ(DMPlexSetConeSize(sdm, hybedge, 2)); 1117 CHKERRQ(DMPlexSetSupportSize(sdm, hybedge, qf)); 1118 CHKERRQ(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 CHKERRQ(DMPlexGetSupport(dm, oldp, &support)); 1123 for (e = 0, qn = 0, qp = 0, qf = 0; e < supportSize; ++e) { 1124 PetscInt val; 1125 1126 CHKERRQ(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 CHKERRQ(DMPlexSetSupportSize(sdm, newp, qn+1)); 1133 /* Split new edge: Faces into split edge and cohesive face (negative side?) */ 1134 CHKERRQ(DMPlexSetSupportSize(sdm, splitp, qp+1)); 1135 /* Add hybrid face */ 1136 CHKERRQ(DMPlexSetConeSize(sdm, hybface, 4)); 1137 CHKERRQ(DMPlexSetSupportSize(sdm, hybface, qf)); 1138 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, oldp, &coneSize)); 1150 CHKERRQ(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1151 CHKERRQ(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 CHKERRQ(PetscFindInt(support[s], numSplitPoints[dep+1], splitPoints[dep+1], &e)); 1158 if (e >= 0) ++qf; 1159 } 1160 CHKERRQ(DMPlexSetSupportSize(sdm, newp, qf+2)); 1161 /* Add hybrid edge */ 1162 CHKERRQ(DMPlexSetConeSize(sdm, hybedge, 2)); 1163 for (e = 0, qf = 0; e < supportSize; ++e) { 1164 PetscInt val; 1165 1166 CHKERRQ(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 CHKERRQ(DMPlexSetSupportSize(sdm, hybedge, qf)); 1172 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMPlexSetSupportSize(sdm, newp, qf+2)); 1184 /* Add hybrid face */ 1185 for (e = 0, qf = 0; e < supportSize; ++e) { 1186 CHKERRQ(DMLabelGetValue(label, support[e], &val)); 1187 if (val == dim-1) ++qf; 1188 } 1189 CHKERRQ(DMPlexSetConeSize(sdm, hybface, 4)); 1190 CHKERRQ(DMPlexSetSupportSize(sdm, hybface, qf)); 1191 CHKERRQ(DMPlexSetCellType(sdm, hybface, DM_POLYTOPE_SEG_PRISM_TENSOR)); 1192 } 1193 } 1194 } 1195 /* Step 4: Setup split DM */ 1196 CHKERRQ(DMSetUp(sdm)); 1197 CHKERRQ(DMPlexShiftPoints_Internal(dm, depthShift, sdm)); 1198 CHKERRQ(DMPlexGetMaxSizes(sdm, &maxConeSizeNew, &maxSupportSizeNew)); 1199 CHKERRQ(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 CHKERRQ(DMPlexGetCellType(dm, oldp, &ct)); 1211 CHKERRQ(DMPlexGetConeSize(dm, oldp, &coneSize)); 1212 CHKERRQ(DMPlexGetCone(dm, oldp, &cone)); 1213 CHKERRQ(DMPlexGetConeOrientation(dm, oldp, &ornt)); 1214 CHKERRQ(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1215 CHKERRQ(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 CHKERRQ(DMPlexGetSupport(sdm, newp, &supportF)); 1223 CHKERRQ(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 CHKERRQ(PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v)); 1228 if (v < 0) { 1229 CHKERRQ(PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v)); 1230 PetscCheckFalse(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 CHKERRQ(DMPlexGetConeSize(dm, cone[q], &econeSize)); 1240 CHKERRQ(DMPlexGetCone(dm, cone[q], &econe)); 1241 for (r = 0; r < econeSize; ++r) { 1242 CHKERRQ(PetscFindInt(econe[r], numSplitPoints[dep-2], splitPoints[dep-2], &vs)); 1243 CHKERRQ(PetscFindInt(econe[r], numUnsplitPoints[dep-2], unsplitPoints[dep-2], &vu)); 1244 if (vs >= 0) continue; 1245 PetscCheckFalse(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 CHKERRQ(DMPlexSetCone(sdm, splitp, &coneNew[2])); 1252 CHKERRQ(DMPlexSetConeOrientation(sdm, splitp, ornt)); 1253 /* Face support */ 1254 for (s = 0; s < supportSize; ++s) { 1255 PetscInt val; 1256 1257 CHKERRQ(DMLabelGetValue(label, support[s], &val)); 1258 if (val < 0) { 1259 /* Split old face: Replace negative side cell with cohesive cell */ 1260 CHKERRQ(DMPlexInsertSupport(sdm, newp, s, hybcell)); 1261 } else { 1262 /* Split new face: Replace positive side cell with cohesive cell */ 1263 CHKERRQ(DMPlexInsertSupport(sdm, splitp, s, hybcell)); 1264 /* Get orientation for cohesive face */ 1265 { 1266 const PetscInt *ncone, *nconeO; 1267 PetscInt nconeSize, nc; 1268 1269 CHKERRQ(DMPlexGetConeSize(dm, support[s], &nconeSize)); 1270 CHKERRQ(DMPlexGetCone(dm, support[s], &ncone)); 1271 CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(PetscFindInt(cone[qa], numSplitPoints[dep-1], splitPoints[dep-1], &v)); 1295 if (v < 0) { 1296 CHKERRQ(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 CHKERRQ(DMPlexSetCone(sdm, hybcell, coneNew)); 1304 CHKERRQ(DMPlexSetConeOrientation(sdm, hybcell, coneONew)); 1305 /* Label the hybrid cells on the boundary of the split */ 1306 if (hasUnsplit) CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, support[e], &val)); 1326 if (val == 1) { 1327 CHKERRQ(PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge)); 1328 PetscCheckFalse(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 CHKERRQ(DMPlexSetSupport(sdm, splitp, supportNew)); 1336 /* Hybrid edge: Old and new split vertex */ 1337 coneNew[0] = newp; 1338 coneNew[1] = splitp; 1339 CHKERRQ(DMPlexSetCone(sdm, hybedge, coneNew)); 1340 for (e = 0, qf = 0; e < supportSize; ++e) { 1341 PetscInt val, edge; 1342 1343 CHKERRQ(DMLabelGetValue(label, support[e], &val)); 1344 if (val == 1) { 1345 CHKERRQ(PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge)); 1346 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(PetscFindInt(cone[q], numSplitPoints[dep-1], splitPoints[dep-1], &v)); 1358 if (v < 0) { 1359 CHKERRQ(PetscFindInt(cone[q], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v)); 1360 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, support[e], &val)); 1385 if (val == dim-1) { 1386 CHKERRQ(PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face)); 1387 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(PetscFindInt(cone[v], numSplitPoints[dep-1], splitPoints[dep-1], &vertex)); 1401 if (vertex < 0) { 1402 CHKERRQ(PetscFindInt(cone[v], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &vertex)); 1403 PetscCheckFalse(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 CHKERRQ(DMPlexSetCone(sdm, hybface, coneNew)); 1410 for (e = 0, qf = 0; e < supportSize; ++e) { 1411 PetscInt val, face; 1412 1413 CHKERRQ(DMLabelGetValue(label, support[e], &val)); 1414 if (val == dim-1) { 1415 CHKERRQ(PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &face)); 1416 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, oldp, &coneSize)); 1432 CHKERRQ(DMPlexGetCone(dm, oldp, &cone)); 1433 CHKERRQ(DMPlexGetSupportSize(dm, oldp, &supportSize)); 1434 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 PetscCheckFalse(q != supportSizeNew,comm, PETSC_ERR_ARG_WRONG, "Support size %d != %d for vertex %d", q, supportSizeNew, newp); 1450 CHKERRQ(DMPlexSetSupport(sdm, newp, supportNew)); 1451 /* Hybrid edge */ 1452 coneNew[0] = newp; 1453 coneNew[1] = newp; 1454 CHKERRQ(DMPlexSetCone(sdm, hybedge, coneNew)); 1455 for (e = 0, qf = 0; e < supportSize; ++e) { 1456 PetscInt val, edge; 1457 1458 CHKERRQ(DMLabelGetValue(label, support[e], &val)); 1459 if (val == 1) { 1460 CHKERRQ(PetscFindInt(support[e], numSplitPoints[dep+1], splitPoints[dep+1], &edge)); 1461 PetscCheckFalse(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 CHKERRQ(PetscFindInt(support[e], numUnsplitPoints[dep+1], unsplitPoints[dep+1], &edge)); 1465 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, support[f], &val)); 1478 if (val == dim-1) { 1479 CHKERRQ(PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face)); 1480 PetscCheckFalse(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 CHKERRQ(DMPlexGetSupportSize(sdm, newp, &supportSizeNew)); 1490 PetscCheckFalse(qf != supportSizeNew,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for unsplit edge %d is %d != %d", newp, qf, supportSizeNew); 1491 CHKERRQ(DMPlexSetSupport(sdm, newp, supportNew)); 1492 /* Add hybrid face */ 1493 coneNew[0] = newp; 1494 coneNew[1] = newp; 1495 CHKERRQ(PetscFindInt(cone[0], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v)); 1496 PetscCheckFalse(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 CHKERRQ(PetscFindInt(cone[1], numUnsplitPoints[dep-1], unsplitPoints[dep-1], &v)); 1499 PetscCheckFalse(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 CHKERRQ(DMPlexSetCone(sdm, hybface, coneNew)); 1502 for (f = 0, qf = 0; f < supportSize; ++f) { 1503 PetscInt val, face; 1504 1505 CHKERRQ(DMLabelGetValue(label, support[f], &val)); 1506 if (val == dim-1) { 1507 CHKERRQ(PetscFindInt(support[f], numSplitPoints[dep+1], splitPoints[dep+1], &face)); 1508 supportNew[qf++] = face + pMaxNew[dep+2] + numSplitPoints[dep+2]; 1509 } 1510 } 1511 CHKERRQ(DMPlexGetSupportSize(sdm, hybface, &supportSizeNew)); 1512 PetscCheckFalse(qf != supportSizeNew,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support size for hybrid face %d is %d != %d", hybface, qf, supportSizeNew); 1513 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, dep, &pIS)); 1526 if (!pIS) continue; 1527 dep = -dep - shift; 1528 CHKERRQ(ISGetLocalSize(pIS, &numPoints)); 1529 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(sdm, newp, &coneSize)); 1540 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, coldp, &val)); 1546 if (val == dep-1) { 1547 CHKERRQ(PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp)); 1548 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(ISRestoreIndices(pIS, &points)); 1558 CHKERRQ(ISDestroy(&pIS)); 1559 } 1560 /* Step 7: Coordinates */ 1561 CHKERRQ(DMPlexShiftCoordinates_Internal(dm, depthShift, sdm)); 1562 CHKERRQ(DMGetCoordinateSection(sdm, &coordSection)); 1563 CHKERRQ(DMGetCoordinatesLocal(sdm, &coordinates)); 1564 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSection, newp, &dof)); 1571 CHKERRQ(PetscSectionGetOffset(coordSection, newp, &off)); 1572 CHKERRQ(PetscSectionGetOffset(coordSection, splitp, &soff)); 1573 for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d]; 1574 } 1575 CHKERRQ(VecRestoreArray(coordinates, &coords)); 1576 /* Step 8: SF, if I can figure this out we can split the mesh in parallel */ 1577 CHKERRQ(DMPlexShiftSF_Internal(dm, depthShift, sdm)); 1578 /* Step 9: Labels */ 1579 CHKERRQ(DMPlexShiftLabels_Internal(dm, depthShift, sdm)); 1580 CHKERRQ(DMPlexCreateVTKLabel_Internal(dm, PETSC_FALSE, sdm)); 1581 CHKERRQ(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 CHKERRQ(DMLabelSetValue(splitLabel, newp, val)); 1592 CHKERRQ(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 CHKERRQ(DMGetLabelName(sdm, l, &lname)); 1601 CHKERRQ(PetscStrcmp(lname, "depth", &isDepth)); 1602 if (isDepth) continue; 1603 CHKERRQ(DMGetLabel(sdm, lname, &mlabel)); 1604 CHKERRQ(DMLabelGetValue(mlabel, newp, &val)); 1605 if (val >= 0) { 1606 CHKERRQ(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]) CHKERRQ(ISRestoreIndices(splitIS[dep], &splitPoints[dep])); 1616 CHKERRQ(ISDestroy(&splitIS[dep])); 1617 if (unsplitIS[dep]) CHKERRQ(ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep])); 1618 CHKERRQ(ISDestroy(&unsplitIS[dep])); 1619 } 1620 if (label) { 1621 CHKERRQ(ISRestoreIndices(valueIS, &values)); 1622 CHKERRQ(ISDestroy(&valueIS)); 1623 } 1624 for (d = 0; d <= depth; ++d) { 1625 CHKERRQ(DMPlexGetDepthStratum(sdm, d, NULL, &pEnd)); 1626 pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d]; 1627 } 1628 CHKERRQ(PetscFree3(coneNew, coneONew, supportNew)); 1629 CHKERRQ(PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld)); 1630 CHKERRQ(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 CHKERRQ(DMCreate(PetscObjectComm((PetscObject)dm), &sdm)); 1660 CHKERRQ(DMSetType(sdm, DMPLEX)); 1661 CHKERRQ(DMGetDimension(dm, &dim)); 1662 CHKERRQ(DMSetDimension(sdm, dim)); 1663 switch (dim) { 1664 case 2: 1665 case 3: 1666 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 1684 CHKERRQ(DMPlexGetConeSize(dm, cell, &coneSize)); 1685 CHKERRQ(DMPlexGetCone(dm, cell, &cone)); 1686 CHKERRQ(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 CHKERRQ(PetscFindInt(cell, numSubpoints, subpoints, &subpoint)); 1696 CHKERRQ(PetscFindInt(face, numSubpoints, subpoints, &subface)); 1697 CHKERRQ(DMPlexGetConeSize(subdm, subpoint, &subconeSize)); 1698 CHKERRQ(DMPlexGetCone(subdm, subpoint, &subcone)); 1699 CHKERRQ(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 PetscCheckFalse(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 PetscCheckFalse(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 1747 CHKERRQ(DMGetDimension(dm, &dim)); 1748 CHKERRQ(DMPlexGetDepthLabel(dm, &depthLabel)); 1749 if (subdm) { 1750 CHKERRQ(DMPlexGetSubpointIS(subdm, &subpointIS)); 1751 if (subpointIS) { 1752 CHKERRQ(ISGetLocalSize(subpointIS, &numSubpoints)); 1753 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, dim-1, &dimIS)); 1758 if (!dimIS) PetscFunctionReturn(0); 1759 CHKERRQ(ISGetLocalSize(dimIS, &numPoints)); 1760 CHKERRQ(ISGetIndices(dimIS, &points)); 1761 for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */ 1762 const PetscInt *support; 1763 PetscInt supportSize, s; 1764 1765 CHKERRQ(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 CHKERRQ(PetscSFGetGraph(dm->sf, NULL, &Nlp, &lp, NULL)); 1774 CHKERRQ(PetscFindInt(points[p], Nlp, lp, &pind)); 1775 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos)); 1785 if (pos) CHKERRQ(DMLabelSetValue(label, support[s], rev*(shift+dim))); 1786 else CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, support[s], &coneSize)); 1790 CHKERRQ(DMPlexGetCone(dm, support[s], &cone)); 1791 for (c = 0; c < coneSize; ++c) { 1792 const PetscInt point = cone[c]; 1793 1794 CHKERRQ(DMLabelGetValue(label, point, &val)); 1795 if (val == -1) { 1796 PetscInt *closure = NULL; 1797 PetscInt closureSize, cl; 1798 1799 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, clp, &val)); 1805 if (blabel) CHKERRQ(DMLabelGetValue(blabel, clp, &bval)); 1806 if ((val >= 0) && (val < dim-1) && (bval < 0)) { 1807 CHKERRQ(DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1))); 1808 break; 1809 } 1810 } 1811 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 1812 } 1813 } 1814 } 1815 } 1816 CHKERRQ(ISRestoreIndices(dimIS, &points)); 1817 CHKERRQ(ISDestroy(&dimIS)); 1818 if (subpointIS) CHKERRQ(ISRestoreIndices(subpointIS, &subpoints)); 1819 /* Mark boundary points as unsplit */ 1820 if (blabel) { 1821 CHKERRQ(DMLabelGetStratumIS(blabel, 1, &dimIS)); 1822 CHKERRQ(ISGetLocalSize(dimIS, &numPoints)); 1823 CHKERRQ(ISGetIndices(dimIS, &points)); 1824 for (p = 0; p < numPoints; ++p) { 1825 const PetscInt point = points[p]; 1826 PetscInt val, bval; 1827 1828 CHKERRQ(DMLabelGetValue(blabel, point, &bval)); 1829 if (bval >= 0) { 1830 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, point, &val)); 1848 PetscCheckFalse((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 CHKERRQ(DMLabelClearValue(label, point, val)); 1850 CHKERRQ(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 CHKERRQ(DMPlexGetSupport(dm, point, &support)); 1855 CHKERRQ(DMPlexGetSupportSize(dm, point, &supportSize)); 1856 for (s = 0; s < supportSize; ++s) { 1857 CHKERRQ(DMPlexGetCone(dm, support[s], &cone)); 1858 CHKERRQ(DMPlexGetConeSize(dm, support[s], &coneSize)); 1859 PetscCheckFalse(coneSize != 2,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize); 1860 CHKERRQ(DMLabelGetValue(blabel, cone[0], &valA)); 1861 CHKERRQ(DMLabelGetValue(blabel, cone[1], &valB)); 1862 CHKERRQ(DMLabelGetValue(blabel, support[s], &valE)); 1863 if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) CHKERRQ(DMLabelSetValue(blabel, support[s], 2)); 1864 } 1865 } 1866 } 1867 CHKERRQ(ISRestoreIndices(dimIS, &points)); 1868 CHKERRQ(ISDestroy(&dimIS)); 1869 } 1870 /* Search for other cells/faces/edges connected to the fault by a vertex */ 1871 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1872 CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 1873 CHKERRQ(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) CHKERRQ(DMLabelGetStratumIS(blabel, 2, &crossEdgeIS)); 1876 if (dimIS && crossEdgeIS) { 1877 IS vertIS = dimIS; 1878 1879 CHKERRQ(ISExpand(vertIS, crossEdgeIS, &dimIS)); 1880 CHKERRQ(ISDestroy(&crossEdgeIS)); 1881 CHKERRQ(ISDestroy(&vertIS)); 1882 } 1883 if (!dimIS) { 1884 PetscFunctionReturn(0); 1885 } 1886 CHKERRQ(ISGetLocalSize(dimIS, &numPoints)); 1887 CHKERRQ(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 CHKERRQ(DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star)); 1895 while (again) { 1896 PetscCheckFalse(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 CHKERRQ(DMLabelGetValue(label, point, &val)); 1905 if (val != -1) continue; 1906 again = again == 1 ? 1 : 2; 1907 CHKERRQ(DMPlexGetConeSize(dm, point, &coneSize)); 1908 CHKERRQ(DMPlexGetCone(dm, point, &cone)); 1909 for (c = 0; c < coneSize; ++c) { 1910 CHKERRQ(DMLabelGetValue(label, cone[c], &val)); 1911 if (val != -1) { 1912 const PetscInt *ccone; 1913 PetscInt cconeSize, cc, side; 1914 1915 PetscCheckFalse(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 CHKERRQ(DMLabelSetValue(label, point, side*(shift+dim))); 1919 /* Mark cell faces which touch the fault */ 1920 CHKERRQ(DMPlexGetConeSize(dm, point, &cconeSize)); 1921 CHKERRQ(DMPlexGetCone(dm, point, &ccone)); 1922 for (cc = 0; cc < cconeSize; ++cc) { 1923 PetscInt *closure = NULL; 1924 PetscInt closureSize, cl; 1925 1926 CHKERRQ(DMLabelGetValue(label, ccone[cc], &val)); 1927 if (val != -1) continue; 1928 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, clp, &val)); 1933 if (val == -1) continue; 1934 CHKERRQ(DMLabelSetValue(label, ccone[cc], side*(shift+dim-1))); 1935 break; 1936 } 1937 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, point, &val)); 1950 if (val == -1) { 1951 PetscInt *sstar = NULL; 1952 PetscInt sstarSize, ss; 1953 PetscBool marked = PETSC_FALSE, isHybrid; 1954 1955 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, spoint, &val)); 1961 PetscCheckFalse(val == -1,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point); 1962 CHKERRQ(DMLabelGetValue(depthLabel, point, &dep)); 1963 if (val > 0) { 1964 CHKERRQ(DMLabelSetValue(label, point, shift+dep)); 1965 } else { 1966 CHKERRQ(DMLabelSetValue(label, point, -(shift+dep))); 1967 } 1968 marked = PETSC_TRUE; 1969 break; 1970 } 1971 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar)); 1972 CHKERRQ(DMPlexCellIsHybrid_Internal(dm, point, &isHybrid)); 1973 PetscCheckFalse(!isHybrid && !marked,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point); 1974 } 1975 } 1976 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star)); 1977 } 1978 CHKERRQ(ISRestoreIndices(dimIS, &points)); 1979 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(label, shift+dim-1, &facePosIS)); 1983 CHKERRQ(DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS)); 1984 CHKERRQ(ISExpand(facePosIS, faceNegIS, &dimIS)); 1985 CHKERRQ(ISDestroy(&facePosIS)); 1986 CHKERRQ(ISDestroy(&faceNegIS)); 1987 CHKERRQ(ISGetLocalSize(dimIS, &numPoints)); 1988 CHKERRQ(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 CHKERRQ(DMPlexGetSupportSize(dm, point, &supportSize)); 1995 if (supportSize != 2) continue; 1996 CHKERRQ(DMPlexGetSupport(dm, point, &support)); 1997 CHKERRQ(DMLabelGetValue(label, support[0], &valA)); 1998 CHKERRQ(DMLabelGetValue(label, support[1], &valB)); 1999 if ((valA == -1) || (valB == -1)) continue; 2000 if (valA*valB > 0) continue; 2001 /* Split the face */ 2002 CHKERRQ(DMLabelGetValue(label, point, &valA)); 2003 CHKERRQ(DMLabelClearValue(label, point, valA)); 2004 CHKERRQ(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 CHKERRQ(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2015 for (cl = 0; cl < closureSize*2; cl += 2) { 2016 CHKERRQ(DMLabelGetValue(label, closure[cl], &valA)); 2017 if (valA == -1) { /* Mark as unsplit */ 2018 CHKERRQ(DMLabelGetValue(depthLabel, closure[cl], &dep)); 2019 CHKERRQ(DMLabelSetValue(label, closure[cl], shift2+dep)); 2020 } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) { 2021 CHKERRQ(DMLabelGetValue(depthLabel, closure[cl], &dep)); 2022 CHKERRQ(DMLabelClearValue(label, closure[cl], valA)); 2023 CHKERRQ(DMLabelSetValue(label, closure[cl], dep)); 2024 } 2025 } 2026 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure)); 2027 } 2028 } 2029 CHKERRQ(ISRestoreIndices(dimIS, &points)); 2030 CHKERRQ(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 CHKERRQ(DMLabelGetDefaultValue(label, &defaultValue)); 2044 CHKERRQ(DMPlexGetSubpointIS(subdm, &subpointIS)); 2045 if (!subpointIS) PetscFunctionReturn(0); 2046 CHKERRQ(DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd)); 2047 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 2048 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, closure[cl], &value)); 2060 if (value == defaultValue) {invalidCell = PETSC_FALSE; break;} 2061 } 2062 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure)); 2063 if (invalidCell) { 2064 CHKERRQ(ISRestoreIndices(subpointIS, &dmpoints)); 2065 CHKERRQ(ISDestroy(&subpointIS)); 2066 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 2121 CHKERRQ(DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm)); 2122 CHKERRQ(DMPlexCheckValidSubmesh_Private(dm, label, idm)); 2123 CHKERRQ(DMPlexOrient(idm)); 2124 CHKERRQ(DMPlexGetSubpointMap(idm, &subpointMap)); 2125 CHKERRQ(DMLabelDuplicate(subpointMap, &hlabel)); 2126 CHKERRQ(DMLabelClearStratum(hlabel, dim)); 2127 if (splitLabel) { 2128 const char *name; 2129 char sname[PETSC_MAX_PATH_LEN]; 2130 2131 CHKERRQ(PetscObjectGetName((PetscObject) hlabel, &name)); 2132 CHKERRQ(PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN)); 2133 CHKERRQ(PetscStrcat(sname, " split")); 2134 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, sname, &slabel)); 2135 } 2136 CHKERRQ(DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm)); 2137 if (dmInterface) {*dmInterface = idm;} 2138 else CHKERRQ(DMDestroy(&idm)); 2139 CHKERRQ(DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid)); 2140 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *dmHybrid)); 2141 if (hybridLabel) *hybridLabel = hlabel; 2142 else CHKERRQ(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 CHKERRQ(DMGetCoordinateDM(*dmHybrid, &cdm)); 2150 CHKERRQ(DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel)); 2151 CHKERRQ(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 CHKERRQ(DMPlexGetDepth(dm, &depth)); 2171 CHKERRQ(DMGetDimension(dm, &dim)); 2172 pSize = PetscMax(depth, dim) + 1; 2173 CHKERRQ(PetscMalloc2(pSize, &pStart, pSize, &pEnd)); 2174 for (d = 0; d <= depth; ++d) { 2175 CHKERRQ(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) CHKERRQ(DMLabelGetStratumIS(vertexLabel, value, &subvertexIS)); 2179 if (subvertexIS) { 2180 CHKERRQ(ISGetSize(subvertexIS, &numSubVerticesInitial)); 2181 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(subpointMap, cell, &cellLoc)); 2200 if (cellLoc == 2) continue; 2201 PetscCheckFalse(cellLoc >= 0,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 2202 CHKERRQ(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 CHKERRQ(DMLabelGetValue(vertexLabel, point, &vertexLoc)); 2210 if (vertexLoc == value) closure[faceSize++] = point; 2211 } 2212 } 2213 if (!(*nFV)) CHKERRQ(DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV)); 2214 PetscCheckFalse(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 CHKERRQ(DMLabelSetValue(subpointMap, closure[cl], 0)); 2222 } 2223 CHKERRQ(DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells)); 2224 for (nc = 0; nc < numCells; ++nc) { 2225 CHKERRQ(DMLabelSetValue(subpointMap, cells[nc], 2)); 2226 } 2227 CHKERRQ(DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells)); 2228 } 2229 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 2230 } 2231 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star)); 2232 } 2233 if (subvertexIS) { 2234 CHKERRQ(ISRestoreIndices(subvertexIS, &subvertices)); 2235 } 2236 CHKERRQ(ISDestroy(&subvertexIS)); 2237 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 2250 CHKERRQ(PetscMalloc2(dim+1, &pStart, dim+1, &pEnd)); 2251 for (d = 0; d <= dim; ++d) { 2252 CHKERRQ(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 CHKERRQ(DMLabelGetStratumIS(vertexLabel, value, &subvertexIS)); 2257 if (subvertexIS) { 2258 CHKERRQ(ISGetSize(subvertexIS, &numSubVerticesInitial)); 2259 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelGetValue(subpointMap, face, &faceLoc)); 2287 if (faceLoc == dim-1) continue; 2288 PetscCheckFalse(faceLoc >= 0,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 2289 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelSetValue(subpointMap, point, d)); 2309 break; 2310 } 2311 } 2312 } 2313 CHKERRQ(DMPlexGetSupportSize(dm, face, &supportSize)); 2314 CHKERRQ(DMPlexGetSupport(dm, face, &support)); 2315 for (s = 0; s < supportSize; ++s) { 2316 CHKERRQ(DMLabelSetValue(subpointMap, support[s], dim)); 2317 } 2318 } 2319 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure)); 2320 } 2321 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star)); 2322 } 2323 if (subvertexIS) CHKERRQ(ISRestoreIndices(subvertexIS, &subvertices)); 2324 CHKERRQ(ISDestroy(&subvertexIS)); 2325 CHKERRQ(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) CHKERRQ(DMGetLabel(dm, labelname, &label)); 2339 *subCells = NULL; 2340 CHKERRQ(DMGetDimension(dm, &dim)); 2341 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, c, &val)); 2348 if (val == value) { 2349 ++(*numFaces); 2350 CHKERRQ(DMPlexGetConeSize(dm, c, &coneSize)); 2351 } 2352 } 2353 } else { 2354 *numFaces = cEnd - cMax; 2355 CHKERRQ(DMPlexGetConeSize(dm, cMax, &coneSize)); 2356 } 2357 CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, c, &val)); 2368 if (val != value) continue; 2369 } 2370 CHKERRQ(DMPlexGetCone(dm, c, &cone)); 2371 for (p = 0; p < *nFV; ++p) { 2372 CHKERRQ(DMLabelSetValue(subpointMap, cone[p], 0)); 2373 } 2374 /* Negative face */ 2375 CHKERRQ(DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells)); 2376 /* Not true in parallel 2377 PetscCheckFalse(numCells != 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2378 for (p = 0; p < numCells; ++p) { 2379 CHKERRQ(DMLabelSetValue(subpointMap, cells[p], 2)); 2380 (*subCells)[subc++] = cells[p]; 2381 } 2382 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 2395 CHKERRQ(DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd)); 2396 if (cMax < 0) PetscFunctionReturn(0); 2397 CHKERRQ(PetscMalloc2(dim+1,&pStart,dim+1,&pEnd)); 2398 for (d = 0; d <= dim; ++d) CHKERRQ(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 CHKERRQ(DMLabelGetValue(label, c, &val)); 2406 if (val != value) continue; 2407 } 2408 CHKERRQ(DMPlexGetConeSize(dm, c, &coneSize)); 2409 CHKERRQ(DMPlexGetCone(dm, c, &cone)); 2410 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelSetValue(subpointMap, point, d)); 2420 break; 2421 } 2422 } 2423 } 2424 CHKERRQ(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 CHKERRQ(DMPlexGetSupportSize(dm, cone[cl], &supportSize)); 2431 /* PetscCheckFalse(supportSize != 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2432 CHKERRQ(DMPlexGetSupport(dm, cone[cl], &support)); 2433 for (s = 0; s < supportSize; ++s) { 2434 CHKERRQ(DMLabelSetValue(subpointMap, support[s], dim)); 2435 } 2436 } 2437 } 2438 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 2451 CHKERRQ(DMGetDimension(dm, &cellDim)); 2452 if (debug) CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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) CHKERRQ(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) CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(dm, cell, &coneSize)); 2767 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 2818 CHKERRQ(DMPlexGetConeSize(subdm, subcell, &coneSize)); 2819 PetscCheckFalse(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 CHKERRQ(DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces)); 2823 #else 2824 { 2825 PetscInt f; 2826 2827 numFaces = 0; 2828 CHKERRQ(DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces)); 2829 for (f = firstFace; f < *newFacePoint; ++f) { 2830 PetscInt dof, off, d; 2831 2832 CHKERRQ(PetscSectionGetDof(submesh->coneSection, f, &dof)); 2833 CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(DMPlexSetCone(subdm, subcell, faces)); 2855 } else { 2856 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2857 PetscBool posOriented; 2858 2859 CHKERRQ(DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices)); 2860 origVertices = &orientedVertices[numFaceVertices]; 2861 indices = &orientedVertices[numFaceVertices*2]; 2862 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2863 CHKERRQ(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 PetscCheckFalse(ov == numFaceVertices,comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2874 } 2875 CHKERRQ(DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices)); 2876 CHKERRQ(DMPlexSetCone(subdm, subcell, newFacePoint)); 2877 CHKERRQ(DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices)); 2878 ++(*newFacePoint); 2879 } 2880 #if 0 2881 CHKERRQ(DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces)); 2882 #else 2883 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 2900 /* Create subpointMap which marks the submesh */ 2901 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 2902 CHKERRQ(DMPlexSetSubpointMap(subdm, subpointMap)); 2903 CHKERRQ(DMLabelDestroy(&subpointMap)); 2904 if (vertexLabel) CHKERRQ(DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm)); 2905 /* Setup chart */ 2906 CHKERRQ(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 2907 CHKERRQ(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 2908 CHKERRQ(DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices)); 2909 CHKERRQ(DMPlexSetVTKCellHeight(subdm, 1)); 2910 /* Set cone sizes */ 2911 firstSubVertex = numSubCells; 2912 firstSubFace = numSubCells+numSubVertices; 2913 newFacePoint = firstSubFace; 2914 CHKERRQ(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 2915 if (subvertexIS) CHKERRQ(ISGetIndices(subvertexIS, &subVertices)); 2916 CHKERRQ(DMLabelGetStratumIS(subpointMap, 2, &subcellIS)); 2917 if (subcellIS) CHKERRQ(ISGetIndices(subcellIS, &subCells)); 2918 for (c = 0; c < numSubCells; ++c) { 2919 CHKERRQ(DMPlexSetConeSize(subdm, c, 1)); 2920 } 2921 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2922 CHKERRQ(DMPlexSetConeSize(subdm, f, nFV)); 2923 } 2924 CHKERRQ(DMSetUp(subdm)); 2925 /* Create face cones */ 2926 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 2927 CHKERRQ(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 2928 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscFindInt(point, numSubVertices, subVertices, &subVertex)); 2943 if (subVertex >= 0) { 2944 closure[faceSize] = point; 2945 subface[faceSize] = firstSubVertex+subVertex; 2946 ++faceSize; 2947 } 2948 } 2949 } 2950 PetscCheckFalse(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 CHKERRQ(DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint)); 2953 } 2954 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 2955 } 2956 CHKERRQ(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 2957 CHKERRQ(DMPlexSymmetrize(subdm)); 2958 CHKERRQ(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 CHKERRQ(DMGetCoordinateSection(dm, &coordSection)); 2968 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates)); 2969 CHKERRQ(DMGetCoordinateSection(subdm, &subCoordSection)); 2970 CHKERRQ(PetscSectionSetNumFields(subCoordSection, 1)); 2971 CHKERRQ(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 2972 CHKERRQ(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 2973 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSection, vertex, &dof)); 2980 CHKERRQ(PetscSectionSetDof(subCoordSection, subvertex, dof)); 2981 CHKERRQ(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 2982 } 2983 CHKERRQ(PetscSectionSetUp(subCoordSection)); 2984 CHKERRQ(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 2985 CHKERRQ(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 2986 CHKERRQ(PetscObjectGetName((PetscObject)coordinates,&name)); 2987 CHKERRQ(PetscObjectSetName((PetscObject)subCoordinates,name)); 2988 CHKERRQ(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 2989 CHKERRQ(VecSetType(subCoordinates,VECSTANDARD)); 2990 if (coordSize) { 2991 CHKERRQ(VecGetArray(coordinates, &coords)); 2992 CHKERRQ(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 CHKERRQ(PetscSectionGetDof(coordSection, vertex, &dof)); 2999 CHKERRQ(PetscSectionGetOffset(coordSection, vertex, &off)); 3000 CHKERRQ(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3001 CHKERRQ(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3002 PetscCheckFalse(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 CHKERRQ(VecRestoreArray(coordinates, &coords)); 3006 CHKERRQ(VecRestoreArray(subCoordinates, &subCoords)); 3007 } 3008 CHKERRQ(DMSetCoordinatesLocal(subdm, subCoordinates)); 3009 CHKERRQ(VecDestroy(&subCoordinates)); 3010 } 3011 /* Cleanup */ 3012 if (subvertexIS) CHKERRQ(ISRestoreIndices(subvertexIS, &subVertices)); 3013 CHKERRQ(ISDestroy(&subvertexIS)); 3014 if (subcellIS) CHKERRQ(ISRestoreIndices(subcellIS, &subCells)); 3015 CHKERRQ(ISDestroy(&subcellIS)); 3016 PetscFunctionReturn(0); 3017 } 3018 3019 static inline PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 3020 { 3021 PetscInt subPoint; 3022 PetscErrorCode ierr; 3023 3024 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 3025 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 3026 } 3027 3028 static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm) 3029 { 3030 PetscInt Nl, l, d; 3031 3032 PetscFunctionBegin; 3033 CHKERRQ(DMGetNumLabels(dm, &Nl)); 3034 for (l = 0; l < Nl; ++l) { 3035 DMLabel label, newlabel; 3036 const char *lname; 3037 PetscBool isDepth, isDim, isCelltype, isVTK; 3038 IS valueIS; 3039 const PetscInt *values; 3040 PetscInt Nv, v; 3041 3042 CHKERRQ(DMGetLabelName(dm, l, &lname)); 3043 CHKERRQ(PetscStrcmp(lname, "depth", &isDepth)); 3044 CHKERRQ(PetscStrcmp(lname, "dim", &isDim)); 3045 CHKERRQ(PetscStrcmp(lname, "celltype", &isCelltype)); 3046 CHKERRQ(PetscStrcmp(lname, "vtk", &isVTK)); 3047 if (isDepth || isDim || isCelltype || isVTK) continue; 3048 CHKERRQ(DMCreateLabel(subdm, lname)); 3049 CHKERRQ(DMGetLabel(dm, lname, &label)); 3050 CHKERRQ(DMGetLabel(subdm, lname, &newlabel)); 3051 CHKERRQ(DMLabelGetDefaultValue(label, &v)); 3052 CHKERRQ(DMLabelSetDefaultValue(newlabel, v)); 3053 CHKERRQ(DMLabelGetValueIS(label, &valueIS)); 3054 CHKERRQ(ISGetLocalSize(valueIS, &Nv)); 3055 CHKERRQ(ISGetIndices(valueIS, &values)); 3056 for (v = 0; v < Nv; ++v) { 3057 IS pointIS; 3058 const PetscInt *points; 3059 PetscInt Np, p; 3060 3061 CHKERRQ(DMLabelGetStratumIS(label, values[v], &pointIS)); 3062 CHKERRQ(ISGetLocalSize(pointIS, &Np)); 3063 CHKERRQ(ISGetIndices(pointIS, &points)); 3064 for (p = 0; p < Np; ++p) { 3065 const PetscInt point = points[p]; 3066 PetscInt subp; 3067 3068 CHKERRQ(DMPlexGetPointDepth(dm, point, &d)); 3069 subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]); 3070 if (subp >= 0) CHKERRQ(DMLabelSetValue(newlabel, subp, values[v])); 3071 } 3072 CHKERRQ(ISRestoreIndices(pointIS, &points)); 3073 CHKERRQ(ISDestroy(&pointIS)); 3074 } 3075 CHKERRQ(ISRestoreIndices(valueIS, &values)); 3076 CHKERRQ(ISDestroy(&valueIS)); 3077 } 3078 PetscFunctionReturn(0); 3079 } 3080 3081 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 3082 { 3083 MPI_Comm comm; 3084 DMLabel subpointMap; 3085 IS *subpointIS; 3086 const PetscInt **subpoints; 3087 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 3088 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 3089 PetscMPIInt rank; 3090 3091 PetscFunctionBegin; 3092 CHKERRQ(PetscObjectGetComm((PetscObject)dm,&comm)); 3093 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 3094 /* Create subpointMap which marks the submesh */ 3095 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3096 CHKERRQ(DMPlexSetSubpointMap(subdm, subpointMap)); 3097 if (cellHeight) { 3098 if (isCohesive) CHKERRQ(DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm)); 3099 else CHKERRQ(DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm)); 3100 } else { 3101 DMLabel depth; 3102 IS pointIS; 3103 const PetscInt *points; 3104 PetscInt numPoints=0; 3105 3106 CHKERRQ(DMPlexGetDepthLabel(dm, &depth)); 3107 CHKERRQ(DMLabelGetStratumIS(label, value, &pointIS)); 3108 if (pointIS) { 3109 CHKERRQ(ISGetIndices(pointIS, &points)); 3110 CHKERRQ(ISGetLocalSize(pointIS, &numPoints)); 3111 } 3112 for (p = 0; p < numPoints; ++p) { 3113 PetscInt *closure = NULL; 3114 PetscInt closureSize, c, pdim; 3115 3116 CHKERRQ(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3117 for (c = 0; c < closureSize*2; c += 2) { 3118 CHKERRQ(DMLabelGetValue(depth, closure[c], &pdim)); 3119 CHKERRQ(DMLabelSetValue(subpointMap, closure[c], pdim)); 3120 } 3121 CHKERRQ(DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure)); 3122 } 3123 if (pointIS) CHKERRQ(ISRestoreIndices(pointIS, &points)); 3124 CHKERRQ(ISDestroy(&pointIS)); 3125 } 3126 /* Setup chart */ 3127 CHKERRQ(DMGetDimension(dm, &dim)); 3128 CHKERRQ(PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints)); 3129 for (d = 0; d <= dim; ++d) { 3130 CHKERRQ(DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d])); 3131 totSubPoints += numSubPoints[d]; 3132 } 3133 CHKERRQ(DMPlexSetChart(subdm, 0, totSubPoints)); 3134 CHKERRQ(DMPlexSetVTKCellHeight(subdm, cellHeight)); 3135 /* Set cone sizes */ 3136 firstSubPoint[dim] = 0; 3137 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 3138 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 3139 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 3140 for (d = 0; d <= dim; ++d) { 3141 CHKERRQ(DMLabelGetStratumIS(subpointMap, d, &subpointIS[d])); 3142 if (subpointIS[d]) CHKERRQ(ISGetIndices(subpointIS[d], &subpoints[d])); 3143 } 3144 /* We do not want this label automatically computed, instead we compute it here */ 3145 CHKERRQ(DMCreateLabel(subdm, "celltype")); 3146 for (d = 0; d <= dim; ++d) { 3147 for (p = 0; p < numSubPoints[d]; ++p) { 3148 const PetscInt point = subpoints[d][p]; 3149 const PetscInt subpoint = firstSubPoint[d] + p; 3150 const PetscInt *cone; 3151 PetscInt coneSize, coneSizeNew, c, val; 3152 DMPolytopeType ct; 3153 3154 CHKERRQ(DMPlexGetConeSize(dm, point, &coneSize)); 3155 CHKERRQ(DMPlexSetConeSize(subdm, subpoint, coneSize)); 3156 CHKERRQ(DMPlexGetCellType(dm, point, &ct)); 3157 CHKERRQ(DMPlexSetCellType(subdm, subpoint, ct)); 3158 if (cellHeight && (d == dim)) { 3159 CHKERRQ(DMPlexGetCone(dm, point, &cone)); 3160 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3161 CHKERRQ(DMLabelGetValue(subpointMap, cone[c], &val)); 3162 if (val >= 0) coneSizeNew++; 3163 } 3164 CHKERRQ(DMPlexSetConeSize(subdm, subpoint, coneSizeNew)); 3165 CHKERRQ(DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST)); 3166 } 3167 } 3168 } 3169 CHKERRQ(DMLabelDestroy(&subpointMap)); 3170 CHKERRQ(DMSetUp(subdm)); 3171 /* Set cones */ 3172 CHKERRQ(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3173 CHKERRQ(PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew)); 3174 for (d = 0; d <= dim; ++d) { 3175 for (p = 0; p < numSubPoints[d]; ++p) { 3176 const PetscInt point = subpoints[d][p]; 3177 const PetscInt subpoint = firstSubPoint[d] + p; 3178 const PetscInt *cone, *ornt; 3179 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 3180 3181 if (d == dim-1) { 3182 const PetscInt *support, *cone, *ornt; 3183 PetscInt supportSize, coneSize, s, subc; 3184 3185 CHKERRQ(DMPlexGetSupport(dm, point, &support)); 3186 CHKERRQ(DMPlexGetSupportSize(dm, point, &supportSize)); 3187 for (s = 0; s < supportSize; ++s) { 3188 PetscBool isHybrid; 3189 3190 CHKERRQ(DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid)); 3191 if (!isHybrid) continue; 3192 CHKERRQ(PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc)); 3193 if (subc >= 0) { 3194 const PetscInt ccell = subpoints[d+1][subc]; 3195 3196 CHKERRQ(DMPlexGetCone(dm, ccell, &cone)); 3197 CHKERRQ(DMPlexGetConeSize(dm, ccell, &coneSize)); 3198 CHKERRQ(DMPlexGetConeOrientation(dm, ccell, &ornt)); 3199 for (c = 0; c < coneSize; ++c) { 3200 if (cone[c] == point) { 3201 fornt = ornt[c]; 3202 break; 3203 } 3204 } 3205 break; 3206 } 3207 } 3208 } 3209 CHKERRQ(DMPlexGetConeSize(dm, point, &coneSize)); 3210 CHKERRQ(DMPlexGetConeSize(subdm, subpoint, &subconeSize)); 3211 CHKERRQ(DMPlexGetCone(dm, point, &cone)); 3212 CHKERRQ(DMPlexGetConeOrientation(dm, point, &ornt)); 3213 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3214 CHKERRQ(PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc)); 3215 if (subc >= 0) { 3216 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 3217 orntNew[coneSizeNew] = ornt[c]; 3218 ++coneSizeNew; 3219 } 3220 } 3221 PetscCheckFalse(coneSizeNew != subconeSize,comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 3222 CHKERRQ(DMPlexSetCone(subdm, subpoint, coneNew)); 3223 CHKERRQ(DMPlexSetConeOrientation(subdm, subpoint, orntNew)); 3224 if (fornt < 0) CHKERRQ(DMPlexOrientPoint(subdm, subpoint, fornt)); 3225 } 3226 } 3227 CHKERRQ(PetscFree2(coneNew,orntNew)); 3228 CHKERRQ(DMPlexSymmetrize(subdm)); 3229 CHKERRQ(DMPlexStratify(subdm)); 3230 /* Build coordinates */ 3231 { 3232 PetscSection coordSection, subCoordSection; 3233 Vec coordinates, subCoordinates; 3234 PetscScalar *coords, *subCoords; 3235 PetscInt cdim, numComp, coordSize; 3236 const char *name; 3237 3238 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 3239 CHKERRQ(DMGetCoordinateSection(dm, &coordSection)); 3240 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates)); 3241 CHKERRQ(DMGetCoordinateSection(subdm, &subCoordSection)); 3242 CHKERRQ(PetscSectionSetNumFields(subCoordSection, 1)); 3243 CHKERRQ(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3244 CHKERRQ(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3245 CHKERRQ(PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0])); 3246 for (v = 0; v < numSubPoints[0]; ++v) { 3247 const PetscInt vertex = subpoints[0][v]; 3248 const PetscInt subvertex = firstSubPoint[0]+v; 3249 PetscInt dof; 3250 3251 CHKERRQ(PetscSectionGetDof(coordSection, vertex, &dof)); 3252 CHKERRQ(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3253 CHKERRQ(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3254 } 3255 CHKERRQ(PetscSectionSetUp(subCoordSection)); 3256 CHKERRQ(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3257 CHKERRQ(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3258 CHKERRQ(PetscObjectGetName((PetscObject)coordinates,&name)); 3259 CHKERRQ(PetscObjectSetName((PetscObject)subCoordinates,name)); 3260 CHKERRQ(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3261 CHKERRQ(VecSetBlockSize(subCoordinates, cdim)); 3262 CHKERRQ(VecSetType(subCoordinates,VECSTANDARD)); 3263 CHKERRQ(VecGetArray(coordinates, &coords)); 3264 CHKERRQ(VecGetArray(subCoordinates, &subCoords)); 3265 for (v = 0; v < numSubPoints[0]; ++v) { 3266 const PetscInt vertex = subpoints[0][v]; 3267 const PetscInt subvertex = firstSubPoint[0]+v; 3268 PetscInt dof, off, sdof, soff, d; 3269 3270 CHKERRQ(PetscSectionGetDof(coordSection, vertex, &dof)); 3271 CHKERRQ(PetscSectionGetOffset(coordSection, vertex, &off)); 3272 CHKERRQ(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3273 CHKERRQ(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3274 PetscCheckFalse(dof != sdof,comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3275 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3276 } 3277 CHKERRQ(VecRestoreArray(coordinates, &coords)); 3278 CHKERRQ(VecRestoreArray(subCoordinates, &subCoords)); 3279 CHKERRQ(DMSetCoordinatesLocal(subdm, subCoordinates)); 3280 CHKERRQ(VecDestroy(&subCoordinates)); 3281 } 3282 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3283 { 3284 PetscSF sfPoint, sfPointSub; 3285 IS subpIS; 3286 const PetscSFNode *remotePoints; 3287 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3288 const PetscInt *localPoints, *subpoints; 3289 PetscInt *slocalPoints; 3290 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3291 PetscMPIInt rank; 3292 3293 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 3294 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 3295 CHKERRQ(DMGetPointSF(subdm, &sfPointSub)); 3296 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 3297 CHKERRQ(DMPlexGetChart(subdm, NULL, &numSubroots)); 3298 CHKERRQ(DMPlexGetSubpointIS(subdm, &subpIS)); 3299 if (subpIS) { 3300 CHKERRQ(ISGetIndices(subpIS, &subpoints)); 3301 CHKERRQ(ISGetLocalSize(subpIS, &numSubpoints)); 3302 } 3303 CHKERRQ(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3304 if (numRoots >= 0) { 3305 CHKERRQ(PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners)); 3306 for (p = 0; p < pEnd-pStart; ++p) { 3307 newLocalPoints[p].rank = -2; 3308 newLocalPoints[p].index = -2; 3309 } 3310 /* Set subleaves */ 3311 for (l = 0; l < numLeaves; ++l) { 3312 const PetscInt point = localPoints[l]; 3313 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3314 3315 if (subpoint < 0) continue; 3316 newLocalPoints[point-pStart].rank = rank; 3317 newLocalPoints[point-pStart].index = subpoint; 3318 ++numSubleaves; 3319 } 3320 /* Must put in owned subpoints */ 3321 for (p = pStart; p < pEnd; ++p) { 3322 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3323 3324 if (subpoint < 0) { 3325 newOwners[p-pStart].rank = -3; 3326 newOwners[p-pStart].index = -3; 3327 } else { 3328 newOwners[p-pStart].rank = rank; 3329 newOwners[p-pStart].index = subpoint; 3330 } 3331 } 3332 CHKERRQ(PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3333 CHKERRQ(PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3334 CHKERRQ(PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3335 CHKERRQ(PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3336 CHKERRQ(PetscMalloc1(numSubleaves, &slocalPoints)); 3337 CHKERRQ(PetscMalloc1(numSubleaves, &sremotePoints)); 3338 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3339 const PetscInt point = localPoints[l]; 3340 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3341 3342 if (subpoint < 0) continue; 3343 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3344 slocalPoints[sl] = subpoint; 3345 sremotePoints[sl].rank = newLocalPoints[point].rank; 3346 sremotePoints[sl].index = newLocalPoints[point].index; 3347 PetscCheckFalse(sremotePoints[sl].rank < 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3348 PetscCheckFalse(sremotePoints[sl].index < 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3349 ++sl; 3350 } 3351 PetscCheckFalse(sl + ll != numSubleaves,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3352 CHKERRQ(PetscFree2(newLocalPoints,newOwners)); 3353 CHKERRQ(PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3354 } 3355 if (subpIS) { 3356 CHKERRQ(ISRestoreIndices(subpIS, &subpoints)); 3357 } 3358 } 3359 /* Filter labels */ 3360 CHKERRQ(DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm)); 3361 /* Cleanup */ 3362 for (d = 0; d <= dim; ++d) { 3363 if (subpointIS[d]) CHKERRQ(ISRestoreIndices(subpointIS[d], &subpoints[d])); 3364 CHKERRQ(ISDestroy(&subpointIS[d])); 3365 } 3366 CHKERRQ(PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints)); 3367 PetscFunctionReturn(0); 3368 } 3369 3370 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm) 3371 { 3372 PetscFunctionBegin; 3373 CHKERRQ(DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm)); 3374 PetscFunctionReturn(0); 3375 } 3376 3377 /*@ 3378 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3379 3380 Input Parameters: 3381 + dm - The original mesh 3382 . vertexLabel - The DMLabel marking points contained in the surface 3383 . value - The label value to use 3384 - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked 3385 3386 Output Parameter: 3387 . subdm - The surface mesh 3388 3389 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3390 3391 Level: developer 3392 3393 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3394 @*/ 3395 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm) 3396 { 3397 DMPlexInterpolatedFlag interpolated; 3398 PetscInt dim, cdim; 3399 3400 PetscFunctionBegin; 3401 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3402 PetscValidPointer(subdm, 5); 3403 CHKERRQ(DMGetDimension(dm, &dim)); 3404 CHKERRQ(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3405 CHKERRQ(DMSetType(*subdm, DMPLEX)); 3406 CHKERRQ(DMSetDimension(*subdm, dim-1)); 3407 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 3408 CHKERRQ(DMSetCoordinateDim(*subdm, cdim)); 3409 CHKERRQ(DMPlexIsInterpolated(dm, &interpolated)); 3410 PetscCheckFalse(interpolated == DMPLEX_INTERPOLATED_PARTIAL,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 3411 if (interpolated) { 3412 CHKERRQ(DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm)); 3413 } else { 3414 CHKERRQ(DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm)); 3415 } 3416 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *subdm)); 3417 PetscFunctionReturn(0); 3418 } 3419 3420 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3421 { 3422 MPI_Comm comm; 3423 DMLabel subpointMap; 3424 IS subvertexIS; 3425 const PetscInt *subVertices; 3426 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3427 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3428 PetscInt c, f; 3429 3430 PetscFunctionBegin; 3431 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 3432 /* Create subpointMap which marks the submesh */ 3433 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap)); 3434 CHKERRQ(DMPlexSetSubpointMap(subdm, subpointMap)); 3435 CHKERRQ(DMLabelDestroy(&subpointMap)); 3436 CHKERRQ(DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm)); 3437 /* Setup chart */ 3438 CHKERRQ(DMLabelGetStratumSize(subpointMap, 0, &numSubVertices)); 3439 CHKERRQ(DMLabelGetStratumSize(subpointMap, 2, &numSubCells)); 3440 CHKERRQ(DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices)); 3441 CHKERRQ(DMPlexSetVTKCellHeight(subdm, 1)); 3442 /* Set cone sizes */ 3443 firstSubVertex = numSubCells; 3444 firstSubFace = numSubCells+numSubVertices; 3445 newFacePoint = firstSubFace; 3446 CHKERRQ(DMLabelGetStratumIS(subpointMap, 0, &subvertexIS)); 3447 if (subvertexIS) CHKERRQ(ISGetIndices(subvertexIS, &subVertices)); 3448 for (c = 0; c < numSubCells; ++c) { 3449 CHKERRQ(DMPlexSetConeSize(subdm, c, 1)); 3450 } 3451 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3452 CHKERRQ(DMPlexSetConeSize(subdm, f, nFV)); 3453 } 3454 CHKERRQ(DMSetUp(subdm)); 3455 /* Create face cones */ 3456 CHKERRQ(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 3457 CHKERRQ(DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 3458 for (c = 0; c < numSubCells; ++c) { 3459 const PetscInt cell = subCells[c]; 3460 const PetscInt subcell = c; 3461 const PetscInt *cone, *cells; 3462 PetscBool isHybrid; 3463 PetscInt numCells, subVertex, p, v; 3464 3465 CHKERRQ(DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid)); 3466 if (!isHybrid) continue; 3467 CHKERRQ(DMPlexGetCone(dm, cell, &cone)); 3468 for (v = 0; v < nFV; ++v) { 3469 CHKERRQ(PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex)); 3470 subface[v] = firstSubVertex+subVertex; 3471 } 3472 CHKERRQ(DMPlexSetCone(subdm, newFacePoint, subface)); 3473 CHKERRQ(DMPlexSetCone(subdm, subcell, &newFacePoint)); 3474 CHKERRQ(DMPlexGetJoin(dm, nFV, cone, &numCells, &cells)); 3475 /* Not true in parallel 3476 PetscCheckFalse(numCells != 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3477 for (p = 0; p < numCells; ++p) { 3478 PetscInt negsubcell; 3479 PetscBool isHybrid; 3480 3481 CHKERRQ(DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid)); 3482 if (isHybrid) continue; 3483 /* I know this is a crap search */ 3484 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3485 if (subCells[negsubcell] == cells[p]) break; 3486 } 3487 PetscCheckFalse(negsubcell == numSubCells,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3488 CHKERRQ(DMPlexSetCone(subdm, negsubcell, &newFacePoint)); 3489 } 3490 CHKERRQ(DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells)); 3491 ++newFacePoint; 3492 } 3493 CHKERRQ(DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface)); 3494 CHKERRQ(DMPlexSymmetrize(subdm)); 3495 CHKERRQ(DMPlexStratify(subdm)); 3496 /* Build coordinates */ 3497 { 3498 PetscSection coordSection, subCoordSection; 3499 Vec coordinates, subCoordinates; 3500 PetscScalar *coords, *subCoords; 3501 PetscInt cdim, numComp, coordSize, v; 3502 const char *name; 3503 3504 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 3505 CHKERRQ(DMGetCoordinateSection(dm, &coordSection)); 3506 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates)); 3507 CHKERRQ(DMGetCoordinateSection(subdm, &subCoordSection)); 3508 CHKERRQ(PetscSectionSetNumFields(subCoordSection, 1)); 3509 CHKERRQ(PetscSectionGetFieldComponents(coordSection, 0, &numComp)); 3510 CHKERRQ(PetscSectionSetFieldComponents(subCoordSection, 0, numComp)); 3511 CHKERRQ(PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices)); 3512 for (v = 0; v < numSubVertices; ++v) { 3513 const PetscInt vertex = subVertices[v]; 3514 const PetscInt subvertex = firstSubVertex+v; 3515 PetscInt dof; 3516 3517 CHKERRQ(PetscSectionGetDof(coordSection, vertex, &dof)); 3518 CHKERRQ(PetscSectionSetDof(subCoordSection, subvertex, dof)); 3519 CHKERRQ(PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof)); 3520 } 3521 CHKERRQ(PetscSectionSetUp(subCoordSection)); 3522 CHKERRQ(PetscSectionGetStorageSize(subCoordSection, &coordSize)); 3523 CHKERRQ(VecCreate(PETSC_COMM_SELF, &subCoordinates)); 3524 CHKERRQ(PetscObjectGetName((PetscObject)coordinates,&name)); 3525 CHKERRQ(PetscObjectSetName((PetscObject)subCoordinates,name)); 3526 CHKERRQ(VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE)); 3527 CHKERRQ(VecSetBlockSize(subCoordinates, cdim)); 3528 CHKERRQ(VecSetType(subCoordinates,VECSTANDARD)); 3529 CHKERRQ(VecGetArray(coordinates, &coords)); 3530 CHKERRQ(VecGetArray(subCoordinates, &subCoords)); 3531 for (v = 0; v < numSubVertices; ++v) { 3532 const PetscInt vertex = subVertices[v]; 3533 const PetscInt subvertex = firstSubVertex+v; 3534 PetscInt dof, off, sdof, soff, d; 3535 3536 CHKERRQ(PetscSectionGetDof(coordSection, vertex, &dof)); 3537 CHKERRQ(PetscSectionGetOffset(coordSection, vertex, &off)); 3538 CHKERRQ(PetscSectionGetDof(subCoordSection, subvertex, &sdof)); 3539 CHKERRQ(PetscSectionGetOffset(subCoordSection, subvertex, &soff)); 3540 PetscCheckFalse(dof != sdof,comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3541 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3542 } 3543 CHKERRQ(VecRestoreArray(coordinates, &coords)); 3544 CHKERRQ(VecRestoreArray(subCoordinates, &subCoords)); 3545 CHKERRQ(DMSetCoordinatesLocal(subdm, subCoordinates)); 3546 CHKERRQ(VecDestroy(&subCoordinates)); 3547 } 3548 /* Build SF */ 3549 CHKMEMQ; 3550 { 3551 PetscSF sfPoint, sfPointSub; 3552 const PetscSFNode *remotePoints; 3553 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3554 const PetscInt *localPoints; 3555 PetscInt *slocalPoints; 3556 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3557 PetscMPIInt rank; 3558 3559 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 3560 CHKERRQ(DMGetPointSF(dm, &sfPoint)); 3561 CHKERRQ(DMGetPointSF(subdm, &sfPointSub)); 3562 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 3563 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 3564 CHKERRQ(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 3565 if (numRoots >= 0) { 3566 /* Only vertices should be shared */ 3567 CHKERRQ(PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners)); 3568 for (p = 0; p < pEnd-pStart; ++p) { 3569 newLocalPoints[p].rank = -2; 3570 newLocalPoints[p].index = -2; 3571 } 3572 /* Set subleaves */ 3573 for (l = 0; l < numLeaves; ++l) { 3574 const PetscInt point = localPoints[l]; 3575 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3576 3577 PetscCheckFalse((point < vStart) && (point >= vEnd),PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3578 if (subPoint < 0) continue; 3579 newLocalPoints[point-pStart].rank = rank; 3580 newLocalPoints[point-pStart].index = subPoint; 3581 ++numSubLeaves; 3582 } 3583 /* Must put in owned subpoints */ 3584 for (p = pStart; p < pEnd; ++p) { 3585 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3586 3587 if (subPoint < 0) { 3588 newOwners[p-pStart].rank = -3; 3589 newOwners[p-pStart].index = -3; 3590 } else { 3591 newOwners[p-pStart].rank = rank; 3592 newOwners[p-pStart].index = subPoint; 3593 } 3594 } 3595 CHKERRQ(PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3596 CHKERRQ(PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC)); 3597 CHKERRQ(PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3598 CHKERRQ(PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE)); 3599 CHKERRQ(PetscMalloc1(numSubLeaves, &slocalPoints)); 3600 CHKERRQ(PetscMalloc1(numSubLeaves, &sremotePoints)); 3601 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3602 const PetscInt point = localPoints[l]; 3603 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3604 3605 if (subPoint < 0) continue; 3606 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3607 slocalPoints[sl] = subPoint; 3608 sremotePoints[sl].rank = newLocalPoints[point].rank; 3609 sremotePoints[sl].index = newLocalPoints[point].index; 3610 PetscCheckFalse(sremotePoints[sl].rank < 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3611 PetscCheckFalse(sremotePoints[sl].index < 0,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3612 ++sl; 3613 } 3614 CHKERRQ(PetscFree2(newLocalPoints,newOwners)); 3615 PetscCheckFalse(sl + ll != numSubLeaves,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3616 CHKERRQ(PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER)); 3617 } 3618 } 3619 CHKMEMQ; 3620 /* Cleanup */ 3621 if (subvertexIS) CHKERRQ(ISRestoreIndices(subvertexIS, &subVertices)); 3622 CHKERRQ(ISDestroy(&subvertexIS)); 3623 CHKERRQ(PetscFree(subCells)); 3624 PetscFunctionReturn(0); 3625 } 3626 3627 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3628 { 3629 DMLabel label = NULL; 3630 3631 PetscFunctionBegin; 3632 if (labelname) CHKERRQ(DMGetLabel(dm, labelname, &label)); 3633 CHKERRQ(DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm)); 3634 PetscFunctionReturn(0); 3635 } 3636 3637 /*@C 3638 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. 3639 3640 Input Parameters: 3641 + dm - The original mesh 3642 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3643 . label - A label name, or NULL 3644 - value - A label value 3645 3646 Output Parameter: 3647 . subdm - The surface mesh 3648 3649 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3650 3651 Level: developer 3652 3653 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3654 @*/ 3655 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3656 { 3657 PetscInt dim, cdim, depth; 3658 3659 PetscFunctionBegin; 3660 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3661 PetscValidPointer(subdm, 5); 3662 CHKERRQ(DMGetDimension(dm, &dim)); 3663 CHKERRQ(DMPlexGetDepth(dm, &depth)); 3664 CHKERRQ(DMCreate(PetscObjectComm((PetscObject)dm), subdm)); 3665 CHKERRQ(DMSetType(*subdm, DMPLEX)); 3666 CHKERRQ(DMSetDimension(*subdm, dim-1)); 3667 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 3668 CHKERRQ(DMSetCoordinateDim(*subdm, cdim)); 3669 if (depth == dim) { 3670 CHKERRQ(DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm)); 3671 } else { 3672 CHKERRQ(DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm)); 3673 } 3674 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *subdm)); 3675 PetscFunctionReturn(0); 3676 } 3677 3678 /*@ 3679 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3680 3681 Input Parameters: 3682 + dm - The original mesh 3683 . cellLabel - The DMLabel marking cells contained in the new mesh 3684 - value - The label value to use 3685 3686 Output Parameter: 3687 . subdm - The new mesh 3688 3689 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3690 3691 Level: developer 3692 3693 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3694 @*/ 3695 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3696 { 3697 PetscInt dim; 3698 3699 PetscFunctionBegin; 3700 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3701 PetscValidPointer(subdm, 4); 3702 CHKERRQ(DMGetDimension(dm, &dim)); 3703 CHKERRQ(DMCreate(PetscObjectComm((PetscObject) dm), subdm)); 3704 CHKERRQ(DMSetType(*subdm, DMPLEX)); 3705 CHKERRQ(DMSetDimension(*subdm, dim)); 3706 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3707 CHKERRQ(DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm)); 3708 CHKERRQ(DMPlexCopy_Internal(dm, PETSC_TRUE, *subdm)); 3709 PetscFunctionReturn(0); 3710 } 3711 3712 /*@ 3713 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3714 3715 Input Parameter: 3716 . dm - The submesh DM 3717 3718 Output Parameter: 3719 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3720 3721 Level: developer 3722 3723 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS() 3724 @*/ 3725 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3726 { 3727 PetscFunctionBegin; 3728 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3729 PetscValidPointer(subpointMap, 2); 3730 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3731 PetscFunctionReturn(0); 3732 } 3733 3734 /*@ 3735 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3736 3737 Input Parameters: 3738 + dm - The submesh DM 3739 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3740 3741 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3742 3743 Level: developer 3744 3745 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS() 3746 @*/ 3747 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3748 { 3749 DM_Plex *mesh = (DM_Plex *) dm->data; 3750 DMLabel tmp; 3751 3752 PetscFunctionBegin; 3753 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3754 tmp = mesh->subpointMap; 3755 mesh->subpointMap = subpointMap; 3756 CHKERRQ(PetscObjectReference((PetscObject) mesh->subpointMap)); 3757 CHKERRQ(DMLabelDestroy(&tmp)); 3758 PetscFunctionReturn(0); 3759 } 3760 3761 static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS) 3762 { 3763 DMLabel spmap; 3764 PetscInt depth, d; 3765 3766 PetscFunctionBegin; 3767 CHKERRQ(DMPlexGetSubpointMap(dm, &spmap)); 3768 CHKERRQ(DMPlexGetDepth(dm, &depth)); 3769 if (spmap && depth >= 0) { 3770 DM_Plex *mesh = (DM_Plex *) dm->data; 3771 PetscInt *points, *depths; 3772 PetscInt pStart, pEnd, p, off; 3773 3774 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 3775 PetscCheck(!pStart,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3776 CHKERRQ(PetscMalloc1(pEnd, &points)); 3777 CHKERRQ(DMGetWorkArray(dm, depth+1, MPIU_INT, &depths)); 3778 depths[0] = depth; 3779 depths[1] = 0; 3780 for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3781 for (d = 0, off = 0; d <= depth; ++d) { 3782 const PetscInt dep = depths[d]; 3783 PetscInt depStart, depEnd, n; 3784 3785 CHKERRQ(DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd)); 3786 CHKERRQ(DMLabelGetStratumSize(spmap, dep, &n)); 3787 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3788 PetscCheckFalse(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); 3789 } else { 3790 if (!n) { 3791 if (d == 0) { 3792 /* Missing cells */ 3793 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3794 } else { 3795 /* Missing faces */ 3796 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3797 } 3798 } 3799 } 3800 if (n) { 3801 IS is; 3802 const PetscInt *opoints; 3803 3804 CHKERRQ(DMLabelGetStratumIS(spmap, dep, &is)); 3805 CHKERRQ(ISGetIndices(is, &opoints)); 3806 for (p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3807 CHKERRQ(ISRestoreIndices(is, &opoints)); 3808 CHKERRQ(ISDestroy(&is)); 3809 } 3810 } 3811 CHKERRQ(DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths)); 3812 PetscCheckFalse(off != pEnd,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3813 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS)); 3814 CHKERRQ(PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState)); 3815 } 3816 PetscFunctionReturn(0); 3817 } 3818 3819 /*@ 3820 DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data 3821 3822 Input Parameter: 3823 . dm - The submesh DM 3824 3825 Output Parameter: 3826 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3827 3828 Note: This IS is guaranteed to be sorted by the construction of the submesh 3829 3830 Level: developer 3831 3832 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3833 @*/ 3834 PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS) 3835 { 3836 DM_Plex *mesh = (DM_Plex *) dm->data; 3837 DMLabel spmap; 3838 PetscObjectState state; 3839 3840 PetscFunctionBegin; 3841 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3842 PetscValidPointer(subpointIS, 2); 3843 CHKERRQ(DMPlexGetSubpointMap(dm, &spmap)); 3844 CHKERRQ(PetscObjectStateGet((PetscObject) spmap, &state)); 3845 if (state != mesh->subpointState || !mesh->subpointIS) CHKERRQ(DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS)); 3846 *subpointIS = mesh->subpointIS; 3847 PetscFunctionReturn(0); 3848 } 3849 3850 /*@ 3851 DMGetEnclosureRelation - Get the relationship between dmA and dmB 3852 3853 Input Parameters: 3854 + dmA - The first DM 3855 - dmB - The second DM 3856 3857 Output Parameter: 3858 . rel - The relation of dmA to dmB 3859 3860 Level: intermediate 3861 3862 .seealso: DMGetEnclosurePoint() 3863 @*/ 3864 PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel) 3865 { 3866 DM plexA, plexB, sdm; 3867 DMLabel spmap; 3868 PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB; 3869 3870 PetscFunctionBegin; 3871 PetscValidPointer(rel, 3); 3872 *rel = DM_ENC_NONE; 3873 if (!dmA || !dmB) PetscFunctionReturn(0); 3874 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 3875 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 3876 if (dmA == dmB) {*rel = DM_ENC_EQUALITY; PetscFunctionReturn(0);} 3877 CHKERRQ(DMConvert(dmA, DMPLEX, &plexA)); 3878 CHKERRQ(DMConvert(dmB, DMPLEX, &plexB)); 3879 CHKERRQ(DMPlexGetChart(plexA, &pStartA, &pEndA)); 3880 CHKERRQ(DMPlexGetChart(plexB, &pStartB, &pEndB)); 3881 /* Assumption 1: subDMs have smaller charts than the DMs that they originate from 3882 - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */ 3883 if ((pStartA == pStartB) && (pEndA == pEndB)) { 3884 *rel = DM_ENC_EQUALITY; 3885 goto end; 3886 } 3887 NpA = pEndA - pStartA; 3888 NpB = pEndB - pStartB; 3889 if (NpA == NpB) goto end; 3890 sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */ 3891 CHKERRQ(DMPlexGetSubpointMap(sdm, &spmap)); 3892 if (!spmap) goto end; 3893 /* TODO Check the space mapped to by subpointMap is same size as dm */ 3894 if (NpA > NpB) { 3895 *rel = DM_ENC_SUPERMESH; 3896 } else { 3897 *rel = DM_ENC_SUBMESH; 3898 } 3899 end: 3900 CHKERRQ(DMDestroy(&plexA)); 3901 CHKERRQ(DMDestroy(&plexB)); 3902 PetscFunctionReturn(0); 3903 } 3904 3905 /*@ 3906 DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB 3907 3908 Input Parameters: 3909 + dmA - The first DM 3910 . dmB - The second DM 3911 . etype - The type of enclosure relation that dmA has to dmB 3912 - pB - A point of dmB 3913 3914 Output Parameter: 3915 . pA - The corresponding point of dmA 3916 3917 Level: intermediate 3918 3919 .seealso: DMGetEnclosureRelation() 3920 @*/ 3921 PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA) 3922 { 3923 DM sdm; 3924 IS subpointIS; 3925 const PetscInt *subpoints; 3926 PetscInt numSubpoints; 3927 3928 PetscFunctionBegin; 3929 /* TODO Cache the IS, making it look like an index */ 3930 switch (etype) { 3931 case DM_ENC_SUPERMESH: 3932 sdm = dmB; 3933 CHKERRQ(DMPlexGetSubpointIS(sdm, &subpointIS)); 3934 CHKERRQ(ISGetIndices(subpointIS, &subpoints)); 3935 *pA = subpoints[pB]; 3936 CHKERRQ(ISRestoreIndices(subpointIS, &subpoints)); 3937 break; 3938 case DM_ENC_SUBMESH: 3939 sdm = dmA; 3940 CHKERRQ(DMPlexGetSubpointIS(sdm, &subpointIS)); 3941 CHKERRQ(ISGetLocalSize(subpointIS, &numSubpoints)); 3942 CHKERRQ(ISGetIndices(subpointIS, &subpoints)); 3943 CHKERRQ(PetscFindInt(pB, numSubpoints, subpoints, pA)); 3944 if (*pA < 0) { 3945 CHKERRQ(DMViewFromOptions(dmA, NULL, "-dm_enc_A_view")); 3946 CHKERRQ(DMViewFromOptions(dmB, NULL, "-dm_enc_B_view")); 3947 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB); 3948 } 3949 CHKERRQ(ISRestoreIndices(subpointIS, &subpoints)); 3950 break; 3951 case DM_ENC_EQUALITY: 3952 case DM_ENC_NONE: 3953 *pA = pB;break; 3954 case DM_ENC_UNKNOWN: 3955 { 3956 DMEnclosureType enc; 3957 3958 CHKERRQ(DMGetEnclosureRelation(dmA, dmB, &enc)); 3959 CHKERRQ(DMGetEnclosurePoint(dmA, dmB, enc, pB, pA)); 3960 } 3961 break; 3962 default: SETERRQ(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype); 3963 } 3964 PetscFunctionReturn(0); 3965 } 3966