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