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