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