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