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