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