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