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