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