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