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 /* TODO Why are we including cross edges here? Shouldn't they be in the star of boundary vertices? */ 1883 if (blabel) {ierr = DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);CHKERRQ(ierr);} 1884 if (dimIS && crossEdgeIS) { 1885 IS vertIS = dimIS; 1886 1887 ierr = ISExpand(vertIS, crossEdgeIS, &dimIS);CHKERRQ(ierr); 1888 ierr = ISDestroy(&crossEdgeIS);CHKERRQ(ierr); 1889 ierr = ISDestroy(&vertIS);CHKERRQ(ierr); 1890 } 1891 if (!dimIS) { 1892 PetscFunctionReturn(0); 1893 } 1894 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1895 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1896 for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */ 1897 PetscInt *star = NULL; 1898 PetscInt starSize, s; 1899 PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */ 1900 1901 /* All points connected to the fault are inside a cell, so at the top level we will only check cells */ 1902 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1903 while (again) { 1904 if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault"); 1905 again = 0; 1906 for (s = 0; s < starSize*2; s += 2) { 1907 const PetscInt point = star[s]; 1908 const PetscInt *cone; 1909 PetscInt coneSize, c; 1910 1911 if ((point < cStart) || (point >= cEnd)) continue; 1912 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1913 if (val != -1) continue; 1914 again = again == 1 ? 1 : 2; 1915 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 1916 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 1917 for (c = 0; c < coneSize; ++c) { 1918 ierr = DMLabelGetValue(label, cone[c], &val);CHKERRQ(ierr); 1919 if (val != -1) { 1920 const PetscInt *ccone; 1921 PetscInt cconeSize, cc, side; 1922 1923 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); 1924 if (val > 0) side = 1; 1925 else side = -1; 1926 ierr = DMLabelSetValue(label, point, side*(shift+dim));CHKERRQ(ierr); 1927 /* Mark cell faces which touch the fault */ 1928 ierr = DMPlexGetConeSize(dm, point, &cconeSize);CHKERRQ(ierr); 1929 ierr = DMPlexGetCone(dm, point, &ccone);CHKERRQ(ierr); 1930 for (cc = 0; cc < cconeSize; ++cc) { 1931 PetscInt *closure = NULL; 1932 PetscInt closureSize, cl; 1933 1934 ierr = DMLabelGetValue(label, ccone[cc], &val);CHKERRQ(ierr); 1935 if (val != -1) continue; 1936 ierr = DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1937 for (cl = 0; cl < closureSize*2; cl += 2) { 1938 const PetscInt clp = closure[cl]; 1939 1940 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1941 if (val == -1) continue; 1942 ierr = DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));CHKERRQ(ierr); 1943 break; 1944 } 1945 ierr = DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1946 } 1947 again = 1; 1948 break; 1949 } 1950 } 1951 } 1952 } 1953 /* Classify the rest by cell membership */ 1954 for (s = 0; s < starSize*2; s += 2) { 1955 const PetscInt point = star[s]; 1956 1957 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1958 if (val == -1) { 1959 PetscInt *sstar = NULL; 1960 PetscInt sstarSize, ss; 1961 PetscBool marked = PETSC_FALSE, isHybrid; 1962 1963 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1964 for (ss = 0; ss < sstarSize*2; ss += 2) { 1965 const PetscInt spoint = sstar[ss]; 1966 1967 if ((spoint < cStart) || (spoint >= cEnd)) continue; 1968 ierr = DMLabelGetValue(label, spoint, &val);CHKERRQ(ierr); 1969 if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point); 1970 ierr = DMLabelGetValue(depthLabel, point, &dep);CHKERRQ(ierr); 1971 if (val > 0) { 1972 ierr = DMLabelSetValue(label, point, shift+dep);CHKERRQ(ierr); 1973 } else { 1974 ierr = DMLabelSetValue(label, point, -(shift+dep));CHKERRQ(ierr); 1975 } 1976 marked = PETSC_TRUE; 1977 break; 1978 } 1979 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1980 ierr = DMPlexCellIsHybrid_Internal(dm, point, &isHybrid);CHKERRQ(ierr); 1981 if (!isHybrid && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point); 1982 } 1983 } 1984 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1985 } 1986 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1987 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1988 /* If any faces touching the fault divide cells on either side, split them 1989 This only occurs without a surface boundary */ 1990 ierr = DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);CHKERRQ(ierr); 1991 ierr = DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);CHKERRQ(ierr); 1992 ierr = ISExpand(facePosIS, faceNegIS, &dimIS);CHKERRQ(ierr); 1993 ierr = ISDestroy(&facePosIS);CHKERRQ(ierr); 1994 ierr = ISDestroy(&faceNegIS);CHKERRQ(ierr); 1995 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1996 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1997 for (p = 0; p < numPoints; ++p) { 1998 const PetscInt point = points[p]; 1999 const PetscInt *support; 2000 PetscInt supportSize, valA, valB; 2001 2002 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 2003 if (supportSize != 2) continue; 2004 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 2005 ierr = DMLabelGetValue(label, support[0], &valA);CHKERRQ(ierr); 2006 ierr = DMLabelGetValue(label, support[1], &valB);CHKERRQ(ierr); 2007 if ((valA == -1) || (valB == -1)) continue; 2008 if (valA*valB > 0) continue; 2009 /* Split the face */ 2010 ierr = DMLabelGetValue(label, point, &valA);CHKERRQ(ierr); 2011 ierr = DMLabelClearValue(label, point, valA);CHKERRQ(ierr); 2012 ierr = DMLabelSetValue(label, point, dim-1);CHKERRQ(ierr); 2013 /* Label its closure: 2014 unmarked: label as unsplit 2015 incident: relabel as split 2016 split: do nothing 2017 */ 2018 { 2019 PetscInt *closure = NULL; 2020 PetscInt closureSize, cl; 2021 2022 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2023 for (cl = 0; cl < closureSize*2; cl += 2) { 2024 ierr = DMLabelGetValue(label, closure[cl], &valA);CHKERRQ(ierr); 2025 if (valA == -1) { /* Mark as unsplit */ 2026 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 2027 ierr = DMLabelSetValue(label, closure[cl], shift2+dep);CHKERRQ(ierr); 2028 } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) { 2029 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 2030 ierr = DMLabelClearValue(label, closure[cl], valA);CHKERRQ(ierr); 2031 ierr = DMLabelSetValue(label, closure[cl], dep);CHKERRQ(ierr); 2032 } 2033 } 2034 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2035 } 2036 } 2037 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 2038 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 2039 PetscFunctionReturn(0); 2040 } 2041 2042 /* Check that no cell have all vertices on the fault */ 2043 PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm) 2044 { 2045 IS subpointIS; 2046 const PetscInt *dmpoints; 2047 PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd; 2048 PetscErrorCode ierr; 2049 2050 PetscFunctionBegin; 2051 if (!label) PetscFunctionReturn(0); 2052 ierr = DMLabelGetDefaultValue(label, &defaultValue);CHKERRQ(ierr); 2053 ierr = DMPlexGetSubpointIS(subdm, &subpointIS);CHKERRQ(ierr); 2054 if (!subpointIS) PetscFunctionReturn(0); 2055 ierr = DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2056 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2057 ierr = ISGetIndices(subpointIS, &dmpoints);CHKERRQ(ierr); 2058 for (c = cStart; c < cEnd; ++c) { 2059 PetscBool invalidCell = PETSC_TRUE; 2060 PetscInt *closure = NULL; 2061 PetscInt closureSize, cl; 2062 2063 ierr = DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2064 for (cl = 0; cl < closureSize*2; cl += 2) { 2065 PetscInt value = 0; 2066 2067 if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue; 2068 ierr = DMLabelGetValue(label, closure[cl], &value);CHKERRQ(ierr); 2069 if (value == defaultValue) {invalidCell = PETSC_FALSE; break;} 2070 } 2071 ierr = DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2072 if (invalidCell) { 2073 ierr = ISRestoreIndices(subpointIS, &dmpoints);CHKERRQ(ierr); 2074 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 2075 ierr = DMDestroy(&subdm);CHKERRQ(ierr); 2076 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]); 2077 } 2078 } 2079 ierr = ISRestoreIndices(subpointIS, &dmpoints);CHKERRQ(ierr); 2080 PetscFunctionReturn(0); 2081 } 2082 2083 /*@ 2084 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 2085 2086 Collective on dm 2087 2088 Input Parameters: 2089 + dm - The original DM 2090 . label - The label specifying the interface vertices 2091 - bdlabel - The optional label specifying the interface boundary vertices 2092 2093 Output Parameters: 2094 + hybridLabel - The label fully marking the interface, or NULL if no output is desired 2095 . splitLabel - The label containing the split points, or NULL if no output is desired 2096 . dmInterface - The new interface DM, or NULL 2097 - dmHybrid - The new DM with cohesive cells 2098 2099 Note: The hybridLabel indicates what parts of the original mesh impinged on the on division surface. For points 2100 directly on the division surface, they are labeled with their dimension, so an edge 7 on the division surface would be 2101 7 (1) in hybridLabel. For points that impinge from the positive side, they are labeled with 100+dim, so an edge 6 with 2102 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 2103 surface also hits vertex 3, it would be 9 (-101) in hybridLabel. 2104 2105 The splitLabel indicates what points in the new hybrid mesh were the result of splitting points in the original 2106 mesh. The label value is +=100+dim for each point. For example, if two edges 10 and 14 in the hybrid resulting from 2107 splitting an edge in the original mesh, you would have 10 (101) and 14 (-101) in the splitLabel. 2108 2109 The dmInterface is a DM built from the original division surface. It has a label which can be retrieved using 2110 DMPlexGetSubpointMap() which maps each point back to the point in the surface of the original mesh. 2111 2112 Level: developer 2113 2114 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMPlexGetSubpointMap(), DMCreate() 2115 @*/ 2116 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DMLabel *splitLabel, DM *dmInterface, DM *dmHybrid) 2117 { 2118 DM idm; 2119 DMLabel subpointMap, hlabel, slabel = NULL; 2120 PetscInt dim; 2121 PetscErrorCode ierr; 2122 2123 PetscFunctionBegin; 2124 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2125 if (bdlabel) PetscValidPointer(bdlabel, 3); 2126 if (hybridLabel) PetscValidPointer(hybridLabel, 4); 2127 if (splitLabel) PetscValidPointer(splitLabel, 5); 2128 if (dmInterface) PetscValidPointer(dmInterface, 6); 2129 PetscValidPointer(dmHybrid, 7); 2130 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2131 ierr = DMPlexCreateSubmesh(dm, label, 1, PETSC_FALSE, &idm);CHKERRQ(ierr); 2132 ierr = DMPlexCheckValidSubmesh_Private(dm, label, idm);CHKERRQ(ierr); 2133 ierr = DMPlexOrient(idm);CHKERRQ(ierr); 2134 ierr = DMPlexGetSubpointMap(idm, &subpointMap);CHKERRQ(ierr); 2135 ierr = DMLabelDuplicate(subpointMap, &hlabel);CHKERRQ(ierr); 2136 ierr = DMLabelClearStratum(hlabel, dim);CHKERRQ(ierr); 2137 if (splitLabel) { 2138 const char *name; 2139 char sname[PETSC_MAX_PATH_LEN]; 2140 2141 ierr = PetscObjectGetName((PetscObject) hlabel, &name);CHKERRQ(ierr); 2142 ierr = PetscStrncpy(sname, name, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 2143 ierr = PetscStrcat(sname, " split");CHKERRQ(ierr); 2144 ierr = DMLabelCreate(PETSC_COMM_SELF, sname, &slabel);CHKERRQ(ierr); 2145 } 2146 ierr = DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);CHKERRQ(ierr); 2147 if (dmInterface) {*dmInterface = idm;} 2148 else {ierr = DMDestroy(&idm);CHKERRQ(ierr);} 2149 ierr = DMPlexConstructCohesiveCells(dm, hlabel, slabel, dmHybrid);CHKERRQ(ierr); 2150 if (hybridLabel) *hybridLabel = hlabel; 2151 else {ierr = DMLabelDestroy(&hlabel);CHKERRQ(ierr);} 2152 if (splitLabel) *splitLabel = slabel; 2153 { 2154 DM cdm; 2155 DMLabel ctLabel; 2156 2157 /* We need to somehow share the celltype label with the coordinate dm */ 2158 ierr = DMGetCoordinateDM(*dmHybrid, &cdm);CHKERRQ(ierr); 2159 ierr = DMPlexGetCellTypeLabel(*dmHybrid, &ctLabel);CHKERRQ(ierr); 2160 ierr = DMSetLabel(cdm, ctLabel);CHKERRQ(ierr); 2161 } 2162 PetscFunctionReturn(0); 2163 } 2164 2165 /* Here we need the explicit assumption that: 2166 2167 For any marked cell, the marked vertices constitute a single face 2168 */ 2169 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 2170 { 2171 IS subvertexIS = NULL; 2172 const PetscInt *subvertices; 2173 PetscInt *pStart, *pEnd, pSize; 2174 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 2175 PetscErrorCode ierr; 2176 2177 PetscFunctionBegin; 2178 *numFaces = 0; 2179 *nFV = 0; 2180 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2181 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2182 pSize = PetscMax(depth, dim) + 1; 2183 ierr = PetscMalloc2(pSize, &pStart, pSize, &pEnd);CHKERRQ(ierr); 2184 for (d = 0; d <= depth; ++d) { 2185 ierr = DMPlexGetSimplexOrBoxCells(dm, depth-d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 2186 } 2187 /* Loop over initial vertices and mark all faces in the collective star() */ 2188 if (vertexLabel) {ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr);} 2189 if (subvertexIS) { 2190 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 2191 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2192 } 2193 for (v = 0; v < numSubVerticesInitial; ++v) { 2194 const PetscInt vertex = subvertices[v]; 2195 PetscInt *star = NULL; 2196 PetscInt starSize, s, numCells = 0, c; 2197 2198 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2199 for (s = 0; s < starSize*2; s += 2) { 2200 const PetscInt point = star[s]; 2201 if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point; 2202 } 2203 for (c = 0; c < numCells; ++c) { 2204 const PetscInt cell = star[c]; 2205 PetscInt *closure = NULL; 2206 PetscInt closureSize, cl; 2207 PetscInt cellLoc, numCorners = 0, faceSize = 0; 2208 2209 ierr = DMLabelGetValue(subpointMap, cell, &cellLoc);CHKERRQ(ierr); 2210 if (cellLoc == 2) continue; 2211 if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 2212 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2213 for (cl = 0; cl < closureSize*2; cl += 2) { 2214 const PetscInt point = closure[cl]; 2215 PetscInt vertexLoc; 2216 2217 if ((point >= pStart[0]) && (point < pEnd[0])) { 2218 ++numCorners; 2219 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2220 if (vertexLoc == value) closure[faceSize++] = point; 2221 } 2222 } 2223 if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);} 2224 if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2225 if (faceSize == *nFV) { 2226 const PetscInt *cells = NULL; 2227 PetscInt numCells, nc; 2228 2229 ++(*numFaces); 2230 for (cl = 0; cl < faceSize; ++cl) { 2231 ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr); 2232 } 2233 ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 2234 for (nc = 0; nc < numCells; ++nc) { 2235 ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr); 2236 } 2237 ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 2238 } 2239 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2240 } 2241 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2242 } 2243 if (subvertexIS) { 2244 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2245 } 2246 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2247 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2248 PetscFunctionReturn(0); 2249 } 2250 2251 static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DMLabel subpointMap, DM subdm) 2252 { 2253 IS subvertexIS = NULL; 2254 const PetscInt *subvertices; 2255 PetscInt *pStart, *pEnd; 2256 PetscInt dim, d, numSubVerticesInitial = 0, v; 2257 PetscErrorCode ierr; 2258 2259 PetscFunctionBegin; 2260 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2261 ierr = PetscMalloc2(dim+1, &pStart, dim+1, &pEnd);CHKERRQ(ierr); 2262 for (d = 0; d <= dim; ++d) { 2263 ierr = DMPlexGetSimplexOrBoxCells(dm, dim-d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 2264 } 2265 /* Loop over initial vertices and mark all faces in the collective star() */ 2266 if (vertexLabel) { 2267 ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr); 2268 if (subvertexIS) { 2269 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 2270 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2271 } 2272 } 2273 for (v = 0; v < numSubVerticesInitial; ++v) { 2274 const PetscInt vertex = subvertices[v]; 2275 PetscInt *star = NULL; 2276 PetscInt starSize, s, numFaces = 0, f; 2277 2278 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2279 for (s = 0; s < starSize*2; s += 2) { 2280 const PetscInt point = star[s]; 2281 PetscInt faceLoc; 2282 2283 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) { 2284 if (markedFaces) { 2285 ierr = DMLabelGetValue(vertexLabel, point, &faceLoc);CHKERRQ(ierr); 2286 if (faceLoc < 0) continue; 2287 } 2288 star[numFaces++] = point; 2289 } 2290 } 2291 for (f = 0; f < numFaces; ++f) { 2292 const PetscInt face = star[f]; 2293 PetscInt *closure = NULL; 2294 PetscInt closureSize, c; 2295 PetscInt faceLoc; 2296 2297 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 2298 if (faceLoc == dim-1) continue; 2299 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 2300 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2301 for (c = 0; c < closureSize*2; c += 2) { 2302 const PetscInt point = closure[c]; 2303 PetscInt vertexLoc; 2304 2305 if ((point >= pStart[0]) && (point < pEnd[0])) { 2306 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2307 if (vertexLoc != value) break; 2308 } 2309 } 2310 if (c == closureSize*2) { 2311 const PetscInt *support; 2312 PetscInt supportSize, s; 2313 2314 for (c = 0; c < closureSize*2; c += 2) { 2315 const PetscInt point = closure[c]; 2316 2317 for (d = 0; d < dim; ++d) { 2318 if ((point >= pStart[d]) && (point < pEnd[d])) { 2319 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2320 break; 2321 } 2322 } 2323 } 2324 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2325 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2326 for (s = 0; s < supportSize; ++s) { 2327 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2328 } 2329 } 2330 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2331 } 2332 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2333 } 2334 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr);} 2335 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2336 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2337 PetscFunctionReturn(0); 2338 } 2339 2340 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 2341 { 2342 DMLabel label = NULL; 2343 const PetscInt *cone; 2344 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1; 2345 PetscErrorCode ierr; 2346 2347 PetscFunctionBegin; 2348 *numFaces = 0; 2349 *nFV = 0; 2350 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 2351 *subCells = NULL; 2352 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2353 ierr = DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);CHKERRQ(ierr); 2354 if (cMax < 0) PetscFunctionReturn(0); 2355 if (label) { 2356 for (c = cMax; c < cEnd; ++c) { 2357 PetscInt val; 2358 2359 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2360 if (val == value) { 2361 ++(*numFaces); 2362 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2363 } 2364 } 2365 } else { 2366 *numFaces = cEnd - cMax; 2367 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 2368 } 2369 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 2370 if (!(*numFaces)) PetscFunctionReturn(0); 2371 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 2372 for (c = cMax; c < cEnd; ++c) { 2373 const PetscInt *cells; 2374 PetscInt numCells; 2375 2376 if (label) { 2377 PetscInt val; 2378 2379 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2380 if (val != value) continue; 2381 } 2382 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2383 for (p = 0; p < *nFV; ++p) { 2384 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 2385 } 2386 /* Negative face */ 2387 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2388 /* Not true in parallel 2389 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2390 for (p = 0; p < numCells; ++p) { 2391 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 2392 (*subCells)[subc++] = cells[p]; 2393 } 2394 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2395 /* Positive face is not included */ 2396 } 2397 PetscFunctionReturn(0); 2398 } 2399 2400 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 2401 { 2402 PetscInt *pStart, *pEnd; 2403 PetscInt dim, cMax, cEnd, c, d; 2404 PetscErrorCode ierr; 2405 2406 PetscFunctionBegin; 2407 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2408 ierr = DMPlexGetTensorPrismBounds_Internal(dm, dim, &cMax, &cEnd);CHKERRQ(ierr); 2409 if (cMax < 0) PetscFunctionReturn(0); 2410 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 2411 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 2412 for (c = cMax; c < cEnd; ++c) { 2413 const PetscInt *cone; 2414 PetscInt *closure = NULL; 2415 PetscInt fconeSize, coneSize, closureSize, cl, val; 2416 2417 if (label) { 2418 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2419 if (val != value) continue; 2420 } 2421 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2422 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2423 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 2424 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2425 /* Negative face */ 2426 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2427 for (cl = 0; cl < closureSize*2; cl += 2) { 2428 const PetscInt point = closure[cl]; 2429 2430 for (d = 0; d <= dim; ++d) { 2431 if ((point >= pStart[d]) && (point < pEnd[d])) { 2432 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2433 break; 2434 } 2435 } 2436 } 2437 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2438 /* Cells -- positive face is not included */ 2439 for (cl = 0; cl < 1; ++cl) { 2440 const PetscInt *support; 2441 PetscInt supportSize, s; 2442 2443 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 2444 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2445 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 2446 for (s = 0; s < supportSize; ++s) { 2447 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2448 } 2449 } 2450 } 2451 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2452 PetscFunctionReturn(0); 2453 } 2454 2455 static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2456 { 2457 MPI_Comm comm; 2458 PetscBool posOrient = PETSC_FALSE; 2459 const PetscInt debug = 0; 2460 PetscInt cellDim, faceSize, f; 2461 PetscErrorCode ierr; 2462 2463 PetscFunctionBegin; 2464 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2465 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 2466 if (debug) {ierr = PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 2467 2468 if (cellDim == 1 && numCorners == 2) { 2469 /* Triangle */ 2470 faceSize = numCorners-1; 2471 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2472 } else if (cellDim == 2 && numCorners == 3) { 2473 /* Triangle */ 2474 faceSize = numCorners-1; 2475 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2476 } else if (cellDim == 3 && numCorners == 4) { 2477 /* Tetrahedron */ 2478 faceSize = numCorners-1; 2479 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2480 } else if (cellDim == 1 && numCorners == 3) { 2481 /* Quadratic line */ 2482 faceSize = 1; 2483 posOrient = PETSC_TRUE; 2484 } else if (cellDim == 2 && numCorners == 4) { 2485 /* Quads */ 2486 faceSize = 2; 2487 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2488 posOrient = PETSC_TRUE; 2489 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2490 posOrient = PETSC_TRUE; 2491 } else { 2492 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2493 posOrient = PETSC_FALSE; 2494 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2495 } 2496 } else if (cellDim == 2 && numCorners == 6) { 2497 /* Quadratic triangle (I hate this) */ 2498 /* Edges are determined by the first 2 vertices (corners of edges) */ 2499 const PetscInt faceSizeTri = 3; 2500 PetscInt sortedIndices[3], i, iFace; 2501 PetscBool found = PETSC_FALSE; 2502 PetscInt faceVerticesTriSorted[9] = { 2503 0, 3, 4, /* bottom */ 2504 1, 4, 5, /* right */ 2505 2, 3, 5, /* left */ 2506 }; 2507 PetscInt faceVerticesTri[9] = { 2508 0, 3, 4, /* bottom */ 2509 1, 4, 5, /* right */ 2510 2, 5, 3, /* left */ 2511 }; 2512 2513 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2514 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2515 for (iFace = 0; iFace < 3; ++iFace) { 2516 const PetscInt ii = iFace*faceSizeTri; 2517 PetscInt fVertex, cVertex; 2518 2519 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2520 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2521 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2522 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2523 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2524 faceVertices[fVertex] = origVertices[cVertex]; 2525 break; 2526 } 2527 } 2528 } 2529 found = PETSC_TRUE; 2530 break; 2531 } 2532 } 2533 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2534 if (posOriented) *posOriented = PETSC_TRUE; 2535 PetscFunctionReturn(0); 2536 } else if (cellDim == 2 && numCorners == 9) { 2537 /* Quadratic quad (I hate this) */ 2538 /* Edges are determined by the first 2 vertices (corners of edges) */ 2539 const PetscInt faceSizeQuad = 3; 2540 PetscInt sortedIndices[3], i, iFace; 2541 PetscBool found = PETSC_FALSE; 2542 PetscInt faceVerticesQuadSorted[12] = { 2543 0, 1, 4, /* bottom */ 2544 1, 2, 5, /* right */ 2545 2, 3, 6, /* top */ 2546 0, 3, 7, /* left */ 2547 }; 2548 PetscInt faceVerticesQuad[12] = { 2549 0, 1, 4, /* bottom */ 2550 1, 2, 5, /* right */ 2551 2, 3, 6, /* top */ 2552 3, 0, 7, /* left */ 2553 }; 2554 2555 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2556 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2557 for (iFace = 0; iFace < 4; ++iFace) { 2558 const PetscInt ii = iFace*faceSizeQuad; 2559 PetscInt fVertex, cVertex; 2560 2561 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2562 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2563 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2564 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2565 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2566 faceVertices[fVertex] = origVertices[cVertex]; 2567 break; 2568 } 2569 } 2570 } 2571 found = PETSC_TRUE; 2572 break; 2573 } 2574 } 2575 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2576 if (posOriented) *posOriented = PETSC_TRUE; 2577 PetscFunctionReturn(0); 2578 } else if (cellDim == 3 && numCorners == 8) { 2579 /* Hexes 2580 A hex is two oriented quads with the normal of the first 2581 pointing up at the second. 2582 2583 7---6 2584 /| /| 2585 4---5 | 2586 | 1-|-2 2587 |/ |/ 2588 0---3 2589 2590 Faces are determined by the first 4 vertices (corners of faces) */ 2591 const PetscInt faceSizeHex = 4; 2592 PetscInt sortedIndices[4], i, iFace; 2593 PetscBool found = PETSC_FALSE; 2594 PetscInt faceVerticesHexSorted[24] = { 2595 0, 1, 2, 3, /* bottom */ 2596 4, 5, 6, 7, /* top */ 2597 0, 3, 4, 5, /* front */ 2598 2, 3, 5, 6, /* right */ 2599 1, 2, 6, 7, /* back */ 2600 0, 1, 4, 7, /* left */ 2601 }; 2602 PetscInt faceVerticesHex[24] = { 2603 1, 2, 3, 0, /* bottom */ 2604 4, 5, 6, 7, /* top */ 2605 0, 3, 5, 4, /* front */ 2606 3, 2, 6, 5, /* right */ 2607 2, 1, 7, 6, /* back */ 2608 1, 0, 4, 7, /* left */ 2609 }; 2610 2611 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2612 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2613 for (iFace = 0; iFace < 6; ++iFace) { 2614 const PetscInt ii = iFace*faceSizeHex; 2615 PetscInt fVertex, cVertex; 2616 2617 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2618 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2619 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2620 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2621 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2622 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2623 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2624 faceVertices[fVertex] = origVertices[cVertex]; 2625 break; 2626 } 2627 } 2628 } 2629 found = PETSC_TRUE; 2630 break; 2631 } 2632 } 2633 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2634 if (posOriented) *posOriented = PETSC_TRUE; 2635 PetscFunctionReturn(0); 2636 } else if (cellDim == 3 && numCorners == 10) { 2637 /* Quadratic tet */ 2638 /* Faces are determined by the first 3 vertices (corners of faces) */ 2639 const PetscInt faceSizeTet = 6; 2640 PetscInt sortedIndices[6], i, iFace; 2641 PetscBool found = PETSC_FALSE; 2642 PetscInt faceVerticesTetSorted[24] = { 2643 0, 1, 2, 6, 7, 8, /* bottom */ 2644 0, 3, 4, 6, 7, 9, /* front */ 2645 1, 4, 5, 7, 8, 9, /* right */ 2646 2, 3, 5, 6, 8, 9, /* left */ 2647 }; 2648 PetscInt faceVerticesTet[24] = { 2649 0, 1, 2, 6, 7, 8, /* bottom */ 2650 0, 4, 3, 6, 7, 9, /* front */ 2651 1, 5, 4, 7, 8, 9, /* right */ 2652 2, 3, 5, 8, 6, 9, /* left */ 2653 }; 2654 2655 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2656 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2657 for (iFace=0; iFace < 4; ++iFace) { 2658 const PetscInt ii = iFace*faceSizeTet; 2659 PetscInt fVertex, cVertex; 2660 2661 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2662 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2663 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2664 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2665 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2666 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2667 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2668 faceVertices[fVertex] = origVertices[cVertex]; 2669 break; 2670 } 2671 } 2672 } 2673 found = PETSC_TRUE; 2674 break; 2675 } 2676 } 2677 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2678 if (posOriented) *posOriented = PETSC_TRUE; 2679 PetscFunctionReturn(0); 2680 } else if (cellDim == 3 && numCorners == 27) { 2681 /* Quadratic hexes (I hate this) 2682 A hex is two oriented quads with the normal of the first 2683 pointing up at the second. 2684 2685 7---6 2686 /| /| 2687 4---5 | 2688 | 3-|-2 2689 |/ |/ 2690 0---1 2691 2692 Faces are determined by the first 4 vertices (corners of faces) */ 2693 const PetscInt faceSizeQuadHex = 9; 2694 PetscInt sortedIndices[9], i, iFace; 2695 PetscBool found = PETSC_FALSE; 2696 PetscInt faceVerticesQuadHexSorted[54] = { 2697 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2698 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2699 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2700 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2701 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2702 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2703 }; 2704 PetscInt faceVerticesQuadHex[54] = { 2705 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2706 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2707 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2708 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2709 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2710 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2711 }; 2712 2713 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2714 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2715 for (iFace = 0; iFace < 6; ++iFace) { 2716 const PetscInt ii = iFace*faceSizeQuadHex; 2717 PetscInt fVertex, cVertex; 2718 2719 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2720 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2721 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2722 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2723 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2724 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2725 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2726 faceVertices[fVertex] = origVertices[cVertex]; 2727 break; 2728 } 2729 } 2730 } 2731 found = PETSC_TRUE; 2732 break; 2733 } 2734 } 2735 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2736 if (posOriented) *posOriented = PETSC_TRUE; 2737 PetscFunctionReturn(0); 2738 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2739 if (!posOrient) { 2740 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2741 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2742 } else { 2743 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2744 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2745 } 2746 if (posOriented) *posOriented = posOrient; 2747 PetscFunctionReturn(0); 2748 } 2749 2750 /*@ 2751 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 2752 in faceVertices. The orientation is such that the face normal points out of the cell 2753 2754 Not collective 2755 2756 Input Parameters: 2757 + dm - The original mesh 2758 . cell - The cell mesh point 2759 . faceSize - The number of vertices on the face 2760 . face - The face vertices 2761 . numCorners - The number of vertices on the cell 2762 . indices - Local numbering of face vertices in cell cone 2763 - origVertices - Original face vertices 2764 2765 Output Parameters: 2766 + faceVertices - The face vertices properly oriented 2767 - posOriented - PETSC_TRUE if the face was oriented with outward normal 2768 2769 Level: developer 2770 2771 .seealso: DMPlexGetCone() 2772 @*/ 2773 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2774 { 2775 const PetscInt *cone = NULL; 2776 PetscInt coneSize, v, f, v2; 2777 PetscInt oppositeVertex = -1; 2778 PetscErrorCode ierr; 2779 2780 PetscFunctionBegin; 2781 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2782 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2783 for (v = 0, v2 = 0; v < coneSize; ++v) { 2784 PetscBool found = PETSC_FALSE; 2785 2786 for (f = 0; f < faceSize; ++f) { 2787 if (face[f] == cone[v]) { 2788 found = PETSC_TRUE; break; 2789 } 2790 } 2791 if (found) { 2792 indices[v2] = v; 2793 origVertices[v2] = cone[v]; 2794 ++v2; 2795 } else { 2796 oppositeVertex = v; 2797 } 2798 } 2799 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2800 PetscFunctionReturn(0); 2801 } 2802 2803 /* 2804 DMPlexInsertFace_Internal - Puts a face into the mesh 2805 2806 Not collective 2807 2808 Input Parameters: 2809 + dm - The DMPlex 2810 . numFaceVertex - The number of vertices in the face 2811 . faceVertices - The vertices in the face for dm 2812 . subfaceVertices - The vertices in the face for subdm 2813 . numCorners - The number of vertices in the cell 2814 . cell - A cell in dm containing the face 2815 . subcell - A cell in subdm containing the face 2816 . firstFace - First face in the mesh 2817 - newFacePoint - Next face in the mesh 2818 2819 Output Parameters: 2820 . newFacePoint - Contains next face point number on input, updated on output 2821 2822 Level: developer 2823 */ 2824 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) 2825 { 2826 MPI_Comm comm; 2827 DM_Plex *submesh = (DM_Plex*) subdm->data; 2828 const PetscInt *faces; 2829 PetscInt numFaces, coneSize; 2830 PetscErrorCode ierr; 2831 2832 PetscFunctionBegin; 2833 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2834 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2835 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2836 #if 0 2837 /* Cannot use this because support() has not been constructed yet */ 2838 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2839 #else 2840 { 2841 PetscInt f; 2842 2843 numFaces = 0; 2844 ierr = DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);CHKERRQ(ierr); 2845 for (f = firstFace; f < *newFacePoint; ++f) { 2846 PetscInt dof, off, d; 2847 2848 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2849 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2850 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2851 for (d = 0; d < dof; ++d) { 2852 const PetscInt p = submesh->cones[off+d]; 2853 PetscInt v; 2854 2855 for (v = 0; v < numFaceVertices; ++v) { 2856 if (subfaceVertices[v] == p) break; 2857 } 2858 if (v == numFaceVertices) break; 2859 } 2860 if (d == dof) { 2861 numFaces = 1; 2862 ((PetscInt*) faces)[0] = f; 2863 } 2864 } 2865 } 2866 #endif 2867 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2868 else if (numFaces == 1) { 2869 /* Add the other cell neighbor for this face */ 2870 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2871 } else { 2872 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2873 PetscBool posOriented; 2874 2875 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);CHKERRQ(ierr); 2876 origVertices = &orientedVertices[numFaceVertices]; 2877 indices = &orientedVertices[numFaceVertices*2]; 2878 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2879 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2880 /* TODO: I know that routine should return a permutation, not the indices */ 2881 for (v = 0; v < numFaceVertices; ++v) { 2882 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2883 for (ov = 0; ov < numFaceVertices; ++ov) { 2884 if (orientedVertices[ov] == vertex) { 2885 orientedSubVertices[ov] = subvertex; 2886 break; 2887 } 2888 } 2889 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2890 } 2891 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2892 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2893 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);CHKERRQ(ierr); 2894 ++(*newFacePoint); 2895 } 2896 #if 0 2897 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2898 #else 2899 ierr = DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);CHKERRQ(ierr); 2900 #endif 2901 PetscFunctionReturn(0); 2902 } 2903 2904 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2905 { 2906 MPI_Comm comm; 2907 DMLabel subpointMap; 2908 IS subvertexIS, subcellIS; 2909 const PetscInt *subVertices, *subCells; 2910 PetscInt numSubVertices, firstSubVertex, numSubCells; 2911 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2912 PetscInt vStart, vEnd, c, f; 2913 PetscErrorCode ierr; 2914 2915 PetscFunctionBegin; 2916 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2917 /* Create subpointMap which marks the submesh */ 2918 ierr = DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);CHKERRQ(ierr); 2919 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2920 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2921 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2922 /* Setup chart */ 2923 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2924 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2925 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2926 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2927 /* Set cone sizes */ 2928 firstSubVertex = numSubCells; 2929 firstSubFace = numSubCells+numSubVertices; 2930 newFacePoint = firstSubFace; 2931 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2932 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2933 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2934 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2935 for (c = 0; c < numSubCells; ++c) { 2936 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2937 } 2938 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2939 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2940 } 2941 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2942 /* Create face cones */ 2943 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2944 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2945 ierr = DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 2946 for (c = 0; c < numSubCells; ++c) { 2947 const PetscInt cell = subCells[c]; 2948 const PetscInt subcell = c; 2949 PetscInt *closure = NULL; 2950 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2951 2952 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2953 for (cl = 0; cl < closureSize*2; cl += 2) { 2954 const PetscInt point = closure[cl]; 2955 PetscInt subVertex; 2956 2957 if ((point >= vStart) && (point < vEnd)) { 2958 ++numCorners; 2959 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2960 if (subVertex >= 0) { 2961 closure[faceSize] = point; 2962 subface[faceSize] = firstSubVertex+subVertex; 2963 ++faceSize; 2964 } 2965 } 2966 } 2967 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2968 if (faceSize == nFV) { 2969 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2970 } 2971 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2972 } 2973 ierr = DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 2974 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2975 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2976 /* Build coordinates */ 2977 { 2978 PetscSection coordSection, subCoordSection; 2979 Vec coordinates, subCoordinates; 2980 PetscScalar *coords, *subCoords; 2981 PetscInt numComp, coordSize, v; 2982 const char *name; 2983 2984 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2985 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2986 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2987 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2988 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2989 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2990 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2991 for (v = 0; v < numSubVertices; ++v) { 2992 const PetscInt vertex = subVertices[v]; 2993 const PetscInt subvertex = firstSubVertex+v; 2994 PetscInt dof; 2995 2996 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2997 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2998 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2999 } 3000 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3001 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3002 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3003 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3004 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3005 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3006 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3007 if (coordSize) { 3008 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3009 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3010 for (v = 0; v < numSubVertices; ++v) { 3011 const PetscInt vertex = subVertices[v]; 3012 const PetscInt subvertex = firstSubVertex+v; 3013 PetscInt dof, off, sdof, soff, d; 3014 3015 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3016 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3017 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3018 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3019 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3020 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3021 } 3022 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3023 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3024 } 3025 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3026 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3027 } 3028 /* Cleanup */ 3029 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3030 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3031 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 3032 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 3033 PetscFunctionReturn(0); 3034 } 3035 3036 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 3037 { 3038 PetscInt subPoint; 3039 PetscErrorCode ierr; 3040 3041 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 3042 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 3043 } 3044 3045 static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPoints[], const PetscInt *subpoints[], const PetscInt firstSubPoint[], DM subdm) 3046 { 3047 PetscInt Nl, l, d; 3048 PetscErrorCode ierr; 3049 3050 PetscFunctionBegin; 3051 ierr = DMGetNumLabels(dm, &Nl);CHKERRQ(ierr); 3052 for (l = 0; l < Nl; ++l) { 3053 DMLabel label, newlabel; 3054 const char *lname; 3055 PetscBool isDepth, isDim, isCelltype; 3056 IS valueIS; 3057 const PetscInt *values; 3058 PetscInt Nv, v; 3059 3060 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr); 3061 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 3062 ierr = PetscStrcmp(lname, "dim", &isDim);CHKERRQ(ierr); 3063 ierr = PetscStrcmp(lname, "celltype", &isCelltype);CHKERRQ(ierr); 3064 if (isDepth || isDim || isCelltype) continue; 3065 ierr = DMCreateLabel(subdm, lname);CHKERRQ(ierr); 3066 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr); 3067 ierr = DMGetLabel(subdm, lname, &newlabel);CHKERRQ(ierr); 3068 ierr = DMLabelGetDefaultValue(label, &v);CHKERRQ(ierr); 3069 ierr = DMLabelSetDefaultValue(newlabel, v);CHKERRQ(ierr); 3070 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 3071 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr); 3072 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 3073 for (v = 0; v < Nv; ++v) { 3074 IS pointIS; 3075 const PetscInt *points; 3076 PetscInt Np, p; 3077 3078 ierr = DMLabelGetStratumIS(label, values[v], &pointIS);CHKERRQ(ierr); 3079 ierr = ISGetLocalSize(pointIS, &Np);CHKERRQ(ierr); 3080 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 3081 for (p = 0; p < Np; ++p) { 3082 const PetscInt point = points[p]; 3083 PetscInt subp; 3084 3085 ierr = DMPlexGetPointDepth(dm, point, &d);CHKERRQ(ierr); 3086 subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]); 3087 if (subp >= 0) {ierr = DMLabelSetValue(newlabel, subp, values[v]);CHKERRQ(ierr);} 3088 } 3089 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 3090 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 3091 } 3092 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 3093 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 3094 } 3095 PetscFunctionReturn(0); 3096 } 3097 3098 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool markedFaces, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 3099 { 3100 MPI_Comm comm; 3101 DMLabel subpointMap; 3102 IS *subpointIS; 3103 const PetscInt **subpoints; 3104 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 3105 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 3106 PetscMPIInt rank; 3107 PetscErrorCode ierr; 3108 3109 PetscFunctionBegin; 3110 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3111 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3112 /* Create subpointMap which marks the submesh */ 3113 ierr = DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);CHKERRQ(ierr); 3114 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3115 if (cellHeight) { 3116 if (isCohesive) {ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 3117 else {ierr = DMPlexMarkSubmesh_Interpolated(dm, label, value, markedFaces, subpointMap, subdm);CHKERRQ(ierr);} 3118 } else { 3119 DMLabel depth; 3120 IS pointIS; 3121 const PetscInt *points; 3122 PetscInt numPoints=0; 3123 3124 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 3125 ierr = DMLabelGetStratumIS(label, value, &pointIS);CHKERRQ(ierr); 3126 if (pointIS) { 3127 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 3128 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 3129 } 3130 for (p = 0; p < numPoints; ++p) { 3131 PetscInt *closure = NULL; 3132 PetscInt closureSize, c, pdim; 3133 3134 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 3135 for (c = 0; c < closureSize*2; c += 2) { 3136 ierr = DMLabelGetValue(depth, closure[c], &pdim);CHKERRQ(ierr); 3137 ierr = DMLabelSetValue(subpointMap, closure[c], pdim);CHKERRQ(ierr); 3138 } 3139 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 3140 } 3141 if (pointIS) {ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);} 3142 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 3143 } 3144 /* Setup chart */ 3145 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3146 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 3147 for (d = 0; d <= dim; ++d) { 3148 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 3149 totSubPoints += numSubPoints[d]; 3150 } 3151 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 3152 ierr = DMPlexSetVTKCellHeight(subdm, cellHeight);CHKERRQ(ierr); 3153 /* Set cone sizes */ 3154 firstSubPoint[dim] = 0; 3155 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 3156 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 3157 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 3158 for (d = 0; d <= dim; ++d) { 3159 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 3160 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3161 } 3162 /* We do not want this label automatically computed, instead we compute it here */ 3163 ierr = DMCreateLabel(subdm, "celltype");CHKERRQ(ierr); 3164 for (d = 0; d <= dim; ++d) { 3165 for (p = 0; p < numSubPoints[d]; ++p) { 3166 const PetscInt point = subpoints[d][p]; 3167 const PetscInt subpoint = firstSubPoint[d] + p; 3168 const PetscInt *cone; 3169 PetscInt coneSize, coneSizeNew, c, val; 3170 DMPolytopeType ct; 3171 3172 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 3173 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 3174 ierr = DMPlexGetCellType(dm, point, &ct);CHKERRQ(ierr); 3175 ierr = DMPlexSetCellType(subdm, subpoint, ct);CHKERRQ(ierr); 3176 if (cellHeight && (d == dim)) { 3177 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3178 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3179 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 3180 if (val >= 0) coneSizeNew++; 3181 } 3182 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 3183 ierr = DMPlexSetCellType(subdm, subpoint, DM_POLYTOPE_FV_GHOST);CHKERRQ(ierr); 3184 } 3185 } 3186 } 3187 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3188 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3189 /* Set cones */ 3190 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3191 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 3192 for (d = 0; d <= dim; ++d) { 3193 for (p = 0; p < numSubPoints[d]; ++p) { 3194 const PetscInt point = subpoints[d][p]; 3195 const PetscInt subpoint = firstSubPoint[d] + p; 3196 const PetscInt *cone, *ornt; 3197 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 3198 3199 if (d == dim-1) { 3200 const PetscInt *support, *cone, *ornt; 3201 PetscInt supportSize, coneSize, s, subc; 3202 3203 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 3204 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 3205 for (s = 0; s < supportSize; ++s) { 3206 PetscBool isHybrid; 3207 3208 ierr = DMPlexCellIsHybrid_Internal(dm, support[s], &isHybrid);CHKERRQ(ierr); 3209 if (!isHybrid) continue; 3210 ierr = PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);CHKERRQ(ierr); 3211 if (subc >= 0) { 3212 const PetscInt ccell = subpoints[d+1][subc]; 3213 3214 ierr = DMPlexGetCone(dm, ccell, &cone);CHKERRQ(ierr); 3215 ierr = DMPlexGetConeSize(dm, ccell, &coneSize);CHKERRQ(ierr); 3216 ierr = DMPlexGetConeOrientation(dm, ccell, &ornt);CHKERRQ(ierr); 3217 for (c = 0; c < coneSize; ++c) { 3218 if (cone[c] == point) { 3219 fornt = ornt[c]; 3220 break; 3221 } 3222 } 3223 break; 3224 } 3225 } 3226 } 3227 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 3228 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 3229 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3230 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 3231 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3232 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 3233 if (subc >= 0) { 3234 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 3235 orntNew[coneSizeNew] = ornt[c]; 3236 ++coneSizeNew; 3237 } 3238 } 3239 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 3240 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 3241 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 3242 if (fornt < 0) {ierr = DMPlexOrientPoint(subdm, subpoint, fornt);CHKERRQ(ierr);} 3243 } 3244 } 3245 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 3246 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3247 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3248 /* Build coordinates */ 3249 { 3250 PetscSection coordSection, subCoordSection; 3251 Vec coordinates, subCoordinates; 3252 PetscScalar *coords, *subCoords; 3253 PetscInt cdim, numComp, coordSize; 3254 const char *name; 3255 3256 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3257 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3258 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3259 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3260 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3261 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3262 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3263 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 3264 for (v = 0; v < numSubPoints[0]; ++v) { 3265 const PetscInt vertex = subpoints[0][v]; 3266 const PetscInt subvertex = firstSubPoint[0]+v; 3267 PetscInt dof; 3268 3269 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3270 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3271 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3272 } 3273 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3274 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3275 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3276 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3277 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3278 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3279 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3280 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3281 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3282 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3283 for (v = 0; v < numSubPoints[0]; ++v) { 3284 const PetscInt vertex = subpoints[0][v]; 3285 const PetscInt subvertex = firstSubPoint[0]+v; 3286 PetscInt dof, off, sdof, soff, d; 3287 3288 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3289 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3290 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3291 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3292 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3293 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3294 } 3295 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3296 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3297 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3298 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3299 } 3300 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3301 { 3302 PetscSF sfPoint, sfPointSub; 3303 IS subpIS; 3304 const PetscSFNode *remotePoints; 3305 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3306 const PetscInt *localPoints, *subpoints; 3307 PetscInt *slocalPoints; 3308 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3309 PetscMPIInt rank; 3310 3311 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 3312 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3313 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3314 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3315 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 3316 ierr = DMPlexGetSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 3317 if (subpIS) { 3318 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 3319 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 3320 } 3321 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3322 if (numRoots >= 0) { 3323 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3324 for (p = 0; p < pEnd-pStart; ++p) { 3325 newLocalPoints[p].rank = -2; 3326 newLocalPoints[p].index = -2; 3327 } 3328 /* Set subleaves */ 3329 for (l = 0; l < numLeaves; ++l) { 3330 const PetscInt point = localPoints[l]; 3331 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3332 3333 if (subpoint < 0) continue; 3334 newLocalPoints[point-pStart].rank = rank; 3335 newLocalPoints[point-pStart].index = subpoint; 3336 ++numSubleaves; 3337 } 3338 /* Must put in owned subpoints */ 3339 for (p = pStart; p < pEnd; ++p) { 3340 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3341 3342 if (subpoint < 0) { 3343 newOwners[p-pStart].rank = -3; 3344 newOwners[p-pStart].index = -3; 3345 } else { 3346 newOwners[p-pStart].rank = rank; 3347 newOwners[p-pStart].index = subpoint; 3348 } 3349 } 3350 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3351 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3352 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);CHKERRQ(ierr); 3353 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);CHKERRQ(ierr); 3354 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 3355 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 3356 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3357 const PetscInt point = localPoints[l]; 3358 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3359 3360 if (subpoint < 0) continue; 3361 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3362 slocalPoints[sl] = subpoint; 3363 sremotePoints[sl].rank = newLocalPoints[point].rank; 3364 sremotePoints[sl].index = newLocalPoints[point].index; 3365 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3366 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3367 ++sl; 3368 } 3369 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3370 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3371 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3372 } 3373 if (subpIS) { 3374 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 3375 } 3376 } 3377 /* Filter labels */ 3378 ierr = DMPlexFilterLabels_Internal(dm, numSubPoints, subpoints, firstSubPoint, subdm);CHKERRQ(ierr); 3379 /* Cleanup */ 3380 for (d = 0; d <= dim; ++d) { 3381 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3382 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 3383 } 3384 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 3385 PetscFunctionReturn(0); 3386 } 3387 3388 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM subdm) 3389 { 3390 PetscErrorCode ierr; 3391 3392 PetscFunctionBegin; 3393 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, markedFaces, PETSC_FALSE, 1, subdm);CHKERRQ(ierr); 3394 PetscFunctionReturn(0); 3395 } 3396 3397 /*@ 3398 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3399 3400 Input Parameters: 3401 + dm - The original mesh 3402 . vertexLabel - The DMLabel marking points contained in the surface 3403 . value - The label value to use 3404 - markedFaces - PETSC_TRUE if surface faces are marked in addition to vertices, PETSC_FALSE if only vertices are marked 3405 3406 Output Parameter: 3407 . subdm - The surface mesh 3408 3409 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3410 3411 Level: developer 3412 3413 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3414 @*/ 3415 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, PetscBool markedFaces, DM *subdm) 3416 { 3417 DMPlexInterpolatedFlag interpolated; 3418 PetscInt dim, cdim; 3419 PetscErrorCode ierr; 3420 3421 PetscFunctionBegin; 3422 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3423 PetscValidPointer(subdm, 5); 3424 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3425 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3426 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3427 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3428 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3429 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3430 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 3431 if (interpolated == DMPLEX_INTERPOLATED_PARTIAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 3432 if (interpolated) { 3433 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, markedFaces, *subdm);CHKERRQ(ierr); 3434 } else { 3435 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3436 } 3437 PetscFunctionReturn(0); 3438 } 3439 3440 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3441 { 3442 MPI_Comm comm; 3443 DMLabel subpointMap; 3444 IS subvertexIS; 3445 const PetscInt *subVertices; 3446 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3447 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3448 PetscInt c, f; 3449 PetscErrorCode ierr; 3450 3451 PetscFunctionBegin; 3452 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 3453 /* Create subpointMap which marks the submesh */ 3454 ierr = DMLabelCreate(PETSC_COMM_SELF, "subpoint_map", &subpointMap);CHKERRQ(ierr); 3455 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3456 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3457 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 3458 /* Setup chart */ 3459 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 3460 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 3461 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 3462 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 3463 /* Set cone sizes */ 3464 firstSubVertex = numSubCells; 3465 firstSubFace = numSubCells+numSubVertices; 3466 newFacePoint = firstSubFace; 3467 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 3468 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3469 for (c = 0; c < numSubCells; ++c) { 3470 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 3471 } 3472 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3473 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 3474 } 3475 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3476 /* Create face cones */ 3477 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3478 ierr = DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 3479 for (c = 0; c < numSubCells; ++c) { 3480 const PetscInt cell = subCells[c]; 3481 const PetscInt subcell = c; 3482 const PetscInt *cone, *cells; 3483 PetscBool isHybrid; 3484 PetscInt numCells, subVertex, p, v; 3485 3486 ierr = DMPlexCellIsHybrid_Internal(dm, cell, &isHybrid);CHKERRQ(ierr); 3487 if (!isHybrid) continue; 3488 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 3489 for (v = 0; v < nFV; ++v) { 3490 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 3491 subface[v] = firstSubVertex+subVertex; 3492 } 3493 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 3494 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 3495 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3496 /* Not true in parallel 3497 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3498 for (p = 0; p < numCells; ++p) { 3499 PetscInt negsubcell; 3500 PetscBool isHybrid; 3501 3502 ierr = DMPlexCellIsHybrid_Internal(dm, cells[p], &isHybrid);CHKERRQ(ierr); 3503 if (isHybrid) continue; 3504 /* I know this is a crap search */ 3505 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3506 if (subCells[negsubcell] == cells[p]) break; 3507 } 3508 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3509 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 3510 } 3511 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3512 ++newFacePoint; 3513 } 3514 ierr = DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 3515 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3516 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3517 /* Build coordinates */ 3518 { 3519 PetscSection coordSection, subCoordSection; 3520 Vec coordinates, subCoordinates; 3521 PetscScalar *coords, *subCoords; 3522 PetscInt cdim, numComp, coordSize, v; 3523 const char *name; 3524 3525 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3526 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3527 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3528 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3529 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3530 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3531 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3532 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3533 for (v = 0; v < numSubVertices; ++v) { 3534 const PetscInt vertex = subVertices[v]; 3535 const PetscInt subvertex = firstSubVertex+v; 3536 PetscInt dof; 3537 3538 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3539 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3540 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3541 } 3542 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3543 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3544 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3545 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3546 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3547 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3548 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3549 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3550 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3551 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3552 for (v = 0; v < numSubVertices; ++v) { 3553 const PetscInt vertex = subVertices[v]; 3554 const PetscInt subvertex = firstSubVertex+v; 3555 PetscInt dof, off, sdof, soff, d; 3556 3557 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3558 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3559 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3560 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3561 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3562 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3563 } 3564 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3565 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3566 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3567 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3568 } 3569 /* Build SF */ 3570 CHKMEMQ; 3571 { 3572 PetscSF sfPoint, sfPointSub; 3573 const PetscSFNode *remotePoints; 3574 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3575 const PetscInt *localPoints; 3576 PetscInt *slocalPoints; 3577 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3578 PetscMPIInt rank; 3579 3580 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 3581 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3582 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3583 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3584 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3585 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3586 if (numRoots >= 0) { 3587 /* Only vertices should be shared */ 3588 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3589 for (p = 0; p < pEnd-pStart; ++p) { 3590 newLocalPoints[p].rank = -2; 3591 newLocalPoints[p].index = -2; 3592 } 3593 /* Set subleaves */ 3594 for (l = 0; l < numLeaves; ++l) { 3595 const PetscInt point = localPoints[l]; 3596 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3597 3598 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3599 if (subPoint < 0) continue; 3600 newLocalPoints[point-pStart].rank = rank; 3601 newLocalPoints[point-pStart].index = subPoint; 3602 ++numSubLeaves; 3603 } 3604 /* Must put in owned subpoints */ 3605 for (p = pStart; p < pEnd; ++p) { 3606 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3607 3608 if (subPoint < 0) { 3609 newOwners[p-pStart].rank = -3; 3610 newOwners[p-pStart].index = -3; 3611 } else { 3612 newOwners[p-pStart].rank = rank; 3613 newOwners[p-pStart].index = subPoint; 3614 } 3615 } 3616 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3617 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3618 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);CHKERRQ(ierr); 3619 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints,MPI_REPLACE);CHKERRQ(ierr); 3620 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3621 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3622 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3623 const PetscInt point = localPoints[l]; 3624 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3625 3626 if (subPoint < 0) continue; 3627 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3628 slocalPoints[sl] = subPoint; 3629 sremotePoints[sl].rank = newLocalPoints[point].rank; 3630 sremotePoints[sl].index = newLocalPoints[point].index; 3631 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3632 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3633 ++sl; 3634 } 3635 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3636 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3637 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3638 } 3639 } 3640 CHKMEMQ; 3641 /* Cleanup */ 3642 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3643 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3644 ierr = PetscFree(subCells);CHKERRQ(ierr); 3645 PetscFunctionReturn(0); 3646 } 3647 3648 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3649 { 3650 DMLabel label = NULL; 3651 PetscErrorCode ierr; 3652 3653 PetscFunctionBegin; 3654 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3655 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_FALSE, PETSC_TRUE, 1, subdm);CHKERRQ(ierr); 3656 PetscFunctionReturn(0); 3657 } 3658 3659 /*@C 3660 DMPlexCreateCohesiveSubmesh - Extract from a mesh with cohesive cells the hypersurface defined by one face of the cells. Optionally, a Label can be given to restrict the cells. 3661 3662 Input Parameters: 3663 + dm - The original mesh 3664 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3665 . label - A label name, or NULL 3666 - value - A label value 3667 3668 Output Parameter: 3669 . subdm - The surface mesh 3670 3671 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3672 3673 Level: developer 3674 3675 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3676 @*/ 3677 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3678 { 3679 PetscInt dim, cdim, depth; 3680 PetscErrorCode ierr; 3681 3682 PetscFunctionBegin; 3683 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3684 PetscValidPointer(subdm, 5); 3685 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3686 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3687 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3688 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3689 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3690 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3691 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3692 if (depth == dim) { 3693 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3694 } else { 3695 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3696 } 3697 PetscFunctionReturn(0); 3698 } 3699 3700 /*@ 3701 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3702 3703 Input Parameters: 3704 + dm - The original mesh 3705 . cellLabel - The DMLabel marking cells contained in the new mesh 3706 - value - The label value to use 3707 3708 Output Parameter: 3709 . subdm - The new mesh 3710 3711 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3712 3713 Level: developer 3714 3715 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3716 @*/ 3717 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3718 { 3719 PetscInt dim; 3720 PetscErrorCode ierr; 3721 3722 PetscFunctionBegin; 3723 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3724 PetscValidPointer(subdm, 4); 3725 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3726 ierr = DMCreate(PetscObjectComm((PetscObject) dm), subdm);CHKERRQ(ierr); 3727 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3728 ierr = DMSetDimension(*subdm, dim);CHKERRQ(ierr); 3729 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3730 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, PETSC_FALSE, 0, *subdm);CHKERRQ(ierr); 3731 PetscFunctionReturn(0); 3732 } 3733 3734 /*@ 3735 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3736 3737 Input Parameter: 3738 . dm - The submesh DM 3739 3740 Output Parameter: 3741 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3742 3743 Level: developer 3744 3745 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS() 3746 @*/ 3747 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3748 { 3749 PetscFunctionBegin; 3750 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3751 PetscValidPointer(subpointMap, 2); 3752 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3753 PetscFunctionReturn(0); 3754 } 3755 3756 /*@ 3757 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3758 3759 Input Parameters: 3760 + dm - The submesh DM 3761 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3762 3763 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3764 3765 Level: developer 3766 3767 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointIS() 3768 @*/ 3769 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3770 { 3771 DM_Plex *mesh = (DM_Plex *) dm->data; 3772 DMLabel tmp; 3773 PetscErrorCode ierr; 3774 3775 PetscFunctionBegin; 3776 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3777 tmp = mesh->subpointMap; 3778 mesh->subpointMap = subpointMap; 3779 ierr = PetscObjectReference((PetscObject) mesh->subpointMap);CHKERRQ(ierr); 3780 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3781 PetscFunctionReturn(0); 3782 } 3783 3784 static PetscErrorCode DMPlexCreateSubpointIS_Internal(DM dm, IS *subpointIS) 3785 { 3786 DMLabel spmap; 3787 PetscInt depth, d; 3788 PetscErrorCode ierr; 3789 3790 PetscFunctionBegin; 3791 ierr = DMPlexGetSubpointMap(dm, &spmap);CHKERRQ(ierr); 3792 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3793 if (spmap && depth >= 0) { 3794 DM_Plex *mesh = (DM_Plex *) dm->data; 3795 PetscInt *points, *depths; 3796 PetscInt pStart, pEnd, p, off; 3797 3798 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3799 if (pStart) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3800 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3801 ierr = DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);CHKERRQ(ierr); 3802 depths[0] = depth; 3803 depths[1] = 0; 3804 for (d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3805 for (d = 0, off = 0; d <= depth; ++d) { 3806 const PetscInt dep = depths[d]; 3807 PetscInt depStart, depEnd, n; 3808 3809 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3810 ierr = DMLabelGetStratumSize(spmap, dep, &n);CHKERRQ(ierr); 3811 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3812 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); 3813 } else { 3814 if (!n) { 3815 if (d == 0) { 3816 /* Missing cells */ 3817 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3818 } else { 3819 /* Missing faces */ 3820 for (p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3821 } 3822 } 3823 } 3824 if (n) { 3825 IS is; 3826 const PetscInt *opoints; 3827 3828 ierr = DMLabelGetStratumIS(spmap, dep, &is);CHKERRQ(ierr); 3829 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3830 for (p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3831 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3832 ierr = ISDestroy(&is);CHKERRQ(ierr); 3833 } 3834 } 3835 ierr = DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);CHKERRQ(ierr); 3836 if (off != pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3837 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3838 ierr = PetscObjectStateGet((PetscObject) spmap, &mesh->subpointState);CHKERRQ(ierr); 3839 } 3840 PetscFunctionReturn(0); 3841 } 3842 3843 /*@ 3844 DMPlexGetSubpointIS - Returns an IS covering the entire subdm chart with the original points as data 3845 3846 Input Parameter: 3847 . dm - The submesh DM 3848 3849 Output Parameter: 3850 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3851 3852 Note: This IS is guaranteed to be sorted by the construction of the submesh 3853 3854 Level: developer 3855 3856 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3857 @*/ 3858 PetscErrorCode DMPlexGetSubpointIS(DM dm, IS *subpointIS) 3859 { 3860 DM_Plex *mesh = (DM_Plex *) dm->data; 3861 DMLabel spmap; 3862 PetscObjectState state; 3863 PetscErrorCode ierr; 3864 3865 PetscFunctionBegin; 3866 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3867 PetscValidPointer(subpointIS, 2); 3868 ierr = DMPlexGetSubpointMap(dm, &spmap);CHKERRQ(ierr); 3869 ierr = PetscObjectStateGet((PetscObject) spmap, &state);CHKERRQ(ierr); 3870 if (state != mesh->subpointState || !mesh->subpointIS) {ierr = DMPlexCreateSubpointIS_Internal(dm, &mesh->subpointIS);CHKERRQ(ierr);} 3871 *subpointIS = mesh->subpointIS; 3872 PetscFunctionReturn(0); 3873 } 3874 3875 /*@ 3876 DMGetEnclosureRelation - Get the relationship between dmA and dmB 3877 3878 Input Parameters: 3879 + dmA - The first DM 3880 - dmB - The second DM 3881 3882 Output Parameter: 3883 . rel - The relation of dmA to dmB 3884 3885 Level: intermediate 3886 3887 .seealso: DMGetEnclosurePoint() 3888 @*/ 3889 PetscErrorCode DMGetEnclosureRelation(DM dmA, DM dmB, DMEnclosureType *rel) 3890 { 3891 DM plexA, plexB, sdm; 3892 DMLabel spmap; 3893 PetscInt pStartA, pEndA, pStartB, pEndB, NpA, NpB; 3894 PetscErrorCode ierr; 3895 3896 PetscFunctionBegin; 3897 PetscValidPointer(rel, 3); 3898 *rel = DM_ENC_NONE; 3899 if (!dmA || !dmB) PetscFunctionReturn(0); 3900 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 3901 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 3902 if (dmA == dmB) {*rel = DM_ENC_EQUALITY; PetscFunctionReturn(0);} 3903 ierr = DMConvert(dmA, DMPLEX, &plexA);CHKERRQ(ierr); 3904 ierr = DMConvert(dmB, DMPLEX, &plexB);CHKERRQ(ierr); 3905 ierr = DMPlexGetChart(plexA, &pStartA, &pEndA);CHKERRQ(ierr); 3906 ierr = DMPlexGetChart(plexB, &pStartB, &pEndB);CHKERRQ(ierr); 3907 /* Assumption 1: subDMs have smaller charts than the DMs that they originate from 3908 - The degenerate case of a subdomain which includes all of the domain on some process can be treated as equality */ 3909 if ((pStartA == pStartB) && (pEndA == pEndB)) { 3910 *rel = DM_ENC_EQUALITY; 3911 goto end; 3912 } 3913 NpA = pEndA - pStartA; 3914 NpB = pEndB - pStartB; 3915 if (NpA == NpB) goto end; 3916 sdm = NpA > NpB ? plexB : plexA; /* The other is the original, enclosing dm */ 3917 ierr = DMPlexGetSubpointMap(sdm, &spmap);CHKERRQ(ierr); 3918 if (!spmap) goto end; 3919 /* TODO Check the space mapped to by subpointMap is same size as dm */ 3920 if (NpA > NpB) { 3921 *rel = DM_ENC_SUPERMESH; 3922 } else { 3923 *rel = DM_ENC_SUBMESH; 3924 } 3925 end: 3926 ierr = DMDestroy(&plexA);CHKERRQ(ierr); 3927 ierr = DMDestroy(&plexB);CHKERRQ(ierr); 3928 PetscFunctionReturn(0); 3929 } 3930 3931 /*@ 3932 DMGetEnclosurePoint - Get the point pA in dmA which corresponds to the point pB in dmB 3933 3934 Input Parameters: 3935 + dmA - The first DM 3936 . dmB - The second DM 3937 . etype - The type of enclosure relation that dmA has to dmB 3938 - pB - A point of dmB 3939 3940 Output Parameter: 3941 . pA - The corresponding point of dmA 3942 3943 Level: intermediate 3944 3945 .seealso: DMGetEnclosureRelation() 3946 @*/ 3947 PetscErrorCode DMGetEnclosurePoint(DM dmA, DM dmB, DMEnclosureType etype, PetscInt pB, PetscInt *pA) 3948 { 3949 DM sdm; 3950 IS subpointIS; 3951 const PetscInt *subpoints; 3952 PetscInt numSubpoints; 3953 PetscErrorCode ierr; 3954 3955 PetscFunctionBegin; 3956 /* TODO Cache the IS, making it look like an index */ 3957 switch (etype) { 3958 case DM_ENC_SUPERMESH: 3959 sdm = dmB; 3960 ierr = DMPlexGetSubpointIS(sdm, &subpointIS);CHKERRQ(ierr); 3961 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 3962 *pA = subpoints[pB]; 3963 ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr); 3964 break; 3965 case DM_ENC_SUBMESH: 3966 sdm = dmA; 3967 ierr = DMPlexGetSubpointIS(sdm, &subpointIS);CHKERRQ(ierr); 3968 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 3969 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 3970 ierr = PetscFindInt(pB, numSubpoints, subpoints, pA);CHKERRQ(ierr); 3971 if (*pA < 0) { 3972 ierr = DMViewFromOptions(dmA, NULL, "-dm_enc_A_view");CHKERRQ(ierr); 3973 ierr = DMViewFromOptions(dmB, NULL, "-dm_enc_B_view");CHKERRQ(ierr); 3974 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", pB); 3975 } 3976 ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr); 3977 break; 3978 case DM_ENC_EQUALITY: 3979 case DM_ENC_NONE: 3980 *pA = pB;break; 3981 case DM_ENC_UNKNOWN: 3982 { 3983 DMEnclosureType enc; 3984 3985 ierr = DMGetEnclosureRelation(dmA, dmB, &enc);CHKERRQ(ierr); 3986 ierr = DMGetEnclosurePoint(dmA, dmB, enc, pB, pA);CHKERRQ(ierr); 3987 } 3988 break; 3989 default: SETERRQ1(PetscObjectComm((PetscObject) dmA), PETSC_ERR_ARG_OUTOFRANGE, "Invalid enclosure type %d", (int) etype); 3990 } 3991 PetscFunctionReturn(0); 3992 } 3993