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