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