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