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 const PetscInt *values; 50 PetscInt numValues, v, cStart, cEnd; 51 PetscMPIInt size; 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 = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);CHKERRQ(ierr); 87 if (size != 1) { 88 PetscSF sfPoint; 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 = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 408 /* Step 8: Convert coordinates */ 409 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 410 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 411 ierr = DMPlexGetDepthStratum(dmNew, 0, &vStartNew, &vEndNew);CHKERRQ(ierr); 412 ierr = DMPlexGetHeightStratum(dmNew, 0, &cStartNew, &cEndNew);CHKERRQ(ierr); 413 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 414 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);CHKERRQ(ierr); 415 ierr = PetscSectionSetNumFields(newCoordSection, 1);CHKERRQ(ierr); 416 ierr = PetscSectionSetFieldComponents(newCoordSection, 0, dim);CHKERRQ(ierr); 417 ierr = PetscSectionGetChart(coordSection, &sStart, &sEnd);CHKERRQ(ierr); 418 hasCells = sStart == cStart ? PETSC_TRUE : PETSC_FALSE; 419 ierr = PetscSectionSetChart(newCoordSection, hasCells ? cStartNew : vStartNew, vEndNew);CHKERRQ(ierr); 420 if (hasCells) { 421 for (c = cStart; c < cEnd; ++c) { 422 PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof; 423 424 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 425 ierr = PetscSectionSetDof(newCoordSection, cNew, dof);CHKERRQ(ierr); 426 ierr = PetscSectionSetFieldDof(newCoordSection, cNew, 0, dof);CHKERRQ(ierr); 427 } 428 } 429 for (v = vStartNew; v < vEndNew; ++v) { 430 ierr = PetscSectionSetDof(newCoordSection, v, dim);CHKERRQ(ierr); 431 ierr = PetscSectionSetFieldDof(newCoordSection, v, 0, dim);CHKERRQ(ierr); 432 } 433 ierr = PetscSectionSetUp(newCoordSection);CHKERRQ(ierr); 434 ierr = DMSetCoordinateSection(dmNew, PETSC_DETERMINE, newCoordSection);CHKERRQ(ierr); 435 ierr = PetscSectionGetStorageSize(newCoordSection, &coordSize);CHKERRQ(ierr); 436 ierr = VecCreate(PETSC_COMM_SELF, &newCoordinates);CHKERRQ(ierr); 437 ierr = PetscObjectSetName((PetscObject) newCoordinates, "coordinates");CHKERRQ(ierr); 438 ierr = VecSetSizes(newCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 439 ierr = VecSetBlockSize(newCoordinates, dim);CHKERRQ(ierr); 440 ierr = VecSetType(newCoordinates,VECSTANDARD);CHKERRQ(ierr); 441 ierr = DMSetCoordinatesLocal(dmNew, newCoordinates);CHKERRQ(ierr); 442 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 443 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 444 ierr = VecGetArray(newCoordinates, &newCoords);CHKERRQ(ierr); 445 if (hasCells) { 446 for (c = cStart; c < cEnd; ++c) { 447 PetscInt cNew = DMPlexShiftPoint_Internal(c, depth, depthShift), dof, off, noff, d; 448 449 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 450 ierr = PetscSectionGetOffset(coordSection, c, &off);CHKERRQ(ierr); 451 ierr = PetscSectionGetOffset(newCoordSection, cNew, &noff);CHKERRQ(ierr); 452 for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d]; 453 } 454 } 455 for (v = vStart; v < vEnd; ++v) { 456 PetscInt dof, off, noff, d; 457 458 ierr = PetscSectionGetDof(coordSection, v, &dof);CHKERRQ(ierr); 459 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 460 ierr = PetscSectionGetOffset(newCoordSection, DMPlexShiftPoint_Internal(v, depth, depthShift), &noff);CHKERRQ(ierr); 461 for (d = 0; d < dof; ++d) newCoords[noff+d] = coords[off+d]; 462 } 463 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 464 ierr = VecRestoreArray(newCoordinates, &newCoords);CHKERRQ(ierr); 465 ierr = VecDestroy(&newCoordinates);CHKERRQ(ierr); 466 ierr = PetscSectionDestroy(&newCoordSection);CHKERRQ(ierr); 467 PetscFunctionReturn(0); 468 } 469 470 static PetscErrorCode DMPlexShiftSF_Internal(DM dm, PetscInt depthShift[], DM dmNew) 471 { 472 PetscInt depth = 0; 473 PetscSF sfPoint, sfPointNew; 474 const PetscSFNode *remotePoints; 475 PetscSFNode *gremotePoints; 476 const PetscInt *localPoints; 477 PetscInt *glocalPoints, *newLocation, *newRemoteLocation; 478 PetscInt numRoots, numLeaves, l, pStart, pEnd, totShift = 0; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 483 /* Step 9: Convert pointSF */ 484 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 485 ierr = DMGetPointSF(dmNew, &sfPointNew);CHKERRQ(ierr); 486 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 487 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 488 totShift = DMPlexShiftPoint_Internal(pEnd,depth,depthShift) - pEnd; 489 if (numRoots >= 0) { 490 ierr = PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);CHKERRQ(ierr); 491 for (l=0; l<numRoots; l++) newLocation[l] = DMPlexShiftPoint_Internal(l, depth, depthShift); 492 ierr = PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);CHKERRQ(ierr); 493 ierr = PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);CHKERRQ(ierr); 494 ierr = PetscMalloc1(numLeaves, &glocalPoints);CHKERRQ(ierr); 495 ierr = PetscMalloc1(numLeaves, &gremotePoints);CHKERRQ(ierr); 496 for (l = 0; l < numLeaves; ++l) { 497 glocalPoints[l] = DMPlexShiftPoint_Internal(localPoints[l], depth, depthShift); 498 gremotePoints[l].rank = remotePoints[l].rank; 499 gremotePoints[l].index = newRemoteLocation[localPoints[l]]; 500 } 501 ierr = PetscFree2(newLocation,newRemoteLocation);CHKERRQ(ierr); 502 ierr = PetscSFSetGraph(sfPointNew, numRoots + totShift, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 503 } 504 PetscFunctionReturn(0); 505 } 506 507 static PetscErrorCode DMPlexShiftLabels_Internal(DM dm, PetscInt depthShift[], DM dmNew) 508 { 509 PetscSF sfPoint; 510 DMLabel vtkLabel, ghostLabel; 511 const PetscSFNode *leafRemote; 512 const PetscInt *leafLocal; 513 PetscInt depth = 0, numLeaves, numLabels, l, cStart, cEnd, c, fStart, fEnd, f; 514 PetscMPIInt rank; 515 PetscErrorCode ierr; 516 517 PetscFunctionBegin; 518 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 519 /* Step 10: Convert labels */ 520 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 521 for (l = 0; l < numLabels; ++l) { 522 DMLabel label, newlabel; 523 const char *lname; 524 PetscBool isDepth; 525 IS valueIS; 526 const PetscInt *values; 527 PetscInt numValues, val; 528 529 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr); 530 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 531 if (isDepth) continue; 532 ierr = DMCreateLabel(dmNew, lname);CHKERRQ(ierr); 533 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr); 534 ierr = DMGetLabel(dmNew, lname, &newlabel);CHKERRQ(ierr); 535 ierr = DMLabelGetDefaultValue(label,&val);CHKERRQ(ierr); 536 ierr = DMLabelSetDefaultValue(newlabel,val);CHKERRQ(ierr); 537 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 538 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 539 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 540 for (val = 0; val < numValues; ++val) { 541 IS pointIS; 542 const PetscInt *points; 543 PetscInt numPoints, p; 544 545 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 546 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 547 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 548 for (p = 0; p < numPoints; ++p) { 549 const PetscInt newpoint = DMPlexShiftPoint_Internal(points[p], depth, depthShift); 550 551 ierr = DMLabelSetValue(newlabel, newpoint, values[val]);CHKERRQ(ierr); 552 } 553 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 554 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 555 } 556 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 557 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 558 } 559 /* Step 11: Make label for output (vtk) and to mark ghost points (ghost) */ 560 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 561 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 562 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 563 ierr = PetscSFGetGraph(sfPoint, NULL, &numLeaves, &leafLocal, &leafRemote);CHKERRQ(ierr); 564 ierr = DMCreateLabel(dmNew, "vtk");CHKERRQ(ierr); 565 ierr = DMCreateLabel(dmNew, "ghost");CHKERRQ(ierr); 566 ierr = DMGetLabel(dmNew, "vtk", &vtkLabel);CHKERRQ(ierr); 567 ierr = DMGetLabel(dmNew, "ghost", &ghostLabel);CHKERRQ(ierr); 568 for (l = 0, c = cStart; l < numLeaves && c < cEnd; ++l, ++c) { 569 for (; c < leafLocal[l] && c < cEnd; ++c) { 570 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 571 } 572 if (leafLocal[l] >= cEnd) break; 573 if (leafRemote[l].rank == rank) { 574 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 575 } else { 576 ierr = DMLabelSetValue(ghostLabel, c, 2);CHKERRQ(ierr); 577 } 578 } 579 for (; c < cEnd; ++c) { 580 ierr = DMLabelSetValue(vtkLabel, c, 1);CHKERRQ(ierr); 581 } 582 if (0) { 583 ierr = DMLabelView(vtkLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 584 } 585 ierr = DMPlexGetHeightStratum(dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 586 for (f = fStart; f < fEnd; ++f) { 587 PetscInt numCells; 588 589 ierr = DMPlexGetSupportSize(dmNew, f, &numCells);CHKERRQ(ierr); 590 if (numCells < 2) { 591 ierr = DMLabelSetValue(ghostLabel, f, 1);CHKERRQ(ierr); 592 } else { 593 const PetscInt *cells = NULL; 594 PetscInt vA, vB; 595 596 ierr = DMPlexGetSupport(dmNew, f, &cells);CHKERRQ(ierr); 597 ierr = DMLabelGetValue(vtkLabel, cells[0], &vA);CHKERRQ(ierr); 598 ierr = DMLabelGetValue(vtkLabel, cells[1], &vB);CHKERRQ(ierr); 599 if (vA != 1 && vB != 1) {ierr = DMLabelSetValue(ghostLabel, f, 1);CHKERRQ(ierr);} 600 } 601 } 602 if (0) { 603 ierr = DMLabelView(ghostLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 604 } 605 PetscFunctionReturn(0); 606 } 607 608 static PetscErrorCode DMPlexShiftTree_Internal(DM dm, PetscInt depthShift[], DM dmNew) 609 { 610 DM refTree; 611 PetscSection pSec; 612 PetscInt *parents, *childIDs; 613 PetscErrorCode ierr; 614 615 PetscFunctionBegin; 616 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 617 ierr = DMPlexSetReferenceTree(dmNew,refTree);CHKERRQ(ierr); 618 ierr = DMPlexGetTree(dm,&pSec,&parents,&childIDs,NULL,NULL);CHKERRQ(ierr); 619 if (pSec) { 620 PetscInt p, pStart, pEnd, *parentsShifted, pStartShifted, pEndShifted, depth; 621 PetscInt *childIDsShifted; 622 PetscSection pSecShifted; 623 624 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 625 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 626 pStartShifted = DMPlexShiftPoint_Internal(pStart,depth,depthShift); 627 pEndShifted = DMPlexShiftPoint_Internal(pEnd,depth,depthShift); 628 ierr = PetscMalloc2(pEndShifted - pStartShifted,&parentsShifted,pEndShifted-pStartShifted,&childIDsShifted);CHKERRQ(ierr); 629 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dmNew),&pSecShifted);CHKERRQ(ierr); 630 ierr = PetscSectionSetChart(pSecShifted,pStartShifted,pEndShifted);CHKERRQ(ierr); 631 for (p = pStartShifted; p < pEndShifted; p++) { 632 /* start off assuming no children */ 633 ierr = PetscSectionSetDof(pSecShifted,p,0);CHKERRQ(ierr); 634 } 635 for (p = pStart; p < pEnd; p++) { 636 PetscInt dof; 637 PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift); 638 639 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 640 ierr = PetscSectionSetDof(pSecShifted,pNew,dof);CHKERRQ(ierr); 641 } 642 ierr = PetscSectionSetUp(pSecShifted);CHKERRQ(ierr); 643 for (p = pStart; p < pEnd; p++) { 644 PetscInt dof; 645 PetscInt pNew = DMPlexShiftPoint_Internal(p,depth,depthShift); 646 647 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 648 if (dof) { 649 PetscInt off, offNew; 650 651 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 652 ierr = PetscSectionGetOffset(pSecShifted,pNew,&offNew);CHKERRQ(ierr); 653 parentsShifted[offNew] = DMPlexShiftPoint_Internal(parents[off],depth,depthShift); 654 childIDsShifted[offNew] = childIDs[off]; 655 } 656 } 657 ierr = DMPlexSetTree(dmNew,pSecShifted,parentsShifted,childIDsShifted);CHKERRQ(ierr); 658 ierr = PetscFree2(parentsShifted,childIDsShifted);CHKERRQ(ierr); 659 ierr = PetscSectionDestroy(&pSecShifted);CHKERRQ(ierr); 660 } 661 PetscFunctionReturn(0); 662 } 663 664 static PetscErrorCode DMPlexConstructGhostCells_Internal(DM dm, DMLabel label, PetscInt *numGhostCells, DM gdm) 665 { 666 PetscSF sf; 667 IS valueIS; 668 const PetscInt *values, *leaves; 669 PetscInt *depthShift; 670 PetscInt d, depth = 0, nleaves, loc, Ng, numFS, fs, fStart, fEnd, ghostCell, cEnd, c; 671 PetscErrorCode ierr; 672 673 PetscFunctionBegin; 674 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 675 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 676 nleaves = PetscMax(0, nleaves); 677 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 678 /* Count ghost cells */ 679 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 680 ierr = ISGetLocalSize(valueIS, &numFS);CHKERRQ(ierr); 681 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 682 Ng = 0; 683 for (fs = 0; fs < numFS; ++fs) { 684 IS faceIS; 685 const PetscInt *faces; 686 PetscInt numFaces, f, numBdFaces = 0; 687 688 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 689 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 690 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 691 for (f = 0; f < numFaces; ++f) { 692 PetscInt numChildren; 693 694 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 695 ierr = DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);CHKERRQ(ierr); 696 /* non-local and ancestors points don't get to register ghosts */ 697 if (loc >= 0 || numChildren) continue; 698 if ((faces[f] >= fStart) && (faces[f] < fEnd)) ++numBdFaces; 699 } 700 Ng += numBdFaces; 701 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 702 } 703 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 704 ierr = PetscMalloc1(2*(depth+1), &depthShift);CHKERRQ(ierr); 705 for (d = 0; d <= depth; d++) { 706 PetscInt dEnd; 707 708 ierr = DMPlexGetDepthStratum(dm,d,NULL,&dEnd);CHKERRQ(ierr); 709 depthShift[2*d] = dEnd; 710 depthShift[2*d+1] = 0; 711 } 712 if (depth >= 0) depthShift[2*depth+1] = Ng; 713 ierr = DMPlexShiftPointSetUp_Internal(depth,depthShift);CHKERRQ(ierr); 714 ierr = DMPlexShiftSizes_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 715 /* Step 3: Set cone/support sizes for new points */ 716 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 717 ierr = DMPlexSetHybridBounds(gdm, cEnd, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 718 for (c = cEnd; c < cEnd + Ng; ++c) { 719 ierr = DMPlexSetConeSize(gdm, c, 1);CHKERRQ(ierr); 720 } 721 for (fs = 0; fs < numFS; ++fs) { 722 IS faceIS; 723 const PetscInt *faces; 724 PetscInt numFaces, f; 725 726 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 727 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 728 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 729 for (f = 0; f < numFaces; ++f) { 730 PetscInt size, numChildren; 731 732 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 733 ierr = DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);CHKERRQ(ierr); 734 if (loc >= 0 || numChildren) continue; 735 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 736 ierr = DMPlexGetSupportSize(dm, faces[f], &size);CHKERRQ(ierr); 737 if (size != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has boundary face %d with %d support cells", faces[f], size); 738 ierr = DMPlexSetSupportSize(gdm, faces[f] + Ng, 2);CHKERRQ(ierr); 739 } 740 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 741 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 742 } 743 /* Step 4: Setup ghosted DM */ 744 ierr = DMSetUp(gdm);CHKERRQ(ierr); 745 ierr = DMPlexShiftPoints_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 746 /* Step 6: Set cones and supports for new points */ 747 ghostCell = cEnd; 748 for (fs = 0; fs < numFS; ++fs) { 749 IS faceIS; 750 const PetscInt *faces; 751 PetscInt numFaces, f; 752 753 ierr = DMLabelGetStratumIS(label, values[fs], &faceIS);CHKERRQ(ierr); 754 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 755 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 756 for (f = 0; f < numFaces; ++f) { 757 PetscInt newFace = faces[f] + Ng, numChildren; 758 759 ierr = PetscFindInt(faces[f], nleaves, leaves, &loc);CHKERRQ(ierr); 760 ierr = DMPlexGetTreeChildren(dm,faces[f],&numChildren,NULL);CHKERRQ(ierr); 761 if (loc >= 0 || numChildren) continue; 762 if ((faces[f] < fStart) || (faces[f] >= fEnd)) continue; 763 ierr = DMPlexSetCone(gdm, ghostCell, &newFace);CHKERRQ(ierr); 764 ierr = DMPlexInsertSupport(gdm, newFace, 1, ghostCell);CHKERRQ(ierr); 765 ++ghostCell; 766 } 767 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 768 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 769 } 770 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 771 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 772 /* Step 7: Stratify */ 773 ierr = DMPlexStratify(gdm);CHKERRQ(ierr); 774 ierr = DMPlexShiftCoordinates_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 775 ierr = DMPlexShiftSF_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 776 ierr = DMPlexShiftLabels_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 777 ierr = DMPlexShiftTree_Internal(dm, depthShift, gdm);CHKERRQ(ierr); 778 ierr = PetscFree(depthShift);CHKERRQ(ierr); 779 /* Step 7: Periodicity */ 780 if (dm->maxCell) { 781 const PetscReal *maxCell, *L; 782 const DMBoundaryType *bd; 783 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 784 ierr = DMSetPeriodicity(gdm, maxCell, L, bd);CHKERRQ(ierr); 785 } 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 #if 0 1444 /* Do not put cohesive edges into the label */ 1445 if (dep == 0) { 1446 const PetscInt cedge = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 1447 ierr = DMLabelSetValue(mlabel, cedge, val);CHKERRQ(ierr); 1448 } else if (dep == dim-2) { 1449 const PetscInt cface = p + pMaxNew[dep+1] + numSplitPoints[dep+1]; 1450 ierr = DMLabelSetValue(mlabel, cface, val);CHKERRQ(ierr); 1451 } 1452 /* Do not put cohesive faces into the label */ 1453 #endif 1454 } 1455 } 1456 } 1457 } 1458 for (sp = 0; sp < numSP; ++sp) { 1459 const PetscInt dep = values[sp]; 1460 1461 if ((dep < 0) || (dep > depth)) continue; 1462 if (splitIS[dep]) {ierr = ISRestoreIndices(splitIS[dep], &splitPoints[dep]);CHKERRQ(ierr);} 1463 ierr = ISDestroy(&splitIS[dep]);CHKERRQ(ierr); 1464 if (unsplitIS[dep]) {ierr = ISRestoreIndices(unsplitIS[dep], &unsplitPoints[dep]);CHKERRQ(ierr);} 1465 ierr = ISDestroy(&unsplitIS[dep]);CHKERRQ(ierr); 1466 } 1467 if (label) { 1468 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 1469 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 1470 } 1471 for (d = 0; d <= depth; ++d) { 1472 ierr = DMPlexGetDepthStratum(sdm, d, NULL, &pEnd);CHKERRQ(ierr); 1473 pMaxNew[d] = pEnd - numHybridPoints[d] - numHybridPointsOld[d]; 1474 } 1475 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); 1476 ierr = PetscFree3(coneNew, coneONew, supportNew);CHKERRQ(ierr); 1477 ierr = PetscFree5(depthMax, depthEnd, depthShift, pMaxNew, numHybridPointsOld);CHKERRQ(ierr); 1478 ierr = PetscFree7(splitIS, unsplitIS, numSplitPoints, numUnsplitPoints, numHybridPoints, splitPoints, unsplitPoints);CHKERRQ(ierr); 1479 PetscFunctionReturn(0); 1480 } 1481 1482 /*@C 1483 DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface 1484 1485 Collective on dm 1486 1487 Input Parameters: 1488 + dm - The original DM 1489 - label - The label specifying the boundary faces (this could be auto-generated) 1490 1491 Output Parameters: 1492 - dmSplit - The new DM 1493 1494 Level: developer 1495 1496 .seealso: DMCreate(), DMPlexLabelCohesiveComplete() 1497 @*/ 1498 PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit) 1499 { 1500 DM sdm; 1501 PetscInt dim; 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1506 PetscValidPointer(dmSplit, 3); 1507 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &sdm);CHKERRQ(ierr); 1508 ierr = DMSetType(sdm, DMPLEX);CHKERRQ(ierr); 1509 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1510 ierr = DMSetDimension(sdm, dim);CHKERRQ(ierr); 1511 switch (dim) { 1512 case 2: 1513 case 3: 1514 ierr = DMPlexConstructCohesiveCells_Internal(dm, label, sdm);CHKERRQ(ierr); 1515 break; 1516 default: 1517 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct cohesive cells for dimension %d", dim); 1518 } 1519 *dmSplit = sdm; 1520 PetscFunctionReturn(0); 1521 } 1522 1523 /* Returns the side of the surface for a given cell with a face on the surface */ 1524 static PetscErrorCode GetSurfaceSide_Static(DM dm, DM subdm, PetscInt numSubpoints, const PetscInt *subpoints, PetscInt cell, PetscInt face, PetscBool *pos) 1525 { 1526 const PetscInt *cone, *ornt; 1527 PetscInt dim, coneSize, c; 1528 PetscErrorCode ierr; 1529 1530 PetscFunctionBegin; 1531 *pos = PETSC_TRUE; 1532 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1533 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 1534 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 1535 ierr = DMPlexGetConeOrientation(dm, cell, &ornt);CHKERRQ(ierr); 1536 for (c = 0; c < coneSize; ++c) { 1537 if (cone[c] == face) { 1538 PetscInt o = ornt[c]; 1539 1540 if (subdm) { 1541 const PetscInt *subcone, *subornt; 1542 PetscInt subpoint, subface, subconeSize, sc; 1543 1544 ierr = PetscFindInt(cell, numSubpoints, subpoints, &subpoint);CHKERRQ(ierr); 1545 ierr = PetscFindInt(face, numSubpoints, subpoints, &subface);CHKERRQ(ierr); 1546 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 1547 ierr = DMPlexGetCone(subdm, subpoint, &subcone);CHKERRQ(ierr); 1548 ierr = DMPlexGetConeOrientation(subdm, subpoint, &subornt);CHKERRQ(ierr); 1549 for (sc = 0; sc < subconeSize; ++sc) { 1550 if (subcone[sc] == subface) { 1551 o = subornt[0]; 1552 break; 1553 } 1554 } 1555 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); 1556 } 1557 if (o >= 0) *pos = PETSC_TRUE; 1558 else *pos = PETSC_FALSE; 1559 break; 1560 } 1561 } 1562 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); 1563 PetscFunctionReturn(0); 1564 } 1565 1566 /*@ 1567 DMPlexLabelCohesiveComplete - Starting with a label marking points on an internal surface, we add all other mesh pieces 1568 to complete the surface 1569 1570 Input Parameters: 1571 + dm - The DM 1572 . label - A DMLabel marking the surface 1573 . blabel - A DMLabel marking the vertices on the boundary which will not be duplicated, or NULL to find them automatically 1574 . flip - Flag to flip the submesh normal and replace points on the other side 1575 - subdm - The subDM associated with the label, or NULL 1576 1577 Output Parameter: 1578 . label - A DMLabel marking all surface points 1579 1580 Note: The vertices in blabel are called "unsplit" in the terminology from hybrid cell creation. 1581 1582 Level: developer 1583 1584 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelComplete() 1585 @*/ 1586 PetscErrorCode DMPlexLabelCohesiveComplete(DM dm, DMLabel label, DMLabel blabel, PetscBool flip, DM subdm) 1587 { 1588 DMLabel depthLabel; 1589 IS dimIS, subpointIS = NULL, facePosIS, faceNegIS, crossEdgeIS = NULL; 1590 const PetscInt *points, *subpoints; 1591 const PetscInt rev = flip ? -1 : 1; 1592 PetscInt *pMax; 1593 PetscInt shift = 100, shift2 = 200, dim, depth, pSize, dep, cStart, cEnd, cMax, fStart, fEnd, vStart, vEnd, numPoints, numSubpoints, p, val; 1594 PetscErrorCode ierr; 1595 1596 PetscFunctionBegin; 1597 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1598 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1599 pSize = PetscMax(depth, dim) + 1; 1600 ierr = PetscMalloc1(pSize,&pMax);CHKERRQ(ierr); 1601 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1602 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1603 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1604 if (subdm) { 1605 ierr = DMPlexCreateSubpointIS(subdm, &subpointIS);CHKERRQ(ierr); 1606 if (subpointIS) { 1607 ierr = ISGetLocalSize(subpointIS, &numSubpoints);CHKERRQ(ierr); 1608 ierr = ISGetIndices(subpointIS, &subpoints);CHKERRQ(ierr); 1609 } 1610 } 1611 /* Mark cell on the fault, and its faces which touch the fault: cell orientation for face gives the side of the fault */ 1612 ierr = DMLabelGetStratumIS(label, dim-1, &dimIS);CHKERRQ(ierr); 1613 if (!dimIS) { 1614 ierr = PetscFree(pMax);CHKERRQ(ierr); 1615 ierr = ISDestroy(&subpointIS);CHKERRQ(ierr); 1616 PetscFunctionReturn(0); 1617 } 1618 ierr = ISGetLocalSize(dimIS, &numPoints);CHKERRQ(ierr); 1619 ierr = ISGetIndices(dimIS, &points);CHKERRQ(ierr); 1620 for (p = 0; p < numPoints; ++p) { /* Loop over fault faces */ 1621 const PetscInt *support; 1622 PetscInt supportSize, s; 1623 1624 ierr = DMPlexGetSupportSize(dm, points[p], &supportSize);CHKERRQ(ierr); 1625 if (supportSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Split face %d has %d != 2 supports", points[p], supportSize); 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 /*@ 1887 DMPlexCreateHybridMesh - Create a mesh with hybrid cells along an internal interface 1888 1889 Collective on dm 1890 1891 Input Parameters: 1892 + dm - The original DM 1893 - labelName - The label specifying the interface vertices 1894 1895 Output Parameters: 1896 + hybridLabel - The label fully marking the interface 1897 - dmHybrid - The new DM 1898 1899 Level: developer 1900 1901 .seealso: DMPlexConstructCohesiveCells(), DMPlexLabelCohesiveComplete(), DMCreate() 1902 @*/ 1903 PetscErrorCode DMPlexCreateHybridMesh(DM dm, DMLabel label, DMLabel *hybridLabel, DM *dmHybrid) 1904 { 1905 DM idm; 1906 DMLabel subpointMap, hlabel; 1907 PetscInt dim; 1908 PetscErrorCode ierr; 1909 1910 PetscFunctionBegin; 1911 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1912 if (hybridLabel) PetscValidPointer(hybridLabel, 3); 1913 PetscValidPointer(dmHybrid, 4); 1914 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1915 ierr = DMPlexCreateSubmesh(dm, label, 1, &idm);CHKERRQ(ierr); 1916 ierr = DMPlexOrient(idm);CHKERRQ(ierr); 1917 ierr = DMPlexGetSubpointMap(idm, &subpointMap);CHKERRQ(ierr); 1918 ierr = DMLabelDuplicate(subpointMap, &hlabel);CHKERRQ(ierr); 1919 ierr = DMLabelClearStratum(hlabel, dim);CHKERRQ(ierr); 1920 ierr = DMPlexLabelCohesiveComplete(dm, hlabel, NULL, PETSC_FALSE, idm);CHKERRQ(ierr); 1921 ierr = DMDestroy(&idm);CHKERRQ(ierr); 1922 ierr = DMPlexConstructCohesiveCells(dm, hlabel, dmHybrid);CHKERRQ(ierr); 1923 if (hybridLabel) *hybridLabel = hlabel; 1924 else {ierr = DMLabelDestroy(&hlabel);CHKERRQ(ierr);} 1925 PetscFunctionReturn(0); 1926 } 1927 1928 /* Here we need the explicit assumption that: 1929 1930 For any marked cell, the marked vertices constitute a single face 1931 */ 1932 static PetscErrorCode DMPlexMarkSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, DM subdm) 1933 { 1934 IS subvertexIS = NULL; 1935 const PetscInt *subvertices; 1936 PetscInt *pStart, *pEnd, *pMax, pSize; 1937 PetscInt depth, dim, d, numSubVerticesInitial = 0, v; 1938 PetscErrorCode ierr; 1939 1940 PetscFunctionBegin; 1941 *numFaces = 0; 1942 *nFV = 0; 1943 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1944 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1945 pSize = PetscMax(depth, dim) + 1; 1946 ierr = PetscMalloc3(pSize,&pStart,pSize,&pEnd,pSize,&pMax);CHKERRQ(ierr); 1947 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 1948 for (d = 0; d <= depth; ++d) { 1949 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 1950 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 1951 } 1952 /* Loop over initial vertices and mark all faces in the collective star() */ 1953 if (vertexLabel) {ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr);} 1954 if (subvertexIS) { 1955 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 1956 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 1957 } 1958 for (v = 0; v < numSubVerticesInitial; ++v) { 1959 const PetscInt vertex = subvertices[v]; 1960 PetscInt *star = NULL; 1961 PetscInt starSize, s, numCells = 0, c; 1962 1963 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1964 for (s = 0; s < starSize*2; s += 2) { 1965 const PetscInt point = star[s]; 1966 if ((point >= pStart[depth]) && (point < pEnd[depth])) star[numCells++] = point; 1967 } 1968 for (c = 0; c < numCells; ++c) { 1969 const PetscInt cell = star[c]; 1970 PetscInt *closure = NULL; 1971 PetscInt closureSize, cl; 1972 PetscInt cellLoc, numCorners = 0, faceSize = 0; 1973 1974 ierr = DMLabelGetValue(subpointMap, cell, &cellLoc);CHKERRQ(ierr); 1975 if (cellLoc == 2) continue; 1976 if (cellLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Cell %d has dimension %d in the surface label", cell, cellLoc); 1977 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1978 for (cl = 0; cl < closureSize*2; cl += 2) { 1979 const PetscInt point = closure[cl]; 1980 PetscInt vertexLoc; 1981 1982 if ((point >= pStart[0]) && (point < pEnd[0])) { 1983 ++numCorners; 1984 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 1985 if (vertexLoc == value) closure[faceSize++] = point; 1986 } 1987 } 1988 if (!(*nFV)) {ierr = DMPlexGetNumFaceVertices(dm, dim, numCorners, nFV);CHKERRQ(ierr);} 1989 if (faceSize > *nFV) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 1990 if (faceSize == *nFV) { 1991 const PetscInt *cells = NULL; 1992 PetscInt numCells, nc; 1993 1994 ++(*numFaces); 1995 for (cl = 0; cl < faceSize; ++cl) { 1996 ierr = DMLabelSetValue(subpointMap, closure[cl], 0);CHKERRQ(ierr); 1997 } 1998 ierr = DMPlexGetJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 1999 for (nc = 0; nc < numCells; ++nc) { 2000 ierr = DMLabelSetValue(subpointMap, cells[nc], 2);CHKERRQ(ierr); 2001 } 2002 ierr = DMPlexRestoreJoin(dm, faceSize, closure, &numCells, &cells);CHKERRQ(ierr); 2003 } 2004 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2005 } 2006 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2007 } 2008 if (subvertexIS) { 2009 ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2010 } 2011 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2012 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 2013 PetscFunctionReturn(0); 2014 } 2015 2016 static PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DMLabel subpointMap, DM subdm) 2017 { 2018 IS subvertexIS = NULL; 2019 const PetscInt *subvertices; 2020 PetscInt *pStart, *pEnd, *pMax; 2021 PetscInt dim, d, numSubVerticesInitial = 0, v; 2022 PetscErrorCode ierr; 2023 2024 PetscFunctionBegin; 2025 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2026 ierr = PetscMalloc3(dim+1,&pStart,dim+1,&pEnd,dim+1,&pMax);CHKERRQ(ierr); 2027 ierr = DMPlexGetHybridBounds(dm, &pMax[dim], dim>1 ? &pMax[dim-1] : NULL, dim > 2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 2028 for (d = 0; d <= dim; ++d) { 2029 ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr); 2030 if (pMax[d] >= 0) pEnd[d] = PetscMin(pEnd[d], pMax[d]); 2031 } 2032 /* Loop over initial vertices and mark all faces in the collective star() */ 2033 if (vertexLabel) { 2034 ierr = DMLabelGetStratumIS(vertexLabel, value, &subvertexIS);CHKERRQ(ierr); 2035 if (subvertexIS) { 2036 ierr = ISGetSize(subvertexIS, &numSubVerticesInitial);CHKERRQ(ierr); 2037 ierr = ISGetIndices(subvertexIS, &subvertices);CHKERRQ(ierr); 2038 } 2039 } 2040 for (v = 0; v < numSubVerticesInitial; ++v) { 2041 const PetscInt vertex = subvertices[v]; 2042 PetscInt *star = NULL; 2043 PetscInt starSize, s, numFaces = 0, f; 2044 2045 ierr = DMPlexGetTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2046 for (s = 0; s < starSize*2; s += 2) { 2047 const PetscInt point = star[s]; 2048 if ((point >= pStart[dim-1]) && (point < pEnd[dim-1])) star[numFaces++] = point; 2049 } 2050 for (f = 0; f < numFaces; ++f) { 2051 const PetscInt face = star[f]; 2052 PetscInt *closure = NULL; 2053 PetscInt closureSize, c; 2054 PetscInt faceLoc; 2055 2056 ierr = DMLabelGetValue(subpointMap, face, &faceLoc);CHKERRQ(ierr); 2057 if (faceLoc == dim-1) continue; 2058 if (faceLoc >= 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Face %d has dimension %d in the surface label", face, faceLoc); 2059 ierr = DMPlexGetTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2060 for (c = 0; c < closureSize*2; c += 2) { 2061 const PetscInt point = closure[c]; 2062 PetscInt vertexLoc; 2063 2064 if ((point >= pStart[0]) && (point < pEnd[0])) { 2065 ierr = DMLabelGetValue(vertexLabel, point, &vertexLoc);CHKERRQ(ierr); 2066 if (vertexLoc != value) break; 2067 } 2068 } 2069 if (c == closureSize*2) { 2070 const PetscInt *support; 2071 PetscInt supportSize, s; 2072 2073 for (c = 0; c < closureSize*2; c += 2) { 2074 const PetscInt point = closure[c]; 2075 2076 for (d = 0; d < dim; ++d) { 2077 if ((point >= pStart[d]) && (point < pEnd[d])) { 2078 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2079 break; 2080 } 2081 } 2082 } 2083 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2084 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2085 for (s = 0; s < supportSize; ++s) { 2086 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2087 } 2088 } 2089 ierr = DMPlexRestoreTransitiveClosure(dm, face, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2090 } 2091 ierr = DMPlexRestoreTransitiveClosure(dm, vertex, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 2092 } 2093 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subvertices);CHKERRQ(ierr);} 2094 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2095 ierr = PetscFree3(pStart,pEnd,pMax);CHKERRQ(ierr); 2096 PetscFunctionReturn(0); 2097 } 2098 2099 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char labelname[], PetscInt value, DMLabel subpointMap, PetscInt *numFaces, PetscInt *nFV, PetscInt *subCells[], DM subdm) 2100 { 2101 DMLabel label = NULL; 2102 const PetscInt *cone; 2103 PetscInt dim, cMax, cEnd, c, subc = 0, p, coneSize = -1; 2104 PetscErrorCode ierr; 2105 2106 PetscFunctionBegin; 2107 *numFaces = 0; 2108 *nFV = 0; 2109 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 2110 *subCells = NULL; 2111 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2112 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2113 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2114 if (cMax < 0) PetscFunctionReturn(0); 2115 if (label) { 2116 for (c = cMax; c < cEnd; ++c) { 2117 PetscInt val; 2118 2119 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2120 if (val == value) { 2121 ++(*numFaces); 2122 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2123 } 2124 } 2125 } else { 2126 *numFaces = cEnd - cMax; 2127 ierr = DMPlexGetConeSize(dm, cMax, &coneSize);CHKERRQ(ierr); 2128 } 2129 ierr = PetscMalloc1(*numFaces *2, subCells);CHKERRQ(ierr); 2130 if (!(*numFaces)) PetscFunctionReturn(0); 2131 *nFV = hasLagrange ? coneSize/3 : coneSize/2; 2132 for (c = cMax; c < cEnd; ++c) { 2133 const PetscInt *cells; 2134 PetscInt numCells; 2135 2136 if (label) { 2137 PetscInt val; 2138 2139 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2140 if (val != value) continue; 2141 } 2142 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2143 for (p = 0; p < *nFV; ++p) { 2144 ierr = DMLabelSetValue(subpointMap, cone[p], 0);CHKERRQ(ierr); 2145 } 2146 /* Negative face */ 2147 ierr = DMPlexGetJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2148 /* Not true in parallel 2149 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 2150 for (p = 0; p < numCells; ++p) { 2151 ierr = DMLabelSetValue(subpointMap, cells[p], 2);CHKERRQ(ierr); 2152 (*subCells)[subc++] = cells[p]; 2153 } 2154 ierr = DMPlexRestoreJoin(dm, *nFV, cone, &numCells, &cells);CHKERRQ(ierr); 2155 /* Positive face is not included */ 2156 } 2157 PetscFunctionReturn(0); 2158 } 2159 2160 static PetscErrorCode DMPlexMarkCohesiveSubmesh_Interpolated(DM dm, DMLabel label, PetscInt value, DMLabel subpointMap, DM subdm) 2161 { 2162 PetscInt *pStart, *pEnd; 2163 PetscInt dim, cMax, cEnd, c, d; 2164 PetscErrorCode ierr; 2165 2166 PetscFunctionBegin; 2167 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2168 ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 2169 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 2170 if (cMax < 0) PetscFunctionReturn(0); 2171 ierr = PetscMalloc2(dim+1,&pStart,dim+1,&pEnd);CHKERRQ(ierr); 2172 for (d = 0; d <= dim; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);CHKERRQ(ierr);} 2173 for (c = cMax; c < cEnd; ++c) { 2174 const PetscInt *cone; 2175 PetscInt *closure = NULL; 2176 PetscInt fconeSize, coneSize, closureSize, cl, val; 2177 2178 if (label) { 2179 ierr = DMLabelGetValue(label, c, &val);CHKERRQ(ierr); 2180 if (val != value) continue; 2181 } 2182 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2183 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2184 ierr = DMPlexGetConeSize(dm, cone[0], &fconeSize);CHKERRQ(ierr); 2185 if (coneSize != (fconeSize ? fconeSize : 1) + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); 2186 /* Negative face */ 2187 ierr = DMPlexGetTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2188 for (cl = 0; cl < closureSize*2; cl += 2) { 2189 const PetscInt point = closure[cl]; 2190 2191 for (d = 0; d <= dim; ++d) { 2192 if ((point >= pStart[d]) && (point < pEnd[d])) { 2193 ierr = DMLabelSetValue(subpointMap, point, d);CHKERRQ(ierr); 2194 break; 2195 } 2196 } 2197 } 2198 ierr = DMPlexRestoreTransitiveClosure(dm, cone[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2199 /* Cells -- positive face is not included */ 2200 for (cl = 0; cl < 1; ++cl) { 2201 const PetscInt *support; 2202 PetscInt supportSize, s; 2203 2204 ierr = DMPlexGetSupportSize(dm, cone[cl], &supportSize);CHKERRQ(ierr); 2205 /* if (supportSize != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive faces should separate two cells"); */ 2206 ierr = DMPlexGetSupport(dm, cone[cl], &support);CHKERRQ(ierr); 2207 for (s = 0; s < supportSize; ++s) { 2208 ierr = DMLabelSetValue(subpointMap, support[s], dim);CHKERRQ(ierr); 2209 } 2210 } 2211 } 2212 ierr = PetscFree2(pStart, pEnd);CHKERRQ(ierr); 2213 PetscFunctionReturn(0); 2214 } 2215 2216 static PetscErrorCode DMPlexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2217 { 2218 MPI_Comm comm; 2219 PetscBool posOrient = PETSC_FALSE; 2220 const PetscInt debug = 0; 2221 PetscInt cellDim, faceSize, f; 2222 PetscErrorCode ierr; 2223 2224 PetscFunctionBegin; 2225 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2226 ierr = DMGetDimension(dm, &cellDim);CHKERRQ(ierr); 2227 if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);CHKERRQ(ierr);} 2228 2229 if (cellDim == 1 && numCorners == 2) { 2230 /* Triangle */ 2231 faceSize = numCorners-1; 2232 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2233 } else if (cellDim == 2 && numCorners == 3) { 2234 /* Triangle */ 2235 faceSize = numCorners-1; 2236 posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2237 } else if (cellDim == 3 && numCorners == 4) { 2238 /* Tetrahedron */ 2239 faceSize = numCorners-1; 2240 posOrient = (oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE; 2241 } else if (cellDim == 1 && numCorners == 3) { 2242 /* Quadratic line */ 2243 faceSize = 1; 2244 posOrient = PETSC_TRUE; 2245 } else if (cellDim == 2 && numCorners == 4) { 2246 /* Quads */ 2247 faceSize = 2; 2248 if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) { 2249 posOrient = PETSC_TRUE; 2250 } else if ((indices[0] == 3) && (indices[1] == 0)) { 2251 posOrient = PETSC_TRUE; 2252 } else { 2253 if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) { 2254 posOrient = PETSC_FALSE; 2255 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge"); 2256 } 2257 } else if (cellDim == 2 && numCorners == 6) { 2258 /* Quadratic triangle (I hate this) */ 2259 /* Edges are determined by the first 2 vertices (corners of edges) */ 2260 const PetscInt faceSizeTri = 3; 2261 PetscInt sortedIndices[3], i, iFace; 2262 PetscBool found = PETSC_FALSE; 2263 PetscInt faceVerticesTriSorted[9] = { 2264 0, 3, 4, /* bottom */ 2265 1, 4, 5, /* right */ 2266 2, 3, 5, /* left */ 2267 }; 2268 PetscInt faceVerticesTri[9] = { 2269 0, 3, 4, /* bottom */ 2270 1, 4, 5, /* right */ 2271 2, 5, 3, /* left */ 2272 }; 2273 2274 for (i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i]; 2275 ierr = PetscSortInt(faceSizeTri, sortedIndices);CHKERRQ(ierr); 2276 for (iFace = 0; iFace < 3; ++iFace) { 2277 const PetscInt ii = iFace*faceSizeTri; 2278 PetscInt fVertex, cVertex; 2279 2280 if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) && 2281 (sortedIndices[1] == faceVerticesTriSorted[ii+1])) { 2282 for (fVertex = 0; fVertex < faceSizeTri; ++fVertex) { 2283 for (cVertex = 0; cVertex < faceSizeTri; ++cVertex) { 2284 if (indices[cVertex] == faceVerticesTri[ii+fVertex]) { 2285 faceVertices[fVertex] = origVertices[cVertex]; 2286 break; 2287 } 2288 } 2289 } 2290 found = PETSC_TRUE; 2291 break; 2292 } 2293 } 2294 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface"); 2295 if (posOriented) *posOriented = PETSC_TRUE; 2296 PetscFunctionReturn(0); 2297 } else if (cellDim == 2 && numCorners == 9) { 2298 /* Quadratic quad (I hate this) */ 2299 /* Edges are determined by the first 2 vertices (corners of edges) */ 2300 const PetscInt faceSizeQuad = 3; 2301 PetscInt sortedIndices[3], i, iFace; 2302 PetscBool found = PETSC_FALSE; 2303 PetscInt faceVerticesQuadSorted[12] = { 2304 0, 1, 4, /* bottom */ 2305 1, 2, 5, /* right */ 2306 2, 3, 6, /* top */ 2307 0, 3, 7, /* left */ 2308 }; 2309 PetscInt faceVerticesQuad[12] = { 2310 0, 1, 4, /* bottom */ 2311 1, 2, 5, /* right */ 2312 2, 3, 6, /* top */ 2313 3, 0, 7, /* left */ 2314 }; 2315 2316 for (i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i]; 2317 ierr = PetscSortInt(faceSizeQuad, sortedIndices);CHKERRQ(ierr); 2318 for (iFace = 0; iFace < 4; ++iFace) { 2319 const PetscInt ii = iFace*faceSizeQuad; 2320 PetscInt fVertex, cVertex; 2321 2322 if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) && 2323 (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) { 2324 for (fVertex = 0; fVertex < faceSizeQuad; ++fVertex) { 2325 for (cVertex = 0; cVertex < faceSizeQuad; ++cVertex) { 2326 if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) { 2327 faceVertices[fVertex] = origVertices[cVertex]; 2328 break; 2329 } 2330 } 2331 } 2332 found = PETSC_TRUE; 2333 break; 2334 } 2335 } 2336 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface"); 2337 if (posOriented) *posOriented = PETSC_TRUE; 2338 PetscFunctionReturn(0); 2339 } else if (cellDim == 3 && numCorners == 8) { 2340 /* Hexes 2341 A hex is two oriented quads with the normal of the first 2342 pointing up at the second. 2343 2344 7---6 2345 /| /| 2346 4---5 | 2347 | 1-|-2 2348 |/ |/ 2349 0---3 2350 2351 Faces are determined by the first 4 vertices (corners of faces) */ 2352 const PetscInt faceSizeHex = 4; 2353 PetscInt sortedIndices[4], i, iFace; 2354 PetscBool found = PETSC_FALSE; 2355 PetscInt faceVerticesHexSorted[24] = { 2356 0, 1, 2, 3, /* bottom */ 2357 4, 5, 6, 7, /* top */ 2358 0, 3, 4, 5, /* front */ 2359 2, 3, 5, 6, /* right */ 2360 1, 2, 6, 7, /* back */ 2361 0, 1, 4, 7, /* left */ 2362 }; 2363 PetscInt faceVerticesHex[24] = { 2364 1, 2, 3, 0, /* bottom */ 2365 4, 5, 6, 7, /* top */ 2366 0, 3, 5, 4, /* front */ 2367 3, 2, 6, 5, /* right */ 2368 2, 1, 7, 6, /* back */ 2369 1, 0, 4, 7, /* left */ 2370 }; 2371 2372 for (i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i]; 2373 ierr = PetscSortInt(faceSizeHex, sortedIndices);CHKERRQ(ierr); 2374 for (iFace = 0; iFace < 6; ++iFace) { 2375 const PetscInt ii = iFace*faceSizeHex; 2376 PetscInt fVertex, cVertex; 2377 2378 if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) && 2379 (sortedIndices[1] == faceVerticesHexSorted[ii+1]) && 2380 (sortedIndices[2] == faceVerticesHexSorted[ii+2]) && 2381 (sortedIndices[3] == faceVerticesHexSorted[ii+3])) { 2382 for (fVertex = 0; fVertex < faceSizeHex; ++fVertex) { 2383 for (cVertex = 0; cVertex < faceSizeHex; ++cVertex) { 2384 if (indices[cVertex] == faceVerticesHex[ii+fVertex]) { 2385 faceVertices[fVertex] = origVertices[cVertex]; 2386 break; 2387 } 2388 } 2389 } 2390 found = PETSC_TRUE; 2391 break; 2392 } 2393 } 2394 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2395 if (posOriented) *posOriented = PETSC_TRUE; 2396 PetscFunctionReturn(0); 2397 } else if (cellDim == 3 && numCorners == 10) { 2398 /* Quadratic tet */ 2399 /* Faces are determined by the first 3 vertices (corners of faces) */ 2400 const PetscInt faceSizeTet = 6; 2401 PetscInt sortedIndices[6], i, iFace; 2402 PetscBool found = PETSC_FALSE; 2403 PetscInt faceVerticesTetSorted[24] = { 2404 0, 1, 2, 6, 7, 8, /* bottom */ 2405 0, 3, 4, 6, 7, 9, /* front */ 2406 1, 4, 5, 7, 8, 9, /* right */ 2407 2, 3, 5, 6, 8, 9, /* left */ 2408 }; 2409 PetscInt faceVerticesTet[24] = { 2410 0, 1, 2, 6, 7, 8, /* bottom */ 2411 0, 4, 3, 6, 7, 9, /* front */ 2412 1, 5, 4, 7, 8, 9, /* right */ 2413 2, 3, 5, 8, 6, 9, /* left */ 2414 }; 2415 2416 for (i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i]; 2417 ierr = PetscSortInt(faceSizeTet, sortedIndices);CHKERRQ(ierr); 2418 for (iFace=0; iFace < 4; ++iFace) { 2419 const PetscInt ii = iFace*faceSizeTet; 2420 PetscInt fVertex, cVertex; 2421 2422 if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) && 2423 (sortedIndices[1] == faceVerticesTetSorted[ii+1]) && 2424 (sortedIndices[2] == faceVerticesTetSorted[ii+2]) && 2425 (sortedIndices[3] == faceVerticesTetSorted[ii+3])) { 2426 for (fVertex = 0; fVertex < faceSizeTet; ++fVertex) { 2427 for (cVertex = 0; cVertex < faceSizeTet; ++cVertex) { 2428 if (indices[cVertex] == faceVerticesTet[ii+fVertex]) { 2429 faceVertices[fVertex] = origVertices[cVertex]; 2430 break; 2431 } 2432 } 2433 } 2434 found = PETSC_TRUE; 2435 break; 2436 } 2437 } 2438 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface"); 2439 if (posOriented) *posOriented = PETSC_TRUE; 2440 PetscFunctionReturn(0); 2441 } else if (cellDim == 3 && numCorners == 27) { 2442 /* Quadratic hexes (I hate this) 2443 A hex is two oriented quads with the normal of the first 2444 pointing up at the second. 2445 2446 7---6 2447 /| /| 2448 4---5 | 2449 | 3-|-2 2450 |/ |/ 2451 0---1 2452 2453 Faces are determined by the first 4 vertices (corners of faces) */ 2454 const PetscInt faceSizeQuadHex = 9; 2455 PetscInt sortedIndices[9], i, iFace; 2456 PetscBool found = PETSC_FALSE; 2457 PetscInt faceVerticesQuadHexSorted[54] = { 2458 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */ 2459 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2460 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */ 2461 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */ 2462 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */ 2463 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */ 2464 }; 2465 PetscInt faceVerticesQuadHex[54] = { 2466 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */ 2467 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */ 2468 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */ 2469 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */ 2470 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */ 2471 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */ 2472 }; 2473 2474 for (i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i]; 2475 ierr = PetscSortInt(faceSizeQuadHex, sortedIndices);CHKERRQ(ierr); 2476 for (iFace = 0; iFace < 6; ++iFace) { 2477 const PetscInt ii = iFace*faceSizeQuadHex; 2478 PetscInt fVertex, cVertex; 2479 2480 if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) && 2481 (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) && 2482 (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) && 2483 (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) { 2484 for (fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) { 2485 for (cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) { 2486 if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) { 2487 faceVertices[fVertex] = origVertices[cVertex]; 2488 break; 2489 } 2490 } 2491 } 2492 found = PETSC_TRUE; 2493 break; 2494 } 2495 } 2496 if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface"); 2497 if (posOriented) *posOriented = PETSC_TRUE; 2498 PetscFunctionReturn(0); 2499 } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation()."); 2500 if (!posOrient) { 2501 if (debug) {ierr = PetscPrintf(comm, " Reversing initial face orientation\n");CHKERRQ(ierr);} 2502 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[faceSize-1 - f]; 2503 } else { 2504 if (debug) {ierr = PetscPrintf(comm, " Keeping initial face orientation\n");CHKERRQ(ierr);} 2505 for (f = 0; f < faceSize; ++f) faceVertices[f] = origVertices[f]; 2506 } 2507 if (posOriented) *posOriented = posOrient; 2508 PetscFunctionReturn(0); 2509 } 2510 2511 /*@ 2512 DMPlexGetOrientedFace - Given a cell and a face, as a set of vertices, return the oriented face, as a set of vertices, 2513 in faceVertices. The orientation is such that the face normal points out of the cell 2514 2515 Not collective 2516 2517 Input Parameters: 2518 + dm - The original mesh 2519 . cell - The cell mesh point 2520 . faceSize - The number of vertices on the face 2521 . face - The face vertices 2522 . numCorners - The number of vertices on the cell 2523 . indices - Local numbering of face vertices in cell cone 2524 - origVertices - Original face vertices 2525 2526 Output Parameter: 2527 + faceVertices - The face vertices properly oriented 2528 - posOriented - PETSC_TRUE if the face was oriented with outward normal 2529 2530 Level: developer 2531 2532 .seealso: DMPlexGetCone() 2533 @*/ 2534 PetscErrorCode DMPlexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) 2535 { 2536 const PetscInt *cone = NULL; 2537 PetscInt coneSize, v, f, v2; 2538 PetscInt oppositeVertex = -1; 2539 PetscErrorCode ierr; 2540 2541 PetscFunctionBegin; 2542 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 2543 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 2544 for (v = 0, v2 = 0; v < coneSize; ++v) { 2545 PetscBool found = PETSC_FALSE; 2546 2547 for (f = 0; f < faceSize; ++f) { 2548 if (face[f] == cone[v]) { 2549 found = PETSC_TRUE; break; 2550 } 2551 } 2552 if (found) { 2553 indices[v2] = v; 2554 origVertices[v2] = cone[v]; 2555 ++v2; 2556 } else { 2557 oppositeVertex = v; 2558 } 2559 } 2560 ierr = DMPlexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);CHKERRQ(ierr); 2561 PetscFunctionReturn(0); 2562 } 2563 2564 /* 2565 DMPlexInsertFace_Internal - Puts a face into the mesh 2566 2567 Not collective 2568 2569 Input Parameters: 2570 + dm - The DMPlex 2571 . numFaceVertex - The number of vertices in the face 2572 . faceVertices - The vertices in the face for dm 2573 . subfaceVertices - The vertices in the face for subdm 2574 . numCorners - The number of vertices in the cell 2575 . cell - A cell in dm containing the face 2576 . subcell - A cell in subdm containing the face 2577 . firstFace - First face in the mesh 2578 - newFacePoint - Next face in the mesh 2579 2580 Output Parameters: 2581 . newFacePoint - Contains next face point number on input, updated on output 2582 2583 Level: developer 2584 */ 2585 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) 2586 { 2587 MPI_Comm comm; 2588 DM_Plex *submesh = (DM_Plex*) subdm->data; 2589 const PetscInt *faces; 2590 PetscInt numFaces, coneSize; 2591 PetscErrorCode ierr; 2592 2593 PetscFunctionBegin; 2594 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2595 ierr = DMPlexGetConeSize(subdm, subcell, &coneSize);CHKERRQ(ierr); 2596 if (coneSize != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size of cell %d is %d != 1", cell, coneSize); 2597 #if 0 2598 /* Cannot use this because support() has not been constructed yet */ 2599 ierr = DMPlexGetJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2600 #else 2601 { 2602 PetscInt f; 2603 2604 numFaces = 0; 2605 ierr = DMGetWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2606 for (f = firstFace; f < *newFacePoint; ++f) { 2607 PetscInt dof, off, d; 2608 2609 ierr = PetscSectionGetDof(submesh->coneSection, f, &dof);CHKERRQ(ierr); 2610 ierr = PetscSectionGetOffset(submesh->coneSection, f, &off);CHKERRQ(ierr); 2611 /* Yes, I know this is quadratic, but I expect the sizes to be <5 */ 2612 for (d = 0; d < dof; ++d) { 2613 const PetscInt p = submesh->cones[off+d]; 2614 PetscInt v; 2615 2616 for (v = 0; v < numFaceVertices; ++v) { 2617 if (subfaceVertices[v] == p) break; 2618 } 2619 if (v == numFaceVertices) break; 2620 } 2621 if (d == dof) { 2622 numFaces = 1; 2623 ((PetscInt*) faces)[0] = f; 2624 } 2625 } 2626 } 2627 #endif 2628 if (numFaces > 1) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Vertex set had %d faces, not one", numFaces); 2629 else if (numFaces == 1) { 2630 /* Add the other cell neighbor for this face */ 2631 ierr = DMPlexSetCone(subdm, subcell, faces);CHKERRQ(ierr); 2632 } else { 2633 PetscInt *indices, *origVertices, *orientedVertices, *orientedSubVertices, v, ov; 2634 PetscBool posOriented; 2635 2636 ierr = DMGetWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2637 origVertices = &orientedVertices[numFaceVertices]; 2638 indices = &orientedVertices[numFaceVertices*2]; 2639 orientedSubVertices = &orientedVertices[numFaceVertices*3]; 2640 ierr = DMPlexGetOrientedFace(dm, cell, numFaceVertices, faceVertices, numCorners, indices, origVertices, orientedVertices, &posOriented);CHKERRQ(ierr); 2641 /* TODO: I know that routine should return a permutation, not the indices */ 2642 for (v = 0; v < numFaceVertices; ++v) { 2643 const PetscInt vertex = faceVertices[v], subvertex = subfaceVertices[v]; 2644 for (ov = 0; ov < numFaceVertices; ++ov) { 2645 if (orientedVertices[ov] == vertex) { 2646 orientedSubVertices[ov] = subvertex; 2647 break; 2648 } 2649 } 2650 if (ov == numFaceVertices) SETERRQ1(comm, PETSC_ERR_PLIB, "Could not find face vertex %d in orientated set", vertex); 2651 } 2652 ierr = DMPlexSetCone(subdm, *newFacePoint, orientedSubVertices);CHKERRQ(ierr); 2653 ierr = DMPlexSetCone(subdm, subcell, newFacePoint);CHKERRQ(ierr); 2654 ierr = DMRestoreWorkArray(subdm, 4*numFaceVertices * sizeof(PetscInt), PETSC_INT, &orientedVertices);CHKERRQ(ierr); 2655 ++(*newFacePoint); 2656 } 2657 #if 0 2658 ierr = DMPlexRestoreJoin(subdm, numFaceVertices, subfaceVertices, &numFaces, &faces);CHKERRQ(ierr); 2659 #else 2660 ierr = DMRestoreWorkArray(subdm, 1, PETSC_INT, (void **) &faces);CHKERRQ(ierr); 2661 #endif 2662 PetscFunctionReturn(0); 2663 } 2664 2665 static PetscErrorCode DMPlexCreateSubmesh_Uninterpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 2666 { 2667 MPI_Comm comm; 2668 DMLabel subpointMap; 2669 IS subvertexIS, subcellIS; 2670 const PetscInt *subVertices, *subCells; 2671 PetscInt numSubVertices, firstSubVertex, numSubCells; 2672 PetscInt *subface, maxConeSize, numSubFaces = 0, firstSubFace, newFacePoint, nFV = 0; 2673 PetscInt vStart, vEnd, c, f; 2674 PetscErrorCode ierr; 2675 2676 PetscFunctionBegin; 2677 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2678 /* Create subpointMap which marks the submesh */ 2679 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2680 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2681 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2682 if (vertexLabel) {ierr = DMPlexMarkSubmesh_Uninterpolated(dm, vertexLabel, value, subpointMap, &numSubFaces, &nFV, subdm);CHKERRQ(ierr);} 2683 /* Setup chart */ 2684 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 2685 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 2686 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 2687 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 2688 /* Set cone sizes */ 2689 firstSubVertex = numSubCells; 2690 firstSubFace = numSubCells+numSubVertices; 2691 newFacePoint = firstSubFace; 2692 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 2693 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2694 ierr = DMLabelGetStratumIS(subpointMap, 2, &subcellIS);CHKERRQ(ierr); 2695 if (subcellIS) {ierr = ISGetIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2696 for (c = 0; c < numSubCells; ++c) { 2697 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 2698 } 2699 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 2700 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 2701 } 2702 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2703 /* Create face cones */ 2704 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2705 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2706 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2707 for (c = 0; c < numSubCells; ++c) { 2708 const PetscInt cell = subCells[c]; 2709 const PetscInt subcell = c; 2710 PetscInt *closure = NULL; 2711 PetscInt closureSize, cl, numCorners = 0, faceSize = 0; 2712 2713 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2714 for (cl = 0; cl < closureSize*2; cl += 2) { 2715 const PetscInt point = closure[cl]; 2716 PetscInt subVertex; 2717 2718 if ((point >= vStart) && (point < vEnd)) { 2719 ++numCorners; 2720 ierr = PetscFindInt(point, numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 2721 if (subVertex >= 0) { 2722 closure[faceSize] = point; 2723 subface[faceSize] = firstSubVertex+subVertex; 2724 ++faceSize; 2725 } 2726 } 2727 } 2728 if (faceSize > nFV) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid submesh: Too many vertices %d of an element on the surface", faceSize); 2729 if (faceSize == nFV) { 2730 ierr = DMPlexInsertFace_Internal(dm, subdm, faceSize, closure, subface, numCorners, cell, subcell, firstSubFace, &newFacePoint);CHKERRQ(ierr); 2731 } 2732 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2733 } 2734 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 2735 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2736 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2737 /* Build coordinates */ 2738 { 2739 PetscSection coordSection, subCoordSection; 2740 Vec coordinates, subCoordinates; 2741 PetscScalar *coords, *subCoords; 2742 PetscInt numComp, coordSize, v; 2743 const char *name; 2744 2745 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2746 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2747 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2748 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2749 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2750 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2751 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 2752 for (v = 0; v < numSubVertices; ++v) { 2753 const PetscInt vertex = subVertices[v]; 2754 const PetscInt subvertex = firstSubVertex+v; 2755 PetscInt dof; 2756 2757 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2758 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2759 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2760 } 2761 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2762 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2763 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 2764 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 2765 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 2766 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2767 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2768 if (coordSize) { 2769 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2770 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2771 for (v = 0; v < numSubVertices; ++v) { 2772 const PetscInt vertex = subVertices[v]; 2773 const PetscInt subvertex = firstSubVertex+v; 2774 PetscInt dof, off, sdof, soff, d; 2775 2776 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2777 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2778 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2779 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2780 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2781 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2782 } 2783 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2784 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2785 } 2786 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2787 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2788 } 2789 /* Cleanup */ 2790 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 2791 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 2792 if (subcellIS) {ierr = ISRestoreIndices(subcellIS, &subCells);CHKERRQ(ierr);} 2793 ierr = ISDestroy(&subcellIS);CHKERRQ(ierr); 2794 PetscFunctionReturn(0); 2795 } 2796 2797 PETSC_STATIC_INLINE PetscInt DMPlexFilterPoint_Internal(PetscInt point, PetscInt firstSubPoint, PetscInt numSubPoints, const PetscInt subPoints[]) 2798 { 2799 PetscInt subPoint; 2800 PetscErrorCode ierr; 2801 2802 ierr = PetscFindInt(point, numSubPoints, subPoints, &subPoint); if (ierr < 0) return ierr; 2803 return subPoint < 0 ? subPoint : firstSubPoint+subPoint; 2804 } 2805 2806 static PetscErrorCode DMPlexCreateSubmeshGeneric_Interpolated(DM dm, DMLabel label, PetscInt value, PetscBool isCohesive, PetscInt cellHeight, DM subdm) 2807 { 2808 MPI_Comm comm; 2809 DMLabel subpointMap; 2810 IS *subpointIS; 2811 const PetscInt **subpoints; 2812 PetscInt *numSubPoints, *firstSubPoint, *coneNew, *orntNew; 2813 PetscInt totSubPoints = 0, maxConeSize, dim, p, d, v; 2814 PetscErrorCode ierr; 2815 2816 PetscFunctionBegin; 2817 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2818 /* Create subpointMap which marks the submesh */ 2819 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 2820 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 2821 if (cellHeight) { 2822 if (isCohesive) {ierr = DMPlexMarkCohesiveSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2823 else {ierr = DMPlexMarkSubmesh_Interpolated(dm, label, value, subpointMap, subdm);CHKERRQ(ierr);} 2824 } else { 2825 DMLabel depth; 2826 IS pointIS; 2827 const PetscInt *points; 2828 PetscInt numPoints; 2829 2830 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 2831 ierr = DMLabelGetStratumSize(label, value, &numPoints);CHKERRQ(ierr); 2832 ierr = DMLabelGetStratumIS(label, value, &pointIS);CHKERRQ(ierr); 2833 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2834 for (p = 0; p < numPoints; ++p) { 2835 PetscInt *closure = NULL; 2836 PetscInt closureSize, c, pdim; 2837 2838 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2839 for (c = 0; c < closureSize*2; c += 2) { 2840 ierr = DMLabelGetValue(depth, closure[c], &pdim);CHKERRQ(ierr); 2841 ierr = DMLabelSetValue(subpointMap, closure[c], pdim);CHKERRQ(ierr); 2842 } 2843 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2844 } 2845 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2846 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2847 } 2848 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 2849 /* Setup chart */ 2850 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2851 ierr = PetscMalloc4(dim+1,&numSubPoints,dim+1,&firstSubPoint,dim+1,&subpointIS,dim+1,&subpoints);CHKERRQ(ierr); 2852 for (d = 0; d <= dim; ++d) { 2853 ierr = DMLabelGetStratumSize(subpointMap, d, &numSubPoints[d]);CHKERRQ(ierr); 2854 totSubPoints += numSubPoints[d]; 2855 } 2856 ierr = DMPlexSetChart(subdm, 0, totSubPoints);CHKERRQ(ierr); 2857 ierr = DMPlexSetVTKCellHeight(subdm, cellHeight);CHKERRQ(ierr); 2858 /* Set cone sizes */ 2859 firstSubPoint[dim] = 0; 2860 firstSubPoint[0] = firstSubPoint[dim] + numSubPoints[dim]; 2861 if (dim > 1) {firstSubPoint[dim-1] = firstSubPoint[0] + numSubPoints[0];} 2862 if (dim > 2) {firstSubPoint[dim-2] = firstSubPoint[dim-1] + numSubPoints[dim-1];} 2863 for (d = 0; d <= dim; ++d) { 2864 ierr = DMLabelGetStratumIS(subpointMap, d, &subpointIS[d]);CHKERRQ(ierr); 2865 if (subpointIS[d]) {ierr = ISGetIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 2866 } 2867 for (d = 0; d <= dim; ++d) { 2868 for (p = 0; p < numSubPoints[d]; ++p) { 2869 const PetscInt point = subpoints[d][p]; 2870 const PetscInt subpoint = firstSubPoint[d] + p; 2871 const PetscInt *cone; 2872 PetscInt coneSize, coneSizeNew, c, val; 2873 2874 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2875 ierr = DMPlexSetConeSize(subdm, subpoint, coneSize);CHKERRQ(ierr); 2876 if (cellHeight && (d == dim)) { 2877 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2878 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2879 ierr = DMLabelGetValue(subpointMap, cone[c], &val);CHKERRQ(ierr); 2880 if (val >= 0) coneSizeNew++; 2881 } 2882 ierr = DMPlexSetConeSize(subdm, subpoint, coneSizeNew);CHKERRQ(ierr); 2883 } 2884 } 2885 } 2886 ierr = DMSetUp(subdm);CHKERRQ(ierr); 2887 /* Set cones */ 2888 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2889 ierr = PetscMalloc2(maxConeSize,&coneNew,maxConeSize,&orntNew);CHKERRQ(ierr); 2890 for (d = 0; d <= dim; ++d) { 2891 for (p = 0; p < numSubPoints[d]; ++p) { 2892 const PetscInt point = subpoints[d][p]; 2893 const PetscInt subpoint = firstSubPoint[d] + p; 2894 const PetscInt *cone, *ornt; 2895 PetscInt coneSize, subconeSize, coneSizeNew, c, subc; 2896 2897 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 2898 ierr = DMPlexGetConeSize(subdm, subpoint, &subconeSize);CHKERRQ(ierr); 2899 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 2900 ierr = DMPlexGetConeOrientation(dm, point, &ornt);CHKERRQ(ierr); 2901 for (c = 0, coneSizeNew = 0; c < coneSize; ++c) { 2902 ierr = PetscFindInt(cone[c], numSubPoints[d-1], subpoints[d-1], &subc);CHKERRQ(ierr); 2903 if (subc >= 0) { 2904 coneNew[coneSizeNew] = firstSubPoint[d-1] + subc; 2905 orntNew[coneSizeNew] = ornt[c]; 2906 ++coneSizeNew; 2907 } 2908 } 2909 if (coneSizeNew != subconeSize) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of cone points located %d does not match subcone size %d", coneSizeNew, subconeSize); 2910 ierr = DMPlexSetCone(subdm, subpoint, coneNew);CHKERRQ(ierr); 2911 ierr = DMPlexSetConeOrientation(subdm, subpoint, orntNew);CHKERRQ(ierr); 2912 } 2913 } 2914 ierr = PetscFree2(coneNew,orntNew);CHKERRQ(ierr); 2915 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 2916 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 2917 /* Build coordinates */ 2918 { 2919 PetscSection coordSection, subCoordSection; 2920 Vec coordinates, subCoordinates; 2921 PetscScalar *coords, *subCoords; 2922 PetscInt numComp, coordSize; 2923 const char *name; 2924 2925 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2926 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2927 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 2928 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 2929 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 2930 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 2931 ierr = PetscSectionSetChart(subCoordSection, firstSubPoint[0], firstSubPoint[0]+numSubPoints[0]);CHKERRQ(ierr); 2932 for (v = 0; v < numSubPoints[0]; ++v) { 2933 const PetscInt vertex = subpoints[0][v]; 2934 const PetscInt subvertex = firstSubPoint[0]+v; 2935 PetscInt dof; 2936 2937 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2938 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 2939 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 2940 } 2941 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 2942 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 2943 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 2944 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 2945 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 2946 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2947 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 2948 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2949 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2950 for (v = 0; v < numSubPoints[0]; ++v) { 2951 const PetscInt vertex = subpoints[0][v]; 2952 const PetscInt subvertex = firstSubPoint[0]+v; 2953 PetscInt dof, off, sdof, soff, d; 2954 2955 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 2956 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 2957 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 2958 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 2959 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 2960 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 2961 } 2962 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2963 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 2964 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 2965 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 2966 } 2967 /* Build SF: We need this complexity because subpoints might not be selected on the owning process */ 2968 { 2969 PetscSF sfPoint, sfPointSub; 2970 IS subpIS; 2971 const PetscSFNode *remotePoints; 2972 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 2973 const PetscInt *localPoints, *subpoints; 2974 PetscInt *slocalPoints; 2975 PetscInt numRoots, numLeaves, numSubpoints = 0, numSubroots, numSubleaves = 0, l, sl, ll, pStart, pEnd, p; 2976 PetscMPIInt rank; 2977 2978 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2979 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2980 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 2981 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2982 ierr = DMPlexGetChart(subdm, NULL, &numSubroots);CHKERRQ(ierr); 2983 ierr = DMPlexCreateSubpointIS(subdm, &subpIS);CHKERRQ(ierr); 2984 if (subpIS) { 2985 ierr = ISGetIndices(subpIS, &subpoints);CHKERRQ(ierr); 2986 ierr = ISGetLocalSize(subpIS, &numSubpoints);CHKERRQ(ierr); 2987 } 2988 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 2989 if (numRoots >= 0) { 2990 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 2991 for (p = 0; p < pEnd-pStart; ++p) { 2992 newLocalPoints[p].rank = -2; 2993 newLocalPoints[p].index = -2; 2994 } 2995 /* Set subleaves */ 2996 for (l = 0; l < numLeaves; ++l) { 2997 const PetscInt point = localPoints[l]; 2998 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 2999 3000 if (subpoint < 0) continue; 3001 newLocalPoints[point-pStart].rank = rank; 3002 newLocalPoints[point-pStart].index = subpoint; 3003 ++numSubleaves; 3004 } 3005 /* Must put in owned subpoints */ 3006 for (p = pStart; p < pEnd; ++p) { 3007 const PetscInt subpoint = DMPlexFilterPoint_Internal(p, 0, numSubpoints, subpoints); 3008 3009 if (subpoint < 0) { 3010 newOwners[p-pStart].rank = -3; 3011 newOwners[p-pStart].index = -3; 3012 } else { 3013 newOwners[p-pStart].rank = rank; 3014 newOwners[p-pStart].index = subpoint; 3015 } 3016 } 3017 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3018 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3019 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3020 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3021 ierr = PetscMalloc1(numSubleaves, &slocalPoints);CHKERRQ(ierr); 3022 ierr = PetscMalloc1(numSubleaves, &sremotePoints);CHKERRQ(ierr); 3023 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3024 const PetscInt point = localPoints[l]; 3025 const PetscInt subpoint = DMPlexFilterPoint_Internal(point, 0, numSubpoints, subpoints); 3026 3027 if (subpoint < 0) continue; 3028 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3029 slocalPoints[sl] = subpoint; 3030 sremotePoints[sl].rank = newLocalPoints[point].rank; 3031 sremotePoints[sl].index = newLocalPoints[point].index; 3032 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3033 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3034 ++sl; 3035 } 3036 if (sl + ll != numSubleaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubleaves); 3037 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3038 ierr = PetscSFSetGraph(sfPointSub, numSubroots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3039 } 3040 if (subpIS) { 3041 ierr = ISRestoreIndices(subpIS, &subpoints);CHKERRQ(ierr); 3042 ierr = ISDestroy(&subpIS);CHKERRQ(ierr); 3043 } 3044 } 3045 /* Cleanup */ 3046 for (d = 0; d <= dim; ++d) { 3047 if (subpointIS[d]) {ierr = ISRestoreIndices(subpointIS[d], &subpoints[d]);CHKERRQ(ierr);} 3048 ierr = ISDestroy(&subpointIS[d]);CHKERRQ(ierr); 3049 } 3050 ierr = PetscFree4(numSubPoints,firstSubPoint,subpointIS,subpoints);CHKERRQ(ierr); 3051 PetscFunctionReturn(0); 3052 } 3053 3054 static PetscErrorCode DMPlexCreateSubmesh_Interpolated(DM dm, DMLabel vertexLabel, PetscInt value, DM subdm) 3055 { 3056 PetscErrorCode ierr; 3057 3058 PetscFunctionBegin; 3059 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, vertexLabel, value, PETSC_FALSE, 1, subdm);CHKERRQ(ierr); 3060 PetscFunctionReturn(0); 3061 } 3062 3063 /*@ 3064 DMPlexCreateSubmesh - Extract a hypersurface from the mesh using vertices defined by a label 3065 3066 Input Parameters: 3067 + dm - The original mesh 3068 . vertexLabel - The DMLabel marking vertices contained in the surface 3069 - value - The label value to use 3070 3071 Output Parameter: 3072 . subdm - The surface mesh 3073 3074 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3075 3076 Level: developer 3077 3078 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3079 @*/ 3080 PetscErrorCode DMPlexCreateSubmesh(DM dm, DMLabel vertexLabel, PetscInt value, DM *subdm) 3081 { 3082 PetscInt dim, depth; 3083 PetscErrorCode ierr; 3084 3085 PetscFunctionBegin; 3086 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3087 PetscValidPointer(subdm, 3); 3088 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3089 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3090 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3091 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3092 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3093 if (depth == dim) { 3094 ierr = DMPlexCreateSubmesh_Interpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3095 } else { 3096 ierr = DMPlexCreateSubmesh_Uninterpolated(dm, vertexLabel, value, *subdm);CHKERRQ(ierr); 3097 } 3098 PetscFunctionReturn(0); 3099 } 3100 3101 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Uninterpolated(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM subdm) 3102 { 3103 MPI_Comm comm; 3104 DMLabel subpointMap; 3105 IS subvertexIS; 3106 const PetscInt *subVertices; 3107 PetscInt numSubVertices, firstSubVertex, numSubCells, *subCells = NULL; 3108 PetscInt *subface, maxConeSize, numSubFaces, firstSubFace, newFacePoint, nFV; 3109 PetscInt cMax, c, f; 3110 PetscErrorCode ierr; 3111 3112 PetscFunctionBegin; 3113 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 3114 /* Create subpointMap which marks the submesh */ 3115 ierr = DMLabelCreate("subpoint_map", &subpointMap);CHKERRQ(ierr); 3116 ierr = DMPlexSetSubpointMap(subdm, subpointMap);CHKERRQ(ierr); 3117 ierr = DMLabelDestroy(&subpointMap);CHKERRQ(ierr); 3118 ierr = DMPlexMarkCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, subpointMap, &numSubFaces, &nFV, &subCells, subdm);CHKERRQ(ierr); 3119 /* Setup chart */ 3120 ierr = DMLabelGetStratumSize(subpointMap, 0, &numSubVertices);CHKERRQ(ierr); 3121 ierr = DMLabelGetStratumSize(subpointMap, 2, &numSubCells);CHKERRQ(ierr); 3122 ierr = DMPlexSetChart(subdm, 0, numSubCells+numSubFaces+numSubVertices);CHKERRQ(ierr); 3123 ierr = DMPlexSetVTKCellHeight(subdm, 1);CHKERRQ(ierr); 3124 /* Set cone sizes */ 3125 firstSubVertex = numSubCells; 3126 firstSubFace = numSubCells+numSubVertices; 3127 newFacePoint = firstSubFace; 3128 ierr = DMLabelGetStratumIS(subpointMap, 0, &subvertexIS);CHKERRQ(ierr); 3129 if (subvertexIS) {ierr = ISGetIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3130 for (c = 0; c < numSubCells; ++c) { 3131 ierr = DMPlexSetConeSize(subdm, c, 1);CHKERRQ(ierr); 3132 } 3133 for (f = firstSubFace; f < firstSubFace+numSubFaces; ++f) { 3134 ierr = DMPlexSetConeSize(subdm, f, nFV);CHKERRQ(ierr); 3135 } 3136 ierr = DMSetUp(subdm);CHKERRQ(ierr); 3137 /* Create face cones */ 3138 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 3139 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 3140 ierr = DMGetWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 3141 for (c = 0; c < numSubCells; ++c) { 3142 const PetscInt cell = subCells[c]; 3143 const PetscInt subcell = c; 3144 const PetscInt *cone, *cells; 3145 PetscInt numCells, subVertex, p, v; 3146 3147 if (cell < cMax) continue; 3148 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr); 3149 for (v = 0; v < nFV; ++v) { 3150 ierr = PetscFindInt(cone[v], numSubVertices, subVertices, &subVertex);CHKERRQ(ierr); 3151 subface[v] = firstSubVertex+subVertex; 3152 } 3153 ierr = DMPlexSetCone(subdm, newFacePoint, subface);CHKERRQ(ierr); 3154 ierr = DMPlexSetCone(subdm, subcell, &newFacePoint);CHKERRQ(ierr); 3155 ierr = DMPlexGetJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3156 /* Not true in parallel 3157 if (numCells != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cohesive cells should separate two cells"); */ 3158 for (p = 0; p < numCells; ++p) { 3159 PetscInt negsubcell; 3160 3161 if (cells[p] >= cMax) continue; 3162 /* I know this is a crap search */ 3163 for (negsubcell = 0; negsubcell < numSubCells; ++negsubcell) { 3164 if (subCells[negsubcell] == cells[p]) break; 3165 } 3166 if (negsubcell == numSubCells) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find negative face neighbor for cohesive cell %d", cell); 3167 ierr = DMPlexSetCone(subdm, negsubcell, &newFacePoint);CHKERRQ(ierr); 3168 } 3169 ierr = DMPlexRestoreJoin(dm, nFV, cone, &numCells, &cells);CHKERRQ(ierr); 3170 ++newFacePoint; 3171 } 3172 ierr = DMRestoreWorkArray(subdm, maxConeSize, PETSC_INT, (void**) &subface);CHKERRQ(ierr); 3173 ierr = DMPlexSymmetrize(subdm);CHKERRQ(ierr); 3174 ierr = DMPlexStratify(subdm);CHKERRQ(ierr); 3175 /* Build coordinates */ 3176 { 3177 PetscSection coordSection, subCoordSection; 3178 Vec coordinates, subCoordinates; 3179 PetscScalar *coords, *subCoords; 3180 PetscInt numComp, coordSize, v; 3181 const char *name; 3182 3183 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3184 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 3185 ierr = DMGetCoordinateSection(subdm, &subCoordSection);CHKERRQ(ierr); 3186 ierr = PetscSectionSetNumFields(subCoordSection, 1);CHKERRQ(ierr); 3187 ierr = PetscSectionGetFieldComponents(coordSection, 0, &numComp);CHKERRQ(ierr); 3188 ierr = PetscSectionSetFieldComponents(subCoordSection, 0, numComp);CHKERRQ(ierr); 3189 ierr = PetscSectionSetChart(subCoordSection, firstSubVertex, firstSubVertex+numSubVertices);CHKERRQ(ierr); 3190 for (v = 0; v < numSubVertices; ++v) { 3191 const PetscInt vertex = subVertices[v]; 3192 const PetscInt subvertex = firstSubVertex+v; 3193 PetscInt dof; 3194 3195 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3196 ierr = PetscSectionSetDof(subCoordSection, subvertex, dof);CHKERRQ(ierr); 3197 ierr = PetscSectionSetFieldDof(subCoordSection, subvertex, 0, dof);CHKERRQ(ierr); 3198 } 3199 ierr = PetscSectionSetUp(subCoordSection);CHKERRQ(ierr); 3200 ierr = PetscSectionGetStorageSize(subCoordSection, &coordSize);CHKERRQ(ierr); 3201 ierr = VecCreate(PETSC_COMM_SELF, &subCoordinates);CHKERRQ(ierr); 3202 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 3203 ierr = PetscObjectSetName((PetscObject)subCoordinates,name);CHKERRQ(ierr); 3204 ierr = VecSetSizes(subCoordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3205 ierr = VecSetType(subCoordinates,VECSTANDARD);CHKERRQ(ierr); 3206 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3207 ierr = VecGetArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3208 for (v = 0; v < numSubVertices; ++v) { 3209 const PetscInt vertex = subVertices[v]; 3210 const PetscInt subvertex = firstSubVertex+v; 3211 PetscInt dof, off, sdof, soff, d; 3212 3213 ierr = PetscSectionGetDof(coordSection, vertex, &dof);CHKERRQ(ierr); 3214 ierr = PetscSectionGetOffset(coordSection, vertex, &off);CHKERRQ(ierr); 3215 ierr = PetscSectionGetDof(subCoordSection, subvertex, &sdof);CHKERRQ(ierr); 3216 ierr = PetscSectionGetOffset(subCoordSection, subvertex, &soff);CHKERRQ(ierr); 3217 if (dof != sdof) SETERRQ4(comm, PETSC_ERR_PLIB, "Coordinate dimension %d on subvertex %d, vertex %d should be %d", sdof, subvertex, vertex, dof); 3218 for (d = 0; d < dof; ++d) subCoords[soff+d] = coords[off+d]; 3219 } 3220 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3221 ierr = VecRestoreArray(subCoordinates, &subCoords);CHKERRQ(ierr); 3222 ierr = DMSetCoordinatesLocal(subdm, subCoordinates);CHKERRQ(ierr); 3223 ierr = VecDestroy(&subCoordinates);CHKERRQ(ierr); 3224 } 3225 /* Build SF */ 3226 CHKMEMQ; 3227 { 3228 PetscSF sfPoint, sfPointSub; 3229 const PetscSFNode *remotePoints; 3230 PetscSFNode *sremotePoints, *newLocalPoints, *newOwners; 3231 const PetscInt *localPoints; 3232 PetscInt *slocalPoints; 3233 PetscInt numRoots, numLeaves, numSubRoots = numSubCells+numSubFaces+numSubVertices, numSubLeaves = 0, l, sl, ll, pStart, pEnd, p, vStart, vEnd; 3234 PetscMPIInt rank; 3235 3236 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 3237 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 3238 ierr = DMGetPointSF(subdm, &sfPointSub);CHKERRQ(ierr); 3239 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3240 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3241 ierr = PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 3242 if (numRoots >= 0) { 3243 /* Only vertices should be shared */ 3244 ierr = PetscMalloc2(pEnd-pStart,&newLocalPoints,numRoots,&newOwners);CHKERRQ(ierr); 3245 for (p = 0; p < pEnd-pStart; ++p) { 3246 newLocalPoints[p].rank = -2; 3247 newLocalPoints[p].index = -2; 3248 } 3249 /* Set subleaves */ 3250 for (l = 0; l < numLeaves; ++l) { 3251 const PetscInt point = localPoints[l]; 3252 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3253 3254 if ((point < vStart) && (point >= vEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Should not be mapping anything but vertices, %d", point); 3255 if (subPoint < 0) continue; 3256 newLocalPoints[point-pStart].rank = rank; 3257 newLocalPoints[point-pStart].index = subPoint; 3258 ++numSubLeaves; 3259 } 3260 /* Must put in owned subpoints */ 3261 for (p = pStart; p < pEnd; ++p) { 3262 const PetscInt subPoint = DMPlexFilterPoint_Internal(p, firstSubVertex, numSubVertices, subVertices); 3263 3264 if (subPoint < 0) { 3265 newOwners[p-pStart].rank = -3; 3266 newOwners[p-pStart].index = -3; 3267 } else { 3268 newOwners[p-pStart].rank = rank; 3269 newOwners[p-pStart].index = subPoint; 3270 } 3271 } 3272 ierr = PetscSFReduceBegin(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3273 ierr = PetscSFReduceEnd(sfPoint, MPIU_2INT, newLocalPoints, newOwners, MPI_MAXLOC);CHKERRQ(ierr); 3274 ierr = PetscSFBcastBegin(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3275 ierr = PetscSFBcastEnd(sfPoint, MPIU_2INT, newOwners, newLocalPoints);CHKERRQ(ierr); 3276 ierr = PetscMalloc1(numSubLeaves, &slocalPoints);CHKERRQ(ierr); 3277 ierr = PetscMalloc1(numSubLeaves, &sremotePoints);CHKERRQ(ierr); 3278 for (l = 0, sl = 0, ll = 0; l < numLeaves; ++l) { 3279 const PetscInt point = localPoints[l]; 3280 const PetscInt subPoint = DMPlexFilterPoint_Internal(point, firstSubVertex, numSubVertices, subVertices); 3281 3282 if (subPoint < 0) continue; 3283 if (newLocalPoints[point].rank == rank) {++ll; continue;} 3284 slocalPoints[sl] = subPoint; 3285 sremotePoints[sl].rank = newLocalPoints[point].rank; 3286 sremotePoints[sl].index = newLocalPoints[point].index; 3287 if (sremotePoints[sl].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote rank for local point %d", point); 3288 if (sremotePoints[sl].index < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid remote subpoint for local point %d", point); 3289 ++sl; 3290 } 3291 ierr = PetscFree2(newLocalPoints,newOwners);CHKERRQ(ierr); 3292 if (sl + ll != numSubLeaves) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch in number of subleaves %d + %d != %d", sl, ll, numSubLeaves); 3293 ierr = PetscSFSetGraph(sfPointSub, numSubRoots, sl, slocalPoints, PETSC_OWN_POINTER, sremotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 3294 } 3295 } 3296 CHKMEMQ; 3297 /* Cleanup */ 3298 if (subvertexIS) {ierr = ISRestoreIndices(subvertexIS, &subVertices);CHKERRQ(ierr);} 3299 ierr = ISDestroy(&subvertexIS);CHKERRQ(ierr); 3300 ierr = PetscFree(subCells);CHKERRQ(ierr); 3301 PetscFunctionReturn(0); 3302 } 3303 3304 static PetscErrorCode DMPlexCreateCohesiveSubmesh_Interpolated(DM dm, const char labelname[], PetscInt value, DM subdm) 3305 { 3306 DMLabel label = NULL; 3307 PetscErrorCode ierr; 3308 3309 PetscFunctionBegin; 3310 if (labelname) {ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr);} 3311 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, label, value, PETSC_TRUE, 1, subdm);CHKERRQ(ierr); 3312 PetscFunctionReturn(0); 3313 } 3314 3315 /*@C 3316 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. 3317 3318 Input Parameters: 3319 + dm - The original mesh 3320 . hasLagrange - The mesh has Lagrange unknowns in the cohesive cells 3321 . label - A label name, or NULL 3322 - value - A label value 3323 3324 Output Parameter: 3325 . subdm - The surface mesh 3326 3327 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3328 3329 Level: developer 3330 3331 .seealso: DMPlexGetSubpointMap(), DMPlexCreateSubmesh() 3332 @*/ 3333 PetscErrorCode DMPlexCreateCohesiveSubmesh(DM dm, PetscBool hasLagrange, const char label[], PetscInt value, DM *subdm) 3334 { 3335 PetscInt dim, depth; 3336 PetscErrorCode ierr; 3337 3338 PetscFunctionBegin; 3339 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3340 PetscValidPointer(subdm, 5); 3341 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3342 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3343 ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr); 3344 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3345 ierr = DMSetDimension(*subdm, dim-1);CHKERRQ(ierr); 3346 if (depth == dim) { 3347 ierr = DMPlexCreateCohesiveSubmesh_Interpolated(dm, label, value, *subdm);CHKERRQ(ierr); 3348 } else { 3349 ierr = DMPlexCreateCohesiveSubmesh_Uninterpolated(dm, hasLagrange, label, value, *subdm);CHKERRQ(ierr); 3350 } 3351 PetscFunctionReturn(0); 3352 } 3353 3354 /*@ 3355 DMPlexFilter - Extract a subset of mesh cells defined by a label as a separate mesh 3356 3357 Input Parameters: 3358 + dm - The original mesh 3359 . cellLabel - The DMLabel marking cells contained in the new mesh 3360 - value - The label value to use 3361 3362 Output Parameter: 3363 . subdm - The new mesh 3364 3365 Note: This function produces a DMLabel mapping original points in the submesh to their depth. This can be obtained using DMPlexGetSubpointMap(). 3366 3367 Level: developer 3368 3369 .seealso: DMPlexGetSubpointMap(), DMGetLabel(), DMLabelSetValue() 3370 @*/ 3371 PetscErrorCode DMPlexFilter(DM dm, DMLabel cellLabel, PetscInt value, DM *subdm) 3372 { 3373 PetscInt dim; 3374 PetscErrorCode ierr; 3375 3376 PetscFunctionBegin; 3377 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3378 PetscValidPointer(subdm, 3); 3379 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3380 ierr = DMCreate(PetscObjectComm((PetscObject) dm), subdm);CHKERRQ(ierr); 3381 ierr = DMSetType(*subdm, DMPLEX);CHKERRQ(ierr); 3382 ierr = DMSetDimension(*subdm, dim);CHKERRQ(ierr); 3383 /* Extract submesh in place, could be empty on some procs, could have inconsistency if procs do not both extract a shared cell */ 3384 ierr = DMPlexCreateSubmeshGeneric_Interpolated(dm, cellLabel, value, PETSC_FALSE, 0, *subdm);CHKERRQ(ierr); 3385 PetscFunctionReturn(0); 3386 } 3387 3388 /*@ 3389 DMPlexGetSubpointMap - Returns a DMLabel with point dimension as values 3390 3391 Input Parameter: 3392 . dm - The submesh DM 3393 3394 Output Parameter: 3395 . subpointMap - The DMLabel of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3396 3397 Level: developer 3398 3399 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3400 @*/ 3401 PetscErrorCode DMPlexGetSubpointMap(DM dm, DMLabel *subpointMap) 3402 { 3403 PetscFunctionBegin; 3404 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3405 PetscValidPointer(subpointMap, 2); 3406 *subpointMap = ((DM_Plex*) dm->data)->subpointMap; 3407 PetscFunctionReturn(0); 3408 } 3409 3410 /*@ 3411 DMPlexSetSubpointMap - Sets the DMLabel with point dimension as values 3412 3413 Input Parameters: 3414 + dm - The submesh DM 3415 - subpointMap - The DMLabel of all the points from the original mesh in this submesh 3416 3417 Note: Should normally not be called by the user, since it is set in DMPlexCreateSubmesh() 3418 3419 Level: developer 3420 3421 .seealso: DMPlexCreateSubmesh(), DMPlexCreateSubpointIS() 3422 @*/ 3423 PetscErrorCode DMPlexSetSubpointMap(DM dm, DMLabel subpointMap) 3424 { 3425 DM_Plex *mesh = (DM_Plex *) dm->data; 3426 DMLabel tmp; 3427 PetscErrorCode ierr; 3428 3429 PetscFunctionBegin; 3430 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3431 tmp = mesh->subpointMap; 3432 mesh->subpointMap = subpointMap; 3433 ++mesh->subpointMap->refct; 3434 ierr = DMLabelDestroy(&tmp);CHKERRQ(ierr); 3435 PetscFunctionReturn(0); 3436 } 3437 3438 /*@ 3439 DMPlexCreateSubpointIS - Creates an IS covering the entire subdm chart with the original points as data 3440 3441 Input Parameter: 3442 . dm - The submesh DM 3443 3444 Output Parameter: 3445 . subpointIS - The IS of all the points from the original mesh in this submesh, or NULL if this is not a submesh 3446 3447 Note: This IS is guaranteed to be sorted by the construction of the submesh 3448 3449 Level: developer 3450 3451 .seealso: DMPlexCreateSubmesh(), DMPlexGetSubpointMap() 3452 @*/ 3453 PetscErrorCode DMPlexCreateSubpointIS(DM dm, IS *subpointIS) 3454 { 3455 MPI_Comm comm; 3456 DMLabel subpointMap; 3457 IS is; 3458 const PetscInt *opoints; 3459 PetscInt *points, *depths; 3460 PetscInt depth, depStart, depEnd, d, pStart, pEnd, p, n, off; 3461 PetscErrorCode ierr; 3462 3463 PetscFunctionBegin; 3464 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3465 PetscValidPointer(subpointIS, 2); 3466 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 3467 *subpointIS = NULL; 3468 ierr = DMPlexGetSubpointMap(dm, &subpointMap);CHKERRQ(ierr); 3469 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3470 if (subpointMap && depth >= 0) { 3471 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 3472 if (pStart) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Submeshes must start the point numbering at 0, not %d", pStart); 3473 ierr = DMGetWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3474 depths[0] = depth; 3475 depths[1] = 0; 3476 for(d = 2; d <= depth; ++d) {depths[d] = depth+1 - d;} 3477 ierr = PetscMalloc1(pEnd, &points);CHKERRQ(ierr); 3478 for(d = 0, off = 0; d <= depth; ++d) { 3479 const PetscInt dep = depths[d]; 3480 3481 ierr = DMPlexGetDepthStratum(dm, dep, &depStart, &depEnd);CHKERRQ(ierr); 3482 ierr = DMLabelGetStratumSize(subpointMap, dep, &n);CHKERRQ(ierr); 3483 if (((d < 2) && (depth > 1)) || (d == 1)) { /* Only check vertices and cells for now since the map is broken for others */ 3484 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); 3485 } else { 3486 if (!n) { 3487 if (d == 0) { 3488 /* Missing cells */ 3489 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = -1; 3490 } else { 3491 /* Missing faces */ 3492 for(p = 0; p < depEnd-depStart; ++p, ++off) points[off] = PETSC_MAX_INT; 3493 } 3494 } 3495 } 3496 if (n) { 3497 ierr = DMLabelGetStratumIS(subpointMap, dep, &is);CHKERRQ(ierr); 3498 ierr = ISGetIndices(is, &opoints);CHKERRQ(ierr); 3499 for(p = 0; p < n; ++p, ++off) points[off] = opoints[p]; 3500 ierr = ISRestoreIndices(is, &opoints);CHKERRQ(ierr); 3501 ierr = ISDestroy(&is);CHKERRQ(ierr); 3502 } 3503 } 3504 ierr = DMRestoreWorkArray(dm, depth+1, PETSC_INT, &depths);CHKERRQ(ierr); 3505 if (off != pEnd) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "The number of mapped submesh points %d should be %d", off, pEnd); 3506 ierr = ISCreateGeneral(PETSC_COMM_SELF, pEnd, points, PETSC_OWN_POINTER, subpointIS);CHKERRQ(ierr); 3507 } 3508 PetscFunctionReturn(0); 3509 } 3510