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