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