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