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 /* Check that no cell have all vertices on the fault */ 1876 PetscErrorCode DMPlexCheckValidSubmesh_Private(DM dm, DMLabel label, DM subdm) 1877 { 1878 IS subpointIS; 1879 const PetscInt *dmpoints; 1880 PetscInt defaultValue, cStart, cEnd, c, vStart, vEnd; 1881 PetscErrorCode ierr; 1882 1883 PetscFunctionBegin; 1884 if (!label) PetscFunctionReturn(0); 1885 ierr = DMLabelGetDefaultValue(label, &defaultValue);CHKERRQ(ierr); 1886 ierr = DMPlexCreateSubpointIS(subdm, &subpointIS);CHKERRQ(ierr); 1887 if (!subpointIS) PetscFunctionReturn(0); 1888 ierr = DMPlexGetHeightStratum(subdm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1889 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1890 ierr = ISGetIndices(subpointIS, &dmpoints);CHKERRQ(ierr); 1891 for (c = cStart; c < cEnd; ++c) { 1892 PetscBool invalidCell = PETSC_TRUE; 1893 PetscInt *closure = NULL; 1894 PetscInt closureSize, cl; 1895 1896 ierr = DMPlexGetTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1897 for (cl = 0; cl < closureSize*2; cl += 2) { 1898 PetscInt value = 0; 1899 1900 if ((closure[cl] < vStart) || (closure[cl] >= vEnd)) continue; 1901 ierr = DMLabelGetValue(label, closure[cl], &value);CHKERRQ(ierr); 1902 if (value == defaultValue) {invalidCell = PETSC_FALSE; break;} 1903 } 1904 ierr = DMPlexRestoreTransitiveClosure(dm, dmpoints[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1905 if (invalidCell) { 1906 ierr = ISRestoreIndices(subpointIS, &dmpoints);CHKERRQ(ierr); 1907 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1908 ierr = DMDestroy(&subdm);CHKERRQ(ierr); 1909 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ambiguous submesh. Cell %D has all of its vertices on the submesh.", dmpoints[c]); 1910 } 1911 } 1912 ierr = ISRestoreIndices(subpointIS, &dmpoints);CHKERRQ(ierr); 1913 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1914 PetscFunctionReturn(0); 1915 } 1916 1917 1918 /*@ 1919 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 1920 1921 Collective on dm 1922 1923 Input Parameters: 1924 + dm - The original DM 1925 . label - The label specifying the interface vertices 1926 - bdlabel - The optional label specifying the interface boundary vertices 1927 1928 Output Parameters: 1929 + hybridLabel - The label fully marking the interface 1930 . dmInterface - The new interface DM, or NULL 1931 - dmHybrid - The new DM with cohesive cells 1932 1933 Level: developer 1934 1935 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate() 1936 @*/ 1937 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel bdlabel, DMLabel *hybridLabel, DM *dmInterface, DM *dmHybrid) 1938 { 1939 DM idm; 1940 DMLabel subpointMap, hlabel; 1941 PetscInt dim; 1942 PetscErrorCode ierr; 1943 1944 PetscFunctionBegin; 1945 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1946 if (bdlabel) PetscValidPointer(bdlabel, 3); 1947 if (hybridLabel) PetscValidPointer(hybridLabel, 4); 1948 if (dmInterface) PetscValidPointer(dmInterface, 5); 1949 PetscValidPointer(dmHybrid, 6); 1950 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1951 ierr = DMPlexCreateSubmesh(dm, label, 1, &idm);CHKERRQ(ierr); 1952 ierr = DMPlexCheckValidSubmesh_Private(dm, label, idm);CHKERRQ(ierr); 1953 ierr = DMPlexOrient(idm);CHKERRQ(ierr); 1954 ierr = DMPlexGetSubpointMap(idm, &subpointMap);CHKERRQ(ierr); 1955 ierr = DMLabelDuplicate(subpointMap, &hlabel);CHKERRQ(ierr); 1956 ierr = DMLabelClearStratum(hlabel, dim);CHKERRQ(ierr); 1957 ierr = DMPlexLabelCohesiveComplete(dm, hlabel, bdlabel, PETSC_FALSE, idm);CHKERRQ(ierr); 1958 if (dmInterface) {*dmInterface = idm;} 1959 else {ierr = DMDestroy(&idm);CHKERRQ(ierr);} 1960 ierr = DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);CHKERRQ(ierr); 1961 if (hybridLabel) *hybridLabel = hlabel; 1962 else {ierr = DMLabelDestroy(&hlabel);CHKERRQ(ierr);} 1963 PetscFunctionReturn(0); 1964 } 1965 1966 /* Here we need the explicit assumption that: 1967 1968 For any marked cell, the marked vertices constitute a single face 1969 */ 1970 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 1971 { 1972 IS subvertexIS = NULL; 1973 const PetscInt *subvertices; 1974 PetscInt *pStart, *pEnd, *pMax, pSize; 1975 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 1976 PetscErrorCode ierr; 1977 1978 PetscFunctionBegin; 1979 *numFaces = 0; 1980 *nFV = 0; 1981 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1982 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1983 pSize = PetscMax(depth, dim) + 1; 1984 ierr = PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);CHKERRQ(ierr); 1985 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1986 for (d = 0; d <= depth; ++d) { 1987 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1988 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1989 } 1990 /* Loop over initial vertices and mark all faces in the collective star() */ 1991 if (vertexLabel) {ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr);} 1992 if (subvertexIS) { 1993 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1994 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1995 } 1996 for (v = 0; v < numSubVerticesInitial; ++v) { 1997 const PetscInt vertex = subvertices[v]; 1998 PetscInt *star = NULL; 1999 PetscInt starSize, s, numCells = 0, c; 2000 2001 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2002 for (s = 0; s < starSize*2; s += 2) { 2003 const PetscInt point = star[s]; 2004 if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point; 2005 } 2006 for (c = 0; c < numCells; ++c) { 2007 const PetscInt cell = star[c]; 2008 PetscInt *closure = NULL; 2009 PetscInt closureSize, cl; 2010 PetscInt cellLoc, numCorners = 0, faceSize = 0; 2011 2012 ierr = DMLabelGetValue(subpointMap, cell, &cellLoc);CHKERRQ(ierr); 2013 if (cellLoc == 2) continue; 2014 if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 2015 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2016 for (cl = 0; cl < closureSize*2; cl += 2) { 2017 const PetscInt point = closure[cl]; 2018 PetscInt vertexLoc; 2019 2020 if ((point >= pStart[0]) && (point < pEnd[0])) { 2021 ++numCorners; 2022 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2023 if (vertexLoc == value) closure[faceSize++] = point; 2024 } 2025 } 2026 if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);} 2027 if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2028 if (faceSize == *nFV) { 2029 const PetscInt *cells = NULL; 2030 PetscInt numCells, nc; 2031 2032 ++(*numFaces); 2033 for (cl = 0; cl < faceSize; ++cl) { 2034 ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr); 2035 } 2036 ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 2037 for (nc = 0; nc < numCells; ++nc) { 2038 ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr); 2039 } 2040 ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 2041 } 2042 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2043 } 2044 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2045 } 2046 if (subvertexIS) { 2047 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2048 } 2049 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2050 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 2051 PetscFunctionReturn(0); 2052 } 2053 2054 static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm) 2055 { 2056 IS subvertexIS = NULL; 2057 const PetscInt *subvertices; 2058 PetscInt *pStart, *pEnd, *pMax; 2059 PetscInt dim, d, numSubVerticesInitial = 0, v; 2060 PetscErrorCode ierr; 2061 2062 PetscFunctionBegin; 2063 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2064 ierr = PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);CHKERRQ(ierr); 2065 ierr = DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 2066 for (d = 0; d <= dim; ++d) { 2067 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 2068 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 2069 } 2070 /* Loop over initial vertices and mark all faces in the collective star() */ 2071 if (vertexLabel) { 2072 ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr); 2073 if (subvertexIS) { 2074 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 2075 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2076 } 2077 } 2078 for (v = 0; v < numSubVerticesInitial; ++v) { 2079 const PetscInt vertex = subvertices[v]; 2080 PetscInt *star = NULL; 2081 PetscInt starSize, s, numFaces = 0, f; 2082 2083 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2084 for (s = 0; s < starSize*2; s += 2) { 2085 const PetscInt point = star[s]; 2086 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point; 2087 } 2088 for (f = 0; f < numFaces; ++f) { 2089 const PetscInt face = star[f]; 2090 PetscInt *closure = NULL; 2091 PetscInt closureSize, c; 2092 PetscInt faceLoc; 2093 2094 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 2095 if (faceLoc == dim-1) continue; 2096 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 2097 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2098 for (c = 0; c < closureSize*2; c += 2) { 2099 const PetscInt point = closure[c]; 2100 PetscInt vertexLoc; 2101 2102 if ((point >= pStart[0]) && (point < pEnd[0])) { 2103 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2104 if (vertexLoc != value) break; 2105 } 2106 } 2107 if (c == closureSize*2) { 2108 const PetscInt *support; 2109 PetscInt supportSize, s; 2110 2111 for (c = 0; c < closureSize*2; c += 2) { 2112 const PetscInt point = closure[c]; 2113 2114 for (d = 0; d < dim; ++d) { 2115 if ((point >= pStart[d]) && (point < pEnd[d])) { 2116 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2117 break; 2118 } 2119 } 2120 } 2121 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2122 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2123 for (s = 0; s < supportSize; ++s) { 2124 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2125 } 2126 } 2127 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2128 } 2129 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2130 } 2131 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr);} 2132 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2133 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 2134 PetscFunctionReturn(0); 2135 } 2136 2137 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 2138 { 2139 DMLabel label = NULL; 2140 const PetscInt *cone; 2141 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1; 2142 PetscErrorCode ierr; 2143 2144 PetscFunctionBegin; 2145 *numFaces = 0; 2146 *nFV = 0; 2147 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 2148 *subCells = NULL; 2149 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2150 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2151 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2152 if (cMax < 0) PetscFunctionReturn(0); 2153 if (label) { 2154 for (c = cMax; c < cEnd; ++c) { 2155 PetscInt val; 2156 2157 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2158 if (val == value) { 2159 ++(*numFaces); 2160 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2161 } 2162 } 2163 } else { 2164 *numFaces = cEnd - cMax; 2165 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 2166 } 2167 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 2168 if (!(*numFaces)) PetscFunctionReturn(0); 2169 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 2170 for (c = cMax; c < cEnd; ++c) { 2171 const PetscInt *cells; 2172 PetscInt numCells; 2173 2174 if (label) { 2175 PetscInt val; 2176 2177 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2178 if (val != value) continue; 2179 } 2180 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2181 for (p = 0; p < *nFV; ++p) { 2182 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 2183 } 2184 /* Negative face */ 2185 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2186 /* Not true in parallel 2187 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2188 for (p = 0; p < numCells; ++p) { 2189 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 2190 (*subCells)[subc++] = cells[p]; 2191 } 2192 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2193 /* Positive face is not included */ 2194 } 2195 PetscFunctionReturn(0); 2196 } 2197 2198 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 2199 { 2200 PetscInt *pStart, *pEnd; 2201 PetscInt dim, cMax, cEnd, c, d; 2202 PetscErrorCode ierr; 2203 2204 PetscFunctionBegin; 2205 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2206 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2207 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2208 if (cMax < 0) PetscFunctionReturn(0); 2209 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 2210 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 2211 for (c = cMax; c < cEnd; ++c) { 2212 const PetscInt *cone; 2213 PetscInt *closure = NULL; 2214 PetscInt fconeSize, coneSize, closureSize, cl, val; 2215 2216 if (label) { 2217 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2218 if (val != value) continue; 2219 } 2220 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2221 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2222 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 2223 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2224 /* Negative face */ 2225 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2226 for (cl = 0; cl < closureSize*2; cl += 2) { 2227 const PetscInt point = closure[cl]; 2228 2229 for (d = 0; d <= dim; ++d) { 2230 if ((point >= pStart[d]) && (point < pEnd[d])) { 2231 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2232 break; 2233 } 2234 } 2235 } 2236 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2237 /* Cells -- positive face is not included */ 2238 for (cl = 0; cl < 1; ++cl) { 2239 const PetscInt *support; 2240 PetscInt supportSize, s; 2241 2242 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 2243 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2244 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 2245 for (s = 0; s < supportSize; ++s) { 2246 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2247 } 2248 } 2249 } 2250 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2251 PetscFunctionReturn(0); 2252 } 2253 2254 static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2255 { 2256 MPI_Comm comm; 2257 PetscBool posOrient = PETSC_FALSE; 2258 const PetscInt debug = 0; 2259 PetscInt cellDim, faceSize, f; 2260 PetscErrorCode ierr; 2261 2262 PetscFunctionBegin; 2263 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2264 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 2265 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 2266 2267 if (cellDim == 1 && numCorners == 2) { 2268 /* Triangle */ 2269 faceSize = numCorners-1; 2270 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2271 } else if (cellDim == 2 && numCorners == 3) { 2272 /* Triangle */ 2273 faceSize = numCorners-1; 2274 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2275 } else if (cellDim == 3 && numCorners == 4) { 2276 /* Tetrahedron */ 2277 faceSize = numCorners-1; 2278 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2279 } else if (cellDim == 1 && numCorners == 3) { 2280 /* Quadratic line */ 2281 faceSize = 1; 2282 posOrient = PETSC_TRUE; 2283 } else if (cellDim == 2 && numCorners == 4) { 2284 /* Quads */ 2285 faceSize = 2; 2286 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2287 posOrient = PETSC_TRUE; 2288 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2289 posOrient = PETSC_TRUE; 2290 } else { 2291 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2292 posOrient = PETSC_FALSE; 2293 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2294 } 2295 } else if (cellDim == 2 && numCorners == 6) { 2296 /* Quadratic triangle (I hate this) */ 2297 /* Edges are determined by the first 2 vertices (corners of edges) */ 2298 const PetscInt faceSizeTri = 3; 2299 PetscInt sortedIndices[3], i, iFace; 2300 PetscBool found = PETSC_FALSE; 2301 PetscInt faceVerticesTriSorted[9] = { 2302 0, 3, 4, /* bottom */ 2303 1, 4, 5, /* right */ 2304 2, 3, 5, /* left */ 2305 }; 2306 PetscInt faceVerticesTri[9] = { 2307 0, 3, 4, /* bottom */ 2308 1, 4, 5, /* right */ 2309 2, 5, 3, /* left */ 2310 }; 2311 2312 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2313 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2314 for (iFace = 0; iFace < 3; ++iFace) { 2315 const PetscInt ii = iFace*faceSizeTri; 2316 PetscInt fVertex, cVertex; 2317 2318 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2319 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2320 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2321 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2322 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2323 faceVertices[fVertex] = origVertices[cVertex]; 2324 break; 2325 } 2326 } 2327 } 2328 found = PETSC_TRUE; 2329 break; 2330 } 2331 } 2332 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2333 if (posOriented) *posOriented = PETSC_TRUE; 2334 PetscFunctionReturn(0); 2335 } else if (cellDim == 2 && numCorners == 9) { 2336 /* Quadratic quad (I hate this) */ 2337 /* Edges are determined by the first 2 vertices (corners of edges) */ 2338 const PetscInt faceSizeQuad = 3; 2339 PetscInt sortedIndices[3], i, iFace; 2340 PetscBool found = PETSC_FALSE; 2341 PetscInt faceVerticesQuadSorted[12] = { 2342 0, 1, 4, /* bottom */ 2343 1, 2, 5, /* right */ 2344 2, 3, 6, /* top */ 2345 0, 3, 7, /* left */ 2346 }; 2347 PetscInt faceVerticesQuad[12] = { 2348 0, 1, 4, /* bottom */ 2349 1, 2, 5, /* right */ 2350 2, 3, 6, /* top */ 2351 3, 0, 7, /* left */ 2352 }; 2353 2354 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2355 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2356 for (iFace = 0; iFace < 4; ++iFace) { 2357 const PetscInt ii = iFace*faceSizeQuad; 2358 PetscInt fVertex, cVertex; 2359 2360 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2361 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2362 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2363 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2364 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2365 faceVertices[fVertex] = origVertices[cVertex]; 2366 break; 2367 } 2368 } 2369 } 2370 found = PETSC_TRUE; 2371 break; 2372 } 2373 } 2374 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2375 if (posOriented) *posOriented = PETSC_TRUE; 2376 PetscFunctionReturn(0); 2377 } else if (cellDim == 3 && numCorners == 8) { 2378 /* Hexes 2379 A hex is two oriented quads with the normal of the first 2380 pointing up at the second. 2381 2382 7---6 2383 /| /| 2384 4---5 | 2385 | 1-|-2 2386 |/ |/ 2387 0---3 2388 2389 Faces are determined by the first 4 vertices (corners of faces) */ 2390 const PetscInt faceSizeHex = 4; 2391 PetscInt sortedIndices[4], i, iFace; 2392 PetscBool found = PETSC_FALSE; 2393 PetscInt faceVerticesHexSorted[24] = { 2394 0, 1, 2, 3, /* bottom */ 2395 4, 5, 6, 7, /* top */ 2396 0, 3, 4, 5, /* front */ 2397 2, 3, 5, 6, /* right */ 2398 1, 2, 6, 7, /* back */ 2399 0, 1, 4, 7, /* left */ 2400 }; 2401 PetscInt faceVerticesHex[24] = { 2402 1, 2, 3, 0, /* bottom */ 2403 4, 5, 6, 7, /* top */ 2404 0, 3, 5, 4, /* front */ 2405 3, 2, 6, 5, /* right */ 2406 2, 1, 7, 6, /* back */ 2407 1, 0, 4, 7, /* left */ 2408 }; 2409 2410 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2411 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2412 for (iFace = 0; iFace < 6; ++iFace) { 2413 const PetscInt ii = iFace*faceSizeHex; 2414 PetscInt fVertex, cVertex; 2415 2416 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2417 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2418 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2419 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2420 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2421 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2422 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2423 faceVertices[fVertex] = origVertices[cVertex]; 2424 break; 2425 } 2426 } 2427 } 2428 found = PETSC_TRUE; 2429 break; 2430 } 2431 } 2432 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2433 if (posOriented) *posOriented = PETSC_TRUE; 2434 PetscFunctionReturn(0); 2435 } else if (cellDim == 3 && numCorners == 10) { 2436 /* Quadratic tet */ 2437 /* Faces are determined by the first 3 vertices (corners of faces) */ 2438 const PetscInt faceSizeTet = 6; 2439 PetscInt sortedIndices[6], i, iFace; 2440 PetscBool found = PETSC_FALSE; 2441 PetscInt faceVerticesTetSorted[24] = { 2442 0, 1, 2, 6, 7, 8, /* bottom */ 2443 0, 3, 4, 6, 7, 9, /* front */ 2444 1, 4, 5, 7, 8, 9, /* right */ 2445 2, 3, 5, 6, 8, 9, /* left */ 2446 }; 2447 PetscInt faceVerticesTet[24] = { 2448 0, 1, 2, 6, 7, 8, /* bottom */ 2449 0, 4, 3, 6, 7, 9, /* front */ 2450 1, 5, 4, 7, 8, 9, /* right */ 2451 2, 3, 5, 8, 6, 9, /* left */ 2452 }; 2453 2454 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2455 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2456 for (iFace=0; iFace < 4; ++iFace) { 2457 const PetscInt ii = iFace*faceSizeTet; 2458 PetscInt fVertex, cVertex; 2459 2460 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2461 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2462 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2463 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2464 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2465 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2466 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2467 faceVertices[fVertex] = origVertices[cVertex]; 2468 break; 2469 } 2470 } 2471 } 2472 found = PETSC_TRUE; 2473 break; 2474 } 2475 } 2476 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2477 if (posOriented) *posOriented = PETSC_TRUE; 2478 PetscFunctionReturn(0); 2479 } else if (cellDim == 3 && numCorners == 27) { 2480 /* Quadratic hexes (I hate this) 2481 A hex is two oriented quads with the normal of the first 2482 pointing up at the second. 2483 2484 7---6 2485 /| /| 2486 4---5 | 2487 | 3-|-2 2488 |/ |/ 2489 0---1 2490 2491 Faces are determined by the first 4 vertices (corners of faces) */ 2492 const PetscInt faceSizeQuadHex = 9; 2493 PetscInt sortedIndices[9], i, iFace; 2494 PetscBool found = PETSC_FALSE; 2495 PetscInt faceVerticesQuadHexSorted[54] = { 2496 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2497 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2498 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2499 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2500 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2501 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2502 }; 2503 PetscInt faceVerticesQuadHex[54] = { 2504 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2505 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2506 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2507 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2508 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2509 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2510 }; 2511 2512 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2513 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2514 for (iFace = 0; iFace < 6; ++iFace) { 2515 const PetscInt ii = iFace*faceSizeQuadHex; 2516 PetscInt fVertex, cVertex; 2517 2518 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2519 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2520 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2521 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2522 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2523 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2524 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2525 faceVertices[fVertex] = origVertices[cVertex]; 2526 break; 2527 } 2528 } 2529 } 2530 found = PETSC_TRUE; 2531 break; 2532 } 2533 } 2534 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2535 if (posOriented) *posOriented = PETSC_TRUE; 2536 PetscFunctionReturn(0); 2537 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2538 if (!posOrient) { 2539 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2540 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2541 } else { 2542 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2543 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2544 } 2545 if (posOriented) *posOriented = posOrient; 2546 PetscFunctionReturn(0); 2547 } 2548 2549 /*@ 2550 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 2551 in faceVertices. The orientation is such that the face normal points out of the cell 2552 2553 Not collective 2554 2555 Input Parameters: 2556 + dm - The original mesh 2557 . cell - The cell mesh point 2558 . faceSize - The number of vertices on the face 2559 . face - The face vertices 2560 . numCorners - The number of vertices on the cell 2561 . indices - Local numbering of face vertices in cell cone 2562 - origVertices - Original face vertices 2563 2564 Output Parameter: 2565 + faceVertices - The face vertices properly oriented 2566 - posOriented - PETSC_TRUE if the face was oriented with outward normal 2567 2568 Level: developer 2569 2570 .seealso: DMPlexGetCone() 2571 @*/ 2572 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2573 { 2574 const PetscInt *cone = NULL; 2575 PetscInt coneSize, v, f, v2; 2576 PetscInt oppositeVertex = -1; 2577 PetscErrorCode ierr; 2578 2579 PetscFunctionBegin; 2580 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2581 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2582 for (v = 0, v2 = 0; v < coneSize; ++v) { 2583 PetscBool found = PETSC_FALSE; 2584 2585 for (f = 0; f < faceSize; ++f) { 2586 if (face[f] == cone[v]) { 2587 found = PETSC_TRUE; break; 2588 } 2589 } 2590 if (found) { 2591 indices[v2] = v; 2592 origVertices[v2] = cone[v]; 2593 ++v2; 2594 } else { 2595 oppositeVertex = v; 2596 } 2597 } 2598 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2599 PetscFunctionReturn(0); 2600 } 2601 2602 /* 2603 DMPlexInsertFace_Internal - Puts a face into the mesh 2604 2605 Not collective 2606 2607 Input Parameters: 2608 + dm - The DMPlex 2609 . numFaceVertex - The number of vertices in the face 2610 . faceVertices - The vertices in the face for dm 2611 . subfaceVertices - The vertices in the face for subdm 2612 . numCorners - The number of vertices in the cell 2613 . cell - A cell in dm containing the face 2614 . subcell - A cell in subdm containing the face 2615 . firstFace - First face in the mesh 2616 - newFacePoint - Next face in the mesh 2617 2618 Output Parameters: 2619 . newFacePoint - Contains next face point number on input, updated on output 2620 2621 Level: developer 2622 */ 2623 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) 2624 { 2625 MPI_Comm comm; 2626 DM_Plex *submesh = (DM_Plex*) subdm->data; 2627 const PetscInt *faces; 2628 PetscInt numFaces, coneSize; 2629 PetscErrorCode ierr; 2630 2631 PetscFunctionBegin; 2632 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2633 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2634 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2635 #if 0 2636 /* Cannot use this because support() has not been constructed yet */ 2637 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2638 #else 2639 { 2640 PetscInt f; 2641 2642 numFaces = 0; 2643 ierr = DMGetWorkArray(subdm, 1, MPIU_INT, (void **) &faces);CHKERRQ(ierr); 2644 for (f = firstFace; f < *newFacePoint; ++f) { 2645 PetscInt dof, off, d; 2646 2647 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2648 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2649 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2650 for (d = 0; d < dof; ++d) { 2651 const PetscInt p = submesh->cones[off+d]; 2652 PetscInt v; 2653 2654 for (v = 0; v < numFaceVertices; ++v) { 2655 if (subfaceVertices[v] == p) break; 2656 } 2657 if (v == numFaceVertices) break; 2658 } 2659 if (d == dof) { 2660 numFaces = 1; 2661 ((PetscInt*) faces)[0] = f; 2662 } 2663 } 2664 } 2665 #endif 2666 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2667 else if (numFaces == 1) { 2668 /* Add the other cell neighbor for this face */ 2669 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2670 } else { 2671 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2672 PetscBool posOriented; 2673 2674 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);CHKERRQ(ierr); 2675 origVertices = &orientedVertices[numFaceVertices]; 2676 indices = &orientedVertices[numFaceVertices*2]; 2677 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2678 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2679 /* TODO: I know that routine should return a permutation, not the indices */ 2680 for (v = 0; v < numFaceVertices; ++v) { 2681 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2682 for (ov = 0; ov < numFaceVertices; ++ov) { 2683 if (orientedVertices[ov] == vertex) { 2684 orientedSubVertices[ov] = subvertex; 2685 break; 2686 } 2687 } 2688 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2689 } 2690 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2691 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2692 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), MPIU_INT, &orientedVertices);CHKERRQ(ierr); 2693 ++(*newFacePoint); 2694 } 2695 #if 0 2696 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2697 #else 2698 ierr = DMRestoreWorkArray(subdm, 1, MPIU_INT, (void **) &faces);CHKERRQ(ierr); 2699 #endif 2700 PetscFunctionReturn(0); 2701 } 2702 2703 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2704 { 2705 MPI_Comm comm; 2706 DMLabel subpointMap; 2707 IS subvertexIS, subcellIS; 2708 const PetscInt *subVertices, *subCells; 2709 PetscInt numSubVertices, firstSubVertex, numSubCells; 2710 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2711 PetscInt vStart, vEnd, c, f; 2712 PetscErrorCode ierr; 2713 2714 PetscFunctionBegin; 2715 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2716 /* Create subpointMap which marks the submesh */ 2717 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2718 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2719 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2720 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2721 /* Setup chart */ 2722 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2723 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2724 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2725 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2726 /* Set cone sizes */ 2727 firstSubVertex = numSubCells; 2728 firstSubFace = numSubCells+numSubVertices; 2729 newFacePoint = firstSubFace; 2730 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2731 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2732 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2733 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2734 for (c = 0; c < numSubCells; ++c) { 2735 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2736 } 2737 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2738 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2739 } 2740 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2741 /* Create face cones */ 2742 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2743 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2744 ierr = DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 2745 for (c = 0; c < numSubCells; ++c) { 2746 const PetscInt cell = subCells[c]; 2747 const PetscInt subcell = c; 2748 PetscInt *closure = NULL; 2749 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2750 2751 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2752 for (cl = 0; cl < closureSize*2; cl += 2) { 2753 const PetscInt point = closure[cl]; 2754 PetscInt subVertex; 2755 2756 if ((point >= vStart) && (point < vEnd)) { 2757 ++numCorners; 2758 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2759 if (subVertex >= 0) { 2760 closure[faceSize] = point; 2761 subface[faceSize] = firstSubVertex+subVertex; 2762 ++faceSize; 2763 } 2764 } 2765 } 2766 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2767 if (faceSize == nFV) { 2768 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2769 } 2770 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2771 } 2772 ierr = DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 2773 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2774 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2775 /* Build coordinates */ 2776 { 2777 PetscSection coordSection, subCoordSection; 2778 Vec coordinates, subCoordinates; 2779 PetscScalar *coords, *subCoords; 2780 PetscInt numComp, coordSize, v; 2781 const char *name; 2782 2783 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2784 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2785 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2786 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2787 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2788 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2789 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2790 for (v = 0; v < numSubVertices; ++v) { 2791 const PetscInt vertex = subVertices[v]; 2792 const PetscInt subvertex = firstSubVertex+v; 2793 PetscInt dof; 2794 2795 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2796 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2797 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2798 } 2799 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2800 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2801 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 2802 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 2803 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 2804 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2805 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2806 if (coordSize) { 2807 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2808 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2809 for (v = 0; v < numSubVertices; ++v) { 2810 const PetscInt vertex = subVertices[v]; 2811 const PetscInt subvertex = firstSubVertex+v; 2812 PetscInt dof, off, sdof, soff, d; 2813 2814 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2815 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2816 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2817 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2818 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2819 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2820 } 2821 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2822 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2823 } 2824 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2825 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2826 } 2827 /* Cleanup */ 2828 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2829 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2830 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2831 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 2832 PetscFunctionReturn(0); 2833 } 2834 2835 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 2836 { 2837 PetscInt subPoint; 2838 PetscErrorCode ierr; 2839 2840 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 2841 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 2842 } 2843 2844 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 2845 { 2846 MPI_Comm comm; 2847 DMLabel subpointMap; 2848 IS *subpointIS; 2849 const PetscInt **subpoints; 2850 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 2851 PetscInt totSubPoints = 0, maxConeSize, cMax, cEnd, dim, p, d, v; 2852 PetscMPIInt rank; 2853 PetscErrorCode ierr; 2854 2855 PetscFunctionBegin; 2856 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2857 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2858 /* Create subpointMap which marks the submesh */ 2859 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2860 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2861 if (cellHeight) { 2862 if (isCohesive) {ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2863 else {ierr = DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2864 } else { 2865 DMLabel depth; 2866 IS pointIS; 2867 const PetscInt *points; 2868 PetscInt numPoints; 2869 2870 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 2871 ierr = DMLabelGetStratumSize(label, value, &numPoints);CHKERRQ(ierr); 2872 ierr = DMLabelGetStratumIS(label, value, &pointIS);CHKERRQ(ierr); 2873 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2874 for (p = 0; p < numPoints; ++p) { 2875 PetscInt *closure = NULL; 2876 PetscInt closureSize, c, pdim; 2877 2878 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2879 for (c = 0; c < closureSize*2; c += 2) { 2880 ierr = DMLabelGetValue(depth, closure[c], &pdim);CHKERRQ(ierr); 2881 ierr = DMLabelSetValue(subpointMap, closure[c], pdim);CHKERRQ(ierr); 2882 } 2883 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2884 } 2885 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2886 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2887 } 2888 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2889 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2890 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2891 cMax = (cMax < 0) ? cEnd : cMax; 2892 /* Setup chart */ 2893 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2894 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2895 for (d = 0; d <= dim; ++d) { 2896 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2897 totSubPoints += numSubPoints[d]; 2898 } 2899 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2900 ierr = DMPlexSetVTKCellHeight(subdm, cellHeight);CHKERRQ(ierr); 2901 /* Set cone sizes */ 2902 firstSubPoint[dim] = 0; 2903 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2904 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2905 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2906 for (d = 0; d <= dim; ++d) { 2907 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2908 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2909 } 2910 for (d = 0; d <= dim; ++d) { 2911 for (p = 0; p < numSubPoints[d]; ++p) { 2912 const PetscInt point = subpoints[d][p]; 2913 const PetscInt subpoint = firstSubPoint[d] + p; 2914 const PetscInt *cone; 2915 PetscInt coneSize, coneSizeNew, c, val; 2916 2917 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2918 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2919 if (cellHeight && (d == dim)) { 2920 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2921 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2922 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2923 if (val >= 0) coneSizeNew++; 2924 } 2925 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2926 } 2927 } 2928 } 2929 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2930 /* Set cones */ 2931 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2932 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 2933 for (d = 0; d <= dim; ++d) { 2934 for (p = 0; p < numSubPoints[d]; ++p) { 2935 const PetscInt point = subpoints[d][p]; 2936 const PetscInt subpoint = firstSubPoint[d] + p; 2937 const PetscInt *cone, *ornt; 2938 PetscInt coneSize, subconeSize, coneSizeNew, c, subc, fornt = 0; 2939 2940 if (d == dim-1) { 2941 const PetscInt *support, *cone, *ornt; 2942 PetscInt supportSize, coneSize, s, subc; 2943 2944 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 2945 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 2946 for (s = 0; s < supportSize; ++s) { 2947 if ((support[s] < cMax) || (support[s] >= cEnd)) continue; 2948 ierr = PetscFindInt(support[s], numSubPoints[d+1], subpoints[d+1], &subc);CHKERRQ(ierr); 2949 if (subc >= 0) { 2950 const PetscInt ccell = subpoints[d+1][subc]; 2951 2952 ierr = DMPlexGetCone(dm, ccell, &cone);CHKERRQ(ierr); 2953 ierr = DMPlexGetConeSize(dm, ccell, &coneSize);CHKERRQ(ierr); 2954 ierr = DMPlexGetConeOrientation(dm, ccell, &ornt);CHKERRQ(ierr); 2955 for (c = 0; c < coneSize; ++c) { 2956 if (cone[c] == point) { 2957 fornt = ornt[c]; 2958 break; 2959 } 2960 } 2961 break; 2962 } 2963 } 2964 } 2965 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2966 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 2967 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2968 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 2969 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2970 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 2971 if (subc >= 0) { 2972 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 2973 orntNew[coneSizeNew] = ornt[c]; 2974 ++coneSizeNew; 2975 } 2976 } 2977 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 2978 if (fornt < 0) { 2979 /* This should be replaced by a call to DMPlexReverseCell() */ 2980 #if 0 2981 ierr = DMPlexReverseCell(subdm, subpoint);CHKERRQ(ierr); 2982 #else 2983 for (c = 0; c < coneSizeNew/2 + coneSizeNew%2; ++c) { 2984 PetscInt faceSize, tmp; 2985 2986 tmp = coneNew[c]; 2987 coneNew[c] = coneNew[coneSizeNew-1-c]; 2988 coneNew[coneSizeNew-1-c] = tmp; 2989 ierr = DMPlexGetConeSize(dm, cone[c], &faceSize);CHKERRQ(ierr); 2990 tmp = orntNew[c] >= 0 ? -(faceSize-orntNew[c]) : faceSize+orntNew[c]; 2991 orntNew[c] = orntNew[coneSizeNew-1-c] >= 0 ? -(faceSize-orntNew[coneSizeNew-1-c]) : faceSize+orntNew[coneSizeNew-1-c]; 2992 orntNew[coneSizeNew-1-c] = tmp; 2993 } 2994 } 2995 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 2996 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 2997 #endif 2998 } 2999 } 3000 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 3001 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3002 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3003 /* Build coordinates */ 3004 { 3005 PetscSection coordSection, subCoordSection; 3006 Vec coordinates, subCoordinates; 3007 PetscScalar *coords, *subCoords; 3008 PetscInt cdim, numComp, coordSize; 3009 const char *name; 3010 3011 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3012 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3013 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3014 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3015 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3016 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3017 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3018 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 3019 for (v = 0; v < numSubPoints[0]; ++v) { 3020 const PetscInt vertex = subpoints[0][v]; 3021 const PetscInt subvertex = firstSubPoint[0]+v; 3022 PetscInt dof; 3023 3024 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3025 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3026 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3027 } 3028 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3029 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3030 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3031 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3032 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3033 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3034 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3035 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3036 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3037 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3038 for (v = 0; v < numSubPoints[0]; ++v) { 3039 const PetscInt vertex = subpoints[0][v]; 3040 const PetscInt subvertex = firstSubPoint[0]+v; 3041 PetscInt dof, off, sdof, soff, d; 3042 3043 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3044 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3045 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3046 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3047 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3048 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3049 } 3050 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3051 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3052 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3053 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3054 } 3055 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 3056 { 3057 PetscSF sfPoint, sfPointSub; 3058 IS subpIS; 3059 const PetscSFNode *remotePoints; 3060 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3061 const PetscInt *localPoints, *subpoints; 3062 PetscInt *slocalPoints; 3063 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 3064 PetscMPIInt rank; 3065 3066 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3067 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3068 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3069 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3070 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 3071 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 3072 if (subpIS) { 3073 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 3074 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 3075 } 3076 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3077 if (numRoots >= 0) { 3078 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3079 for (p = 0; p < pEnd-pStart; ++p) { 3080 newLocalPoints[p].rank = -2; 3081 newLocalPoints[p].index = -2; 3082 } 3083 /* Set subleaves */ 3084 for (l = 0; l < numLeaves; ++l) { 3085 const PetscInt point = localPoints[l]; 3086 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3087 3088 if (subpoint < 0) continue; 3089 newLocalPoints[point-pStart].rank = rank; 3090 newLocalPoints[point-pStart].index = subpoint; 3091 ++numSubleaves; 3092 } 3093 /* Must put in owned subpoints */ 3094 for (p = pStart; p < pEnd; ++p) { 3095 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3096 3097 if (subpoint < 0) { 3098 newOwners[p-pStart].rank = -3; 3099 newOwners[p-pStart].index = -3; 3100 } else { 3101 newOwners[p-pStart].rank = rank; 3102 newOwners[p-pStart].index = subpoint; 3103 } 3104 } 3105 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3106 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3107 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3108 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3109 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 3110 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 3111 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3112 const PetscInt point = localPoints[l]; 3113 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3114 3115 if (subpoint < 0) continue; 3116 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3117 slocalPoints[sl] = subpoint; 3118 sremotePoints[sl].rank = newLocalPoints[point].rank; 3119 sremotePoints[sl].index = newLocalPoints[point].index; 3120 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3121 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3122 ++sl; 3123 } 3124 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3125 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3126 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3127 } 3128 if (subpIS) { 3129 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 3130 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 3131 } 3132 } 3133 /* Cleanup */ 3134 for (d = 0; d <= dim; ++d) { 3135 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3136 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 3137 } 3138 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 3139 PetscFunctionReturn(0); 3140 } 3141 3142 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 3143 { 3144 PetscErrorCode ierr; 3145 3146 PetscFunctionBegin; 3147 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);CHKERRQ(ierr); 3148 PetscFunctionReturn(0); 3149 } 3150 3151 /*@ 3152 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3153 3154 Input Parameters: 3155 + dm - The original mesh 3156 . vertexLabel - The DMLabel marking vertices contained in the surface 3157 - value - The label value to use 3158 3159 Output Parameter: 3160 . subdm - The surface mesh 3161 3162 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3163 3164 Level: developer 3165 3166 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3167 @*/ 3168 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 3169 { 3170 PetscInt dim, cdim, depth; 3171 PetscErrorCode ierr; 3172 3173 PetscFunctionBegin; 3174 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3175 PetscValidPointer(subdm, 3); 3176 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3177 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3178 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3179 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3180 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3181 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3182 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3183 if (depth == dim) { 3184 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3185 } else { 3186 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3187 } 3188 PetscFunctionReturn(0); 3189 } 3190 3191 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3192 { 3193 MPI_Comm comm; 3194 DMLabel subpointMap; 3195 IS subvertexIS; 3196 const PetscInt *subVertices; 3197 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3198 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3199 PetscInt cMax, c, f; 3200 PetscErrorCode ierr; 3201 3202 PetscFunctionBegin; 3203 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 3204 /* Create subpointMap which marks the submesh */ 3205 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 3206 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3207 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3208 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 3209 /* Setup chart */ 3210 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 3211 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 3212 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 3213 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 3214 /* Set cone sizes */ 3215 firstSubVertex = numSubCells; 3216 firstSubFace = numSubCells+numSubVertices; 3217 newFacePoint = firstSubFace; 3218 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 3219 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3220 for (c = 0; c < numSubCells; ++c) { 3221 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 3222 } 3223 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3224 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 3225 } 3226 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3227 /* Create face cones */ 3228 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3229 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3230 ierr = DMGetWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 3231 for (c = 0; c < numSubCells; ++c) { 3232 const PetscInt cell = subCells[c]; 3233 const PetscInt subcell = c; 3234 const PetscInt *cone, *cells; 3235 PetscInt numCells, subVertex, p, v; 3236 3237 if (cell < cMax) continue; 3238 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 3239 for (v = 0; v < nFV; ++v) { 3240 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 3241 subface[v] = firstSubVertex+subVertex; 3242 } 3243 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 3244 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 3245 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3246 /* Not true in parallel 3247 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3248 for (p = 0; p < numCells; ++p) { 3249 PetscInt negsubcell; 3250 3251 if (cells[p] >= cMax) continue; 3252 /* I know this is a crap search */ 3253 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3254 if (subCells[negsubcell] == cells[p]) break; 3255 } 3256 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3257 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 3258 } 3259 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3260 ++newFacePoint; 3261 } 3262 ierr = DMRestoreWorkArray(subdm, maxConeSize, MPIU_INT, (void**) &subface);CHKERRQ(ierr); 3263 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3264 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3265 /* Build coordinates */ 3266 { 3267 PetscSection coordSection, subCoordSection; 3268 Vec coordinates, subCoordinates; 3269 PetscScalar *coords, *subCoords; 3270 PetscInt cdim, numComp, coordSize, v; 3271 const char *name; 3272 3273 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3274 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3275 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3276 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3277 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3278 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3279 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3280 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3281 for (v = 0; v < numSubVertices; ++v) { 3282 const PetscInt vertex = subVertices[v]; 3283 const PetscInt subvertex = firstSubVertex+v; 3284 PetscInt dof; 3285 3286 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3287 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3288 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3289 } 3290 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3291 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3292 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3293 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3294 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3295 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3296 ierr = VecSetBlockSize(subCoordinates, cdim);CHKERRQ(ierr); 3297 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3298 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3299 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3300 for (v = 0; v < numSubVertices; ++v) { 3301 const PetscInt vertex = subVertices[v]; 3302 const PetscInt subvertex = firstSubVertex+v; 3303 PetscInt dof, off, sdof, soff, d; 3304 3305 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3306 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3307 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3308 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3309 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3310 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3311 } 3312 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3313 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3314 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3315 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3316 } 3317 /* Build SF */ 3318 CHKMEMQ; 3319 { 3320 PetscSF sfPoint, sfPointSub; 3321 const PetscSFNode *remotePoints; 3322 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3323 const PetscInt *localPoints; 3324 PetscInt *slocalPoints; 3325 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3326 PetscMPIInt rank; 3327 3328 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3329 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3330 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3331 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3332 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3333 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3334 if (numRoots >= 0) { 3335 /* Only vertices should be shared */ 3336 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3337 for (p = 0; p < pEnd-pStart; ++p) { 3338 newLocalPoints[p].rank = -2; 3339 newLocalPoints[p].index = -2; 3340 } 3341 /* Set subleaves */ 3342 for (l = 0; l < numLeaves; ++l) { 3343 const PetscInt point = localPoints[l]; 3344 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3345 3346 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3347 if (subPoint < 0) continue; 3348 newLocalPoints[point-pStart].rank = rank; 3349 newLocalPoints[point-pStart].index = subPoint; 3350 ++numSubLeaves; 3351 } 3352 /* Must put in owned subpoints */ 3353 for (p = pStart; p < pEnd; ++p) { 3354 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3355 3356 if (subPoint < 0) { 3357 newOwners[p-pStart].rank = -3; 3358 newOwners[p-pStart].index = -3; 3359 } else { 3360 newOwners[p-pStart].rank = rank; 3361 newOwners[p-pStart].index = subPoint; 3362 } 3363 } 3364 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3365 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3366 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3367 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3368 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3369 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3370 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3371 const PetscInt point = localPoints[l]; 3372 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3373 3374 if (subPoint < 0) continue; 3375 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3376 slocalPoints[sl] = subPoint; 3377 sremotePoints[sl].rank = newLocalPoints[point].rank; 3378 sremotePoints[sl].index = newLocalPoints[point].index; 3379 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3380 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3381 ++sl; 3382 } 3383 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3384 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3385 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3386 } 3387 } 3388 CHKMEMQ; 3389 /* Cleanup */ 3390 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3391 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3392 ierr = PetscFree(subCells);CHKERRQ(ierr); 3393 PetscFunctionReturn(0); 3394 } 3395 3396 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3397 { 3398 DMLabel label = NULL; 3399 PetscErrorCode ierr; 3400 3401 PetscFunctionBegin; 3402 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3403 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);CHKERRQ(ierr); 3404 PetscFunctionReturn(0); 3405 } 3406 3407 /*@C 3408 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. 3409 3410 Input Parameters: 3411 + dm - The original mesh 3412 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3413 . label - A label name, or NULL 3414 - value - A label value 3415 3416 Output Parameter: 3417 . subdm - The surface mesh 3418 3419 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3420 3421 Level: developer 3422 3423 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3424 @*/ 3425 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3426 { 3427 PetscInt dim, cdim, depth; 3428 PetscErrorCode ierr; 3429 3430 PetscFunctionBegin; 3431 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3432 PetscValidPointer(subdm, 5); 3433 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3434 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3435 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3436 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3437 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3438 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 3439 ierr = DMSetCoordinateDim(*subdm, cdim);CHKERRQ(ierr); 3440 if (depth == dim) { 3441 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3442 } else { 3443 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3444 } 3445 PetscFunctionReturn(0); 3446 } 3447 3448 /*@ 3449 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3450 3451 Input Parameters: 3452 + dm - The original mesh 3453 . cellLabel - The DMLabel marking cells contained in the new mesh 3454 - value - The label value to use 3455 3456 Output Parameter: 3457 . subdm - The new mesh 3458 3459 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3460 3461 Level: developer 3462 3463 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3464 @*/ 3465 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3466 { 3467 PetscInt dim; 3468 PetscErrorCode ierr; 3469 3470 PetscFunctionBegin; 3471 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3472 PetscValidPointer(subdm, 3); 3473 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3474 ierr = DMCreate(PetscObjectComm((PetscObject) dm), subdm);CHKERRQ(ierr); 3475 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3476 ierr = DMSetDimension(*subdm, dim);CHKERRQ(ierr); 3477 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3478 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);CHKERRQ(ierr); 3479 PetscFunctionReturn(0); 3480 } 3481 3482 /*@ 3483 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3484 3485 Input Parameter: 3486 . dm - The submesh DM 3487 3488 Output Parameter: 3489 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3490 3491 Level: developer 3492 3493 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3494 @*/ 3495 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3496 { 3497 PetscFunctionBegin; 3498 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3499 PetscValidPointer(subpointMap, 2); 3500 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3501 PetscFunctionReturn(0); 3502 } 3503 3504 /*@ 3505 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3506 3507 Input Parameters: 3508 + dm - The submesh DM 3509 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3510 3511 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3512 3513 Level: developer 3514 3515 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3516 @*/ 3517 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3518 { 3519 DM_Plex *mesh = (DM_Plex *) dm->data; 3520 DMLabel tmp; 3521 PetscErrorCode ierr; 3522 3523 PetscFunctionBegin; 3524 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3525 tmp = mesh->subpointMap; 3526 mesh->subpointMap = subpointMap; 3527 ++mesh->subpointMap->refct; 3528 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3529 PetscFunctionReturn(0); 3530 } 3531 3532 /*@ 3533 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3534 3535 Input Parameter: 3536 . dm - The submesh DM 3537 3538 Output Parameter: 3539 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3540 3541 Note: This IS is guaranteed to be sorted by the construction of the submesh 3542 3543 Level: developer 3544 3545 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3546 @*/ 3547 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3548 { 3549 MPI_Comm comm; 3550 DMLabel subpointMap; 3551 IS is; 3552 const PetscInt *opoints; 3553 PetscInt *points, *depths; 3554 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3555 PetscErrorCode ierr; 3556 3557 PetscFunctionBegin; 3558 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3559 PetscValidPointer(subpointIS, 2); 3560 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3561 *subpointIS = NULL; 3562 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3563 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3564 if (subpointMap && depth >= 0) { 3565 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3566 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3567 ierr = DMGetWorkArray(dm, depth+1, MPIU_INT, &depths);CHKERRQ(ierr); 3568 depths[0] = depth; 3569 depths[1] = 0; 3570 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3571 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3572 for(d = 0, off = 0; d <= depth; ++d) { 3573 const PetscInt dep = depths[d]; 3574 3575 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3576 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3577 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3578 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); 3579 } else { 3580 if (!n) { 3581 if (d == 0) { 3582 /* Missing cells */ 3583 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3584 } else { 3585 /* Missing faces */ 3586 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3587 } 3588 } 3589 } 3590 if (n) { 3591 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3592 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3593 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3594 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3595 ierr = ISDestroy(&is);CHKERRQ(ierr); 3596 } 3597 } 3598 ierr = DMRestoreWorkArray(dm, depth+1, MPIU_INT, &depths);CHKERRQ(ierr); 3599 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3600 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3601 } 3602 PetscFunctionReturn(0); 3603 } 3604