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