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