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