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