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