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