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