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