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, &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, 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 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point; 2125 } 2126 for (f = 0; f < numFaces; ++f) { 2127 const PetscInt face = star[f]; 2128 PetscInt *closure = NULL; 2129 PetscInt closureSize, c; 2130 PetscInt faceLoc; 2131 2132 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 2133 if (faceLoc == dim-1) continue; 2134 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 2135 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2136 for (c = 0; c < closureSize*2; c += 2) { 2137 const PetscInt point = closure[c]; 2138 PetscInt vertexLoc; 2139 2140 if ((point >= pStart[0]) && (point < pEnd[0])) { 2141 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2142 if (vertexLoc != value) break; 2143 } 2144 } 2145 if (c == closureSize*2) { 2146 const PetscInt *support; 2147 PetscInt supportSize, s; 2148 2149 for (c = 0; c < closureSize*2; c += 2) { 2150 const PetscInt point = closure[c]; 2151 2152 for (d = 0; d < dim; ++d) { 2153 if ((point >= pStart[d]) && (point < pEnd[d])) { 2154 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2155 break; 2156 } 2157 } 2158 } 2159 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2160 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2161 for (s = 0; s < supportSize; ++s) { 2162 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2163 } 2164 } 2165 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2166 } 2167 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2168 } 2169 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr);} 2170 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2171 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 2172 PetscFunctionReturn(0); 2173 } 2174 2175 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 2176 { 2177 DMLabel label = NULL; 2178 const PetscInt *cone; 2179 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1; 2180 PetscErrorCode ierr; 2181 2182 PetscFunctionBegin; 2183 *numFaces = 0; 2184 *nFV = 0; 2185 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 2186 *subCells = NULL; 2187 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2188 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2189 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2190 if (cMax < 0) PetscFunctionReturn(0); 2191 if (label) { 2192 for (c = cMax; c < cEnd; ++c) { 2193 PetscInt val; 2194 2195 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2196 if (val == value) { 2197 ++(*numFaces); 2198 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2199 } 2200 } 2201 } else { 2202 *numFaces = cEnd - cMax; 2203 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 2204 } 2205 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 2206 if (!(*numFaces)) PetscFunctionReturn(0); 2207 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 2208 for (c = cMax; c < cEnd; ++c) { 2209 const PetscInt *cells; 2210 PetscInt numCells; 2211 2212 if (label) { 2213 PetscInt val; 2214 2215 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2216 if (val != value) continue; 2217 } 2218 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2219 for (p = 0; p < *nFV; ++p) { 2220 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 2221 } 2222 /* Negative face */ 2223 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2224 /* Not true in parallel 2225 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2226 for (p = 0; p < numCells; ++p) { 2227 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 2228 (*subCells)[subc++] = cells[p]; 2229 } 2230 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2231 /* Positive face is not included */ 2232 } 2233 PetscFunctionReturn(0); 2234 } 2235 2236 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 2237 { 2238 PetscInt *pStart, *pEnd; 2239 PetscInt dim, cMax, cEnd, c, d; 2240 PetscErrorCode ierr; 2241 2242 PetscFunctionBegin; 2243 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2244 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2245 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2246 if (cMax < 0) PetscFunctionReturn(0); 2247 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 2248 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 2249 for (c = cMax; c < cEnd; ++c) { 2250 const PetscInt *cone; 2251 PetscInt *closure = NULL; 2252 PetscInt fconeSize, coneSize, closureSize, cl, val; 2253 2254 if (label) { 2255 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2256 if (val != value) continue; 2257 } 2258 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2259 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2260 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 2261 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2262 /* Negative face */ 2263 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2264 for (cl = 0; cl < closureSize*2; cl += 2) { 2265 const PetscInt point = closure[cl]; 2266 2267 for (d = 0; d <= dim; ++d) { 2268 if ((point >= pStart[d]) && (point < pEnd[d])) { 2269 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2270 break; 2271 } 2272 } 2273 } 2274 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2275 /* Cells -- positive face is not included */ 2276 for (cl = 0; cl < 1; ++cl) { 2277 const PetscInt *support; 2278 PetscInt supportSize, s; 2279 2280 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 2281 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2282 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 2283 for (s = 0; s < supportSize; ++s) { 2284 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2285 } 2286 } 2287 } 2288 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2289 PetscFunctionReturn(0); 2290 } 2291 2292 static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2293 { 2294 MPI_Comm comm; 2295 PetscBool posOrient = PETSC_FALSE; 2296 const PetscInt debug = 0; 2297 PetscInt cellDim, faceSize, f; 2298 PetscErrorCode ierr; 2299 2300 PetscFunctionBegin; 2301 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2302 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 2303 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 2304 2305 if (cellDim == 1 && numCorners == 2) { 2306 /* Triangle */ 2307 faceSize = numCorners-1; 2308 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2309 } else if (cellDim == 2 && numCorners == 3) { 2310 /* Triangle */ 2311 faceSize = numCorners-1; 2312 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2313 } else if (cellDim == 3 && numCorners == 4) { 2314 /* Tetrahedron */ 2315 faceSize = numCorners-1; 2316 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2317 } else if (cellDim == 1 && numCorners == 3) { 2318 /* Quadratic line */ 2319 faceSize = 1; 2320 posOrient = PETSC_TRUE; 2321 } else if (cellDim == 2 && numCorners == 4) { 2322 /* Quads */ 2323 faceSize = 2; 2324 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2325 posOrient = PETSC_TRUE; 2326 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2327 posOrient = PETSC_TRUE; 2328 } else { 2329 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2330 posOrient = PETSC_FALSE; 2331 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2332 } 2333 } else if (cellDim == 2 && numCorners == 6) { 2334 /* Quadratic triangle (I hate this) */ 2335 /* Edges are determined by the first 2 vertices (corners of edges) */ 2336 const PetscInt faceSizeTri = 3; 2337 PetscInt sortedIndices[3], i, iFace; 2338 PetscBool found = PETSC_FALSE; 2339 PetscInt faceVerticesTriSorted[9] = { 2340 0, 3, 4, /* bottom */ 2341 1, 4, 5, /* right */ 2342 2, 3, 5, /* left */ 2343 }; 2344 PetscInt faceVerticesTri[9] = { 2345 0, 3, 4, /* bottom */ 2346 1, 4, 5, /* right */ 2347 2, 5, 3, /* left */ 2348 }; 2349 2350 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2351 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2352 for (iFace = 0; iFace < 3; ++iFace) { 2353 const PetscInt ii = iFace*faceSizeTri; 2354 PetscInt fVertex, cVertex; 2355 2356 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2357 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2358 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2359 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2360 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2361 faceVertices[fVertex] = origVertices[cVertex]; 2362 break; 2363 } 2364 } 2365 } 2366 found = PETSC_TRUE; 2367 break; 2368 } 2369 } 2370 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2371 if (posOriented) *posOriented = PETSC_TRUE; 2372 PetscFunctionReturn(0); 2373 } else if (cellDim == 2 && numCorners == 9) { 2374 /* Quadratic quad (I hate this) */ 2375 /* Edges are determined by the first 2 vertices (corners of edges) */ 2376 const PetscInt faceSizeQuad = 3; 2377 PetscInt sortedIndices[3], i, iFace; 2378 PetscBool found = PETSC_FALSE; 2379 PetscInt faceVerticesQuadSorted[12] = { 2380 0, 1, 4, /* bottom */ 2381 1, 2, 5, /* right */ 2382 2, 3, 6, /* top */ 2383 0, 3, 7, /* left */ 2384 }; 2385 PetscInt faceVerticesQuad[12] = { 2386 0, 1, 4, /* bottom */ 2387 1, 2, 5, /* right */ 2388 2, 3, 6, /* top */ 2389 3, 0, 7, /* left */ 2390 }; 2391 2392 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2393 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2394 for (iFace = 0; iFace < 4; ++iFace) { 2395 const PetscInt ii = iFace*faceSizeQuad; 2396 PetscInt fVertex, cVertex; 2397 2398 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2399 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2400 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2401 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2402 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2403 faceVertices[fVertex] = origVertices[cVertex]; 2404 break; 2405 } 2406 } 2407 } 2408 found = PETSC_TRUE; 2409 break; 2410 } 2411 } 2412 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2413 if (posOriented) *posOriented = PETSC_TRUE; 2414 PetscFunctionReturn(0); 2415 } else if (cellDim == 3 && numCorners == 8) { 2416 /* Hexes 2417 A hex is two oriented quads with the normal of the first 2418 pointing up at the second. 2419 2420 7---6 2421 /| /| 2422 4---5 | 2423 | 1-|-2 2424 |/ |/ 2425 0---3 2426 2427 Faces are determined by the first 4 vertices (corners of faces) */ 2428 const PetscInt faceSizeHex = 4; 2429 PetscInt sortedIndices[4], i, iFace; 2430 PetscBool found = PETSC_FALSE; 2431 PetscInt faceVerticesHexSorted[24] = { 2432 0, 1, 2, 3, /* bottom */ 2433 4, 5, 6, 7, /* top */ 2434 0, 3, 4, 5, /* front */ 2435 2, 3, 5, 6, /* right */ 2436 1, 2, 6, 7, /* back */ 2437 0, 1, 4, 7, /* left */ 2438 }; 2439 PetscInt faceVerticesHex[24] = { 2440 1, 2, 3, 0, /* bottom */ 2441 4, 5, 6, 7, /* top */ 2442 0, 3, 5, 4, /* front */ 2443 3, 2, 6, 5, /* right */ 2444 2, 1, 7, 6, /* back */ 2445 1, 0, 4, 7, /* left */ 2446 }; 2447 2448 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2449 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2450 for (iFace = 0; iFace < 6; ++iFace) { 2451 const PetscInt ii = iFace*faceSizeHex; 2452 PetscInt fVertex, cVertex; 2453 2454 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2455 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2456 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2457 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2458 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2459 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2460 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2461 faceVertices[fVertex] = origVertices[cVertex]; 2462 break; 2463 } 2464 } 2465 } 2466 found = PETSC_TRUE; 2467 break; 2468 } 2469 } 2470 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2471 if (posOriented) *posOriented = PETSC_TRUE; 2472 PetscFunctionReturn(0); 2473 } else if (cellDim == 3 && numCorners == 10) { 2474 /* Quadratic tet */ 2475 /* Faces are determined by the first 3 vertices (corners of faces) */ 2476 const PetscInt faceSizeTet = 6; 2477 PetscInt sortedIndices[6], i, iFace; 2478 PetscBool found = PETSC_FALSE; 2479 PetscInt faceVerticesTetSorted[24] = { 2480 0, 1, 2, 6, 7, 8, /* bottom */ 2481 0, 3, 4, 6, 7, 9, /* front */ 2482 1, 4, 5, 7, 8, 9, /* right */ 2483 2, 3, 5, 6, 8, 9, /* left */ 2484 }; 2485 PetscInt faceVerticesTet[24] = { 2486 0, 1, 2, 6, 7, 8, /* bottom */ 2487 0, 4, 3, 6, 7, 9, /* front */ 2488 1, 5, 4, 7, 8, 9, /* right */ 2489 2, 3, 5, 8, 6, 9, /* left */ 2490 }; 2491 2492 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2493 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2494 for (iFace=0; iFace < 4; ++iFace) { 2495 const PetscInt ii = iFace*faceSizeTet; 2496 PetscInt fVertex, cVertex; 2497 2498 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2499 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2500 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2501 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2502 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2503 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2504 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2505 faceVertices[fVertex] = origVertices[cVertex]; 2506 break; 2507 } 2508 } 2509 } 2510 found = PETSC_TRUE; 2511 break; 2512 } 2513 } 2514 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2515 if (posOriented) *posOriented = PETSC_TRUE; 2516 PetscFunctionReturn(0); 2517 } else if (cellDim == 3 && numCorners == 27) { 2518 /* Quadratic hexes (I hate this) 2519 A hex is two oriented quads with the normal of the first 2520 pointing up at the second. 2521 2522 7---6 2523 /| /| 2524 4---5 | 2525 | 3-|-2 2526 |/ |/ 2527 0---1 2528 2529 Faces are determined by the first 4 vertices (corners of faces) */ 2530 const PetscInt faceSizeQuadHex = 9; 2531 PetscInt sortedIndices[9], i, iFace; 2532 PetscBool found = PETSC_FALSE; 2533 PetscInt faceVerticesQuadHexSorted[54] = { 2534 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2535 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2536 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2537 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2538 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2539 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2540 }; 2541 PetscInt faceVerticesQuadHex[54] = { 2542 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2543 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2544 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2545 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2546 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2547 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2548 }; 2549 2550 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2551 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2552 for (iFace = 0; iFace < 6; ++iFace) { 2553 const PetscInt ii = iFace*faceSizeQuadHex; 2554 PetscInt fVertex, cVertex; 2555 2556 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2557 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2558 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2559 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2560 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2561 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2562 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2563 faceVertices[fVertex] = origVertices[cVertex]; 2564 break; 2565 } 2566 } 2567 } 2568 found = PETSC_TRUE; 2569 break; 2570 } 2571 } 2572 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2573 if (posOriented) *posOriented = PETSC_TRUE; 2574 PetscFunctionReturn(0); 2575 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2576 if (!posOrient) { 2577 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2578 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2579 } else { 2580 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2581 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2582 } 2583 if (posOriented) *posOriented = posOrient; 2584 PetscFunctionReturn(0); 2585 } 2586 2587 /*@ 2588 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 2589 in faceVertices. The orientation is such that the face normal points out of the cell 2590 2591 Not collective 2592 2593 Input Parameters: 2594 + dm - The original mesh 2595 . cell - The cell mesh point 2596 . faceSize - The number of vertices on the face 2597 . face - The face vertices 2598 . numCorners - The number of vertices on the cell 2599 . indices - Local numbering of face vertices in cell cone 2600 - origVertices - Original face vertices 2601 2602 Output Parameter: 2603 + faceVertices - The face vertices properly oriented 2604 - posOriented - PETSC_TRUE if the face was oriented with outward normal 2605 2606 Level: developer 2607 2608 .seealso: DMPlexGetCone() 2609 @*/ 2610 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2611 { 2612 const PetscInt *cone = NULL; 2613 PetscInt coneSize, v, f, v2; 2614 PetscInt oppositeVertex = -1; 2615 PetscErrorCode ierr; 2616 2617 PetscFunctionBegin; 2618 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2619 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2620 for (v = 0, v2 = 0; v < coneSize; ++v) { 2621 PetscBool found = PETSC_FALSE; 2622 2623 for (f = 0; f < faceSize; ++f) { 2624 if (face[f] == cone[v]) { 2625 found = PETSC_TRUE; break; 2626 } 2627 } 2628 if (found) { 2629 indices[v2] = v; 2630 origVertices[v2] = cone[v]; 2631 ++v2; 2632 } else { 2633 oppositeVertex = v; 2634 } 2635 } 2636 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 /* 2641 DMPlexInsertFace_Internal - Puts a face into the mesh 2642 2643 Not collective 2644 2645 Input Parameters: 2646 + dm - The DMPlex 2647 . numFaceVertex - The number of vertices in the face 2648 . faceVertices - The vertices in the face for dm 2649 . subfaceVertices - The vertices in the face for subdm 2650 . numCorners - The number of vertices in the cell 2651 . cell - A cell in dm containing the face 2652 . subcell - A cell in subdm containing the face 2653 . firstFace - First face in the mesh 2654 - newFacePoint - Next face in the mesh 2655 2656 Output Parameters: 2657 . newFacePoint - Contains next face point number on input, updated on output 2658 2659 Level: developer 2660 */ 2661 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) 2662 { 2663 MPI_Comm comm; 2664 DM_Plex *submesh = (DM_Plex*) subdm->data; 2665 const PetscInt *faces; 2666 PetscInt numFaces, coneSize; 2667 PetscErrorCode ierr; 2668 2669 PetscFunctionBegin; 2670 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2671 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2672 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2673 #if 0 2674 /* Cannot use this because support() has not been constructed yet */ 2675 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2676 #else 2677 { 2678 PetscInt f; 2679 2680 numFaces = 0; 2681 ierr = DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);CHKERRQ(ierr); 2682 for (f = firstFace; f < *newFacePoint; ++f) { 2683 PetscInt dof, off, d; 2684 2685 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2686 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2687 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2688 for (d = 0; d < dof; ++d) { 2689 const PetscInt p = submesh->cones[off+d]; 2690 PetscInt v; 2691 2692 for (v = 0; v < numFaceVertices; ++v) { 2693 if (subfaceVertices[v] == p) break; 2694 } 2695 if (v == numFaceVertices) break; 2696 } 2697 if (d == dof) { 2698 numFaces = 1; 2699 ((PetscInt*) faces)[0] = f; 2700 } 2701 } 2702 } 2703 #endif 2704 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2705 else if (numFaces == 1) { 2706 /* Add the other cell neighbor for this face */ 2707 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2708 } else { 2709 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2710 PetscBool posOriented; 2711 2712 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);CHKERRQ(ierr); 2713 origVertices = &orientedVertices[numFaceVertices]; 2714 indices = &orientedVertices[numFaceVertices*2]; 2715 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2716 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2717 /* TODO: I know that routine should return a permutation, not the indices */ 2718 for (v = 0; v < numFaceVertices; ++v) { 2719 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2720 for (ov = 0; ov < numFaceVertices; ++ov) { 2721 if (orientedVertices[ov] == vertex) { 2722 orientedSubVertices[ov] = subvertex; 2723 break; 2724 } 2725 } 2726 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2727 } 2728 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2729 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2730 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);CHKERRQ(ierr); 2731 ++(*newFacePoint); 2732 } 2733 #if 0 2734 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2735 #else 2736 ierr = DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);CHKERRQ(ierr); 2737 #endif 2738 PetscFunctionReturn(0); 2739 } 2740 2741 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2742 { 2743 MPI_Comm comm; 2744 DMLabel subpointMap; 2745 IS subvertexIS, subcellIS; 2746 const PetscInt *subVertices, *subCells; 2747 PetscInt numSubVertices, firstSubVertex, numSubCells; 2748 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2749 PetscInt vStart, vEnd, c, f; 2750 PetscErrorCode ierr; 2751 2752 PetscFunctionBegin; 2753 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2754 /* Create subpointMap which marks the submesh */ 2755 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2756 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2757 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2758 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2759 /* Setup chart */ 2760 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2761 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2762 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2763 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2764 /* Set cone sizes */ 2765 firstSubVertex = numSubCells; 2766 firstSubFace = numSubCells+numSubVertices; 2767 newFacePoint = firstSubFace; 2768 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2769 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2770 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2771 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2772 for (c = 0; c < numSubCells; ++c) { 2773 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2774 } 2775 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2776 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2777 } 2778 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2779 /* Create face cones */ 2780 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2781 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2782 ierr = DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 2783 for (c = 0; c < numSubCells; ++c) { 2784 const PetscInt cell = subCells[c]; 2785 const PetscInt subcell = c; 2786 PetscInt *closure = NULL; 2787 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2788 2789 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2790 for (cl = 0; cl < closureSize*2; cl += 2) { 2791 const PetscInt point = closure[cl]; 2792 PetscInt subVertex; 2793 2794 if ((point >= vStart) && (point < vEnd)) { 2795 ++numCorners; 2796 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2797 if (subVertex >= 0) { 2798 closure[faceSize] = point; 2799 subface[faceSize] = firstSubVertex+subVertex; 2800 ++faceSize; 2801 } 2802 } 2803 } 2804 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2805 if (faceSize == nFV) { 2806 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2807 } 2808 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2809 } 2810 ierr = DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 2811 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2812 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2813 /* Build coordinates */ 2814 { 2815 PetscSection coordSection, subCoordSection; 2816 Vec coordinates, subCoordinates; 2817 PetscScalar *coords, *subCoords; 2818 PetscInt numComp, coordSize, v; 2819 const char *name; 2820 2821 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2822 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2823 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2824 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2825 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2826 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2827 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2828 for (v = 0; v < numSubVertices; ++v) { 2829 const PetscInt vertex = subVertices[v]; 2830 const PetscInt subvertex = firstSubVertex+v; 2831 PetscInt dof; 2832 2833 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2834 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2835 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2836 } 2837 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2838 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2839 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 2840 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 2841 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 2842 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2843 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2844 if (coordSize) { 2845 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2846 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2847 for (v = 0; v < numSubVertices; ++v) { 2848 const PetscInt vertex = subVertices[v]; 2849 const PetscInt subvertex = firstSubVertex+v; 2850 PetscInt dof, off, sdof, soff, d; 2851 2852 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2853 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2854 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2855 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2856 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2857 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2858 } 2859 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2860 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2861 } 2862 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2863 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2864 } 2865 /* Cleanup */ 2866 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2867 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2868 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2869 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 2870 PetscFunctionReturn(0); 2871 } 2872 2873 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 2874 { 2875 PetscInt subPoint; 2876 PetscErrorCode ierr; 2877 2878 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 2879 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 2880 } 2881 2882 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 2883 { 2884 MPI_Comm comm; 2885 DMLabel subpointMap; 2886 IS *subpointIS; 2887 const PetscInt **subpoints; 2888 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 2889 PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v; 2890 PetscMPIInt rank; 2891 PetscErrorCode ierr; 2892 2893 PetscFunctionBegin; 2894 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2895 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2896 /* Create subpointMap which marks the submesh */ 2897 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2898 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2899 if (cellHeight) { 2900 if (isCohesive) {ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2901 else {ierr = DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2902 } else { 2903 DMLabel depth; 2904 IS pointIS; 2905 const PetscInt *points; 2906 PetscInt numPoints; 2907 2908 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 2909 ierr = DMLabelGetStratumSize(label, value, &numPoints);CHKERRQ(ierr); 2910 ierr = DMLabelGetStratumIS(label, value, &pointIS);CHKERRQ(ierr); 2911 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2912 for (p = 0; p < numPoints; ++p) { 2913 PetscInt *closure = NULL; 2914 PetscInt closureSize, c, pdim; 2915 2916 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2917 for (c = 0; c < closureSize*2; c += 2) { 2918 ierr = DMLabelGetValue(depth, closure[c], &pdim);CHKERRQ(ierr); 2919 ierr = DMLabelSetValue(subpointMap, closure[c], pdim);CHKERRQ(ierr); 2920 } 2921 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2922 } 2923 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2924 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2925 } 2926 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2927 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2928 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2929 cMax = (cMax < 0) ? cEnd : cMax; 2930 /* Setup chart */ 2931 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2932 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2933 for (d = 0; d <= dim; ++d) { 2934 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2935 totSubPoints += numSubPoints[d]; 2936 } 2937 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2938 ierr = DMPlexSetVTKCellHeight(subdm, cellHeight);CHKERRQ(ierr); 2939 /* Set cone sizes */ 2940 firstSubPoint[dim] = 0; 2941 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2942 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2943 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2944 for (d = 0; d <= dim; ++d) { 2945 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2946 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2947 } 2948 for (d = 0; d <= dim; ++d) { 2949 for (p = 0; p < numSubPoints[d]; ++p) { 2950 const PetscInt point = subpoints[d][p]; 2951 const PetscInt subpoint = firstSubPoint[d] + p; 2952 const PetscInt *cone; 2953 PetscInt coneSize, coneSizeNew, c, val; 2954 2955 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2956 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2957 if (cellHeight && (d == dim)) { 2958 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2959 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2960 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2961 if (val >= 0) coneSizeNew++; 2962 } 2963 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2964 } 2965 } 2966 } 2967 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2968 /* Set cones */ 2969 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2970 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 2971 for (d = 0; d <= dim; ++d) { 2972 for (p = 0; p < numSubPoints[d]; ++p) { 2973 const PetscInt point = subpoints[d][p]; 2974 const PetscInt subpoint = firstSubPoint[d] + p; 2975 const PetscInt *cone, *ornt; 2976 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 2977 2978 if (d == dim-1) { 2979 const PetscInt *support, *cone, *ornt; 2980 PetscInt supportSize, coneSize, s, subc; 2981 2982 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 2983 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 2984 for (s = 0; s < supportSize; ++s) { 2985 if ((support[s] < cMax) || (support[s] >= cEnd)) continue; 2986 ierr = PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);CHKERRQ(ierr); 2987 if (subc >= 0) { 2988 const PetscInt ccell = subpoints[d+1][subc]; 2989 2990 ierr = DMPlexGetCone(dm, ccell, &cone);CHKERRQ(ierr); 2991 ierr = DMPlexGetConeSize(dm, ccell, &coneSize);CHKERRQ(ierr); 2992 ierr = DMPlexGetConeOrientation(dm, ccell, &ornt);CHKERRQ(ierr); 2993 for (c = 0; c < coneSize; ++c) { 2994 if (cone[c] == point) { 2995 fornt = ornt[c]; 2996 break; 2997 } 2998 } 2999 break; 3000 } 3001 } 3002 } 3003 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 3004 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 3005 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3006 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 3007 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 3008 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 3009 if (subc >= 0) { 3010 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 3011 orntNew[coneSizeNew] = ornt[c]; 3012 ++coneSizeNew; 3013 } 3014 } 3015 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 3016 if (fornt < 0) { 3017 /* This should be replaced by a call to DMPlexReverseCell() */ 3018 #if 0 3019 ierr = DMPlexReverseCell(subdm, subpoint);CHKERRQ(ierr); 3020 #else 3021 for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) { 3022 PetscInt faceSize, tmp; 3023 3024 tmp = coneNew[c]; 3025 coneNew[c] = coneNew[coneSizeNew-1-c]; 3026 coneNew[coneSizeNew-1-c] = tmp; 3027 ierr = DMPlexGetConeSize(dm, cone[c], &faceSize);CHKERRQ(ierr); 3028 tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c]; 3029 orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c]; 3030 orntNew[coneSizeNew-1-c] = tmp; 3031 } 3032 } 3033 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 3034 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 3035 #endif 3036 } 3037 } 3038 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 3039 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3040 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3041 /* Build coordinates */ 3042 { 3043 PetscSection coordSection, subCoordSection; 3044 Vec coordinates, subCoordinates; 3045 PetscScalar *coords, *subCoords; 3046 PetscInt cdim, numComp, coordSize; 3047 const char *name; 3048 3049 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3050 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3051 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3052 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3053 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3054 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3055 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3056 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 3057 for (v = 0; v < numSubPoints[0]; ++v) { 3058 const PetscInt vertex = subpoints[0][v]; 3059 const PetscInt subvertex = firstSubPoint[0]+v; 3060 PetscInt dof; 3061 3062 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3063 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3064 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3065 } 3066 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3067 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3068 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3069 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3070 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3071 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3072 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3073 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3074 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3075 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3076 for (v = 0; v < numSubPoints[0]; ++v) { 3077 const PetscInt vertex = subpoints[0][v]; 3078 const PetscInt subvertex = firstSubPoint[0]+v; 3079 PetscInt dof, off, sdof, soff, d; 3080 3081 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3082 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3083 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3084 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3085 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3086 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3087 } 3088 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3089 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3090 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3091 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3092 } 3093 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3094 { 3095 PetscSF sfPoint, sfPointSub; 3096 IS subpIS; 3097 const PetscSFNode *remotePoints; 3098 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3099 const PetscInt *localPoints, *subpoints; 3100 PetscInt *slocalPoints; 3101 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3102 PetscMPIInt rank; 3103 3104 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3105 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3106 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3107 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3108 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 3109 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 3110 if (subpIS) { 3111 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 3112 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 3113 } 3114 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3115 if (numRoots >= 0) { 3116 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3117 for (p = 0; p < pEnd-pStart; ++p) { 3118 newLocalPoints[p].rank = -2; 3119 newLocalPoints[p].index = -2; 3120 } 3121 /* Set subleaves */ 3122 for (l = 0; l < numLeaves; ++l) { 3123 const PetscInt point = localPoints[l]; 3124 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3125 3126 if (subpoint < 0) continue; 3127 newLocalPoints[point-pStart].rank = rank; 3128 newLocalPoints[point-pStart].index = subpoint; 3129 ++numSubleaves; 3130 } 3131 /* Must put in owned subpoints */ 3132 for (p = pStart; p < pEnd; ++p) { 3133 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3134 3135 if (subpoint < 0) { 3136 newOwners[p-pStart].rank = -3; 3137 newOwners[p-pStart].index = -3; 3138 } else { 3139 newOwners[p-pStart].rank = rank; 3140 newOwners[p-pStart].index = subpoint; 3141 } 3142 } 3143 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3144 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3145 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3146 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3147 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 3148 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 3149 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3150 const PetscInt point = localPoints[l]; 3151 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3152 3153 if (subpoint < 0) continue; 3154 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3155 slocalPoints[sl] = subpoint; 3156 sremotePoints[sl].rank = newLocalPoints[point].rank; 3157 sremotePoints[sl].index = newLocalPoints[point].index; 3158 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3159 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3160 ++sl; 3161 } 3162 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3163 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3164 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3165 } 3166 if (subpIS) { 3167 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 3168 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 3169 } 3170 } 3171 /* Cleanup */ 3172 for (d = 0; d <= dim; ++d) { 3173 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3174 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 3175 } 3176 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 3177 PetscFunctionReturn(0); 3178 } 3179 3180 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 3181 { 3182 PetscErrorCode ierr; 3183 3184 PetscFunctionBegin; 3185 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);CHKERRQ(ierr); 3186 PetscFunctionReturn(0); 3187 } 3188 3189 /*@ 3190 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3191 3192 Input Parameters: 3193 + dm - The original mesh 3194 . vertexLabel - The DMLabel marking vertices contained in the surface 3195 - value - The label value to use 3196 3197 Output Parameter: 3198 . subdm - The surface mesh 3199 3200 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3201 3202 Level: developer 3203 3204 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3205 @*/ 3206 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 3207 { 3208 PetscInt dim, cdim, depth; 3209 PetscErrorCode ierr; 3210 3211 PetscFunctionBegin; 3212 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3213 PetscValidPointer(subdm, 3); 3214 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3215 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3216 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3217 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3218 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3219 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3220 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3221 if (depth == dim) { 3222 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3223 } else { 3224 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3225 } 3226 PetscFunctionReturn(0); 3227 } 3228 3229 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3230 { 3231 MPI_Comm comm; 3232 DMLabel subpointMap; 3233 IS subvertexIS; 3234 const PetscInt *subVertices; 3235 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3236 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3237 PetscInt cMax, c, f; 3238 PetscErrorCode ierr; 3239 3240 PetscFunctionBegin; 3241 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 3242 /* Create subpointMap which marks the submesh */ 3243 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 3244 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3245 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3246 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 3247 /* Setup chart */ 3248 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 3249 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 3250 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 3251 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 3252 /* Set cone sizes */ 3253 firstSubVertex = numSubCells; 3254 firstSubFace = numSubCells+numSubVertices; 3255 newFacePoint = firstSubFace; 3256 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 3257 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3258 for (c = 0; c < numSubCells; ++c) { 3259 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 3260 } 3261 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3262 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 3263 } 3264 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3265 /* Create face cones */ 3266 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3267 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3268 ierr = DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 3269 for (c = 0; c < numSubCells; ++c) { 3270 const PetscInt cell = subCells[c]; 3271 const PetscInt subcell = c; 3272 const PetscInt *cone, *cells; 3273 PetscInt numCells, subVertex, p, v; 3274 3275 if (cell < cMax) continue; 3276 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 3277 for (v = 0; v < nFV; ++v) { 3278 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 3279 subface[v] = firstSubVertex+subVertex; 3280 } 3281 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 3282 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 3283 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3284 /* Not true in parallel 3285 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3286 for (p = 0; p < numCells; ++p) { 3287 PetscInt negsubcell; 3288 3289 if (cells[p] >= cMax) continue; 3290 /* I know this is a crap search */ 3291 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3292 if (subCells[negsubcell] == cells[p]) break; 3293 } 3294 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3295 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 3296 } 3297 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3298 ++newFacePoint; 3299 } 3300 ierr = DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 3301 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3302 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3303 /* Build coordinates */ 3304 { 3305 PetscSection coordSection, subCoordSection; 3306 Vec coordinates, subCoordinates; 3307 PetscScalar *coords, *subCoords; 3308 PetscInt cdim, numComp, coordSize, v; 3309 const char *name; 3310 3311 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3312 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3313 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3314 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3315 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3316 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3317 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3318 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3319 for (v = 0; v < numSubVertices; ++v) { 3320 const PetscInt vertex = subVertices[v]; 3321 const PetscInt subvertex = firstSubVertex+v; 3322 PetscInt dof; 3323 3324 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3325 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3326 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3327 } 3328 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3329 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3330 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3331 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3332 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3333 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3334 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3335 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3336 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3337 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3338 for (v = 0; v < numSubVertices; ++v) { 3339 const PetscInt vertex = subVertices[v]; 3340 const PetscInt subvertex = firstSubVertex+v; 3341 PetscInt dof, off, sdof, soff, d; 3342 3343 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3344 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3345 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3346 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3347 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3348 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3349 } 3350 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3351 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3352 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3353 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3354 } 3355 /* Build SF */ 3356 CHKMEMQ; 3357 { 3358 PetscSF sfPoint, sfPointSub; 3359 const PetscSFNode *remotePoints; 3360 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3361 const PetscInt *localPoints; 3362 PetscInt *slocalPoints; 3363 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3364 PetscMPIInt rank; 3365 3366 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3367 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3368 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3369 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3370 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3371 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3372 if (numRoots >= 0) { 3373 /* Only vertices should be shared */ 3374 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3375 for (p = 0; p < pEnd-pStart; ++p) { 3376 newLocalPoints[p].rank = -2; 3377 newLocalPoints[p].index = -2; 3378 } 3379 /* Set subleaves */ 3380 for (l = 0; l < numLeaves; ++l) { 3381 const PetscInt point = localPoints[l]; 3382 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3383 3384 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3385 if (subPoint < 0) continue; 3386 newLocalPoints[point-pStart].rank = rank; 3387 newLocalPoints[point-pStart].index = subPoint; 3388 ++numSubLeaves; 3389 } 3390 /* Must put in owned subpoints */ 3391 for (p = pStart; p < pEnd; ++p) { 3392 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3393 3394 if (subPoint < 0) { 3395 newOwners[p-pStart].rank = -3; 3396 newOwners[p-pStart].index = -3; 3397 } else { 3398 newOwners[p-pStart].rank = rank; 3399 newOwners[p-pStart].index = subPoint; 3400 } 3401 } 3402 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3403 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3404 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3405 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3406 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3407 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3408 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3409 const PetscInt point = localPoints[l]; 3410 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3411 3412 if (subPoint < 0) continue; 3413 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3414 slocalPoints[sl] = subPoint; 3415 sremotePoints[sl].rank = newLocalPoints[point].rank; 3416 sremotePoints[sl].index = newLocalPoints[point].index; 3417 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3418 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3419 ++sl; 3420 } 3421 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3422 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3423 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3424 } 3425 } 3426 CHKMEMQ; 3427 /* Cleanup */ 3428 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3429 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3430 ierr = PetscFree(subCells);CHKERRQ(ierr); 3431 PetscFunctionReturn(0); 3432 } 3433 3434 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3435 { 3436 DMLabel label = NULL; 3437 PetscErrorCode ierr; 3438 3439 PetscFunctionBegin; 3440 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3441 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);CHKERRQ(ierr); 3442 PetscFunctionReturn(0); 3443 } 3444 3445 /*@C 3446 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. 3447 3448 Input Parameters: 3449 + dm - The original mesh 3450 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3451 . label - A label name, or NULL 3452 - value - A label value 3453 3454 Output Parameter: 3455 . subdm - The surface mesh 3456 3457 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3458 3459 Level: developer 3460 3461 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3462 @*/ 3463 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3464 { 3465 PetscInt dim, cdim, depth; 3466 PetscErrorCode ierr; 3467 3468 PetscFunctionBegin; 3469 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3470 PetscValidPointer(subdm, 5); 3471 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3472 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3473 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3474 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3475 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3476 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3477 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3478 if (depth == dim) { 3479 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3480 } else { 3481 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3482 } 3483 PetscFunctionReturn(0); 3484 } 3485 3486 /*@ 3487 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3488 3489 Input Parameters: 3490 + dm - The original mesh 3491 . cellLabel - The DMLabel marking cells contained in the new mesh 3492 - value - The label value to use 3493 3494 Output Parameter: 3495 . subdm - The new mesh 3496 3497 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3498 3499 Level: developer 3500 3501 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3502 @*/ 3503 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3504 { 3505 PetscInt dim; 3506 PetscErrorCode ierr; 3507 3508 PetscFunctionBegin; 3509 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3510 PetscValidPointer(subdm, 3); 3511 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3512 ierr = DMCreate(PetscObjectComm((PetscObject) dm), subdm);CHKERRQ(ierr); 3513 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3514 ierr = DMSetDimension(*subdm, dim);CHKERRQ(ierr); 3515 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3516 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);CHKERRQ(ierr); 3517 PetscFunctionReturn(0); 3518 } 3519 3520 /*@ 3521 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3522 3523 Input Parameter: 3524 . dm - The submesh DM 3525 3526 Output Parameter: 3527 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3528 3529 Level: developer 3530 3531 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3532 @*/ 3533 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3534 { 3535 PetscFunctionBegin; 3536 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3537 PetscValidPointer(subpointMap, 2); 3538 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3539 PetscFunctionReturn(0); 3540 } 3541 3542 /*@ 3543 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3544 3545 Input Parameters: 3546 + dm - The submesh DM 3547 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3548 3549 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3550 3551 Level: developer 3552 3553 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3554 @*/ 3555 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3556 { 3557 DM_Plex *mesh = (DM_Plex *) dm->data; 3558 DMLabel tmp; 3559 PetscErrorCode ierr; 3560 3561 PetscFunctionBegin; 3562 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3563 tmp = mesh->subpointMap; 3564 mesh->subpointMap = subpointMap; 3565 ++mesh->subpointMap->refct; 3566 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3567 PetscFunctionReturn(0); 3568 } 3569 3570 /*@ 3571 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3572 3573 Input Parameter: 3574 . dm - The submesh DM 3575 3576 Output Parameter: 3577 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3578 3579 Note: This IS is guaranteed to be sorted by the construction of the submesh 3580 3581 Level: developer 3582 3583 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3584 @*/ 3585 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3586 { 3587 MPI_Comm comm; 3588 DMLabel subpointMap; 3589 IS is; 3590 const PetscInt *opoints; 3591 PetscInt *points, *depths; 3592 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3593 PetscErrorCode ierr; 3594 3595 PetscFunctionBegin; 3596 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3597 PetscValidPointer(subpointIS, 2); 3598 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3599 *subpointIS = NULL; 3600 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3601 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3602 if (subpointMap && depth >= 0) { 3603 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3604 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3605 ierr = DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);CHKERRQ(ierr); 3606 depths[0] = depth; 3607 depths[1] = 0; 3608 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3609 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3610 for(d = 0, off = 0; d <= depth; ++d) { 3611 const PetscInt dep = depths[d]; 3612 3613 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3614 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3615 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3616 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); 3617 } else { 3618 if (!n) { 3619 if (d == 0) { 3620 /* Missing cells */ 3621 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3622 } else { 3623 /* Missing faces */ 3624 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3625 } 3626 } 3627 } 3628 if (n) { 3629 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3630 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3631 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3632 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3633 ierr = ISDestroy(&is);CHKERRQ(ierr); 3634 } 3635 } 3636 ierr = DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);CHKERRQ(ierr); 3637 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3638 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3639 } 3640 PetscFunctionReturn(0); 3641 } 3642 3643 /*@ 3644 DMPlexGetSubpoint - Return the subpoint corresponding to a point in the original mesh. If the DM 3645 is not a submesh, just return the input point. 3646 3647 Note collective 3648 3649 Input Parameters: 3650 + dm - The submesh DM 3651 - p - The point in the original, from which the submesh was created 3652 3653 Output Parameter: 3654 . subp - The point in the submesh 3655 3656 Level: developer 3657 3658 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap(), DMPlexCreateSubpointIS() 3659 @*/ 3660 PetscErrorCode DMPlexGetSubpoint(DM dm, PetscInt p, PetscInt *subp) 3661 { 3662 DMLabel spmap; 3663 PetscErrorCode ierr; 3664 3665 *subp = p; 3666 ierr = DMPlexGetSubpointMap(dm, &spmap);CHKERRQ(ierr); 3667 if (spmap) { 3668 IS subpointIS; 3669 const PetscInt *subpoints; 3670 PetscInt numSubpoints; 3671 3672 /* TODO Cache the IS, making it look like an index */ 3673 ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr); 3674 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 3675 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 3676 ierr = PetscFindInt(p, numSubpoints, subpoints, subp);CHKERRQ(ierr); 3677 if (*subp < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d not found in submesh", p); 3678 ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr); 3679 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 3680 } 3681 PetscFunctionReturn(0); 3682 } 3683