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, 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 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 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 /* Setup chart */ 2739 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2740 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2741 for (d = 0; d <= dim; ++d) { 2742 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2743 totSubPoints += numSubPoints[d]; 2744 } 2745 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2746 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2747 /* Set cone sizes */ 2748 firstSubPoint[dim] = 0; 2749 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2750 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2751 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2752 for (d = 0; d <= dim; ++d) { 2753 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2754 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2755 } 2756 for (d = 0; d <= dim; ++d) { 2757 for (p = 0; p < numSubPoints[d]; ++p) { 2758 const PetscInt point = subpoints[d][p]; 2759 const PetscInt subpoint = firstSubPoint[d] + p; 2760 const PetscInt *cone; 2761 PetscInt coneSize, coneSizeNew, c, val; 2762 2763 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2764 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2765 if (d == dim) { 2766 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2767 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2768 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2769 if (val >= 0) coneSizeNew++; 2770 } 2771 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2772 } 2773 } 2774 } 2775 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2776 /* Set cones */ 2777 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2778 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 2779 for (d = 0; d <= dim; ++d) { 2780 for (p = 0; p < numSubPoints[d]; ++p) { 2781 const PetscInt point = subpoints[d][p]; 2782 const PetscInt subpoint = firstSubPoint[d] + p; 2783 const PetscInt *cone, *ornt; 2784 PetscInt coneSize, subconeSize, coneSizeNew, c, subc; 2785 2786 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2787 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 2788 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2789 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 2790 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2791 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 2792 if (subc >= 0) { 2793 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 2794 orntNew[coneSizeNew] = ornt[c]; 2795 ++coneSizeNew; 2796 } 2797 } 2798 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 2799 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 2800 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 2801 } 2802 } 2803 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 2804 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2805 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2806 /* Build coordinates */ 2807 { 2808 PetscSection coordSection, subCoordSection; 2809 Vec coordinates, subCoordinates; 2810 PetscScalar *coords, *subCoords; 2811 PetscInt numComp, coordSize; 2812 2813 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2814 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2815 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2816 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2817 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2818 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2819 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 2820 for (v = 0; v < numSubPoints[0]; ++v) { 2821 const PetscInt vertex = subpoints[0][v]; 2822 const PetscInt subvertex = firstSubPoint[0]+v; 2823 PetscInt dof; 2824 2825 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2826 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2827 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2828 } 2829 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2830 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2831 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 2832 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2833 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2834 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2835 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2836 for (v = 0; v < numSubPoints[0]; ++v) { 2837 const PetscInt vertex = subpoints[0][v]; 2838 const PetscInt subvertex = firstSubPoint[0]+v; 2839 PetscInt dof, off, sdof, soff, d; 2840 2841 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2842 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2843 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2844 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2845 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2846 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2847 } 2848 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2849 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2850 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2851 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2852 } 2853 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 2854 { 2855 PetscSF sfPoint, sfPointSub; 2856 IS subpIS; 2857 const PetscSFNode *remotePoints; 2858 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 2859 const PetscInt *localPoints, *subpoints; 2860 PetscInt *slocalPoints; 2861 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 2862 PetscMPIInt rank; 2863 2864 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2865 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2866 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 2867 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2868 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 2869 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 2870 if (subpIS) { 2871 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 2872 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 2873 } 2874 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2875 if (numRoots >= 0) { 2876 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 2877 for (p = 0; p < pEnd-pStart; ++p) { 2878 newLocalPoints[p].rank = -2; 2879 newLocalPoints[p].index = -2; 2880 } 2881 /* Set subleaves */ 2882 for (l = 0; l < numLeaves; ++l) { 2883 const PetscInt point = localPoints[l]; 2884 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 2885 2886 if (subpoint < 0) continue; 2887 newLocalPoints[point-pStart].rank = rank; 2888 newLocalPoints[point-pStart].index = subpoint; 2889 ++numSubleaves; 2890 } 2891 /* Must put in owned subpoints */ 2892 for (p = pStart; p < pEnd; ++p) { 2893 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 2894 2895 if (subpoint < 0) { 2896 newOwners[p-pStart].rank = -3; 2897 newOwners[p-pStart].index = -3; 2898 } else { 2899 newOwners[p-pStart].rank = rank; 2900 newOwners[p-pStart].index = subpoint; 2901 } 2902 } 2903 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 2904 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 2905 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 2906 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 2907 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 2908 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 2909 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 2910 const PetscInt point = localPoints[l]; 2911 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 2912 2913 if (subpoint < 0) continue; 2914 if (newLocalPoints[point].rank == rank) {++ll; continue;} 2915 slocalPoints[sl] = subpoint; 2916 sremotePoints[sl].rank = newLocalPoints[point].rank; 2917 sremotePoints[sl].index = newLocalPoints[point].index; 2918 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 2919 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 2920 ++sl; 2921 } 2922 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 2923 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 2924 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 2925 } 2926 if (subpIS) { 2927 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 2928 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 2929 } 2930 } 2931 /* Cleanup */ 2932 for (d = 0; d <= dim; ++d) { 2933 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2934 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 2935 } 2936 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 2937 PetscFunctionReturn(0); 2938 } 2939 2940 #undef __FUNCT__ 2941 #define __FUNCT__ "DMPlexCreateSubmesh_Interpolated" 2942 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2943 { 2944 PetscErrorCode ierr; 2945 2946 PetscFunctionBegin; 2947 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, subdm);CHKERRQ(ierr); 2948 PetscFunctionReturn(0); 2949 } 2950 2951 #undef __FUNCT__ 2952 #define __FUNCT__ "DMPlexCreateSubmesh" 2953 /*@ 2954 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 2955 2956 Input Parameters: 2957 + dm - The original mesh 2958 . vertexLabel - The DMLabel marking vertices contained in the surface 2959 - value - The label value to use 2960 2961 Output Parameter: 2962 . subdm - The surface mesh 2963 2964 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 2965 2966 Level: developer 2967 2968 .seealso: DMPlexGetSubpointMap(), DMPlexGetLabel(), DMLabelSetValue() 2969 @*/ 2970 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 2971 { 2972 PetscInt dim, depth; 2973 PetscErrorCode ierr; 2974 2975 PetscFunctionBegin; 2976 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2977 PetscValidPointer(subdm, 3); 2978 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2979 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2980 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 2981 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 2982 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 2983 if (depth == dim) { 2984 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 2985 } else { 2986 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 2987 } 2988 PetscFunctionReturn(0); 2989 } 2990 2991 #undef __FUNCT__ 2992 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Uninterpolated" 2993 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 2994 { 2995 MPI_Comm comm; 2996 DMLabel subpointMap; 2997 IS subvertexIS; 2998 const PetscInt *subVertices; 2999 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3000 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3001 PetscInt cMax, c, f; 3002 PetscErrorCode ierr; 3003 3004 PetscFunctionBegin; 3005 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 3006 /* Create subpointMap which marks the submesh */ 3007 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 3008 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3009 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3010 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 3011 /* Setup chart */ 3012 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 3013 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 3014 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 3015 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 3016 /* Set cone sizes */ 3017 firstSubVertex = numSubCells; 3018 firstSubFace = numSubCells+numSubVertices; 3019 newFacePoint = firstSubFace; 3020 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 3021 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3022 for (c = 0; c < numSubCells; ++c) { 3023 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 3024 } 3025 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3026 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 3027 } 3028 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3029 /* Create face cones */ 3030 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3031 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3032 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 3033 for (c = 0; c < numSubCells; ++c) { 3034 const PetscInt cell = subCells[c]; 3035 const PetscInt subcell = c; 3036 const PetscInt *cone, *cells; 3037 PetscInt numCells, subVertex, p, v; 3038 3039 if (cell < cMax) continue; 3040 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 3041 for (v = 0; v < nFV; ++v) { 3042 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 3043 subface[v] = firstSubVertex+subVertex; 3044 } 3045 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 3046 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 3047 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3048 /* Not true in parallel 3049 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3050 for (p = 0; p < numCells; ++p) { 3051 PetscInt negsubcell; 3052 3053 if (cells[p] >= cMax) continue; 3054 /* I know this is a crap search */ 3055 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3056 if (subCells[negsubcell] == cells[p]) break; 3057 } 3058 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3059 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 3060 } 3061 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3062 ++newFacePoint; 3063 } 3064 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 3065 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3066 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3067 /* Build coordinates */ 3068 { 3069 PetscSection coordSection, subCoordSection; 3070 Vec coordinates, subCoordinates; 3071 PetscScalar *coords, *subCoords; 3072 PetscInt numComp, coordSize, v; 3073 3074 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3075 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3076 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3077 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3078 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3079 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3080 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3081 for (v = 0; v < numSubVertices; ++v) { 3082 const PetscInt vertex = subVertices[v]; 3083 const PetscInt subvertex = firstSubVertex+v; 3084 PetscInt dof; 3085 3086 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3087 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3088 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3089 } 3090 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3091 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3092 ierr = VecCreate(comm, &subCoordinates);CHKERRQ(ierr); 3093 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3094 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3095 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3096 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3097 for (v = 0; v < numSubVertices; ++v) { 3098 const PetscInt vertex = subVertices[v]; 3099 const PetscInt subvertex = firstSubVertex+v; 3100 PetscInt dof, off, sdof, soff, d; 3101 3102 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3103 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3104 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3105 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3106 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3107 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3108 } 3109 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3110 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3111 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3112 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3113 } 3114 /* Build SF */ 3115 CHKMEMQ; 3116 { 3117 PetscSF sfPoint, sfPointSub; 3118 const PetscSFNode *remotePoints; 3119 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3120 const PetscInt *localPoints; 3121 PetscInt *slocalPoints; 3122 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3123 PetscMPIInt rank; 3124 3125 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3126 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3127 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3128 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3129 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3130 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3131 if (numRoots >= 0) { 3132 /* Only vertices should be shared */ 3133 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3134 for (p = 0; p < pEnd-pStart; ++p) { 3135 newLocalPoints[p].rank = -2; 3136 newLocalPoints[p].index = -2; 3137 } 3138 /* Set subleaves */ 3139 for (l = 0; l < numLeaves; ++l) { 3140 const PetscInt point = localPoints[l]; 3141 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3142 3143 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3144 if (subPoint < 0) continue; 3145 newLocalPoints[point-pStart].rank = rank; 3146 newLocalPoints[point-pStart].index = subPoint; 3147 ++numSubLeaves; 3148 } 3149 /* Must put in owned subpoints */ 3150 for (p = pStart; p < pEnd; ++p) { 3151 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3152 3153 if (subPoint < 0) { 3154 newOwners[p-pStart].rank = -3; 3155 newOwners[p-pStart].index = -3; 3156 } else { 3157 newOwners[p-pStart].rank = rank; 3158 newOwners[p-pStart].index = subPoint; 3159 } 3160 } 3161 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3162 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3163 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3164 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3165 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3166 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3167 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3168 const PetscInt point = localPoints[l]; 3169 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3170 3171 if (subPoint < 0) continue; 3172 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3173 slocalPoints[sl] = subPoint; 3174 sremotePoints[sl].rank = newLocalPoints[point].rank; 3175 sremotePoints[sl].index = newLocalPoints[point].index; 3176 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3177 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3178 ++sl; 3179 } 3180 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3181 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3182 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3183 } 3184 } 3185 CHKMEMQ; 3186 /* Cleanup */ 3187 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3188 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3189 ierr = PetscFree(subCells);CHKERRQ(ierr); 3190 PetscFunctionReturn(0); 3191 } 3192 3193 #undef __FUNCT__ 3194 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh_Interpolated" 3195 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3196 { 3197 DMLabel label = NULL; 3198 PetscErrorCode ierr; 3199 3200 PetscFunctionBegin; 3201 if (labelname) {ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3202 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, subdm);CHKERRQ(ierr); 3203 PetscFunctionReturn(0); 3204 } 3205 3206 #undef __FUNCT__ 3207 #define __FUNCT__ "DMPlexCreateCohesiveSubmesh" 3208 /* 3209 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. 3210 3211 Input Parameters: 3212 + dm - The original mesh 3213 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3214 . label - A label name, or NULL 3215 - value - A label value 3216 3217 Output Parameter: 3218 . subdm - The surface mesh 3219 3220 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3221 3222 Level: developer 3223 3224 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3225 */ 3226 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3227 { 3228 PetscInt dim, depth; 3229 PetscErrorCode ierr; 3230 3231 PetscFunctionBegin; 3232 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3233 PetscValidPointer(subdm, 5); 3234 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3235 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3236 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3237 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3238 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3239 if (depth == dim) { 3240 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3241 } else { 3242 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3243 } 3244 PetscFunctionReturn(0); 3245 } 3246 3247 #undef __FUNCT__ 3248 #define __FUNCT__ "DMPlexGetSubpointMap" 3249 /*@ 3250 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3251 3252 Input Parameter: 3253 . dm - The submesh DM 3254 3255 Output Parameter: 3256 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3257 3258 Level: developer 3259 3260 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3261 @*/ 3262 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3263 { 3264 DM_Plex *mesh = (DM_Plex*) dm->data; 3265 3266 PetscFunctionBegin; 3267 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3268 PetscValidPointer(subpointMap, 2); 3269 *subpointMap = mesh->subpointMap; 3270 PetscFunctionReturn(0); 3271 } 3272 3273 #undef __FUNCT__ 3274 #define __FUNCT__ "DMPlexSetSubpointMap" 3275 /* Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() */ 3276 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3277 { 3278 DM_Plex *mesh = (DM_Plex *) dm->data; 3279 DMLabel tmp; 3280 PetscErrorCode ierr; 3281 3282 PetscFunctionBegin; 3283 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3284 tmp = mesh->subpointMap; 3285 mesh->subpointMap = subpointMap; 3286 ++mesh->subpointMap->refct; 3287 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3288 PetscFunctionReturn(0); 3289 } 3290 3291 #undef __FUNCT__ 3292 #define __FUNCT__ "DMPlexCreateSubpointIS" 3293 /*@ 3294 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3295 3296 Input Parameter: 3297 . dm - The submesh DM 3298 3299 Output Parameter: 3300 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3301 3302 Note: This IS is guaranteed to be sorted by the construction of the submesh 3303 3304 Level: developer 3305 3306 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3307 @*/ 3308 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3309 { 3310 MPI_Comm comm; 3311 DMLabel subpointMap; 3312 IS is; 3313 const PetscInt *opoints; 3314 PetscInt *points, *depths; 3315 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3316 PetscErrorCode ierr; 3317 3318 PetscFunctionBegin; 3319 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3320 PetscValidPointer(subpointIS, 2); 3321 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3322 *subpointIS = NULL; 3323 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3324 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3325 if (subpointMap && depth >= 0) { 3326 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3327 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3328 ierr = DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3329 depths[0] = depth; 3330 depths[1] = 0; 3331 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3332 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3333 for(d = 0, off = 0; d <= depth; ++d) { 3334 const PetscInt dep = depths[d]; 3335 3336 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3337 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3338 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3339 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); 3340 } else { 3341 if (!n) { 3342 if (d == 0) { 3343 /* Missing cells */ 3344 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3345 } else { 3346 /* Missing faces */ 3347 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3348 } 3349 } 3350 } 3351 if (n) { 3352 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3353 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3354 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3355 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3356 ierr = ISDestroy(&is);CHKERRQ(ierr); 3357 } 3358 } 3359 ierr = DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3360 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3361 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3362 } 3363 PetscFunctionReturn(0); 3364 } 3365