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