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