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