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