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