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