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