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 = DMSetCoordinateDim(dmNew, dim);CHKERRQ(ierr); 366 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 367 /* Step 8: Convert coordinates */ 368 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 369 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 370 ierr = DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);CHKERRQ(ierr); 371 ierr = DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);CHKERRQ(ierr); 372 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 373 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);CHKERRQ(ierr); 374 ierr = PetscSectionSetNumFields(newCoordSection, 1);CHKERRQ(ierr); 375 ierr = PetscSectionSetFieldComponents(newCoordSection, 0, dim);CHKERRQ(ierr); 376 ierr = PetscSectionGetChart(coordSection, &sStart, &sEnd);CHKERRQ(ierr); 377 hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE; 378 ierr = PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);CHKERRQ(ierr); 379 if (hasCells) { 380 for (c = cStart; c < cEnd; ++c) { 381 PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof; 382 383 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 384 ierr = PetscSectionSetDof(newCoordSection, cNew, dof);CHKERRQ(ierr); 385 ierr = PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);CHKERRQ(ierr); 386 } 387 } 388 for (v = vStartNew; v < vEndNew; ++v) { 389 ierr = PetscSectionSetDof(newCoordSection, v, dim);CHKERRQ(ierr); 390 ierr = PetscSectionSetFieldDof(newCoordSection, v, 0, dim);CHKERRQ(ierr); 391 } 392 ierr = PetscSectionSetUp(newCoordSection);CHKERRQ(ierr); 393 ierr = DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);CHKERRQ(ierr); 394 ierr = PetscSectionGetStorageSize(newCoordSection, &coordSize);CHKERRQ(ierr); 395 ierr = VecCreate(PETSC_COMM_SELF, &newCoordinates);CHKERRQ(ierr); 396 ierr = PetscObjectSetName((PetscObject) newCoordinates, "coordinates");CHKERRQ(ierr); 397 ierr = VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 398 ierr = VecSetBlockSize(newCoordinates, dim);CHKERRQ(ierr); 399 ierr = VecSetType(newCoordinates,VECSTANDARD);CHKERRQ(ierr); 400 ierr = DMSetCoordinatesLocal(dmNew, newCoordinates);CHKERRQ(ierr); 401 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 402 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 403 ierr = VecGetArray(newCoordinates, &newCoords);CHKERRQ(ierr); 404 if (hasCells) { 405 for (c = cStart; c < cEnd; ++c) { 406 PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d; 407 408 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 409 ierr = PetscSectionGetOffset(coordSection, c, &off);CHKERRQ(ierr); 410 ierr = PetscSectionGetOffset(newCoordSection, cNew, &noff);CHKERRQ(ierr); 411 for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d]; 412 } 413 } 414 for (v = vStart; v < vEnd; ++v) { 415 PetscInt dof, off, noff, d; 416 417 ierr = PetscSectionGetDof(coordSection, v, &dof);CHKERRQ(ierr); 418 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 419 ierr = PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);CHKERRQ(ierr); 420 for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d]; 421 } 422 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 423 ierr = VecRestoreArray(newCoordinates, &newCoords);CHKERRQ(ierr); 424 ierr = VecDestroy(&newCoordinates);CHKERRQ(ierr); 425 ierr = PetscSectionDestroy(&newCoordSection);CHKERRQ(ierr); 426 PetscFunctionReturn(0); 427 } 428 429 static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew) 430 { 431 PetscInt depth = 0; 432 PetscSF sfPoint, sfPointNew; 433 const PetscSFNode *remotePoints; 434 PetscSFNode *gremotePoints; 435 const PetscInt *localPoints; 436 PetscInt *glocalPoints, *newLocation, *newRemoteLocation; 437 PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0; 438 PetscErrorCode ierr; 439 440 PetscFunctionBegin; 441 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 442 /* Step 9: Convert pointSF */ 443 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 444 ierr = DMGetPointSF(dmNew, &sfPointNew);CHKERRQ(ierr); 445 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 446 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 447 totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd; 448 if (numRoots >= 0) { 449 ierr = PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);CHKERRQ(ierr); 450 for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift); 451 ierr = PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);CHKERRQ(ierr); 452 ierr = PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);CHKERRQ(ierr); 453 ierr = PetscMalloc1(numLeaves, &glocalPoints);CHKERRQ(ierr); 454 ierr = PetscMalloc1(numLeaves, &gremotePoints);CHKERRQ(ierr); 455 for (l = 0; l < numLeaves; ++l) { 456 glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift); 457 gremotePoints[l].rank = remotePoints[l].rank; 458 gremotePoints[l].index = newRemoteLocation[localPoints[l]]; 459 } 460 ierr = PetscFree2(newLocation,newRemoteLocation);CHKERRQ(ierr); 461 ierr = PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 462 } 463 PetscFunctionReturn(0); 464 } 465 466 static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew) 467 { 468 PetscSF sfPoint; 469 DMLabel vtkLabel, ghostLabel; 470 const PetscSFNode *leafRemote; 471 const PetscInt *leafLocal; 472 PetscInt depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f; 473 PetscMPIInt rank; 474 PetscErrorCode ierr; 475 476 PetscFunctionBegin; 477 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 478 /* Step 10: Convert labels */ 479 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 480 for (l = 0; l < numLabels; ++l) { 481 DMLabel label, newlabel; 482 const char *lname; 483 PetscBool isDepth; 484 IS valueIS; 485 const PetscInt *values; 486 PetscInt numValues, val; 487 488 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr); 489 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 490 if (isDepth) continue; 491 ierr = DMCreateLabel(dmNew, lname);CHKERRQ(ierr); 492 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr); 493 ierr = DMGetLabel(dmNew, lname, &newlabel);CHKERRQ(ierr); 494 ierr = DMLabelGetDefaultValue(label,&val);CHKERRQ(ierr); 495 ierr = DMLabelSetDefaultValue(newlabel,val);CHKERRQ(ierr); 496 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 497 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 498 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 499 for (val = 0; val < numValues; ++val) { 500 IS pointIS; 501 const PetscInt *points; 502 PetscInt numPoints, p; 503 504 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 505 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 506 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 507 for (p = 0; p < numPoints; ++p) { 508 const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift); 509 510 ierr = DMLabelSetValue(newlabel, newpoint, values[val]);CHKERRQ(ierr); 511 } 512 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 513 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 514 } 515 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 516 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 517 } 518 /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */ 519 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 520 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 521 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 522 ierr = PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);CHKERRQ(ierr); 523 ierr = DMCreateLabel(dmNew, "vtk");CHKERRQ(ierr); 524 ierr = DMCreateLabel(dmNew, "ghost");CHKERRQ(ierr); 525 ierr = DMGetLabel(dmNew, "vtk", &vtkLabel);CHKERRQ(ierr); 526 ierr = DMGetLabel(dmNew, "ghost", &ghostLabel);CHKERRQ(ierr); 527 for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) { 528 for (; c < leafLocal[l] && c < cEnd; ++c) { 529 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 530 } 531 if (leafLocal[l] >= cEnd) break; 532 if (leafRemote[l].rank == rank) { 533 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 534 } else { 535 ierr = DMLabelSetValue(ghostLabel, c, 2);CHKERRQ(ierr); 536 } 537 } 538 for (; c < cEnd; ++c) { 539 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 540 } 541 if (0) { 542 ierr = DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 543 } 544 ierr = DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 545 for (f = fStart; f < fEnd; ++f) { 546 PetscInt numCells; 547 548 ierr = DMPlexGetSupportSize(dmNew, f, &numCells);CHKERRQ(ierr); 549 if (numCells < 2) { 550 ierr = DMLabelSetValue(ghostLabel, f, 1);CHKERRQ(ierr); 551 } else { 552 const PetscInt *cells = NULL; 553 PetscInt vA, vB; 554 555 ierr = DMPlexGetSupport(dmNew, f, &cells);CHKERRQ(ierr); 556 ierr = DMLabelGetValue(vtkLabel, cells[0], &vA);CHKERRQ(ierr); 557 ierr = DMLabelGetValue(vtkLabel, cells[1], &vB);CHKERRQ(ierr); 558 if (vA != 1 && vB != 1) {ierr = DMLabelSetValue(ghostLabel, f, 1);CHKERRQ(ierr);} 559 } 560 } 561 if (0) { 562 ierr = DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 563 } 564 PetscFunctionReturn(0); 565 } 566 567 static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew) 568 { 569 DM refTree; 570 PetscSection pSec; 571 PetscInt *parents, *childIDs; 572 PetscErrorCode ierr; 573 574 PetscFunctionBegin; 575 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 576 ierr = DMPlexSetReferenceTree(dmNew,refTree);CHKERRQ(ierr); 577 ierr = DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);CHKERRQ(ierr); 578 if (pSec) { 579 PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth; 580 PetscInt *childIDsShifted; 581 PetscSection pSecShifted; 582 583 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 584 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 585 pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift); 586 pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift); 587 ierr = PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);CHKERRQ(ierr); 588 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);CHKERRQ(ierr); 589 ierr = PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);CHKERRQ(ierr); 590 for (p = pStartShifted; p < pEndShifted; p++) { 591 /* start off assuming no children */ 592 ierr = PetscSectionSetDof(pSecShifted,p,0);CHKERRQ(ierr); 593 } 594 for (p = pStart; p < pEnd; p++) { 595 PetscInt dof; 596 PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift); 597 598 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 599 ierr = PetscSectionSetDof(pSecShifted,pNew,dof);CHKERRQ(ierr); 600 } 601 ierr = PetscSectionSetUp(pSecShifted);CHKERRQ(ierr); 602 for (p = pStart; p < pEnd; p++) { 603 PetscInt dof; 604 PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift); 605 606 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 607 if (dof) { 608 PetscInt off, offNew; 609 610 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 611 ierr = PetscSectionGetOffset(pSecShifted,pNew,&offNew);CHKERRQ(ierr); 612 parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift); 613 childIDsShifted[offNew] = childIDs[off]; 614 } 615 } 616 ierr = DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);CHKERRQ(ierr); 617 ierr = PetscFree2(parentsShifted,childIDsShifted);CHKERRQ(ierr); 618 ierr = PetscSectionDestroy(&pSecShifted);CHKERRQ(ierr); 619 } 620 PetscFunctionReturn(0); 621 } 622 623 static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm) 624 { 625 PetscSF sf; 626 IS valueIS; 627 const PetscInt *values, *leaves; 628 PetscInt *depthShift; 629 PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c; 630 PetscBool isper; 631 const PetscReal *maxCell, *L; 632 const DMBoundaryType *bd; 633 PetscErrorCode ierr; 634 635 PetscFunctionBegin; 636 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 637 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 638 nleaves = PetscMax(0, nleaves); 639 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 640 /* Count ghost cells */ 641 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 642 ierr = ISGetLocalSize(valueIS, &numFS);CHKERRQ(ierr); 643 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 644 Ng = 0; 645 for (fs = 0; fs < numFS; ++fs) { 646 IS faceIS; 647 const PetscInt *faces; 648 PetscInt numFaces, f, numBdFaces = 0; 649 650 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 651 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 652 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 653 for (f = 0; f < numFaces; ++f) { 654 PetscInt numChildren; 655 656 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 657 ierr = DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);CHKERRQ(ierr); 658 /* non-local and ancestors points don't get to register ghosts */ 659 if (loc >= 0 || numChildren) continue; 660 if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces; 661 } 662 Ng += numBdFaces; 663 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 664 } 665 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 666 ierr = PetscMalloc1(2*(depth+1), &depthShift);CHKERRQ(ierr); 667 for (d = 0; d <= depth; d++) { 668 PetscInt dEnd; 669 670 ierr = DMPlexGetDepthStratum(dm,d,NULL,&dEnd);CHKERRQ(ierr); 671 depthShift[2*d] = dEnd; 672 depthShift[2*d+1] = 0; 673 } 674 if (depth >= 0) depthShift[2*depth+1] = Ng; 675 ierr = DMPlexShiftPointSetUp_Internal(depth,depthShift);CHKERRQ(ierr); 676 ierr = DMPlexShiftSizes_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 677 /* Step 3: Set cone/support sizes for new points */ 678 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 679 ierr = DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 680 for (c = cEnd; c < cEnd + Ng; ++c) { 681 ierr = DMPlexSetConeSize(gdm, c, 1);CHKERRQ(ierr); 682 } 683 for (fs = 0; fs < numFS; ++fs) { 684 IS faceIS; 685 const PetscInt *faces; 686 PetscInt numFaces, f; 687 688 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 689 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 690 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 691 for (f = 0; f < numFaces; ++f) { 692 PetscInt size, numChildren; 693 694 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 695 ierr = DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);CHKERRQ(ierr); 696 if (loc >= 0 || numChildren) continue; 697 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 698 ierr = DMPlexGetSupportSize(dm, faces[f], &size);CHKERRQ(ierr); 699 if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size); 700 ierr = DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);CHKERRQ(ierr); 701 } 702 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 703 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 704 } 705 /* Step 4: Setup ghosted DM */ 706 ierr = DMSetUp(gdm);CHKERRQ(ierr); 707 ierr = DMPlexShiftPoints_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 708 /* Step 6: Set cones and supports for new points */ 709 ghostCell = cEnd; 710 for (fs = 0; fs < numFS; ++fs) { 711 IS faceIS; 712 const PetscInt *faces; 713 PetscInt numFaces, f; 714 715 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 716 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 717 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 718 for (f = 0; f < numFaces; ++f) { 719 PetscInt newFace = faces[f] + Ng, numChildren; 720 721 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 722 ierr = DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);CHKERRQ(ierr); 723 if (loc >= 0 || numChildren) continue; 724 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 725 ierr = DMPlexSetCone(gdm, ghostCell, &newFace);CHKERRQ(ierr); 726 ierr = DMPlexInsertSupport(gdm, newFace, 1, ghostCell);CHKERRQ(ierr); 727 ++ghostCell; 728 } 729 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 730 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 731 } 732 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 733 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 734 /* Step 7: Stratify */ 735 ierr = DMPlexStratify(gdm);CHKERRQ(ierr); 736 ierr = DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 737 ierr = DMPlexShiftSF_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 738 ierr = DMPlexShiftLabels_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 739 ierr = DMPlexShiftTree_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 740 ierr = PetscFree(depthShift);CHKERRQ(ierr); 741 /* Step 7: Periodicity */ 742 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 743 ierr = DMSetPeriodicity(gdm, isper, maxCell, L, bd);CHKERRQ(ierr); 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 } 1402 } 1403 } 1404 } 1405 for (sp = 0; sp < numSP; ++sp) { 1406 const PetscInt dep = values[sp]; 1407 1408 if ((dep < 0) || (dep > depth)) continue; 1409 if (splitIS[dep]) {ierr = ISRestoreIndices(splitIS[dep], &splitPoints[dep]);CHKERRQ(ierr);} 1410 ierr = ISDestroy(&splitIS[dep]);CHKERRQ(ierr); 1411 if (unsplitIS[dep]) {ierr = ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);CHKERRQ(ierr);} 1412 ierr = ISDestroy(&unsplitIS[dep]);CHKERRQ(ierr); 1413 } 1414 if (label) { 1415 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 1416 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 1417 } 1418 for (d = 0; d <= depth; ++d) { 1419 ierr = DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);CHKERRQ(ierr); 1420 pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d]; 1421 } 1422 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); 1423 ierr = PetscFree3(coneNew, coneONew, supportNew);CHKERRQ(ierr); 1424 ierr = PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);CHKERRQ(ierr); 1425 ierr = PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);CHKERRQ(ierr); 1426 PetscFunctionReturn(0); 1427 } 1428 1429 /*@C 1430 DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface 1431 1432 Collective on dm 1433 1434 Input Parameters: 1435 + dm - The original DM 1436 - label - The label specifying the boundary faces (this could be auto-generated) 1437 1438 Output Parameters: 1439 - dmSplit - The new DM 1440 1441 Level: developer 1442 1443 .seealso: DMCreate(), DMPlexLabelCohesiveComplete() 1444 @*/ 1445 PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit) 1446 { 1447 DM sdm; 1448 PetscInt dim; 1449 PetscErrorCode ierr; 1450 1451 PetscFunctionBegin; 1452 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1453 PetscValidPointer(dmSplit, 3); 1454 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &sdm);CHKERRQ(ierr); 1455 ierr = DMSetType(sdm, DMPLEX);CHKERRQ(ierr); 1456 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1457 ierr = DMSetDimension(sdm, dim);CHKERRQ(ierr); 1458 switch (dim) { 1459 case 2: 1460 case 3: 1461 ierr = DMPlexConstructCohesiveCells_Internal(dm, label, sdm);CHKERRQ(ierr); 1462 break; 1463 default: 1464 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim); 1465 } 1466 *dmSplit = sdm; 1467 PetscFunctionReturn(0); 1468 } 1469 1470 /* Returns the side of the surface for a given cell with a face on the surface */ 1471 static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos) 1472 { 1473 const PetscInt *cone, *ornt; 1474 PetscInt dim, coneSize, c; 1475 PetscErrorCode ierr; 1476 1477 PetscFunctionBegin; 1478 *pos = PETSC_TRUE; 1479 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1480 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1481 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 1482 ierr = DMPlexGetConeOrientation(dm, cell, &ornt);CHKERRQ(ierr); 1483 for (c = 0; c < coneSize; ++c) { 1484 if (cone[c] == face) { 1485 PetscInt o = ornt[c]; 1486 1487 if (subdm) { 1488 const PetscInt *subcone, *subornt; 1489 PetscInt subpoint, subface, subconeSize, sc; 1490 1491 ierr = PetscFindInt(cell, numSubpoints, subpoints, &subpoint);CHKERRQ(ierr); 1492 ierr = PetscFindInt(face, numSubpoints, subpoints, &subface);CHKERRQ(ierr); 1493 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 1494 ierr = DMPlexGetCone(subdm, subpoint, &subcone);CHKERRQ(ierr); 1495 ierr = DMPlexGetConeOrientation(subdm, subpoint, &subornt);CHKERRQ(ierr); 1496 for (sc = 0; sc < subconeSize; ++sc) { 1497 if (subcone[sc] == subface) { 1498 o = subornt[0]; 1499 break; 1500 } 1501 } 1502 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); 1503 } 1504 if (o >= 0) *pos = PETSC_TRUE; 1505 else *pos = PETSC_FALSE; 1506 break; 1507 } 1508 } 1509 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); 1510 PetscFunctionReturn(0); 1511 } 1512 1513 /*@ 1514 DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces 1515 to complete the surface 1516 1517 Input Parameters: 1518 + dm - The DM 1519 . label - A DMLabel marking the surface 1520 . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically 1521 . flip - Flag to flip the submesh normal and replace points on the other side 1522 - subdm - The subDM associated with the label, or NULL 1523 1524 Output Parameter: 1525 . label - A DMLabel marking all surface points 1526 1527 Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation. 1528 1529 Level: developer 1530 1531 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete() 1532 @*/ 1533 PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm) 1534 { 1535 DMLabel depthLabel; 1536 IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL; 1537 const PetscInt *points, *subpoints; 1538 const PetscInt rev = flip ? -1 : 1; 1539 PetscInt *pMax; 1540 PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val; 1541 PetscErrorCode ierr; 1542 1543 PetscFunctionBegin; 1544 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1545 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1546 pSize = PetscMax(depth, dim) + 1; 1547 ierr = PetscMalloc1(pSize,&pMax);CHKERRQ(ierr); 1548 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1549 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1550 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1551 if (subdm) { 1552 ierr = DMPlexCreateSubpointIS(subdm, &subpointIS);CHKERRQ(ierr); 1553 if (subpointIS) { 1554 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 1555 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 1556 } 1557 } 1558 /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */ 1559 ierr = DMLabelGetStratumIS(label, dim-1, &dimIS);CHKERRQ(ierr); 1560 if (!dimIS) { 1561 ierr = PetscFree(pMax);CHKERRQ(ierr); 1562 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1563 PetscFunctionReturn(0); 1564 } 1565 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1566 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1567 for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */ 1568 const PetscInt *support; 1569 PetscInt supportSize, s; 1570 1571 ierr = DMPlexGetSupportSize(dm, points[p], &supportSize);CHKERRQ(ierr); 1572 if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize); 1573 ierr = DMPlexGetSupport(dm, points[p], &support);CHKERRQ(ierr); 1574 for (s = 0; s < supportSize; ++s) { 1575 const PetscInt *cone; 1576 PetscInt coneSize, c; 1577 PetscBool pos; 1578 1579 ierr = GetSurfaceSide_Static(dm, subdm, numSubpoints, subpoints, support[s], points[p], &pos);CHKERRQ(ierr); 1580 if (pos) {ierr = DMLabelSetValue(label, support[s], rev*(shift+dim));CHKERRQ(ierr);} 1581 else {ierr = DMLabelSetValue(label, support[s], -rev*(shift+dim));CHKERRQ(ierr);} 1582 if (rev < 0) pos = !pos ? PETSC_TRUE : PETSC_FALSE; 1583 /* Put faces touching the fault in the label */ 1584 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1585 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1586 for (c = 0; c < coneSize; ++c) { 1587 const PetscInt point = cone[c]; 1588 1589 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1590 if (val == -1) { 1591 PetscInt *closure = NULL; 1592 PetscInt closureSize, cl; 1593 1594 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1595 for (cl = 0; cl < closureSize*2; cl += 2) { 1596 const PetscInt clp = closure[cl]; 1597 PetscInt bval = -1; 1598 1599 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1600 if (blabel) {ierr = DMLabelGetValue(blabel, clp, &bval);CHKERRQ(ierr);} 1601 if ((val >= 0) && (val < dim-1) && (bval < 0)) { 1602 ierr = DMLabelSetValue(label, point, pos == PETSC_TRUE ? shift+dim-1 : -(shift+dim-1));CHKERRQ(ierr); 1603 break; 1604 } 1605 } 1606 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1607 } 1608 } 1609 } 1610 } 1611 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1612 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1613 if (subpointIS) {ierr = ISRestoreIndices(subpointIS, &subpoints);CHKERRQ(ierr);} 1614 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1615 /* Mark boundary points as unsplit */ 1616 if (blabel) { 1617 ierr = DMLabelGetStratumIS(blabel, 1, &dimIS);CHKERRQ(ierr); 1618 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1619 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1620 for (p = 0; p < numPoints; ++p) { 1621 const PetscInt point = points[p]; 1622 PetscInt val, bval; 1623 1624 ierr = DMLabelGetValue(blabel, point, &bval);CHKERRQ(ierr); 1625 if (bval >= 0) { 1626 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1627 if ((val < 0) || (val > dim)) { 1628 /* This could be a point added from splitting a vertex on an adjacent fault, otherwise its just wrong */ 1629 ierr = DMLabelClearValue(blabel, point, bval);CHKERRQ(ierr); 1630 } 1631 } 1632 } 1633 for (p = 0; p < numPoints; ++p) { 1634 const PetscInt point = points[p]; 1635 PetscInt val, bval; 1636 1637 ierr = DMLabelGetValue(blabel, point, &bval);CHKERRQ(ierr); 1638 if (bval >= 0) { 1639 const PetscInt *cone, *support; 1640 PetscInt coneSize, supportSize, s, valA, valB, valE; 1641 1642 /* Mark as unsplit */ 1643 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1644 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); 1645 ierr = DMLabelClearValue(label, point, val);CHKERRQ(ierr); 1646 ierr = DMLabelSetValue(label, point, shift2+val);CHKERRQ(ierr); 1647 /* Check for cross-edge 1648 A cross-edge has endpoints which are both on the boundary of the surface, but the edge itself is not. */ 1649 if (val != 0) continue; 1650 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1651 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1652 for (s = 0; s < supportSize; ++s) { 1653 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 1654 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 1655 if (coneSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Edge %D has %D vertices != 2", support[s], coneSize); 1656 ierr = DMLabelGetValue(blabel, cone[0], &valA);CHKERRQ(ierr); 1657 ierr = DMLabelGetValue(blabel, cone[1], &valB);CHKERRQ(ierr); 1658 ierr = DMLabelGetValue(blabel, support[s], &valE);CHKERRQ(ierr); 1659 if ((valE < 0) && (valA >= 0) && (valB >= 0) && (cone[0] != cone[1])) {ierr = DMLabelSetValue(blabel, support[s], 2);CHKERRQ(ierr);} 1660 } 1661 } 1662 } 1663 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1664 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1665 } 1666 /* Search for other cells/faces/edges connected to the fault by a vertex */ 1667 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1668 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1669 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 1670 cMax = cMax < 0 ? cEnd : cMax; 1671 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1672 ierr = DMLabelGetStratumIS(label, 0, &dimIS);CHKERRQ(ierr); 1673 if (blabel) {ierr = DMLabelGetStratumIS(blabel, 2, &crossEdgeIS);CHKERRQ(ierr);} 1674 if (dimIS && crossEdgeIS) { 1675 IS vertIS = dimIS; 1676 1677 ierr = ISExpand(vertIS, crossEdgeIS, &dimIS);CHKERRQ(ierr); 1678 ierr = ISDestroy(&crossEdgeIS);CHKERRQ(ierr); 1679 ierr = ISDestroy(&vertIS);CHKERRQ(ierr); 1680 } 1681 if (!dimIS) { 1682 ierr = PetscFree(pMax);CHKERRQ(ierr); 1683 PetscFunctionReturn(0); 1684 } 1685 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1686 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1687 for (p = 0; p < numPoints; ++p) { /* Loop over fault vertices */ 1688 PetscInt *star = NULL; 1689 PetscInt starSize, s; 1690 PetscInt again = 1; /* 0: Finished 1: Keep iterating after a change 2: No change */ 1691 1692 /* All points connected to the fault are inside a cell, so at the top level we will only check cells */ 1693 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1694 while (again) { 1695 if (again > 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Could not classify all cells connected to the fault"); 1696 again = 0; 1697 for (s = 0; s < starSize*2; s += 2) { 1698 const PetscInt point = star[s]; 1699 const PetscInt *cone; 1700 PetscInt coneSize, c; 1701 1702 if ((point < cStart) || (point >= cMax)) continue; 1703 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1704 if (val != -1) continue; 1705 again = again == 1 ? 1 : 2; 1706 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 1707 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 1708 for (c = 0; c < coneSize; ++c) { 1709 ierr = DMLabelGetValue(label, cone[c], &val);CHKERRQ(ierr); 1710 if (val != -1) { 1711 const PetscInt *ccone; 1712 PetscInt cconeSize, cc, side; 1713 1714 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); 1715 if (val > 0) side = 1; 1716 else side = -1; 1717 ierr = DMLabelSetValue(label, point, side*(shift+dim));CHKERRQ(ierr); 1718 /* Mark cell faces which touch the fault */ 1719 ierr = DMPlexGetConeSize(dm, point, &cconeSize);CHKERRQ(ierr); 1720 ierr = DMPlexGetCone(dm, point, &ccone);CHKERRQ(ierr); 1721 for (cc = 0; cc < cconeSize; ++cc) { 1722 PetscInt *closure = NULL; 1723 PetscInt closureSize, cl; 1724 1725 ierr = DMLabelGetValue(label, ccone[cc], &val);CHKERRQ(ierr); 1726 if (val != -1) continue; 1727 ierr = DMPlexGetTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1728 for (cl = 0; cl < closureSize*2; cl += 2) { 1729 const PetscInt clp = closure[cl]; 1730 1731 ierr = DMLabelGetValue(label, clp, &val);CHKERRQ(ierr); 1732 if (val == -1) continue; 1733 ierr = DMLabelSetValue(label, ccone[cc], side*(shift+dim-1));CHKERRQ(ierr); 1734 break; 1735 } 1736 ierr = DMPlexRestoreTransitiveClosure(dm, ccone[cc], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1737 } 1738 again = 1; 1739 break; 1740 } 1741 } 1742 } 1743 } 1744 /* Classify the rest by cell membership */ 1745 for (s = 0; s < starSize*2; s += 2) { 1746 const PetscInt point = star[s]; 1747 1748 ierr = DMLabelGetValue(label, point, &val);CHKERRQ(ierr); 1749 if (val == -1) { 1750 PetscInt *sstar = NULL; 1751 PetscInt sstarSize, ss; 1752 PetscBool marked = PETSC_FALSE; 1753 1754 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1755 for (ss = 0; ss < sstarSize*2; ss += 2) { 1756 const PetscInt spoint = sstar[ss]; 1757 1758 if ((spoint < cStart) || (spoint >= cMax)) continue; 1759 ierr = DMLabelGetValue(label, spoint, &val);CHKERRQ(ierr); 1760 if (val == -1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d in star of %d does not have a valid label", spoint, point); 1761 ierr = DMLabelGetValue(depthLabel, point, &dep);CHKERRQ(ierr); 1762 if (val > 0) { 1763 ierr = DMLabelSetValue(label, point, shift+dep);CHKERRQ(ierr); 1764 } else { 1765 ierr = DMLabelSetValue(label, point, -(shift+dep));CHKERRQ(ierr); 1766 } 1767 marked = PETSC_TRUE; 1768 break; 1769 } 1770 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &sstarSize, &sstar);CHKERRQ(ierr); 1771 ierr = DMLabelGetValue(depthLabel, point, &dep);CHKERRQ(ierr); 1772 if (point < pMax[dep] && !marked) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Point %d could not be classified", point); 1773 } 1774 } 1775 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1776 } 1777 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1778 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1779 /* If any faces touching the fault divide cells on either side, split them */ 1780 ierr = DMLabelGetStratumIS(label, shift+dim-1, &facePosIS);CHKERRQ(ierr); 1781 ierr = DMLabelGetStratumIS(label, -(shift+dim-1), &faceNegIS);CHKERRQ(ierr); 1782 ierr = ISExpand(facePosIS, faceNegIS, &dimIS);CHKERRQ(ierr); 1783 ierr = ISDestroy(&facePosIS);CHKERRQ(ierr); 1784 ierr = ISDestroy(&faceNegIS);CHKERRQ(ierr); 1785 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1786 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1787 for (p = 0; p < numPoints; ++p) { 1788 const PetscInt point = points[p]; 1789 const PetscInt *support; 1790 PetscInt supportSize, valA, valB; 1791 1792 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1793 if (supportSize != 2) continue; 1794 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1795 ierr = DMLabelGetValue(label, support[0], &valA);CHKERRQ(ierr); 1796 ierr = DMLabelGetValue(label, support[1], &valB);CHKERRQ(ierr); 1797 if ((valA == -1) || (valB == -1)) continue; 1798 if (valA*valB > 0) continue; 1799 /* Split the face */ 1800 ierr = DMLabelGetValue(label, point, &valA);CHKERRQ(ierr); 1801 ierr = DMLabelClearValue(label, point, valA);CHKERRQ(ierr); 1802 ierr = DMLabelSetValue(label, point, dim-1);CHKERRQ(ierr); 1803 /* Label its closure: 1804 unmarked: label as unsplit 1805 incident: relabel as split 1806 split: do nothing 1807 */ 1808 { 1809 PetscInt *closure = NULL; 1810 PetscInt closureSize, cl; 1811 1812 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1813 for (cl = 0; cl < closureSize*2; cl += 2) { 1814 ierr = DMLabelGetValue(label, closure[cl], &valA);CHKERRQ(ierr); 1815 if (valA == -1) { /* Mark as unsplit */ 1816 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 1817 ierr = DMLabelSetValue(label, closure[cl], shift2+dep);CHKERRQ(ierr); 1818 } else if (((valA >= shift) && (valA < shift2)) || ((valA <= -shift) && (valA > -shift2))) { 1819 ierr = DMLabelGetValue(depthLabel, closure[cl], &dep);CHKERRQ(ierr); 1820 ierr = DMLabelClearValue(label, closure[cl], valA);CHKERRQ(ierr); 1821 ierr = DMLabelSetValue(label, closure[cl], dep);CHKERRQ(ierr); 1822 } 1823 } 1824 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1825 } 1826 } 1827 ierr = ISRestoreIndices(dimIS, &points);CHKERRQ(ierr); 1828 ierr = ISDestroy(&dimIS);CHKERRQ(ierr); 1829 ierr = PetscFree(pMax);CHKERRQ(ierr); 1830 PetscFunctionReturn(0); 1831 } 1832 1833 /*@ 1834 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 1835 1836 Collective on dm 1837 1838 Input Parameters: 1839 + dm - The original DM 1840 - labelName - The label specifying the interface vertices 1841 1842 Output Parameters: 1843 + hybridLabel - The label fully marking the interface 1844 - dmHybrid - The new DM 1845 1846 Level: developer 1847 1848 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate() 1849 @*/ 1850 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid) 1851 { 1852 DM idm; 1853 DMLabel subpointMap, hlabel; 1854 PetscInt dim; 1855 PetscErrorCode ierr; 1856 1857 PetscFunctionBegin; 1858 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1859 if (hybridLabel) PetscValidPointer(hybridLabel, 3); 1860 PetscValidPointer(dmHybrid, 4); 1861 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1862 ierr = DMPlexCreateSubmesh(dm, label, 1, &idm);CHKERRQ(ierr); 1863 ierr = DMPlexOrient(idm);CHKERRQ(ierr); 1864 ierr = DMPlexGetSubpointMap(idm, &subpointMap);CHKERRQ(ierr); 1865 ierr = DMLabelDuplicate(subpointMap, &hlabel);CHKERRQ(ierr); 1866 ierr = DMLabelClearStratum(hlabel, dim);CHKERRQ(ierr); 1867 ierr = DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);CHKERRQ(ierr); 1868 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1869 ierr = DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);CHKERRQ(ierr); 1870 if (hybridLabel) *hybridLabel = hlabel; 1871 else {ierr = DMLabelDestroy(&hlabel);CHKERRQ(ierr);} 1872 PetscFunctionReturn(0); 1873 } 1874 1875 /* Here we need the explicit assumption that: 1876 1877 For any marked cell, the marked vertices constitute a single face 1878 */ 1879 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 1880 { 1881 IS subvertexIS = NULL; 1882 const PetscInt *subvertices; 1883 PetscInt *pStart, *pEnd, *pMax, pSize; 1884 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 1885 PetscErrorCode ierr; 1886 1887 PetscFunctionBegin; 1888 *numFaces = 0; 1889 *nFV = 0; 1890 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1891 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1892 pSize = PetscMax(depth, dim) + 1; 1893 ierr = PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);CHKERRQ(ierr); 1894 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1895 for (d = 0; d <= depth; ++d) { 1896 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1897 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1898 } 1899 /* Loop over initial vertices and mark all faces in the collective star() */ 1900 if (vertexLabel) {ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr);} 1901 if (subvertexIS) { 1902 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1903 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1904 } 1905 for (v = 0; v < numSubVerticesInitial; ++v) { 1906 const PetscInt vertex = subvertices[v]; 1907 PetscInt *star = NULL; 1908 PetscInt starSize, s, numCells = 0, c; 1909 1910 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1911 for (s = 0; s < starSize*2; s += 2) { 1912 const PetscInt point = star[s]; 1913 if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point; 1914 } 1915 for (c = 0; c < numCells; ++c) { 1916 const PetscInt cell = star[c]; 1917 PetscInt *closure = NULL; 1918 PetscInt closureSize, cl; 1919 PetscInt cellLoc, numCorners = 0, faceSize = 0; 1920 1921 ierr = DMLabelGetValue(subpointMap, cell, &cellLoc);CHKERRQ(ierr); 1922 if (cellLoc == 2) continue; 1923 if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 1924 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1925 for (cl = 0; cl < closureSize*2; cl += 2) { 1926 const PetscInt point = closure[cl]; 1927 PetscInt vertexLoc; 1928 1929 if ((point >= pStart[0]) && (point < pEnd[0])) { 1930 ++numCorners; 1931 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 1932 if (vertexLoc == value) closure[faceSize++] = point; 1933 } 1934 } 1935 if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);} 1936 if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 1937 if (faceSize == *nFV) { 1938 const PetscInt *cells = NULL; 1939 PetscInt numCells, nc; 1940 1941 ++(*numFaces); 1942 for (cl = 0; cl < faceSize; ++cl) { 1943 ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr); 1944 } 1945 ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1946 for (nc = 0; nc < numCells; ++nc) { 1947 ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr); 1948 } 1949 ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1950 } 1951 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1952 } 1953 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1954 } 1955 if (subvertexIS) { 1956 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1957 } 1958 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 1959 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 1960 PetscFunctionReturn(0); 1961 } 1962 1963 static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm) 1964 { 1965 IS subvertexIS = NULL; 1966 const PetscInt *subvertices; 1967 PetscInt *pStart, *pEnd, *pMax; 1968 PetscInt dim, d, numSubVerticesInitial = 0, v; 1969 PetscErrorCode ierr; 1970 1971 PetscFunctionBegin; 1972 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1973 ierr = PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);CHKERRQ(ierr); 1974 ierr = DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1975 for (d = 0; d <= dim; ++d) { 1976 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1977 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1978 } 1979 /* Loop over initial vertices and mark all faces in the collective star() */ 1980 if (vertexLabel) { 1981 ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr); 1982 if (subvertexIS) { 1983 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1984 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1985 } 1986 } 1987 for (v = 0; v < numSubVerticesInitial; ++v) { 1988 const PetscInt vertex = subvertices[v]; 1989 PetscInt *star = NULL; 1990 PetscInt starSize, s, numFaces = 0, f; 1991 1992 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1993 for (s = 0; s < starSize*2; s += 2) { 1994 const PetscInt point = star[s]; 1995 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point; 1996 } 1997 for (f = 0; f < numFaces; ++f) { 1998 const PetscInt face = star[f]; 1999 PetscInt *closure = NULL; 2000 PetscInt closureSize, c; 2001 PetscInt faceLoc; 2002 2003 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 2004 if (faceLoc == dim-1) continue; 2005 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 2006 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2007 for (c = 0; c < closureSize*2; c += 2) { 2008 const PetscInt point = closure[c]; 2009 PetscInt vertexLoc; 2010 2011 if ((point >= pStart[0]) && (point < pEnd[0])) { 2012 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2013 if (vertexLoc != value) break; 2014 } 2015 } 2016 if (c == closureSize*2) { 2017 const PetscInt *support; 2018 PetscInt supportSize, s; 2019 2020 for (c = 0; c < closureSize*2; c += 2) { 2021 const PetscInt point = closure[c]; 2022 2023 for (d = 0; d < dim; ++d) { 2024 if ((point >= pStart[d]) && (point < pEnd[d])) { 2025 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2026 break; 2027 } 2028 } 2029 } 2030 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2031 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2032 for (s = 0; s < supportSize; ++s) { 2033 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2034 } 2035 } 2036 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2037 } 2038 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2039 } 2040 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr);} 2041 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2042 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 2047 { 2048 DMLabel label = NULL; 2049 const PetscInt *cone; 2050 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1; 2051 PetscErrorCode ierr; 2052 2053 PetscFunctionBegin; 2054 *numFaces = 0; 2055 *nFV = 0; 2056 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 2057 *subCells = NULL; 2058 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2059 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2060 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2061 if (cMax < 0) PetscFunctionReturn(0); 2062 if (label) { 2063 for (c = cMax; c < cEnd; ++c) { 2064 PetscInt val; 2065 2066 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2067 if (val == value) { 2068 ++(*numFaces); 2069 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2070 } 2071 } 2072 } else { 2073 *numFaces = cEnd - cMax; 2074 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 2075 } 2076 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 2077 if (!(*numFaces)) PetscFunctionReturn(0); 2078 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 2079 for (c = cMax; c < cEnd; ++c) { 2080 const PetscInt *cells; 2081 PetscInt numCells; 2082 2083 if (label) { 2084 PetscInt val; 2085 2086 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2087 if (val != value) continue; 2088 } 2089 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2090 for (p = 0; p < *nFV; ++p) { 2091 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 2092 } 2093 /* Negative face */ 2094 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2095 /* Not true in parallel 2096 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2097 for (p = 0; p < numCells; ++p) { 2098 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 2099 (*subCells)[subc++] = cells[p]; 2100 } 2101 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2102 /* Positive face is not included */ 2103 } 2104 PetscFunctionReturn(0); 2105 } 2106 2107 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 2108 { 2109 PetscInt *pStart, *pEnd; 2110 PetscInt dim, cMax, cEnd, c, d; 2111 PetscErrorCode ierr; 2112 2113 PetscFunctionBegin; 2114 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2115 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2116 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2117 if (cMax < 0) PetscFunctionReturn(0); 2118 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 2119 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 2120 for (c = cMax; c < cEnd; ++c) { 2121 const PetscInt *cone; 2122 PetscInt *closure = NULL; 2123 PetscInt fconeSize, coneSize, closureSize, cl, val; 2124 2125 if (label) { 2126 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2127 if (val != value) continue; 2128 } 2129 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2130 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2131 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 2132 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2133 /* Negative face */ 2134 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2135 for (cl = 0; cl < closureSize*2; cl += 2) { 2136 const PetscInt point = closure[cl]; 2137 2138 for (d = 0; d <= dim; ++d) { 2139 if ((point >= pStart[d]) && (point < pEnd[d])) { 2140 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2141 break; 2142 } 2143 } 2144 } 2145 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2146 /* Cells -- positive face is not included */ 2147 for (cl = 0; cl < 1; ++cl) { 2148 const PetscInt *support; 2149 PetscInt supportSize, s; 2150 2151 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 2152 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2153 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 2154 for (s = 0; s < supportSize; ++s) { 2155 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2156 } 2157 } 2158 } 2159 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2160 PetscFunctionReturn(0); 2161 } 2162 2163 static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2164 { 2165 MPI_Comm comm; 2166 PetscBool posOrient = PETSC_FALSE; 2167 const PetscInt debug = 0; 2168 PetscInt cellDim, faceSize, f; 2169 PetscErrorCode ierr; 2170 2171 PetscFunctionBegin; 2172 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2173 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 2174 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 2175 2176 if (cellDim == 1 && numCorners == 2) { 2177 /* Triangle */ 2178 faceSize = numCorners-1; 2179 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2180 } else if (cellDim == 2 && numCorners == 3) { 2181 /* Triangle */ 2182 faceSize = numCorners-1; 2183 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2184 } else if (cellDim == 3 && numCorners == 4) { 2185 /* Tetrahedron */ 2186 faceSize = numCorners-1; 2187 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2188 } else if (cellDim == 1 && numCorners == 3) { 2189 /* Quadratic line */ 2190 faceSize = 1; 2191 posOrient = PETSC_TRUE; 2192 } else if (cellDim == 2 && numCorners == 4) { 2193 /* Quads */ 2194 faceSize = 2; 2195 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2196 posOrient = PETSC_TRUE; 2197 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2198 posOrient = PETSC_TRUE; 2199 } else { 2200 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2201 posOrient = PETSC_FALSE; 2202 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2203 } 2204 } else if (cellDim == 2 && numCorners == 6) { 2205 /* Quadratic triangle (I hate this) */ 2206 /* Edges are determined by the first 2 vertices (corners of edges) */ 2207 const PetscInt faceSizeTri = 3; 2208 PetscInt sortedIndices[3], i, iFace; 2209 PetscBool found = PETSC_FALSE; 2210 PetscInt faceVerticesTriSorted[9] = { 2211 0, 3, 4, /* bottom */ 2212 1, 4, 5, /* right */ 2213 2, 3, 5, /* left */ 2214 }; 2215 PetscInt faceVerticesTri[9] = { 2216 0, 3, 4, /* bottom */ 2217 1, 4, 5, /* right */ 2218 2, 5, 3, /* left */ 2219 }; 2220 2221 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2222 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2223 for (iFace = 0; iFace < 3; ++iFace) { 2224 const PetscInt ii = iFace*faceSizeTri; 2225 PetscInt fVertex, cVertex; 2226 2227 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2228 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2229 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2230 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2231 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2232 faceVertices[fVertex] = origVertices[cVertex]; 2233 break; 2234 } 2235 } 2236 } 2237 found = PETSC_TRUE; 2238 break; 2239 } 2240 } 2241 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2242 if (posOriented) *posOriented = PETSC_TRUE; 2243 PetscFunctionReturn(0); 2244 } else if (cellDim == 2 && numCorners == 9) { 2245 /* Quadratic quad (I hate this) */ 2246 /* Edges are determined by the first 2 vertices (corners of edges) */ 2247 const PetscInt faceSizeQuad = 3; 2248 PetscInt sortedIndices[3], i, iFace; 2249 PetscBool found = PETSC_FALSE; 2250 PetscInt faceVerticesQuadSorted[12] = { 2251 0, 1, 4, /* bottom */ 2252 1, 2, 5, /* right */ 2253 2, 3, 6, /* top */ 2254 0, 3, 7, /* left */ 2255 }; 2256 PetscInt faceVerticesQuad[12] = { 2257 0, 1, 4, /* bottom */ 2258 1, 2, 5, /* right */ 2259 2, 3, 6, /* top */ 2260 3, 0, 7, /* left */ 2261 }; 2262 2263 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2264 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2265 for (iFace = 0; iFace < 4; ++iFace) { 2266 const PetscInt ii = iFace*faceSizeQuad; 2267 PetscInt fVertex, cVertex; 2268 2269 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2270 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2271 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2272 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2273 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2274 faceVertices[fVertex] = origVertices[cVertex]; 2275 break; 2276 } 2277 } 2278 } 2279 found = PETSC_TRUE; 2280 break; 2281 } 2282 } 2283 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2284 if (posOriented) *posOriented = PETSC_TRUE; 2285 PetscFunctionReturn(0); 2286 } else if (cellDim == 3 && numCorners == 8) { 2287 /* Hexes 2288 A hex is two oriented quads with the normal of the first 2289 pointing up at the second. 2290 2291 7---6 2292 /| /| 2293 4---5 | 2294 | 1-|-2 2295 |/ |/ 2296 0---3 2297 2298 Faces are determined by the first 4 vertices (corners of faces) */ 2299 const PetscInt faceSizeHex = 4; 2300 PetscInt sortedIndices[4], i, iFace; 2301 PetscBool found = PETSC_FALSE; 2302 PetscInt faceVerticesHexSorted[24] = { 2303 0, 1, 2, 3, /* bottom */ 2304 4, 5, 6, 7, /* top */ 2305 0, 3, 4, 5, /* front */ 2306 2, 3, 5, 6, /* right */ 2307 1, 2, 6, 7, /* back */ 2308 0, 1, 4, 7, /* left */ 2309 }; 2310 PetscInt faceVerticesHex[24] = { 2311 1, 2, 3, 0, /* bottom */ 2312 4, 5, 6, 7, /* top */ 2313 0, 3, 5, 4, /* front */ 2314 3, 2, 6, 5, /* right */ 2315 2, 1, 7, 6, /* back */ 2316 1, 0, 4, 7, /* left */ 2317 }; 2318 2319 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2320 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2321 for (iFace = 0; iFace < 6; ++iFace) { 2322 const PetscInt ii = iFace*faceSizeHex; 2323 PetscInt fVertex, cVertex; 2324 2325 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2326 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2327 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2328 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2329 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2330 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2331 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2332 faceVertices[fVertex] = origVertices[cVertex]; 2333 break; 2334 } 2335 } 2336 } 2337 found = PETSC_TRUE; 2338 break; 2339 } 2340 } 2341 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2342 if (posOriented) *posOriented = PETSC_TRUE; 2343 PetscFunctionReturn(0); 2344 } else if (cellDim == 3 && numCorners == 10) { 2345 /* Quadratic tet */ 2346 /* Faces are determined by the first 3 vertices (corners of faces) */ 2347 const PetscInt faceSizeTet = 6; 2348 PetscInt sortedIndices[6], i, iFace; 2349 PetscBool found = PETSC_FALSE; 2350 PetscInt faceVerticesTetSorted[24] = { 2351 0, 1, 2, 6, 7, 8, /* bottom */ 2352 0, 3, 4, 6, 7, 9, /* front */ 2353 1, 4, 5, 7, 8, 9, /* right */ 2354 2, 3, 5, 6, 8, 9, /* left */ 2355 }; 2356 PetscInt faceVerticesTet[24] = { 2357 0, 1, 2, 6, 7, 8, /* bottom */ 2358 0, 4, 3, 6, 7, 9, /* front */ 2359 1, 5, 4, 7, 8, 9, /* right */ 2360 2, 3, 5, 8, 6, 9, /* left */ 2361 }; 2362 2363 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2364 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2365 for (iFace=0; iFace < 4; ++iFace) { 2366 const PetscInt ii = iFace*faceSizeTet; 2367 PetscInt fVertex, cVertex; 2368 2369 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2370 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2371 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2372 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2373 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2374 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2375 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2376 faceVertices[fVertex] = origVertices[cVertex]; 2377 break; 2378 } 2379 } 2380 } 2381 found = PETSC_TRUE; 2382 break; 2383 } 2384 } 2385 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2386 if (posOriented) *posOriented = PETSC_TRUE; 2387 PetscFunctionReturn(0); 2388 } else if (cellDim == 3 && numCorners == 27) { 2389 /* Quadratic hexes (I hate this) 2390 A hex is two oriented quads with the normal of the first 2391 pointing up at the second. 2392 2393 7---6 2394 /| /| 2395 4---5 | 2396 | 3-|-2 2397 |/ |/ 2398 0---1 2399 2400 Faces are determined by the first 4 vertices (corners of faces) */ 2401 const PetscInt faceSizeQuadHex = 9; 2402 PetscInt sortedIndices[9], i, iFace; 2403 PetscBool found = PETSC_FALSE; 2404 PetscInt faceVerticesQuadHexSorted[54] = { 2405 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2406 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2407 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2408 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2409 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2410 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2411 }; 2412 PetscInt faceVerticesQuadHex[54] = { 2413 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2414 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2415 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2416 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2417 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2418 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2419 }; 2420 2421 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2422 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2423 for (iFace = 0; iFace < 6; ++iFace) { 2424 const PetscInt ii = iFace*faceSizeQuadHex; 2425 PetscInt fVertex, cVertex; 2426 2427 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2428 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2429 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2430 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2431 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2432 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2433 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2434 faceVertices[fVertex] = origVertices[cVertex]; 2435 break; 2436 } 2437 } 2438 } 2439 found = PETSC_TRUE; 2440 break; 2441 } 2442 } 2443 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2444 if (posOriented) *posOriented = PETSC_TRUE; 2445 PetscFunctionReturn(0); 2446 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2447 if (!posOrient) { 2448 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2449 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2450 } else { 2451 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2452 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2453 } 2454 if (posOriented) *posOriented = posOrient; 2455 PetscFunctionReturn(0); 2456 } 2457 2458 /*@ 2459 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 2460 in faceVertices. The orientation is such that the face normal points out of the cell 2461 2462 Not collective 2463 2464 Input Parameters: 2465 + dm - The original mesh 2466 . cell - The cell mesh point 2467 . faceSize - The number of vertices on the face 2468 . face - The face vertices 2469 . numCorners - The number of vertices on the cell 2470 . indices - Local numbering of face vertices in cell cone 2471 - origVertices - Original face vertices 2472 2473 Output Parameter: 2474 + faceVertices - The face vertices properly oriented 2475 - posOriented - PETSC_TRUE if the face was oriented with outward normal 2476 2477 Level: developer 2478 2479 .seealso: DMPlexGetCone() 2480 @*/ 2481 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2482 { 2483 const PetscInt *cone = NULL; 2484 PetscInt coneSize, v, f, v2; 2485 PetscInt oppositeVertex = -1; 2486 PetscErrorCode ierr; 2487 2488 PetscFunctionBegin; 2489 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2490 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2491 for (v = 0, v2 = 0; v < coneSize; ++v) { 2492 PetscBool found = PETSC_FALSE; 2493 2494 for (f = 0; f < faceSize; ++f) { 2495 if (face[f] == cone[v]) { 2496 found = PETSC_TRUE; break; 2497 } 2498 } 2499 if (found) { 2500 indices[v2] = v; 2501 origVertices[v2] = cone[v]; 2502 ++v2; 2503 } else { 2504 oppositeVertex = v; 2505 } 2506 } 2507 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2508 PetscFunctionReturn(0); 2509 } 2510 2511 /* 2512 DMPlexInsertFace_Internal - Puts a face into the mesh 2513 2514 Not collective 2515 2516 Input Parameters: 2517 + dm - The DMPlex 2518 . numFaceVertex - The number of vertices in the face 2519 . faceVertices - The vertices in the face for dm 2520 . subfaceVertices - The vertices in the face for subdm 2521 . numCorners - The number of vertices in the cell 2522 . cell - A cell in dm containing the face 2523 . subcell - A cell in subdm containing the face 2524 . firstFace - First face in the mesh 2525 - newFacePoint - Next face in the mesh 2526 2527 Output Parameters: 2528 . newFacePoint - Contains next face point number on input, updated on output 2529 2530 Level: developer 2531 */ 2532 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) 2533 { 2534 MPI_Comm comm; 2535 DM_Plex *submesh = (DM_Plex*) subdm->data; 2536 const PetscInt *faces; 2537 PetscInt numFaces, coneSize; 2538 PetscErrorCode ierr; 2539 2540 PetscFunctionBegin; 2541 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2542 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2543 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2544 #if 0 2545 /* Cannot use this because support() has not been constructed yet */ 2546 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2547 #else 2548 { 2549 PetscInt f; 2550 2551 numFaces = 0; 2552 ierr = DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2553 for (f = firstFace; f < *newFacePoint; ++f) { 2554 PetscInt dof, off, d; 2555 2556 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2557 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2558 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2559 for (d = 0; d < dof; ++d) { 2560 const PetscInt p = submesh->cones[off+d]; 2561 PetscInt v; 2562 2563 for (v = 0; v < numFaceVertices; ++v) { 2564 if (subfaceVertices[v] == p) break; 2565 } 2566 if (v == numFaceVertices) break; 2567 } 2568 if (d == dof) { 2569 numFaces = 1; 2570 ((PetscInt*) faces)[0] = f; 2571 } 2572 } 2573 } 2574 #endif 2575 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2576 else if (numFaces == 1) { 2577 /* Add the other cell neighbor for this face */ 2578 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2579 } else { 2580 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2581 PetscBool posOriented; 2582 2583 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2584 origVertices = &orientedVertices[numFaceVertices]; 2585 indices = &orientedVertices[numFaceVertices*2]; 2586 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2587 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2588 /* TODO: I know that routine should return a permutation, not the indices */ 2589 for (v = 0; v < numFaceVertices; ++v) { 2590 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2591 for (ov = 0; ov < numFaceVertices; ++ov) { 2592 if (orientedVertices[ov] == vertex) { 2593 orientedSubVertices[ov] = subvertex; 2594 break; 2595 } 2596 } 2597 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2598 } 2599 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2600 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2601 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2602 ++(*newFacePoint); 2603 } 2604 #if 0 2605 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2606 #else 2607 ierr = DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2608 #endif 2609 PetscFunctionReturn(0); 2610 } 2611 2612 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2613 { 2614 MPI_Comm comm; 2615 DMLabel subpointMap; 2616 IS subvertexIS, subcellIS; 2617 const PetscInt *subVertices, *subCells; 2618 PetscInt numSubVertices, firstSubVertex, numSubCells; 2619 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2620 PetscInt vStart, vEnd, c, f; 2621 PetscErrorCode ierr; 2622 2623 PetscFunctionBegin; 2624 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2625 /* Create subpointMap which marks the submesh */ 2626 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2627 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2628 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2629 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2630 /* Setup chart */ 2631 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2632 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2633 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2634 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2635 /* Set cone sizes */ 2636 firstSubVertex = numSubCells; 2637 firstSubFace = numSubCells+numSubVertices; 2638 newFacePoint = firstSubFace; 2639 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2640 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2641 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2642 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2643 for (c = 0; c < numSubCells; ++c) { 2644 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2645 } 2646 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2647 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2648 } 2649 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2650 /* Create face cones */ 2651 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2652 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2653 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2654 for (c = 0; c < numSubCells; ++c) { 2655 const PetscInt cell = subCells[c]; 2656 const PetscInt subcell = c; 2657 PetscInt *closure = NULL; 2658 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2659 2660 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2661 for (cl = 0; cl < closureSize*2; cl += 2) { 2662 const PetscInt point = closure[cl]; 2663 PetscInt subVertex; 2664 2665 if ((point >= vStart) && (point < vEnd)) { 2666 ++numCorners; 2667 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2668 if (subVertex >= 0) { 2669 closure[faceSize] = point; 2670 subface[faceSize] = firstSubVertex+subVertex; 2671 ++faceSize; 2672 } 2673 } 2674 } 2675 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2676 if (faceSize == nFV) { 2677 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2678 } 2679 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2680 } 2681 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2682 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2683 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2684 /* Build coordinates */ 2685 { 2686 PetscSection coordSection, subCoordSection; 2687 Vec coordinates, subCoordinates; 2688 PetscScalar *coords, *subCoords; 2689 PetscInt numComp, coordSize, v; 2690 const char *name; 2691 2692 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2693 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2694 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2695 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2696 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2697 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2698 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2699 for (v = 0; v < numSubVertices; ++v) { 2700 const PetscInt vertex = subVertices[v]; 2701 const PetscInt subvertex = firstSubVertex+v; 2702 PetscInt dof; 2703 2704 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2705 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2706 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2707 } 2708 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2709 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2710 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 2711 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 2712 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 2713 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2714 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2715 if (coordSize) { 2716 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2717 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2718 for (v = 0; v < numSubVertices; ++v) { 2719 const PetscInt vertex = subVertices[v]; 2720 const PetscInt subvertex = firstSubVertex+v; 2721 PetscInt dof, off, sdof, soff, d; 2722 2723 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2724 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2725 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2726 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2727 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2728 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2729 } 2730 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2731 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2732 } 2733 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2734 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2735 } 2736 /* Cleanup */ 2737 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2738 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2739 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2740 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 2741 PetscFunctionReturn(0); 2742 } 2743 2744 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 2745 { 2746 PetscInt subPoint; 2747 PetscErrorCode ierr; 2748 2749 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 2750 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 2751 } 2752 2753 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 2754 { 2755 MPI_Comm comm; 2756 DMLabel subpointMap; 2757 IS *subpointIS; 2758 const PetscInt **subpoints; 2759 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 2760 PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v; 2761 PetscMPIInt rank; 2762 PetscErrorCode ierr; 2763 2764 PetscFunctionBegin; 2765 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2766 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2767 /* Create subpointMap which marks the submesh */ 2768 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2769 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2770 if (cellHeight) { 2771 if (isCohesive) {ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2772 else {ierr = DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2773 } else { 2774 DMLabel depth; 2775 IS pointIS; 2776 const PetscInt *points; 2777 PetscInt numPoints; 2778 2779 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 2780 ierr = DMLabelGetStratumSize(label, value, &numPoints);CHKERRQ(ierr); 2781 ierr = DMLabelGetStratumIS(label, value, &pointIS);CHKERRQ(ierr); 2782 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2783 for (p = 0; p < numPoints; ++p) { 2784 PetscInt *closure = NULL; 2785 PetscInt closureSize, c, pdim; 2786 2787 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2788 for (c = 0; c < closureSize*2; c += 2) { 2789 ierr = DMLabelGetValue(depth, closure[c], &pdim);CHKERRQ(ierr); 2790 ierr = DMLabelSetValue(subpointMap, closure[c], pdim);CHKERRQ(ierr); 2791 } 2792 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2793 } 2794 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2795 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2796 } 2797 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2798 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2799 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2800 cMax = (cMax < 0) ? cEnd : cMax; 2801 /* Setup chart */ 2802 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2803 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2804 for (d = 0; d <= dim; ++d) { 2805 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2806 totSubPoints += numSubPoints[d]; 2807 } 2808 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2809 ierr = DMPlexSetVTKCellHeight(subdm, cellHeight);CHKERRQ(ierr); 2810 /* Set cone sizes */ 2811 firstSubPoint[dim] = 0; 2812 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2813 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2814 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2815 for (d = 0; d <= dim; ++d) { 2816 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2817 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2818 } 2819 for (d = 0; d <= dim; ++d) { 2820 for (p = 0; p < numSubPoints[d]; ++p) { 2821 const PetscInt point = subpoints[d][p]; 2822 const PetscInt subpoint = firstSubPoint[d] + p; 2823 const PetscInt *cone; 2824 PetscInt coneSize, coneSizeNew, c, val; 2825 2826 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2827 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2828 if (cellHeight && (d == dim)) { 2829 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2830 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2831 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2832 if (val >= 0) coneSizeNew++; 2833 } 2834 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2835 } 2836 } 2837 } 2838 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2839 /* Set cones */ 2840 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2841 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 2842 for (d = 0; d <= dim; ++d) { 2843 for (p = 0; p < numSubPoints[d]; ++p) { 2844 const PetscInt point = subpoints[d][p]; 2845 const PetscInt subpoint = firstSubPoint[d] + p; 2846 const PetscInt *cone, *ornt; 2847 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 2848 2849 if (d == dim-1) { 2850 const PetscInt *support, *cone, *ornt; 2851 PetscInt supportSize, coneSize, s, subc; 2852 2853 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 2854 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 2855 for (s = 0; s < supportSize; ++s) { 2856 if ((support[s] < cMax) || (support[s] >= cEnd)) continue; 2857 ierr = PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);CHKERRQ(ierr); 2858 if (subc >= 0) { 2859 const PetscInt ccell = subpoints[d+1][subc]; 2860 2861 ierr = DMPlexGetCone(dm, ccell, &cone);CHKERRQ(ierr); 2862 ierr = DMPlexGetConeSize(dm, ccell, &coneSize);CHKERRQ(ierr); 2863 ierr = DMPlexGetConeOrientation(dm, ccell, &ornt);CHKERRQ(ierr); 2864 for (c = 0; c < coneSize; ++c) { 2865 if (cone[c] == point) { 2866 fornt = ornt[c]; 2867 break; 2868 } 2869 } 2870 break; 2871 } 2872 } 2873 } 2874 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2875 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 2876 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2877 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 2878 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2879 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 2880 if (subc >= 0) { 2881 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 2882 orntNew[coneSizeNew] = ornt[c]; 2883 ++coneSizeNew; 2884 } 2885 } 2886 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 2887 if (fornt < 0) { 2888 /* This should be replaced by a call to DMPlexReverseCell() */ 2889 #if 0 2890 ierr = DMPlexReverseCell(subdm, subpoint);CHKERRQ(ierr); 2891 #else 2892 for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) { 2893 PetscInt faceSize, tmp; 2894 2895 tmp = coneNew[c]; 2896 coneNew[c] = coneNew[coneSizeNew-1-c]; 2897 coneNew[coneSizeNew-1-c] = tmp; 2898 ierr = DMPlexGetConeSize(dm, cone[c], &faceSize);CHKERRQ(ierr); 2899 tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c]; 2900 orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c]; 2901 orntNew[coneSizeNew-1-c] = tmp; 2902 } 2903 } 2904 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 2905 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 2906 #endif 2907 } 2908 } 2909 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 2910 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2911 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2912 /* Build coordinates */ 2913 { 2914 PetscSection coordSection, subCoordSection; 2915 Vec coordinates, subCoordinates; 2916 PetscScalar *coords, *subCoords; 2917 PetscInt cdim, numComp, coordSize; 2918 const char *name; 2919 2920 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 2921 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2922 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2923 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2924 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2925 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2926 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2927 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 2928 for (v = 0; v < numSubPoints[0]; ++v) { 2929 const PetscInt vertex = subpoints[0][v]; 2930 const PetscInt subvertex = firstSubPoint[0]+v; 2931 PetscInt dof; 2932 2933 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2934 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2935 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2936 } 2937 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2938 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2939 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 2940 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 2941 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 2942 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2943 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 2944 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2945 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2946 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2947 for (v = 0; v < numSubPoints[0]; ++v) { 2948 const PetscInt vertex = subpoints[0][v]; 2949 const PetscInt subvertex = firstSubPoint[0]+v; 2950 PetscInt dof, off, sdof, soff, d; 2951 2952 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2953 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2954 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2955 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2956 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2957 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2958 } 2959 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2960 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2961 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2962 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2963 } 2964 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 2965 { 2966 PetscSF sfPoint, sfPointSub; 2967 IS subpIS; 2968 const PetscSFNode *remotePoints; 2969 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 2970 const PetscInt *localPoints, *subpoints; 2971 PetscInt *slocalPoints; 2972 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 2973 PetscMPIInt rank; 2974 2975 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2976 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2977 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 2978 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2979 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 2980 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 2981 if (subpIS) { 2982 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 2983 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 2984 } 2985 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2986 if (numRoots >= 0) { 2987 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 2988 for (p = 0; p < pEnd-pStart; ++p) { 2989 newLocalPoints[p].rank = -2; 2990 newLocalPoints[p].index = -2; 2991 } 2992 /* Set subleaves */ 2993 for (l = 0; l < numLeaves; ++l) { 2994 const PetscInt point = localPoints[l]; 2995 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 2996 2997 if (subpoint < 0) continue; 2998 newLocalPoints[point-pStart].rank = rank; 2999 newLocalPoints[point-pStart].index = subpoint; 3000 ++numSubleaves; 3001 } 3002 /* Must put in owned subpoints */ 3003 for (p = pStart; p < pEnd; ++p) { 3004 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3005 3006 if (subpoint < 0) { 3007 newOwners[p-pStart].rank = -3; 3008 newOwners[p-pStart].index = -3; 3009 } else { 3010 newOwners[p-pStart].rank = rank; 3011 newOwners[p-pStart].index = subpoint; 3012 } 3013 } 3014 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3015 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3016 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3017 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3018 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 3019 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 3020 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3021 const PetscInt point = localPoints[l]; 3022 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3023 3024 if (subpoint < 0) continue; 3025 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3026 slocalPoints[sl] = subpoint; 3027 sremotePoints[sl].rank = newLocalPoints[point].rank; 3028 sremotePoints[sl].index = newLocalPoints[point].index; 3029 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3030 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3031 ++sl; 3032 } 3033 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3034 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3035 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3036 } 3037 if (subpIS) { 3038 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 3039 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 3040 } 3041 } 3042 /* Cleanup */ 3043 for (d = 0; d <= dim; ++d) { 3044 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3045 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 3046 } 3047 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 3048 PetscFunctionReturn(0); 3049 } 3050 3051 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 3052 { 3053 PetscErrorCode ierr; 3054 3055 PetscFunctionBegin; 3056 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);CHKERRQ(ierr); 3057 PetscFunctionReturn(0); 3058 } 3059 3060 /*@ 3061 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3062 3063 Input Parameters: 3064 + dm - The original mesh 3065 . vertexLabel - The DMLabel marking vertices contained in the surface 3066 - value - The label value to use 3067 3068 Output Parameter: 3069 . subdm - The surface mesh 3070 3071 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3072 3073 Level: developer 3074 3075 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3076 @*/ 3077 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 3078 { 3079 PetscInt dim, cdim, depth; 3080 PetscErrorCode ierr; 3081 3082 PetscFunctionBegin; 3083 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3084 PetscValidPointer(subdm, 3); 3085 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3086 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3087 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3088 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3089 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3090 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3091 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3092 if (depth == dim) { 3093 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3094 } else { 3095 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3096 } 3097 PetscFunctionReturn(0); 3098 } 3099 3100 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3101 { 3102 MPI_Comm comm; 3103 DMLabel subpointMap; 3104 IS subvertexIS; 3105 const PetscInt *subVertices; 3106 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3107 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3108 PetscInt cMax, c, f; 3109 PetscErrorCode ierr; 3110 3111 PetscFunctionBegin; 3112 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 3113 /* Create subpointMap which marks the submesh */ 3114 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 3115 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3116 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3117 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 3118 /* Setup chart */ 3119 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 3120 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 3121 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 3122 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 3123 /* Set cone sizes */ 3124 firstSubVertex = numSubCells; 3125 firstSubFace = numSubCells+numSubVertices; 3126 newFacePoint = firstSubFace; 3127 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 3128 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3129 for (c = 0; c < numSubCells; ++c) { 3130 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 3131 } 3132 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3133 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 3134 } 3135 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3136 /* Create face cones */ 3137 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3138 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3139 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 3140 for (c = 0; c < numSubCells; ++c) { 3141 const PetscInt cell = subCells[c]; 3142 const PetscInt subcell = c; 3143 const PetscInt *cone, *cells; 3144 PetscInt numCells, subVertex, p, v; 3145 3146 if (cell < cMax) continue; 3147 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 3148 for (v = 0; v < nFV; ++v) { 3149 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 3150 subface[v] = firstSubVertex+subVertex; 3151 } 3152 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 3153 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 3154 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3155 /* Not true in parallel 3156 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3157 for (p = 0; p < numCells; ++p) { 3158 PetscInt negsubcell; 3159 3160 if (cells[p] >= cMax) continue; 3161 /* I know this is a crap search */ 3162 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3163 if (subCells[negsubcell] == cells[p]) break; 3164 } 3165 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3166 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 3167 } 3168 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3169 ++newFacePoint; 3170 } 3171 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 3172 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3173 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3174 /* Build coordinates */ 3175 { 3176 PetscSection coordSection, subCoordSection; 3177 Vec coordinates, subCoordinates; 3178 PetscScalar *coords, *subCoords; 3179 PetscInt cdim, numComp, coordSize, v; 3180 const char *name; 3181 3182 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3183 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3184 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3185 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3186 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3187 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3188 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3189 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3190 for (v = 0; v < numSubVertices; ++v) { 3191 const PetscInt vertex = subVertices[v]; 3192 const PetscInt subvertex = firstSubVertex+v; 3193 PetscInt dof; 3194 3195 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3196 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3197 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3198 } 3199 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3200 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3201 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3202 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3203 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3204 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3205 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3206 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3207 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3208 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3209 for (v = 0; v < numSubVertices; ++v) { 3210 const PetscInt vertex = subVertices[v]; 3211 const PetscInt subvertex = firstSubVertex+v; 3212 PetscInt dof, off, sdof, soff, d; 3213 3214 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3215 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3216 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3217 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3218 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3219 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3220 } 3221 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3222 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3223 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3224 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3225 } 3226 /* Build SF */ 3227 CHKMEMQ; 3228 { 3229 PetscSF sfPoint, sfPointSub; 3230 const PetscSFNode *remotePoints; 3231 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3232 const PetscInt *localPoints; 3233 PetscInt *slocalPoints; 3234 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3235 PetscMPIInt rank; 3236 3237 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3238 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3239 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3240 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3241 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3242 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3243 if (numRoots >= 0) { 3244 /* Only vertices should be shared */ 3245 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3246 for (p = 0; p < pEnd-pStart; ++p) { 3247 newLocalPoints[p].rank = -2; 3248 newLocalPoints[p].index = -2; 3249 } 3250 /* Set subleaves */ 3251 for (l = 0; l < numLeaves; ++l) { 3252 const PetscInt point = localPoints[l]; 3253 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3254 3255 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3256 if (subPoint < 0) continue; 3257 newLocalPoints[point-pStart].rank = rank; 3258 newLocalPoints[point-pStart].index = subPoint; 3259 ++numSubLeaves; 3260 } 3261 /* Must put in owned subpoints */ 3262 for (p = pStart; p < pEnd; ++p) { 3263 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3264 3265 if (subPoint < 0) { 3266 newOwners[p-pStart].rank = -3; 3267 newOwners[p-pStart].index = -3; 3268 } else { 3269 newOwners[p-pStart].rank = rank; 3270 newOwners[p-pStart].index = subPoint; 3271 } 3272 } 3273 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3274 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3275 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3276 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3277 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3278 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3279 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3280 const PetscInt point = localPoints[l]; 3281 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3282 3283 if (subPoint < 0) continue; 3284 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3285 slocalPoints[sl] = subPoint; 3286 sremotePoints[sl].rank = newLocalPoints[point].rank; 3287 sremotePoints[sl].index = newLocalPoints[point].index; 3288 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3289 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3290 ++sl; 3291 } 3292 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3293 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3294 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3295 } 3296 } 3297 CHKMEMQ; 3298 /* Cleanup */ 3299 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3300 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3301 ierr = PetscFree(subCells);CHKERRQ(ierr); 3302 PetscFunctionReturn(0); 3303 } 3304 3305 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3306 { 3307 DMLabel label = NULL; 3308 PetscErrorCode ierr; 3309 3310 PetscFunctionBegin; 3311 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3312 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);CHKERRQ(ierr); 3313 PetscFunctionReturn(0); 3314 } 3315 3316 /*@C 3317 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. 3318 3319 Input Parameters: 3320 + dm - The original mesh 3321 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3322 . label - A label name, or NULL 3323 - value - A label value 3324 3325 Output Parameter: 3326 . subdm - The surface mesh 3327 3328 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3329 3330 Level: developer 3331 3332 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3333 @*/ 3334 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3335 { 3336 PetscInt dim, cdim, depth; 3337 PetscErrorCode ierr; 3338 3339 PetscFunctionBegin; 3340 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3341 PetscValidPointer(subdm, 5); 3342 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3343 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3344 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3345 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3346 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3347 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3348 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3349 if (depth == dim) { 3350 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3351 } else { 3352 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3353 } 3354 PetscFunctionReturn(0); 3355 } 3356 3357 /*@ 3358 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3359 3360 Input Parameters: 3361 + dm - The original mesh 3362 . cellLabel - The DMLabel marking cells contained in the new mesh 3363 - value - The label value to use 3364 3365 Output Parameter: 3366 . subdm - The new mesh 3367 3368 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3369 3370 Level: developer 3371 3372 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3373 @*/ 3374 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3375 { 3376 PetscInt dim; 3377 PetscErrorCode ierr; 3378 3379 PetscFunctionBegin; 3380 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3381 PetscValidPointer(subdm, 3); 3382 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3383 ierr = DMCreate(PetscObjectComm((PetscObject) dm), subdm);CHKERRQ(ierr); 3384 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3385 ierr = DMSetDimension(*subdm, dim);CHKERRQ(ierr); 3386 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3387 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);CHKERRQ(ierr); 3388 PetscFunctionReturn(0); 3389 } 3390 3391 /*@ 3392 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3393 3394 Input Parameter: 3395 . dm - The submesh DM 3396 3397 Output Parameter: 3398 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3399 3400 Level: developer 3401 3402 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3403 @*/ 3404 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3405 { 3406 PetscFunctionBegin; 3407 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3408 PetscValidPointer(subpointMap, 2); 3409 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3410 PetscFunctionReturn(0); 3411 } 3412 3413 /*@ 3414 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3415 3416 Input Parameters: 3417 + dm - The submesh DM 3418 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3419 3420 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3421 3422 Level: developer 3423 3424 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3425 @*/ 3426 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3427 { 3428 DM_Plex *mesh = (DM_Plex *) dm->data; 3429 DMLabel tmp; 3430 PetscErrorCode ierr; 3431 3432 PetscFunctionBegin; 3433 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3434 tmp = mesh->subpointMap; 3435 mesh->subpointMap = subpointMap; 3436 ++mesh->subpointMap->refct; 3437 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3438 PetscFunctionReturn(0); 3439 } 3440 3441 /*@ 3442 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3443 3444 Input Parameter: 3445 . dm - The submesh DM 3446 3447 Output Parameter: 3448 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3449 3450 Note: This IS is guaranteed to be sorted by the construction of the submesh 3451 3452 Level: developer 3453 3454 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3455 @*/ 3456 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3457 { 3458 MPI_Comm comm; 3459 DMLabel subpointMap; 3460 IS is; 3461 const PetscInt *opoints; 3462 PetscInt *points, *depths; 3463 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3464 PetscErrorCode ierr; 3465 3466 PetscFunctionBegin; 3467 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3468 PetscValidPointer(subpointIS, 2); 3469 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3470 *subpointIS = NULL; 3471 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3472 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3473 if (subpointMap && depth >= 0) { 3474 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3475 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3476 ierr = DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3477 depths[0] = depth; 3478 depths[1] = 0; 3479 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3480 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3481 for(d = 0, off = 0; d <= depth; ++d) { 3482 const PetscInt dep = depths[d]; 3483 3484 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3485 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3486 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3487 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); 3488 } else { 3489 if (!n) { 3490 if (d == 0) { 3491 /* Missing cells */ 3492 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3493 } else { 3494 /* Missing faces */ 3495 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3496 } 3497 } 3498 } 3499 if (n) { 3500 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3501 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3502 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3503 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3504 ierr = ISDestroy(&is);CHKERRQ(ierr); 3505 } 3506 } 3507 ierr = DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3508 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3509 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3510 } 3511 PetscFunctionReturn(0); 3512 } 3513