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