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