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