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