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