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