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