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