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