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