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