1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/hashmapi.h> 3 #include <petsc/private/hashmapij.h> 4 5 const char *const DMPlexInterpolatedFlags[] = {"none", "partial", "mixed", "full", "DMPlexInterpolatedFlag", "DMPLEX_INTERPOLATED_", NULL}; 6 7 /* HMapIJKL */ 8 9 #include <petsc/private/hashmapijkl.h> 10 11 static PetscSFNode _PetscInvalidSFNode = {-1, -1}; 12 13 typedef struct _PetscHMapIJKLRemoteKey { 14 PetscSFNode i, j, k, l; 15 } PetscHMapIJKLRemoteKey; 16 17 #define PetscHMapIJKLRemoteKeyHash(key) \ 18 PetscHashCombine(PetscHashCombine(PetscHashInt((key).i.rank + (key).i.index), PetscHashInt((key).j.rank + (key).j.index)), PetscHashCombine(PetscHashInt((key).k.rank + (key).k.index), PetscHashInt((key).l.rank + (key).l.index))) 19 20 #define PetscHMapIJKLRemoteKeyEqual(k1, k2) \ 21 (((k1).i.rank == (k2).i.rank) ? ((k1).i.index == (k2).i.index) ? ((k1).j.rank == (k2).j.rank) ? ((k1).j.index == (k2).j.index) ? ((k1).k.rank == (k2).k.rank) ? ((k1).k.index == (k2).k.index) ? ((k1).l.rank == (k2).l.rank) ? ((k1).l.index == (k2).l.index) : 0 : 0 : 0 : 0 : 0 : 0 : 0) 22 23 PetscDisableStaticAnalyzerForExpressionUnderstandingThatThisIsDangerousAndBugprone(PETSC_HASH_MAP(HMapIJKLRemote, PetscHMapIJKLRemoteKey, PetscSFNode, PetscHMapIJKLRemoteKeyHash, PetscHMapIJKLRemoteKeyEqual, _PetscInvalidSFNode)) 24 25 static PetscErrorCode PetscSortSFNode(PetscInt n, PetscSFNode A[]) 26 { 27 PetscInt i; 28 29 PetscFunctionBegin; 30 for (i = 1; i < n; ++i) { 31 PetscSFNode x = A[i]; 32 PetscInt j; 33 34 for (j = i - 1; j >= 0; --j) { 35 if ((A[j].rank > x.rank) || (A[j].rank == x.rank && A[j].index > x.index)) break; 36 A[j + 1] = A[j]; 37 } 38 A[j + 1] = x; 39 } 40 PetscFunctionReturn(PETSC_SUCCESS); 41 } 42 43 /* 44 DMPlexGetRawFaces_Internal - Gets groups of vertices that correspond to faces for the given cone 45 */ 46 PetscErrorCode DMPlexGetRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[]) 47 { 48 DMPolytopeType *typesTmp = NULL; 49 PetscInt *sizesTmp = NULL, *facesTmp = NULL; 50 PetscInt *tmp; 51 PetscInt maxConeSize, maxSupportSize, maxSize; 52 PetscInt getSize = 0; 53 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 56 if (cone) PetscAssertPointer(cone, 3); 57 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 58 maxSize = PetscMax(maxConeSize, maxSupportSize); 59 if (faceTypes) getSize += maxSize; 60 if (faceSizes) getSize += maxSize; 61 if (faces) getSize += PetscSqr(maxSize); 62 PetscCall(DMGetWorkArray(dm, getSize, MPIU_INT, &tmp)); 63 if (faceTypes) { 64 typesTmp = (DMPolytopeType *)tmp; 65 tmp += maxSize; 66 } 67 if (faceSizes) { 68 sizesTmp = tmp; 69 tmp += maxSize; 70 } 71 if (faces) facesTmp = tmp; 72 switch (ct) { 73 case DM_POLYTOPE_POINT: 74 if (numFaces) *numFaces = 0; 75 if (faceTypes) *faceTypes = typesTmp; 76 if (faceSizes) *faceSizes = sizesTmp; 77 if (faces) *faces = facesTmp; 78 break; 79 case DM_POLYTOPE_SEGMENT: 80 if (numFaces) *numFaces = 2; 81 if (faceTypes) { 82 typesTmp[0] = DM_POLYTOPE_POINT; 83 typesTmp[1] = DM_POLYTOPE_POINT; 84 *faceTypes = typesTmp; 85 } 86 if (faceSizes) { 87 sizesTmp[0] = 1; 88 sizesTmp[1] = 1; 89 *faceSizes = sizesTmp; 90 } 91 if (faces) { 92 facesTmp[0] = cone[0]; 93 facesTmp[1] = cone[1]; 94 *faces = facesTmp; 95 } 96 break; 97 case DM_POLYTOPE_POINT_PRISM_TENSOR: 98 if (numFaces) *numFaces = 2; 99 if (faceTypes) { 100 typesTmp[0] = DM_POLYTOPE_POINT; 101 typesTmp[1] = DM_POLYTOPE_POINT; 102 *faceTypes = typesTmp; 103 } 104 if (faceSizes) { 105 sizesTmp[0] = 1; 106 sizesTmp[1] = 1; 107 *faceSizes = sizesTmp; 108 } 109 if (faces) { 110 facesTmp[0] = cone[0]; 111 facesTmp[1] = cone[1]; 112 *faces = facesTmp; 113 } 114 break; 115 case DM_POLYTOPE_TRIANGLE: 116 if (numFaces) *numFaces = 3; 117 if (faceTypes) { 118 typesTmp[0] = DM_POLYTOPE_SEGMENT; 119 typesTmp[1] = DM_POLYTOPE_SEGMENT; 120 typesTmp[2] = DM_POLYTOPE_SEGMENT; 121 *faceTypes = typesTmp; 122 } 123 if (faceSizes) { 124 sizesTmp[0] = 2; 125 sizesTmp[1] = 2; 126 sizesTmp[2] = 2; 127 *faceSizes = sizesTmp; 128 } 129 if (faces) { 130 facesTmp[0] = cone[0]; 131 facesTmp[1] = cone[1]; 132 facesTmp[2] = cone[1]; 133 facesTmp[3] = cone[2]; 134 facesTmp[4] = cone[2]; 135 facesTmp[5] = cone[0]; 136 *faces = facesTmp; 137 } 138 break; 139 case DM_POLYTOPE_QUADRILATERAL: 140 /* Vertices follow right hand rule */ 141 if (numFaces) *numFaces = 4; 142 if (faceTypes) { 143 typesTmp[0] = DM_POLYTOPE_SEGMENT; 144 typesTmp[1] = DM_POLYTOPE_SEGMENT; 145 typesTmp[2] = DM_POLYTOPE_SEGMENT; 146 typesTmp[3] = DM_POLYTOPE_SEGMENT; 147 *faceTypes = typesTmp; 148 } 149 if (faceSizes) { 150 sizesTmp[0] = 2; 151 sizesTmp[1] = 2; 152 sizesTmp[2] = 2; 153 sizesTmp[3] = 2; 154 *faceSizes = sizesTmp; 155 } 156 if (faces) { 157 facesTmp[0] = cone[0]; 158 facesTmp[1] = cone[1]; 159 facesTmp[2] = cone[1]; 160 facesTmp[3] = cone[2]; 161 facesTmp[4] = cone[2]; 162 facesTmp[5] = cone[3]; 163 facesTmp[6] = cone[3]; 164 facesTmp[7] = cone[0]; 165 *faces = facesTmp; 166 } 167 break; 168 case DM_POLYTOPE_SEG_PRISM_TENSOR: 169 if (numFaces) *numFaces = 4; 170 if (faceTypes) { 171 typesTmp[0] = DM_POLYTOPE_SEGMENT; 172 typesTmp[1] = DM_POLYTOPE_SEGMENT; 173 typesTmp[2] = DM_POLYTOPE_POINT_PRISM_TENSOR; 174 typesTmp[3] = DM_POLYTOPE_POINT_PRISM_TENSOR; 175 *faceTypes = typesTmp; 176 } 177 if (faceSizes) { 178 sizesTmp[0] = 2; 179 sizesTmp[1] = 2; 180 sizesTmp[2] = 2; 181 sizesTmp[3] = 2; 182 *faceSizes = sizesTmp; 183 } 184 if (faces) { 185 facesTmp[0] = cone[0]; 186 facesTmp[1] = cone[1]; 187 facesTmp[2] = cone[2]; 188 facesTmp[3] = cone[3]; 189 facesTmp[4] = cone[0]; 190 facesTmp[5] = cone[2]; 191 facesTmp[6] = cone[1]; 192 facesTmp[7] = cone[3]; 193 *faces = facesTmp; 194 } 195 break; 196 case DM_POLYTOPE_TETRAHEDRON: 197 /* Vertices of first face follow right hand rule and normal points away from last vertex */ 198 if (numFaces) *numFaces = 4; 199 if (faceTypes) { 200 typesTmp[0] = DM_POLYTOPE_TRIANGLE; 201 typesTmp[1] = DM_POLYTOPE_TRIANGLE; 202 typesTmp[2] = DM_POLYTOPE_TRIANGLE; 203 typesTmp[3] = DM_POLYTOPE_TRIANGLE; 204 *faceTypes = typesTmp; 205 } 206 if (faceSizes) { 207 sizesTmp[0] = 3; 208 sizesTmp[1] = 3; 209 sizesTmp[2] = 3; 210 sizesTmp[3] = 3; 211 *faceSizes = sizesTmp; 212 } 213 if (faces) { 214 facesTmp[0] = cone[0]; 215 facesTmp[1] = cone[1]; 216 facesTmp[2] = cone[2]; 217 facesTmp[3] = cone[0]; 218 facesTmp[4] = cone[3]; 219 facesTmp[5] = cone[1]; 220 facesTmp[6] = cone[0]; 221 facesTmp[7] = cone[2]; 222 facesTmp[8] = cone[3]; 223 facesTmp[9] = cone[2]; 224 facesTmp[10] = cone[1]; 225 facesTmp[11] = cone[3]; 226 *faces = facesTmp; 227 } 228 break; 229 case DM_POLYTOPE_HEXAHEDRON: 230 /* 7--------6 231 /| /| 232 / | / | 233 4--------5 | 234 | | | | 235 | | | | 236 | 1--------2 237 | / | / 238 |/ |/ 239 0--------3 240 */ 241 if (numFaces) *numFaces = 6; 242 if (faceTypes) { 243 typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; 244 typesTmp[1] = DM_POLYTOPE_QUADRILATERAL; 245 typesTmp[2] = DM_POLYTOPE_QUADRILATERAL; 246 typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; 247 typesTmp[4] = DM_POLYTOPE_QUADRILATERAL; 248 typesTmp[5] = DM_POLYTOPE_QUADRILATERAL; 249 *faceTypes = typesTmp; 250 } 251 if (faceSizes) { 252 sizesTmp[0] = 4; 253 sizesTmp[1] = 4; 254 sizesTmp[2] = 4; 255 sizesTmp[3] = 4; 256 sizesTmp[4] = 4; 257 sizesTmp[5] = 4; 258 *faceSizes = sizesTmp; 259 } 260 if (faces) { 261 facesTmp[0] = cone[0]; 262 facesTmp[1] = cone[1]; 263 facesTmp[2] = cone[2]; 264 facesTmp[3] = cone[3]; /* Bottom */ 265 facesTmp[4] = cone[4]; 266 facesTmp[5] = cone[5]; 267 facesTmp[6] = cone[6]; 268 facesTmp[7] = cone[7]; /* Top */ 269 facesTmp[8] = cone[0]; 270 facesTmp[9] = cone[3]; 271 facesTmp[10] = cone[5]; 272 facesTmp[11] = cone[4]; /* Front */ 273 facesTmp[12] = cone[2]; 274 facesTmp[13] = cone[1]; 275 facesTmp[14] = cone[7]; 276 facesTmp[15] = cone[6]; /* Back */ 277 facesTmp[16] = cone[3]; 278 facesTmp[17] = cone[2]; 279 facesTmp[18] = cone[6]; 280 facesTmp[19] = cone[5]; /* Right */ 281 facesTmp[20] = cone[0]; 282 facesTmp[21] = cone[4]; 283 facesTmp[22] = cone[7]; 284 facesTmp[23] = cone[1]; /* Left */ 285 *faces = facesTmp; 286 } 287 break; 288 case DM_POLYTOPE_TRI_PRISM: 289 if (numFaces) *numFaces = 5; 290 if (faceTypes) { 291 typesTmp[0] = DM_POLYTOPE_TRIANGLE; 292 typesTmp[1] = DM_POLYTOPE_TRIANGLE; 293 typesTmp[2] = DM_POLYTOPE_QUADRILATERAL; 294 typesTmp[3] = DM_POLYTOPE_QUADRILATERAL; 295 typesTmp[4] = DM_POLYTOPE_QUADRILATERAL; 296 *faceTypes = typesTmp; 297 } 298 if (faceSizes) { 299 sizesTmp[0] = 3; 300 sizesTmp[1] = 3; 301 sizesTmp[2] = 4; 302 sizesTmp[3] = 4; 303 sizesTmp[4] = 4; 304 *faceSizes = sizesTmp; 305 } 306 if (faces) { 307 facesTmp[0] = cone[0]; 308 facesTmp[1] = cone[1]; 309 facesTmp[2] = cone[2]; /* Bottom */ 310 facesTmp[3] = cone[3]; 311 facesTmp[4] = cone[4]; 312 facesTmp[5] = cone[5]; /* Top */ 313 facesTmp[6] = cone[0]; 314 facesTmp[7] = cone[2]; 315 facesTmp[8] = cone[4]; 316 facesTmp[9] = cone[3]; /* Back left */ 317 facesTmp[10] = cone[2]; 318 facesTmp[11] = cone[1]; 319 facesTmp[12] = cone[5]; 320 facesTmp[13] = cone[4]; /* Front */ 321 facesTmp[14] = cone[1]; 322 facesTmp[15] = cone[0]; 323 facesTmp[16] = cone[3]; 324 facesTmp[17] = cone[5]; /* Back right */ 325 *faces = facesTmp; 326 } 327 break; 328 case DM_POLYTOPE_TRI_PRISM_TENSOR: 329 if (numFaces) *numFaces = 5; 330 if (faceTypes) { 331 typesTmp[0] = DM_POLYTOPE_TRIANGLE; 332 typesTmp[1] = DM_POLYTOPE_TRIANGLE; 333 typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR; 334 typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; 335 typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR; 336 *faceTypes = typesTmp; 337 } 338 if (faceSizes) { 339 sizesTmp[0] = 3; 340 sizesTmp[1] = 3; 341 sizesTmp[2] = 4; 342 sizesTmp[3] = 4; 343 sizesTmp[4] = 4; 344 *faceSizes = sizesTmp; 345 } 346 if (faces) { 347 facesTmp[0] = cone[0]; 348 facesTmp[1] = cone[1]; 349 facesTmp[2] = cone[2]; /* Bottom */ 350 facesTmp[3] = cone[3]; 351 facesTmp[4] = cone[4]; 352 facesTmp[5] = cone[5]; /* Top */ 353 facesTmp[6] = cone[0]; 354 facesTmp[7] = cone[1]; 355 facesTmp[8] = cone[3]; 356 facesTmp[9] = cone[4]; /* Back left */ 357 facesTmp[10] = cone[1]; 358 facesTmp[11] = cone[2]; 359 facesTmp[12] = cone[4]; 360 facesTmp[13] = cone[5]; /* Back right */ 361 facesTmp[14] = cone[2]; 362 facesTmp[15] = cone[0]; 363 facesTmp[16] = cone[5]; 364 facesTmp[17] = cone[3]; /* Front */ 365 *faces = facesTmp; 366 } 367 break; 368 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 369 /* 7--------6 370 /| /| 371 / | / | 372 4--------5 | 373 | | | | 374 | | | | 375 | 3--------2 376 | / | / 377 |/ |/ 378 0--------1 379 */ 380 if (numFaces) *numFaces = 6; 381 if (faceTypes) { 382 typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; 383 typesTmp[1] = DM_POLYTOPE_QUADRILATERAL; 384 typesTmp[2] = DM_POLYTOPE_SEG_PRISM_TENSOR; 385 typesTmp[3] = DM_POLYTOPE_SEG_PRISM_TENSOR; 386 typesTmp[4] = DM_POLYTOPE_SEG_PRISM_TENSOR; 387 typesTmp[5] = DM_POLYTOPE_SEG_PRISM_TENSOR; 388 *faceTypes = typesTmp; 389 } 390 if (faceSizes) { 391 sizesTmp[0] = 4; 392 sizesTmp[1] = 4; 393 sizesTmp[2] = 4; 394 sizesTmp[3] = 4; 395 sizesTmp[4] = 4; 396 sizesTmp[5] = 4; 397 *faceSizes = sizesTmp; 398 } 399 if (faces) { 400 facesTmp[0] = cone[0]; 401 facesTmp[1] = cone[1]; 402 facesTmp[2] = cone[2]; 403 facesTmp[3] = cone[3]; /* Bottom */ 404 facesTmp[4] = cone[4]; 405 facesTmp[5] = cone[5]; 406 facesTmp[6] = cone[6]; 407 facesTmp[7] = cone[7]; /* Top */ 408 facesTmp[8] = cone[0]; 409 facesTmp[9] = cone[1]; 410 facesTmp[10] = cone[4]; 411 facesTmp[11] = cone[5]; /* Front */ 412 facesTmp[12] = cone[1]; 413 facesTmp[13] = cone[2]; 414 facesTmp[14] = cone[5]; 415 facesTmp[15] = cone[6]; /* Right */ 416 facesTmp[16] = cone[2]; 417 facesTmp[17] = cone[3]; 418 facesTmp[18] = cone[6]; 419 facesTmp[19] = cone[7]; /* Back */ 420 facesTmp[20] = cone[3]; 421 facesTmp[21] = cone[0]; 422 facesTmp[22] = cone[7]; 423 facesTmp[23] = cone[4]; /* Left */ 424 *faces = facesTmp; 425 } 426 break; 427 case DM_POLYTOPE_PYRAMID: 428 /* 429 4---- 430 |\-\ \----- 431 | \ -\ \ 432 | 1--\-----2 433 | / \ / 434 |/ \ / 435 0--------3 436 */ 437 if (numFaces) *numFaces = 5; 438 if (faceTypes) { 439 typesTmp[0] = DM_POLYTOPE_QUADRILATERAL; 440 typesTmp[1] = DM_POLYTOPE_TRIANGLE; 441 typesTmp[2] = DM_POLYTOPE_TRIANGLE; 442 typesTmp[3] = DM_POLYTOPE_TRIANGLE; 443 typesTmp[4] = DM_POLYTOPE_TRIANGLE; 444 *faceTypes = typesTmp; 445 } 446 if (faceSizes) { 447 sizesTmp[0] = 4; 448 sizesTmp[1] = 3; 449 sizesTmp[2] = 3; 450 sizesTmp[3] = 3; 451 sizesTmp[4] = 3; 452 *faceSizes = sizesTmp; 453 } 454 if (faces) { 455 facesTmp[0] = cone[0]; 456 facesTmp[1] = cone[1]; 457 facesTmp[2] = cone[2]; 458 facesTmp[3] = cone[3]; /* Bottom */ 459 facesTmp[4] = cone[0]; 460 facesTmp[5] = cone[3]; 461 facesTmp[6] = cone[4]; /* Front */ 462 facesTmp[7] = cone[3]; 463 facesTmp[8] = cone[2]; 464 facesTmp[9] = cone[4]; /* Right */ 465 facesTmp[10] = cone[2]; 466 facesTmp[11] = cone[1]; 467 facesTmp[12] = cone[4]; /* Back */ 468 facesTmp[13] = cone[1]; 469 facesTmp[14] = cone[0]; 470 facesTmp[15] = cone[4]; /* Left */ 471 *faces = facesTmp; 472 } 473 break; 474 default: 475 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No face description for cell type %s", DMPolytopeTypes[ct]); 476 } 477 PetscFunctionReturn(PETSC_SUCCESS); 478 } 479 480 PetscErrorCode DMPlexRestoreRawFaces_Internal(DM dm, DMPolytopeType ct, const PetscInt cone[], PetscInt *numFaces, const DMPolytopeType *faceTypes[], const PetscInt *faceSizes[], const PetscInt *faces[]) 481 { 482 PetscFunctionBegin; 483 if (faceTypes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceTypes)); 484 else if (faceSizes) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faceSizes)); 485 else if (faces) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, (void *)faces)); 486 if (faceTypes) *faceTypes = NULL; 487 if (faceSizes) *faceSizes = NULL; 488 if (faces) *faces = NULL; 489 PetscFunctionReturn(PETSC_SUCCESS); 490 } 491 492 /* This interpolates faces for cells at some stratum */ 493 PetscErrorCode DMPlexInterpolateFaces_Internal(DM dm, PetscInt cellDepth, DM idm) 494 { 495 DMLabel ctLabel; 496 PetscHMapIJKL faceTable; 497 PetscInt faceTypeNum[DM_NUM_POLYTOPES]; 498 PetscInt depth, pStart, Np, cStart, cEnd, fStart, fEnd, vStart, vEnd; 499 PetscInt cntFaces, *facesId, minCone; 500 501 PetscFunctionBegin; 502 PetscCall(DMPlexGetDepth(dm, &depth)); 503 PetscCall(PetscHMapIJKLCreate(&faceTable)); 504 PetscCall(PetscArrayzero(faceTypeNum, DM_NUM_POLYTOPES)); 505 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 506 PetscCall(DMPlexGetDepthStratum(dm, cellDepth, &cStart, &cEnd)); 507 // If the range incorporates the vertices, it means we have a non-manifold topology, so choose just cells 508 if (cStart <= vStart && cEnd >= vEnd) cEnd = vStart; 509 // Number new faces and save face vertices in hash table 510 // If depth > cellDepth, meaning we are interpolating faces, put the new (d-1)-faces after them 511 // otherwise, we are interpolating cells, so put the faces after the vertices 512 PetscCall(DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart)); 513 fEnd = fStart; 514 515 minCone = PETSC_INT_MAX; 516 cntFaces = 0; 517 for (PetscInt c = cStart; c < cEnd; ++c) { 518 const PetscInt *cone; 519 DMPolytopeType ct; 520 PetscInt numFaces = 0, coneSize; 521 522 PetscCall(DMPlexGetCellType(dm, c, &ct)); 523 PetscCall(DMPlexGetCone(dm, c, &cone)); 524 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 525 for (PetscInt j = 0; j < coneSize; j++) minCone = PetscMin(cone[j], minCone); 526 // Ignore faces since they are interpolated 527 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, NULL, NULL, NULL)); 528 cntFaces += numFaces; 529 } 530 // Encode so that we can use 0 as an excluded value, instead of PETSC_INT_MAX 531 minCone = -(minCone - 1); 532 533 PetscCall(PetscMalloc1(cntFaces, &facesId)); 534 535 cntFaces = 0; 536 for (PetscInt c = cStart; c < cEnd; ++c) { 537 const PetscInt *cone, *faceSizes, *faces; 538 const DMPolytopeType *faceTypes; 539 DMPolytopeType ct; 540 PetscInt numFaces, foff = 0; 541 542 PetscCall(DMPlexGetCellType(dm, c, &ct)); 543 PetscCall(DMPlexGetCone(dm, c, &cone)); 544 // Ignore faces since they are interpolated 545 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) { 546 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 547 } else { 548 numFaces = 0; 549 } 550 for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) { 551 const PetscInt faceSize = faceSizes[cf]; 552 const DMPolytopeType faceType = faceTypes[cf]; 553 const PetscInt *face = &faces[foff]; 554 PetscHashIJKLKey key; 555 PetscHashIter iter; 556 PetscBool missing; 557 558 PetscCheck(faceSize <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Do not support faces of size %" PetscInt_FMT " > 4", faceSize); 559 key.i = face[0] + minCone; 560 key.j = faceSize > 1 ? face[1] + minCone : 0; 561 key.k = faceSize > 2 ? face[2] + minCone : 0; 562 key.l = faceSize > 3 ? face[3] + minCone : 0; 563 PetscCall(PetscSortInt(faceSize, (PetscInt *)&key)); 564 PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing)); 565 if (missing) { 566 facesId[cntFaces] = fEnd; 567 PetscCall(PetscHMapIJKLIterSet(faceTable, iter, fEnd++)); 568 ++faceTypeNum[faceType]; 569 } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces])); 570 cntFaces++; 571 } 572 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 573 } 574 /* We need to number faces contiguously among types */ 575 { 576 PetscInt faceTypeStart[DM_NUM_POLYTOPES], ct, numFT = 0; 577 578 for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) { 579 if (faceTypeNum[ct]) ++numFT; 580 faceTypeStart[ct] = 0; 581 } 582 if (numFT > 1) { 583 PetscCall(PetscHMapIJKLClear(faceTable)); 584 faceTypeStart[0] = fStart; 585 for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) faceTypeStart[ct] = faceTypeStart[ct - 1] + faceTypeNum[ct - 1]; 586 cntFaces = 0; 587 for (PetscInt c = cStart; c < cEnd; ++c) { 588 const PetscInt *cone, *faceSizes, *faces; 589 const DMPolytopeType *faceTypes; 590 DMPolytopeType ct; 591 PetscInt numFaces, foff = 0; 592 593 PetscCall(DMPlexGetCellType(dm, c, &ct)); 594 PetscCall(DMPlexGetCone(dm, c, &cone)); 595 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) { 596 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 597 } else { 598 numFaces = 0; 599 } 600 for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) { 601 const PetscInt faceSize = faceSizes[cf]; 602 const DMPolytopeType faceType = faceTypes[cf]; 603 const PetscInt *face = &faces[foff]; 604 PetscHashIJKLKey key; 605 PetscHashIter iter; 606 PetscBool missing; 607 608 key.i = face[0] + minCone; 609 key.j = faceSize > 1 ? face[1] + minCone : 0; 610 key.k = faceSize > 2 ? face[2] + minCone : 0; 611 key.l = faceSize > 3 ? face[3] + minCone : 0; 612 PetscCall(PetscSortInt(faceSize, (PetscInt *)&key)); 613 PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing)); 614 if (missing) { 615 facesId[cntFaces] = faceTypeStart[faceType]; 616 PetscCall(PetscHMapIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++)); 617 } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces])); 618 cntFaces++; 619 } 620 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 621 } 622 for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) { 623 PetscCheck(faceTypeStart[ct] == faceTypeStart[ct - 1] + faceTypeNum[ct], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent numbering for cell type %s, %" PetscInt_FMT " != %" PetscInt_FMT " + %" PetscInt_FMT, DMPolytopeTypes[ct], faceTypeStart[ct], faceTypeStart[ct - 1], faceTypeNum[ct]); 624 } 625 } 626 } 627 PetscCall(PetscHMapIJKLDestroy(&faceTable)); 628 629 // Add new points, perhaps inserting into the numbering 630 PetscCall(DMPlexGetChart(dm, &pStart, &Np)); 631 PetscCall(DMPlexSetChart(idm, pStart, Np + (fEnd - fStart))); 632 // Set cone sizes 633 // Must create the celltype label here so that we do not automatically try to compute the types 634 PetscCall(DMCreateLabel(idm, "celltype")); 635 PetscCall(DMPlexGetCellTypeLabel(idm, &ctLabel)); 636 for (PetscInt d = 0; d <= depth; ++d) { 637 DMPolytopeType ct; 638 PetscInt coneSize, pStart, pEnd, poff = 0; 639 640 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 641 // Check for non-manifold condition 642 if (d == cellDepth) { 643 if (pEnd == cEnd) continue; 644 else pStart = vEnd; 645 } 646 // Account for insertion 647 if (pStart >= fStart) poff = fEnd - fStart; 648 for (PetscInt p = pStart; p < pEnd; ++p) { 649 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 650 PetscCall(DMPlexSetConeSize(idm, p + poff, coneSize)); 651 PetscCall(DMPlexGetCellType(dm, p, &ct)); 652 PetscCall(DMPlexSetCellType(idm, p + poff, ct)); 653 } 654 } 655 cntFaces = 0; 656 for (PetscInt c = cStart; c < cEnd; ++c) { 657 const PetscInt *cone, *faceSizes; 658 const DMPolytopeType *faceTypes; 659 DMPolytopeType ct; 660 PetscInt numFaces, poff = 0; 661 662 PetscCall(DMPlexGetCellType(dm, c, &ct)); 663 PetscCall(DMPlexGetCone(dm, c, &cone)); 664 if (c >= fStart) poff = fEnd - fStart; 665 if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) { 666 PetscCall(DMPlexSetCellType(idm, c + poff, ct)); 667 PetscCall(DMPlexSetConeSize(idm, c + poff, 2)); 668 continue; 669 } 670 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL)); 671 PetscCall(DMPlexSetCellType(idm, c + poff, ct)); 672 PetscCall(DMPlexSetConeSize(idm, c + poff, numFaces)); 673 for (PetscInt cf = 0; cf < numFaces; ++cf) { 674 const PetscInt f = facesId[cntFaces]; 675 DMPolytopeType faceType = faceTypes[cf]; 676 const PetscInt faceSize = faceSizes[cf]; 677 PetscCall(DMPlexSetConeSize(idm, f, faceSize)); 678 PetscCall(DMPlexSetCellType(idm, f, faceType)); 679 cntFaces++; 680 } 681 PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL)); 682 } 683 PetscCall(DMSetUp(idm)); 684 // Initialize cones so we do not need the bash table to tell us that a cone has been set 685 { 686 PetscSection cs; 687 PetscInt *cones, csize; 688 689 PetscCall(DMPlexGetConeSection(idm, &cs)); 690 PetscCall(DMPlexGetCones(idm, &cones)); 691 PetscCall(PetscSectionGetStorageSize(cs, &csize)); 692 for (PetscInt c = 0; c < csize; ++c) cones[c] = -1; 693 } 694 // Set cones 695 { 696 PetscInt *icone; 697 PetscInt maxConeSize; 698 699 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL)); 700 PetscCall(PetscMalloc1(maxConeSize, &icone)); 701 for (PetscInt d = 0; d <= depth; ++d) { 702 const PetscInt *cone; 703 PetscInt pStart, pEnd, poff = 0, coneSize; 704 705 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 706 // Check for non-manifold condition 707 if (d == cellDepth) { 708 if (pEnd == cEnd) continue; 709 else pStart = vEnd; 710 } 711 // Account for insertion 712 if (pStart >= fStart) poff = fEnd - fStart; 713 for (PetscInt p = pStart; p < pEnd; ++p) { 714 PetscCall(DMPlexGetCone(dm, p, &cone)); 715 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 716 for (PetscInt cp = 0; cp < coneSize; ++cp) icone[cp] = cone[cp] + (cone[cp] >= fStart ? fEnd - fStart : 0); 717 PetscCall(DMPlexSetCone(idm, p + poff, icone)); 718 PetscCall(DMPlexGetConeOrientation(dm, p, &cone)); 719 PetscCall(DMPlexSetConeOrientation(idm, p + poff, cone)); 720 } 721 } 722 cntFaces = 0; 723 for (PetscInt c = cStart; c < cEnd; ++c) { 724 const PetscInt *cone, *faceSizes, *faces; 725 const DMPolytopeType *faceTypes; 726 DMPolytopeType ct; 727 PetscInt coneSize, numFaces, foff = 0, poff = 0; 728 729 PetscCall(DMPlexGetCellType(dm, c, &ct)); 730 PetscCall(DMPlexGetCone(dm, c, &cone)); 731 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 732 if (c >= fStart) poff = fEnd - fStart; 733 if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) { 734 for (PetscInt cp = 0; cp < coneSize; ++cp) icone[cp] = cone[cp] + (cone[cp] >= fStart ? fEnd - fStart : 0); 735 PetscCall(DMPlexSetCone(idm, c + poff, icone)); 736 PetscCall(DMPlexGetConeOrientation(dm, c, &cone)); 737 PetscCall(DMPlexSetConeOrientation(idm, c + poff, cone)); 738 continue; 739 } 740 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 741 for (PetscInt cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) { 742 DMPolytopeType faceType = faceTypes[cf]; 743 const PetscInt faceSize = faceSizes[cf]; 744 const PetscInt f = facesId[cntFaces]; 745 const PetscInt *face = &faces[foff]; 746 const PetscInt *fcone; 747 748 PetscCall(DMPlexInsertCone(idm, c, cf, f)); 749 PetscCall(DMPlexGetCone(idm, f, &fcone)); 750 if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face)); 751 { 752 const PetscInt *fcone2; 753 PetscInt ornt; 754 755 PetscCall(DMPlexGetConeSize(idm, f, &coneSize)); 756 PetscCall(DMPlexGetCone(idm, f, &fcone2)); 757 PetscCheck(coneSize == faceSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of face vertices %" PetscInt_FMT " for face %" PetscInt_FMT " should be %" PetscInt_FMT, coneSize, f, faceSize); 758 /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */ 759 PetscCall(DMPolytopeGetVertexOrientation(faceType, fcone2, face, &ornt)); 760 PetscCall(DMPlexInsertConeOrientation(idm, c + poff, cf, ornt)); 761 } 762 cntFaces++; 763 } 764 PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 765 } 766 PetscCall(PetscFree(icone)); 767 } 768 PetscCall(PetscFree(facesId)); 769 PetscCall(DMPlexSymmetrize(idm)); 770 PetscCall(DMPlexStratify(idm)); 771 PetscFunctionReturn(PETSC_SUCCESS); 772 } 773 774 static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 775 { 776 PetscInt nleaves; 777 PetscMPIInt nranks; 778 const PetscMPIInt *ranks = NULL; 779 const PetscInt *roffset = NULL, *rmine = NULL, *rremote = NULL; 780 PetscInt n, o; 781 782 PetscFunctionBegin; 783 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote)); 784 nleaves = roffset[nranks]; 785 PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1)); 786 for (PetscMPIInt r = 0; r < nranks; r++) { 787 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 788 - to unify order with the other side */ 789 o = roffset[r]; 790 n = roffset[r + 1] - o; 791 PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n)); 792 PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n)); 793 PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o])); 794 } 795 PetscFunctionReturn(PETSC_SUCCESS); 796 } 797 798 PetscErrorCode DMPlexOrientInterface_Internal(DM dm) 799 { 800 PetscSF sf; 801 const PetscInt *locals; 802 const PetscSFNode *remotes; 803 const PetscMPIInt *ranks; 804 const PetscInt *roffset; 805 PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 806 PetscInt nroots, p, nleaves, maxConeSize = 0; 807 PetscMPIInt nranks, r; 808 PetscInt (*roots)[4], (*leaves)[4], mainCone[4]; 809 PetscMPIInt (*rootsRanks)[4], (*leavesRanks)[4]; 810 MPI_Comm comm; 811 PetscMPIInt rank, size; 812 PetscInt debug = 0; 813 814 PetscFunctionBegin; 815 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 816 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 817 PetscCallMPI(MPI_Comm_size(comm, &size)); 818 PetscCall(DMGetPointSF(dm, &sf)); 819 PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view")); 820 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE)); 821 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes)); 822 if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS); 823 PetscCall(PetscSFSetUp(sf)); 824 PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1)); 825 for (p = 0; p < nleaves; ++p) { 826 PetscInt coneSize; 827 PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize)); 828 maxConeSize = PetscMax(maxConeSize, coneSize); 829 } 830 PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize); 831 PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks)); 832 for (p = 0; p < nroots; ++p) { 833 const PetscInt *cone; 834 PetscInt coneSize, c, ind0; 835 836 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 837 PetscCall(DMPlexGetCone(dm, p, &cone)); 838 /* Ignore vertices */ 839 if (coneSize < 2) { 840 for (c = 0; c < 4; c++) { 841 roots[p][c] = -1; 842 rootsRanks[p][c] = -1; 843 } 844 continue; 845 } 846 /* Translate all points to root numbering */ 847 for (c = 0; c < PetscMin(coneSize, 4); c++) { 848 PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0)); 849 if (ind0 < 0) { 850 roots[p][c] = cone[c]; 851 rootsRanks[p][c] = rank; 852 } else { 853 roots[p][c] = remotes[ind0].index; 854 PetscCall(PetscMPIIntCast(remotes[ind0].rank, &rootsRanks[p][c])); 855 } 856 } 857 for (c = coneSize; c < 4; c++) { 858 roots[p][c] = -1; 859 rootsRanks[p][c] = -1; 860 } 861 } 862 for (p = 0; p < nroots; ++p) { 863 PetscInt c; 864 for (c = 0; c < 4; c++) { 865 leaves[p][c] = -2; 866 leavesRanks[p][c] = -2; 867 } 868 } 869 PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE)); 870 PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE)); 871 PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE)); 872 PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE)); 873 if (debug) { 874 PetscCall(PetscSynchronizedFlush(comm, NULL)); 875 if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n")); 876 } 877 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL)); 878 for (p = 0; p < nroots; ++p) { 879 DMPolytopeType ct; 880 const PetscInt *cone; 881 PetscInt coneSize, c, ind0, o, ir; 882 883 if (leaves[p][0] < 0) continue; /* Ignore vertices */ 884 PetscCall(DMPlexGetCellType(dm, p, &ct)); 885 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 886 PetscCall(DMPlexGetCone(dm, p, &cone)); 887 if (debug) { 888 PetscCall(PetscSynchronizedPrintf(comm, "[%d] %4" PetscInt_FMT ": cone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "] roots=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")] leaves=[(%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ") (%d,%4" PetscInt_FMT ")]", rank, p, cone[0], cone[1], cone[2], cone[3], rootsRanks[p][0], roots[p][0], rootsRanks[p][1], roots[p][1], rootsRanks[p][2], roots[p][2], rootsRanks[p][3], roots[p][3], leavesRanks[p][0], leaves[p][0], leavesRanks[p][1], leaves[p][1], leavesRanks[p][2], leaves[p][2], leavesRanks[p][3], leaves[p][3])); 889 } 890 if (leavesRanks[p][0] != rootsRanks[p][0] || leaves[p][0] != roots[p][0] || leavesRanks[p][1] != rootsRanks[p][1] || leaves[p][1] != roots[p][1] || leavesRanks[p][2] != rootsRanks[p][2] || leaves[p][2] != roots[p][2] || leavesRanks[p][3] != rootsRanks[p][3] || leaves[p][3] != roots[p][3]) { 891 /* Translate these leaves to my cone points; mainCone means desired order p's cone points */ 892 for (c = 0; c < PetscMin(coneSize, 4); ++c) { 893 PetscInt rS, rN; 894 895 if (leavesRanks[p][c] == rank) { 896 /* A local leaf is just taken as it is */ 897 mainCone[c] = leaves[p][c]; 898 continue; 899 } 900 /* Find index of rank leavesRanks[p][c] among remote ranks */ 901 /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */ 902 PetscCall(PetscFindMPIInt(leavesRanks[p][c], nranks, ranks, &ir)); 903 PetscCall(PetscMPIIntCast(ir, &r)); 904 PetscCheck(r >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leaf (%d,%" PetscInt_FMT "): leaf rank not found among remote ranks", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]); 905 PetscCheck(ranks[r] >= 0 && ranks[r] < size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%" PetscInt_FMT " c=%" PetscInt_FMT " commsize=%d: ranks[%d] = %d makes no sense", p, c, size, r, ranks[r]); 906 /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */ 907 rS = roffset[r]; 908 rN = roffset[r + 1] - rS; 909 PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0)); 910 PetscCheck(ind0 >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " cone[%" PetscInt_FMT "]=%" PetscInt_FMT " root (%d,%" PetscInt_FMT ") leave (%d,%" PetscInt_FMT "): corresponding remote point not found - it seems there is missing connection in point SF!", p, c, cone[c], rootsRanks[p][c], roots[p][c], leavesRanks[p][c], leaves[p][c]); 911 /* Get the corresponding local point */ 912 mainCone[c] = rmine1[rS + ind0]; 913 } 914 if (debug) PetscCall(PetscSynchronizedPrintf(comm, " mainCone=[%4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT " %4" PetscInt_FMT "]\n", mainCone[0], mainCone[1], mainCone[2], mainCone[3])); 915 /* Set the desired order of p's cone points and fix orientations accordingly */ 916 PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o)); 917 PetscCall(DMPlexOrientPoint(dm, p, o)); 918 } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n")); 919 } 920 if (debug) { 921 PetscCall(PetscSynchronizedFlush(comm, NULL)); 922 PetscCallMPI(MPI_Barrier(comm)); 923 } 924 PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view")); 925 PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks)); 926 PetscCall(PetscFree2(rmine1, rremote1)); 927 PetscFunctionReturn(PETSC_SUCCESS); 928 } 929 930 static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[]) 931 { 932 PetscInt idx; 933 PetscMPIInt rank; 934 PetscBool flg; 935 936 PetscFunctionBegin; 937 PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg)); 938 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 939 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 940 PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name)); 941 for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx])); 942 PetscCall(PetscSynchronizedFlush(comm, NULL)); 943 PetscFunctionReturn(PETSC_SUCCESS); 944 } 945 946 static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[]) 947 { 948 PetscInt idx; 949 PetscMPIInt rank; 950 PetscBool flg; 951 952 PetscFunctionBegin; 953 PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg)); 954 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 955 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 956 PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name)); 957 if (idxname) { 958 for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, idxname, idx, a[idx].rank, a[idx].index)); 959 } else { 960 for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index)); 961 } 962 PetscCall(PetscSynchronizedFlush(comm, NULL)); 963 PetscFunctionReturn(PETSC_SUCCESS); 964 } 965 966 static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed) 967 { 968 PetscSF sf; 969 const PetscInt *locals; 970 PetscMPIInt rank; 971 972 PetscFunctionBegin; 973 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 974 PetscCall(DMGetPointSF(dm, &sf)); 975 PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL)); 976 if (mapFailed) *mapFailed = PETSC_FALSE; 977 if (remotePoint.rank == rank) { 978 *localPoint = remotePoint.index; 979 } else { 980 PetscHashIJKey key; 981 PetscInt l; 982 983 key.i = remotePoint.index; 984 key.j = remotePoint.rank; 985 PetscCall(PetscHMapIJGet(remotehash, key, &l)); 986 if (l >= 0) { 987 *localPoint = locals[l]; 988 } else if (mapFailed) *mapFailed = PETSC_TRUE; 989 } 990 PetscFunctionReturn(PETSC_SUCCESS); 991 } 992 993 static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed) 994 { 995 PetscSF sf; 996 const PetscInt *locals, *rootdegree; 997 const PetscSFNode *remotes; 998 PetscInt Nl, l; 999 PetscMPIInt rank; 1000 1001 PetscFunctionBegin; 1002 if (mapFailed) *mapFailed = PETSC_FALSE; 1003 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 1004 PetscCall(DMGetPointSF(dm, &sf)); 1005 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes)); 1006 if (Nl < 0) goto owned; 1007 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 1008 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 1009 if (rootdegree[localPoint]) goto owned; 1010 PetscCall(PetscFindInt(localPoint, Nl, locals, &l)); 1011 if (l < 0) { 1012 if (mapFailed) *mapFailed = PETSC_TRUE; 1013 } else *remotePoint = remotes[l]; 1014 PetscFunctionReturn(PETSC_SUCCESS); 1015 owned: 1016 remotePoint->rank = rank; 1017 remotePoint->index = localPoint; 1018 PetscFunctionReturn(PETSC_SUCCESS); 1019 } 1020 1021 static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared) 1022 { 1023 PetscSF sf; 1024 const PetscInt *locals, *rootdegree; 1025 PetscInt Nl, idx; 1026 1027 PetscFunctionBegin; 1028 *isShared = PETSC_FALSE; 1029 PetscCall(DMGetPointSF(dm, &sf)); 1030 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL)); 1031 if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS); 1032 PetscCall(PetscFindInt(p, Nl, locals, &idx)); 1033 if (idx >= 0) { 1034 *isShared = PETSC_TRUE; 1035 PetscFunctionReturn(PETSC_SUCCESS); 1036 } 1037 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 1038 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 1039 if (rootdegree[p] > 0) *isShared = PETSC_TRUE; 1040 PetscFunctionReturn(PETSC_SUCCESS); 1041 } 1042 1043 static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared) 1044 { 1045 const PetscInt *cone; 1046 PetscInt coneSize, c; 1047 PetscBool cShared = PETSC_TRUE; 1048 1049 PetscFunctionBegin; 1050 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1051 PetscCall(DMPlexGetCone(dm, p, &cone)); 1052 for (c = 0; c < coneSize; ++c) { 1053 PetscBool pointShared; 1054 1055 PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared)); 1056 cShared = (PetscBool)(cShared && pointShared); 1057 } 1058 *isShared = coneSize ? cShared : PETSC_FALSE; 1059 PetscFunctionReturn(PETSC_SUCCESS); 1060 } 1061 1062 static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin) 1063 { 1064 const PetscInt *cone; 1065 PetscInt coneSize, c; 1066 PetscSFNode cmin = {PETSC_INT_MAX, PETSC_MPI_INT_MAX}, missing = {-1, -1}; 1067 1068 PetscFunctionBegin; 1069 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1070 PetscCall(DMPlexGetCone(dm, p, &cone)); 1071 for (c = 0; c < coneSize; ++c) { 1072 PetscSFNode rcp; 1073 PetscBool mapFailed; 1074 1075 PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed)); 1076 if (mapFailed) { 1077 cmin = missing; 1078 } else { 1079 cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin; 1080 } 1081 } 1082 *cpmin = coneSize ? cmin : missing; 1083 PetscFunctionReturn(PETSC_SUCCESS); 1084 } 1085 1086 /* 1087 Each shared face has an entry in the candidates array: 1088 (-1, coneSize-1), {(global cone point)} 1089 where the set is missing the point p which we use as the key for the face 1090 */ 1091 static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug) 1092 { 1093 MPI_Comm comm; 1094 const PetscInt *support; 1095 PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height; 1096 PetscMPIInt rank; 1097 1098 PetscFunctionBegin; 1099 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1100 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1101 PetscCall(DMPlexGetOverlap(dm, &overlap)); 1102 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 1103 PetscCall(DMPlexGetPointHeight(dm, p, &height)); 1104 if (!overlap && height <= cellHeight + 1) { 1105 /* cells can't be shared for non-overlapping meshes */ 1106 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p)); 1107 PetscFunctionReturn(PETSC_SUCCESS); 1108 } 1109 PetscCall(DMPlexGetSupportSize(dm, p, &supportSize)); 1110 PetscCall(DMPlexGetSupport(dm, p, &support)); 1111 if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off)); 1112 for (s = 0; s < supportSize; ++s) { 1113 const PetscInt face = support[s]; 1114 const PetscInt *cone; 1115 PetscSFNode cpmin = {-1, -1}, rp = {-1, -1}; 1116 PetscInt coneSize, c, f; 1117 PetscBool isShared = PETSC_FALSE; 1118 PetscHashIJKey key; 1119 1120 /* Only add point once */ 1121 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Support face %" PetscInt_FMT "\n", rank, face)); 1122 key.i = p; 1123 key.j = face; 1124 PetscCall(PetscHMapIJGet(faceHash, key, &f)); 1125 if (f >= 0) continue; 1126 PetscCall(DMPlexConeIsShared(dm, face, &isShared)); 1127 PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin)); 1128 PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL)); 1129 if (debug) { 1130 PetscCall(PetscSynchronizedPrintf(comm, "[%d] Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared)); 1131 PetscCall(PetscSynchronizedPrintf(comm, "[%d] Global point (%" PetscInt_FMT ", %" PetscInt_FMT ") Min Cone Point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rp.rank, rp.index, cpmin.rank, cpmin.index)); 1132 } 1133 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) { 1134 PetscCall(PetscHMapIJSet(faceHash, key, p)); 1135 if (candidates) { 1136 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d] ", rank, face, idx, rank)); 1137 PetscCall(DMPlexGetConeSize(dm, face, &coneSize)); 1138 PetscCall(DMPlexGetCone(dm, face, &cone)); 1139 candidates[off + idx].rank = -1; 1140 candidates[off + idx++].index = coneSize - 1; 1141 candidates[off + idx].rank = rank; 1142 candidates[off + idx++].index = face; 1143 for (c = 0; c < coneSize; ++c) { 1144 const PetscInt cp = cone[c]; 1145 1146 if (cp == p) continue; 1147 PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL)); 1148 if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT "),%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index)); 1149 ++idx; 1150 } 1151 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n")); 1152 } else { 1153 /* Add cone size to section */ 1154 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %" PetscInt_FMT "\n", rank, face)); 1155 PetscCall(DMPlexGetConeSize(dm, face, &coneSize)); 1156 PetscCall(PetscHMapIJSet(faceHash, key, p)); 1157 PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1)); 1158 } 1159 } 1160 } 1161 PetscFunctionReturn(PETSC_SUCCESS); 1162 } 1163 1164 /*@ 1165 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation 1166 1167 Collective 1168 1169 Input Parameters: 1170 + dm - The interpolated `DMPLEX` 1171 - pointSF - The initial `PetscSF` without interpolated points 1172 1173 Level: developer 1174 1175 Note: 1176 Debugging for this process can be turned on with the options: `-dm_interp_pre_view` `-petscsf_interp_pre_view` `-petscsection_interp_candidate_view` `-petscsection_interp_candidate_remote_view` `-petscsection_interp_claim_view` `-petscsf_interp_pre_view` `-dmplex_interp_debug` 1177 1178 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()` 1179 @*/ 1180 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 1181 { 1182 MPI_Comm comm; 1183 PetscHMapIJ remoteHash; 1184 PetscHMapI claimshash; 1185 PetscSection candidateSection, candidateRemoteSection, claimSection; 1186 PetscSFNode *candidates, *candidatesRemote, *claims; 1187 const PetscInt *localPoints, *rootdegree; 1188 const PetscSFNode *remotePoints; 1189 PetscInt ov, Nr, r, Nl, l; 1190 PetscInt candidatesSize, candidatesRemoteSize, claimsSize; 1191 PetscBool flg, debug = PETSC_FALSE; 1192 PetscMPIInt rank; 1193 1194 PetscFunctionBegin; 1195 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1196 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 2); 1197 PetscCall(DMPlexIsDistributed(dm, &flg)); 1198 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 1199 /* Set initial SF so that lower level queries work */ 1200 PetscCall(DMSetPointSF(dm, pointSF)); 1201 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1202 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1203 PetscCall(DMPlexGetOverlap(dm, &ov)); 1204 PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1205 PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug)); 1206 PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view")); 1207 PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view")); 1208 PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0)); 1209 /* Step 0: Precalculations */ 1210 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints)); 1211 PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 1212 PetscCall(PetscHMapIJCreate(&remoteHash)); 1213 for (l = 0; l < Nl; ++l) { 1214 PetscHashIJKey key; 1215 1216 key.i = remotePoints[l].index; 1217 key.j = remotePoints[l].rank; 1218 PetscCall(PetscHMapIJSet(remoteHash, key, l)); 1219 } 1220 /* Compute root degree to identify shared points */ 1221 PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree)); 1222 PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree)); 1223 PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree)); 1224 /* 1225 1) Loop over each leaf point $p$ at depth $d$ in the SF 1226 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which 1227 \begin{itemize} 1228 \item all cone points of $f$ are shared 1229 \item $p$ is the cone point with smallest canonical number 1230 \end{itemize} 1231 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$ 1232 \item At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared \label{alg:rootStep} and choose the root face 1233 \item Send the root face from the root back to all leaf process 1234 \item Leaf processes add the shared face to the SF 1235 */ 1236 /* Step 1: Construct section+SFNode array 1237 The section has entries for all shared faces for which we have a leaf point in the cone 1238 The array holds candidate shared faces, each face is referred to by the leaf point */ 1239 PetscCall(PetscSectionCreate(comm, &candidateSection)); 1240 PetscCall(PetscSectionSetChart(candidateSection, 0, Nr)); 1241 { 1242 PetscHMapIJ faceHash; 1243 1244 PetscCall(PetscHMapIJCreate(&faceHash)); 1245 for (l = 0; l < Nl; ++l) { 1246 const PetscInt p = localPoints[l]; 1247 1248 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %" PetscInt_FMT "\n", rank, p)); 1249 PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug)); 1250 } 1251 PetscCall(PetscHMapIJClear(faceHash)); 1252 PetscCall(PetscSectionSetUp(candidateSection)); 1253 PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize)); 1254 PetscCall(PetscMalloc1(candidatesSize, &candidates)); 1255 for (l = 0; l < Nl; ++l) { 1256 const PetscInt p = localPoints[l]; 1257 1258 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %" PetscInt_FMT "\n", rank, p)); 1259 PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug)); 1260 } 1261 PetscCall(PetscHMapIJDestroy(&faceHash)); 1262 if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL)); 1263 } 1264 PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section")); 1265 PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view")); 1266 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates)); 1267 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1268 /* Note that this section is indexed by offsets into leaves, not by point number */ 1269 { 1270 PetscSF sfMulti, sfInverse, sfCandidates; 1271 PetscInt *remoteOffsets; 1272 1273 PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti)); 1274 PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse)); 1275 PetscCall(PetscSectionCreate(comm, &candidateRemoteSection)); 1276 PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection)); 1277 PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates)); 1278 PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize)); 1279 PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote)); 1280 PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_SF_NODE, candidates, candidatesRemote, MPI_REPLACE)); 1281 PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_SF_NODE, candidates, candidatesRemote, MPI_REPLACE)); 1282 PetscCall(PetscSFDestroy(&sfInverse)); 1283 PetscCall(PetscSFDestroy(&sfCandidates)); 1284 PetscCall(PetscFree(remoteOffsets)); 1285 1286 PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section")); 1287 PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view")); 1288 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote)); 1289 } 1290 /* Step 3: At the root, if at least two faces with a given cone are present, including a local face, mark the face as shared and choose the root face */ 1291 { 1292 PetscHMapIJKLRemote faceTable; 1293 PetscInt idx, idx2; 1294 1295 PetscCall(PetscHMapIJKLRemoteCreate(&faceTable)); 1296 /* There is a section point for every leaf attached to a given root point */ 1297 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) { 1298 PetscInt deg; 1299 1300 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1301 PetscInt offset, dof, d; 1302 1303 PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof)); 1304 PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset)); 1305 /* dof may include many faces from the remote process */ 1306 for (d = 0; d < dof; ++d) { 1307 const PetscInt hidx = offset + d; 1308 const PetscInt Np = candidatesRemote[hidx].index + 1; 1309 const PetscSFNode rface = candidatesRemote[hidx + 1]; 1310 const PetscSFNode *fcone = &candidatesRemote[hidx + 2]; 1311 PetscSFNode fcp0; 1312 const PetscSFNode pmax = {-1, -1}; 1313 const PetscInt *join = NULL; 1314 PetscHMapIJKLRemoteKey key; 1315 PetscHashIter iter; 1316 PetscBool missing, mapToLocalPointFailed = PETSC_FALSE; 1317 PetscInt points[1024], p, joinSize; 1318 1319 if (debug) 1320 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking face (%" PetscInt_FMT "), %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with cone size %" PetscInt_FMT "\n", rank, rface.rank, 1321 rface.index, r, idx, d, Np)); 1322 1323 PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle face (%" PetscInt_FMT "), %" PetscInt_FMT ") at (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") with %" PetscInt_FMT " cone points", rface.rank, rface.index, r, idx, d, Np); 1324 fcp0.rank = rank; 1325 fcp0.index = r; 1326 d += Np; 1327 /* Put remote face in hash table */ 1328 key.i = fcp0; 1329 key.j = fcone[0]; 1330 key.k = Np > 2 ? fcone[1] : pmax; 1331 key.l = Np > 3 ? fcone[2] : pmax; 1332 PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key)); 1333 PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing)); 1334 if (missing) { 1335 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT "))\n", rank, rface.index, rface.rank)); 1336 PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface)); 1337 } else { 1338 PetscSFNode oface; 1339 1340 PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface)); 1341 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1342 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank)); 1343 PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface)); 1344 } 1345 } 1346 /* Check for local face */ 1347 points[0] = r; 1348 for (p = 1; p < Np; ++p) { 1349 PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed)); 1350 if (mapToLocalPointFailed) break; /* We got a point not in our overlap */ 1351 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking local candidate %" PetscInt_FMT "\n", rank, points[p])); 1352 } 1353 if (mapToLocalPointFailed) continue; 1354 PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join)); 1355 if (joinSize == 1) { 1356 PetscSFNode lface; 1357 PetscSFNode oface; 1358 1359 /* Always replace with local face */ 1360 lface.rank = rank; 1361 lface.index = join[0]; 1362 PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface)); 1363 if (debug) 1364 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing (%" PetscInt_FMT ", %" PetscInt_FMT ") with local face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, oface.index, oface.rank, lface.index, lface.rank)); 1365 PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface)); 1366 } 1367 PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join)); 1368 } 1369 } 1370 /* Put back faces for this root */ 1371 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) { 1372 PetscInt offset, dof, d; 1373 1374 PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof)); 1375 PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset)); 1376 /* dof may include many faces from the remote process */ 1377 for (d = 0; d < dof; ++d) { 1378 const PetscInt hidx = offset + d; 1379 const PetscInt Np = candidatesRemote[hidx].index + 1; 1380 const PetscSFNode *fcone = &candidatesRemote[hidx + 2]; 1381 PetscSFNode fcp0; 1382 const PetscSFNode pmax = {-1, -1}; 1383 PetscHMapIJKLRemoteKey key; 1384 PetscHashIter iter; 1385 PetscBool missing; 1386 1387 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx)); 1388 PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np); 1389 fcp0.rank = rank; 1390 fcp0.index = r; 1391 d += Np; 1392 /* Find remote face in hash table */ 1393 key.i = fcp0; 1394 key.j = fcone[0]; 1395 key.k = Np > 2 ? fcone[1] : pmax; 1396 key.l = Np > 3 ? fcone[2] : pmax; 1397 PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key)); 1398 if (debug) 1399 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] key (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ") (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, 1400 key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index)); 1401 PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing)); 1402 PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2); 1403 PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx])); 1404 } 1405 } 1406 } 1407 if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL)); 1408 PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable)); 1409 } 1410 /* Step 4: Push back owned faces */ 1411 { 1412 PetscSF sfMulti, sfClaims, sfPointNew; 1413 PetscSFNode *remotePointsNew; 1414 PetscInt *remoteOffsets, *localPointsNew; 1415 PetscInt pStart, pEnd, r, NlNew, p; 1416 1417 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1418 PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti)); 1419 PetscCall(PetscSectionCreate(comm, &claimSection)); 1420 PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection)); 1421 PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims)); 1422 PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize)); 1423 PetscCall(PetscMalloc1(claimsSize, &claims)); 1424 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1425 PetscCall(PetscSFBcastBegin(sfClaims, MPIU_SF_NODE, candidatesRemote, claims, MPI_REPLACE)); 1426 PetscCall(PetscSFBcastEnd(sfClaims, MPIU_SF_NODE, candidatesRemote, claims, MPI_REPLACE)); 1427 PetscCall(PetscSFDestroy(&sfClaims)); 1428 PetscCall(PetscFree(remoteOffsets)); 1429 PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section")); 1430 PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view")); 1431 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims)); 1432 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */ 1433 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */ 1434 PetscCall(PetscHMapICreate(&claimshash)); 1435 for (r = 0; r < Nr; ++r) { 1436 PetscInt dof, off, d; 1437 1438 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %" PetscInt_FMT "\n", rank, r)); 1439 PetscCall(PetscSectionGetDof(candidateSection, r, &dof)); 1440 PetscCall(PetscSectionGetOffset(candidateSection, r, &off)); 1441 for (d = 0; d < dof;) { 1442 if (claims[off + d].rank >= 0) { 1443 const PetscInt faceInd = off + d; 1444 const PetscInt Np = candidates[off + d].index; 1445 const PetscInt *join = NULL; 1446 PetscInt joinSize, points[1024], c; 1447 1448 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index)); 1449 points[0] = r; 1450 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[0])); 1451 for (c = 0, d += 2; c < Np; ++c, ++d) { 1452 PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL)); 1453 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[c + 1])); 1454 } 1455 PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join)); 1456 if (joinSize == 1) { 1457 if (claims[faceInd].rank == rank) { 1458 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0])); 1459 } else { 1460 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found local face %" PetscInt_FMT "\n", rank, join[0])); 1461 PetscCall(PetscHMapISet(claimshash, join[0], faceInd)); 1462 } 1463 } else { 1464 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank)); 1465 } 1466 PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join)); 1467 } else { 1468 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] No claim for point %" PetscInt_FMT "\n", rank, r)); 1469 d += claims[off + d].index + 1; 1470 } 1471 } 1472 } 1473 if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL)); 1474 /* Step 6) Create new pointSF from hashed claims */ 1475 PetscCall(PetscHMapIGetSize(claimshash, &NlNew)); 1476 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1477 PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew)); 1478 PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew)); 1479 for (l = 0; l < Nl; ++l) { 1480 localPointsNew[l] = localPoints[l]; 1481 remotePointsNew[l].index = remotePoints[l].index; 1482 remotePointsNew[l].rank = remotePoints[l].rank; 1483 } 1484 p = Nl; 1485 PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew)); 1486 /* We sort new points, and assume they are numbered after all existing points */ 1487 PetscCall(PetscSortInt(NlNew, PetscSafePointerPlusOffset(localPointsNew, Nl))); 1488 for (p = Nl; p < Nl + NlNew; ++p) { 1489 PetscInt off; 1490 PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off)); 1491 PetscCheck(claims[off].rank >= 0 && claims[off].index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid claim for local point %" PetscInt_FMT ", (%" PetscInt_FMT ", %" PetscInt_FMT ")", localPointsNew[p], claims[off].rank, claims[off].index); 1492 remotePointsNew[p] = claims[off]; 1493 } 1494 PetscCall(PetscSFCreate(comm, &sfPointNew)); 1495 PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1496 PetscCall(PetscSFSetUp(sfPointNew)); 1497 PetscCall(DMSetPointSF(dm, sfPointNew)); 1498 PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view")); 1499 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE)); 1500 PetscCall(PetscSFDestroy(&sfPointNew)); 1501 PetscCall(PetscHMapIDestroy(&claimshash)); 1502 } 1503 PetscCall(PetscHMapIJDestroy(&remoteHash)); 1504 PetscCall(PetscSectionDestroy(&candidateSection)); 1505 PetscCall(PetscSectionDestroy(&candidateRemoteSection)); 1506 PetscCall(PetscSectionDestroy(&claimSection)); 1507 PetscCall(PetscFree(candidates)); 1508 PetscCall(PetscFree(candidatesRemote)); 1509 PetscCall(PetscFree(claims)); 1510 PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0)); 1511 PetscFunctionReturn(PETSC_SUCCESS); 1512 } 1513 1514 /*@ 1515 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1516 1517 Collective 1518 1519 Input Parameter: 1520 . dm - The `DMPLEX` object with only cells and vertices 1521 1522 Output Parameter: 1523 . dmInt - The complete `DMPLEX` object 1524 1525 Level: intermediate 1526 1527 Note: 1528 Labels and coordinates are copied. 1529 1530 Developer Notes: 1531 It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`. 1532 1533 .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()` 1534 @*/ 1535 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1536 { 1537 DMPlexInterpolatedFlag interpolated; 1538 DM idm, odm = dm; 1539 PetscSF sfPoint; 1540 PetscInt depth, dim, d; 1541 const char *name; 1542 PetscBool flg = PETSC_TRUE; 1543 1544 PetscFunctionBegin; 1545 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1546 PetscAssertPointer(dmInt, 2); 1547 PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0)); 1548 PetscCall(DMPlexGetDepth(dm, &depth)); 1549 PetscCall(DMGetDimension(dm, &dim)); 1550 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1551 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1552 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1553 PetscCall(PetscObjectReference((PetscObject)dm)); 1554 idm = dm; 1555 } else { 1556 PetscBool nonmanifold = PETSC_FALSE; 1557 1558 PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_stratify_celltype", &nonmanifold, NULL)); 1559 if (nonmanifold) { 1560 do { 1561 const char *prefix; 1562 PetscInt pStart, pEnd, pdepth; 1563 PetscBool done = PETSC_TRUE; 1564 1565 // Find a point which is not correctly interpolated 1566 PetscCall(DMPlexGetChart(odm, &pStart, &pEnd)); 1567 for (PetscInt p = pStart; p < pEnd; ++p) { 1568 DMPolytopeType ct; 1569 const PetscInt *cone; 1570 PetscInt coneSize, cdepth; 1571 1572 PetscCall(DMPlexGetPointDepth(odm, p, &pdepth)); 1573 PetscCall(DMPlexGetCellType(odm, p, &ct)); 1574 // Check against celltype 1575 if (pdepth != DMPolytopeTypeGetDim(ct)) { 1576 done = PETSC_FALSE; 1577 break; 1578 } 1579 // Check against boundary 1580 PetscCall(DMPlexGetCone(odm, p, &cone)); 1581 PetscCall(DMPlexGetConeSize(odm, p, &coneSize)); 1582 for (PetscInt c = 0; c < coneSize; ++c) { 1583 PetscCall(DMPlexGetPointDepth(odm, cone[c], &cdepth)); 1584 if (cdepth != pdepth - 1) { 1585 done = PETSC_FALSE; 1586 p = pEnd; 1587 break; 1588 } 1589 } 1590 } 1591 if (done) break; 1592 /* Create interpolated mesh */ 1593 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm)); 1594 PetscCall(DMSetType(idm, DMPLEX)); 1595 PetscCall(DMSetDimension(idm, dim)); 1596 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 1597 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix)); 1598 if (depth > 0) { 1599 PetscCall(DMPlexInterpolateFaces_Internal(odm, pdepth, idm)); 1600 PetscCall(DMGetPointSF(odm, &sfPoint)); 1601 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE)); 1602 { 1603 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1604 PetscInt nroots; 1605 PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1606 if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint)); 1607 } 1608 } 1609 if (odm != dm) PetscCall(DMDestroy(&odm)); 1610 odm = idm; 1611 } while (1); 1612 } else { 1613 for (d = 1; d < dim; ++d) { 1614 const char *prefix; 1615 1616 /* Create interpolated mesh */ 1617 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm)); 1618 PetscCall(DMSetType(idm, DMPLEX)); 1619 PetscCall(DMSetDimension(idm, dim)); 1620 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 1621 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix)); 1622 if (depth > 0) { 1623 PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm)); 1624 PetscCall(DMGetPointSF(odm, &sfPoint)); 1625 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE)); 1626 { 1627 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1628 PetscInt nroots; 1629 PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1630 if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint)); 1631 } 1632 } 1633 if (odm != dm) PetscCall(DMDestroy(&odm)); 1634 odm = idm; 1635 } 1636 } 1637 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 1638 PetscCall(PetscObjectSetName((PetscObject)idm, name)); 1639 PetscCall(DMPlexCopyCoordinates(dm, idm)); 1640 PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1641 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL)); 1642 if (flg) PetscCall(DMPlexOrientInterface_Internal(idm)); 1643 } 1644 /* This function makes the mesh fully interpolated on all ranks */ 1645 { 1646 DM_Plex *plex = (DM_Plex *)idm->data; 1647 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1648 } 1649 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm)); 1650 *dmInt = idm; 1651 PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0)); 1652 PetscFunctionReturn(PETSC_SUCCESS); 1653 } 1654 1655 /*@ 1656 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1657 1658 Collective 1659 1660 Input Parameter: 1661 . dmA - The `DMPLEX` object with initial coordinates 1662 1663 Output Parameter: 1664 . dmB - The `DMPLEX` object with copied coordinates 1665 1666 Level: intermediate 1667 1668 Notes: 1669 This is typically used when adding pieces other than vertices to a mesh 1670 1671 This function does not copy localized coordinates. 1672 1673 .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()` 1674 @*/ 1675 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1676 { 1677 Vec coordinatesA, coordinatesB; 1678 VecType vtype; 1679 PetscSection coordSectionA, coordSectionB; 1680 PetscScalar *coordsA, *coordsB; 1681 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1682 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim; 1683 1684 PetscFunctionBegin; 1685 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1686 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1687 if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS); 1688 PetscCall(DMGetCoordinateDim(dmA, &cdim)); 1689 PetscCall(DMSetCoordinateDim(dmB, cdim)); 1690 PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA)); 1691 PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB)); 1692 PetscCheck((vEndA - vStartA) == (vEndB - vStartB), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of vertices in first DM %" PetscInt_FMT " != %" PetscInt_FMT " in the second DM", vEndA - vStartA, vEndB - vStartB); 1693 /* Copy over discretization if it exists */ 1694 { 1695 DM cdmA, cdmB; 1696 PetscDS dsA, dsB; 1697 PetscObject objA, objB; 1698 PetscClassId idA, idB; 1699 const PetscScalar *constants; 1700 PetscInt cdim, Nc; 1701 1702 PetscCall(DMGetCoordinateDM(dmA, &cdmA)); 1703 PetscCall(DMGetCoordinateDM(dmB, &cdmB)); 1704 PetscCall(DMGetField(cdmA, 0, NULL, &objA)); 1705 PetscCall(DMGetField(cdmB, 0, NULL, &objB)); 1706 PetscCall(PetscObjectGetClassId(objA, &idA)); 1707 PetscCall(PetscObjectGetClassId(objB, &idB)); 1708 if ((idA == PETSCFE_CLASSID) && (idA != idB)) { 1709 PetscCall(DMSetField(cdmB, 0, NULL, objA)); 1710 PetscCall(DMCreateDS(cdmB)); 1711 PetscCall(DMGetDS(cdmA, &dsA)); 1712 PetscCall(DMGetDS(cdmB, &dsB)); 1713 PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim)); 1714 PetscCall(PetscDSSetCoordinateDimension(dsB, cdim)); 1715 PetscCall(PetscDSGetConstants(dsA, &Nc, &constants)); 1716 PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants)); 1717 } 1718 } 1719 PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA)); 1720 PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB)); 1721 PetscCall(DMGetCoordinateSection(dmA, &coordSectionA)); 1722 PetscCall(DMGetCoordinateSection(dmB, &coordSectionB)); 1723 if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS); 1724 PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf)); 1725 if (!Nf) PetscFunctionReturn(PETSC_SUCCESS); 1726 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf); 1727 if (!coordSectionB) { 1728 PetscInt dim; 1729 1730 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB)); 1731 PetscCall(DMGetCoordinateDim(dmA, &dim)); 1732 PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB)); 1733 PetscCall(PetscObjectDereference((PetscObject)coordSectionB)); 1734 } 1735 PetscCall(PetscSectionSetNumFields(coordSectionB, 1)); 1736 PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim)); 1737 PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim)); 1738 PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE)); 1739 cS = vStartB; 1740 cE = vEndB; 1741 PetscCall(PetscSectionSetChart(coordSectionB, cS, cE)); 1742 for (v = vStartB; v < vEndB; ++v) { 1743 PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim)); 1744 PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim)); 1745 } 1746 PetscCall(PetscSectionSetUp(coordSectionB)); 1747 PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB)); 1748 PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA)); 1749 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB)); 1750 PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates")); 1751 PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE)); 1752 PetscCall(VecGetBlockSize(coordinatesA, &d)); 1753 PetscCall(VecSetBlockSize(coordinatesB, d)); 1754 PetscCall(VecGetType(coordinatesA, &vtype)); 1755 PetscCall(VecSetType(coordinatesB, vtype)); 1756 PetscCall(VecGetArray(coordinatesA, &coordsA)); 1757 PetscCall(VecGetArray(coordinatesB, &coordsB)); 1758 for (v = 0; v < vEndB - vStartB; ++v) { 1759 PetscInt offA, offB; 1760 1761 PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA)); 1762 PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB)); 1763 for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d]; 1764 } 1765 PetscCall(VecRestoreArray(coordinatesA, &coordsA)); 1766 PetscCall(VecRestoreArray(coordinatesB, &coordsB)); 1767 PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB)); 1768 PetscCall(VecDestroy(&coordinatesB)); 1769 PetscFunctionReturn(PETSC_SUCCESS); 1770 } 1771 1772 /*@ 1773 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1774 1775 Collective 1776 1777 Input Parameter: 1778 . dm - The complete `DMPLEX` object 1779 1780 Output Parameter: 1781 . dmUnint - The `DMPLEX` object with only cells and vertices 1782 1783 Level: intermediate 1784 1785 Note: 1786 It does not copy over the coordinates. 1787 1788 Developer Notes: 1789 Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`. 1790 1791 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()` 1792 @*/ 1793 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1794 { 1795 DMPlexInterpolatedFlag interpolated; 1796 DM udm; 1797 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 1798 1799 PetscFunctionBegin; 1800 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1801 PetscAssertPointer(dmUnint, 2); 1802 PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0)); 1803 PetscCall(DMGetDimension(dm, &dim)); 1804 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1805 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1806 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1807 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1808 PetscCall(PetscObjectReference((PetscObject)dm)); 1809 *dmUnint = dm; 1810 PetscFunctionReturn(PETSC_SUCCESS); 1811 } 1812 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1813 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1814 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm)); 1815 PetscCall(DMSetType(udm, DMPLEX)); 1816 PetscCall(DMSetDimension(udm, dim)); 1817 PetscCall(DMPlexSetChart(udm, cStart, vEnd)); 1818 for (c = cStart; c < cEnd; ++c) { 1819 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1820 1821 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1822 for (cl = 0; cl < closureSize * 2; cl += 2) { 1823 const PetscInt p = closure[cl]; 1824 1825 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1826 } 1827 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1828 PetscCall(DMPlexSetConeSize(udm, c, coneSize)); 1829 maxConeSize = PetscMax(maxConeSize, coneSize); 1830 } 1831 PetscCall(DMSetUp(udm)); 1832 PetscCall(PetscMalloc1(maxConeSize, &cone)); 1833 for (c = cStart; c < cEnd; ++c) { 1834 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1835 1836 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1837 for (cl = 0; cl < closureSize * 2; cl += 2) { 1838 const PetscInt p = closure[cl]; 1839 1840 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1841 } 1842 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1843 PetscCall(DMPlexSetCone(udm, c, cone)); 1844 } 1845 PetscCall(PetscFree(cone)); 1846 PetscCall(DMPlexSymmetrize(udm)); 1847 PetscCall(DMPlexStratify(udm)); 1848 /* Reduce SF */ 1849 { 1850 PetscSF sfPoint, sfPointUn; 1851 const PetscSFNode *remotePoints; 1852 const PetscInt *localPoints; 1853 PetscSFNode *remotePointsUn; 1854 PetscInt *localPointsUn; 1855 PetscInt numRoots, numLeaves, l; 1856 PetscInt numLeavesUn = 0, n = 0; 1857 1858 /* Get original SF information */ 1859 PetscCall(DMGetPointSF(dm, &sfPoint)); 1860 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE)); 1861 PetscCall(DMGetPointSF(udm, &sfPointUn)); 1862 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1863 if (numRoots >= 0) { 1864 /* Allocate space for cells and vertices */ 1865 for (l = 0; l < numLeaves; ++l) { 1866 const PetscInt p = localPoints[l]; 1867 1868 if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++; 1869 } 1870 /* Fill in leaves */ 1871 PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn)); 1872 PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn)); 1873 for (l = 0; l < numLeaves; l++) { 1874 const PetscInt p = localPoints[l]; 1875 1876 if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) { 1877 localPointsUn[n] = p; 1878 remotePointsUn[n].rank = remotePoints[l].rank; 1879 remotePointsUn[n].index = remotePoints[l].index; 1880 ++n; 1881 } 1882 } 1883 PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn); 1884 PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER)); 1885 } 1886 } 1887 /* This function makes the mesh fully uninterpolated on all ranks */ 1888 { 1889 DM_Plex *plex = (DM_Plex *)udm->data; 1890 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1891 } 1892 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm)); 1893 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE)); 1894 *dmUnint = udm; 1895 PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0)); 1896 PetscFunctionReturn(PETSC_SUCCESS); 1897 } 1898 1899 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1900 { 1901 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1902 MPI_Comm comm; 1903 1904 PetscFunctionBegin; 1905 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1906 PetscCall(DMPlexGetDepth(dm, &depth)); 1907 PetscCall(DMGetDimension(dm, &dim)); 1908 1909 if (depth == dim) { 1910 *interpolated = DMPLEX_INTERPOLATED_FULL; 1911 if (!dim) goto finish; 1912 1913 /* Check points at height = dim are vertices (have no cones) */ 1914 PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd)); 1915 for (p = pStart; p < pEnd; p++) { 1916 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1917 if (coneSize) { 1918 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1919 goto finish; 1920 } 1921 } 1922 1923 /* Check points at height < dim have cones */ 1924 for (h = 0; h < dim; h++) { 1925 PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 1926 for (p = pStart; p < pEnd; p++) { 1927 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1928 if (!coneSize) { 1929 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1930 goto finish; 1931 } 1932 } 1933 } 1934 } else if (depth == 1) { 1935 *interpolated = DMPLEX_INTERPOLATED_NONE; 1936 } else { 1937 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1938 } 1939 finish: 1940 PetscFunctionReturn(PETSC_SUCCESS); 1941 } 1942 1943 /*@ 1944 DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated. 1945 1946 Not Collective 1947 1948 Input Parameter: 1949 . dm - The `DMPLEX` object 1950 1951 Output Parameter: 1952 . interpolated - Flag whether the `DM` is interpolated 1953 1954 Level: intermediate 1955 1956 Notes: 1957 Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective 1958 so the results can be different on different ranks in special cases. 1959 However, `DMPlexInterpolate()` guarantees the result is the same on all. 1960 1961 Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`. 1962 1963 Developer Notes: 1964 Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`. 1965 1966 If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called. 1967 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1968 `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`. 1969 1970 If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated. 1971 1972 `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`, 1973 and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`. 1974 1975 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()` 1976 @*/ 1977 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1978 { 1979 DM_Plex *plex = (DM_Plex *)dm->data; 1980 1981 PetscFunctionBegin; 1982 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1983 PetscAssertPointer(interpolated, 2); 1984 if (plex->interpolated < 0) { 1985 PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated)); 1986 } else if (PetscDefined(USE_DEBUG)) { 1987 DMPlexInterpolatedFlag flg; 1988 1989 PetscCall(DMPlexIsInterpolated_Internal(dm, &flg)); 1990 PetscCheck(plex->tr || flg == plex->interpolated, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]); 1991 } 1992 *interpolated = plex->interpolated; 1993 PetscFunctionReturn(PETSC_SUCCESS); 1994 } 1995 1996 /*@ 1997 DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner). 1998 1999 Collective 2000 2001 Input Parameter: 2002 . dm - The `DMPLEX` object 2003 2004 Output Parameter: 2005 . interpolated - Flag whether the `DM` is interpolated 2006 2007 Level: intermediate 2008 2009 Notes: 2010 Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks. 2011 2012 This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks. 2013 2014 Developer Notes: 2015 Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`. 2016 2017 If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated. 2018 `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 2019 if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`, 2020 otherwise sets plex->interpolatedCollective = plex->interpolated. 2021 2022 If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective. 2023 2024 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()` 2025 @*/ 2026 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 2027 { 2028 DM_Plex *plex = (DM_Plex *)dm->data; 2029 PetscBool debug = PETSC_FALSE; 2030 2031 PetscFunctionBegin; 2032 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2033 PetscAssertPointer(interpolated, 2); 2034 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL)); 2035 if (plex->interpolatedCollective < 0) { 2036 DMPlexInterpolatedFlag min, max; 2037 MPI_Comm comm; 2038 2039 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2040 PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective)); 2041 PetscCallMPI(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm)); 2042 PetscCallMPI(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm)); 2043 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 2044 if (debug) { 2045 PetscMPIInt rank; 2046 2047 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2048 PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective])); 2049 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 2050 } 2051 } 2052 *interpolated = plex->interpolatedCollective; 2053 PetscFunctionReturn(PETSC_SUCCESS); 2054 } 2055 2056 /*@ 2057 DMPlexGetInterpolatePreferTensor - Get the flag to prefer tensor order when interpolating a cell 2058 2059 Not Collective 2060 2061 Input Parameter: 2062 . dm - The `DMPLEX` object 2063 2064 Output Parameter: 2065 . preferTensor - Flag to prefer tensor order 2066 2067 Level: intermediate 2068 2069 Notes: 2070 Currently, this just affects triangular prisms, 2071 2072 .seealso: `DMPlexSetInterpolatePreferTensor()`, `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()` 2073 @*/ 2074 PetscErrorCode DMPlexGetInterpolatePreferTensor(DM dm, PetscBool *preferTensor) 2075 { 2076 DM_Plex *plex = (DM_Plex *)dm->data; 2077 2078 PetscFunctionBegin; 2079 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2080 PetscAssertPointer(preferTensor, 2); 2081 *preferTensor = plex->interpolatePreferTensor; 2082 PetscFunctionReturn(PETSC_SUCCESS); 2083 } 2084 2085 /*@ 2086 DMPlexSetInterpolatePreferTensor - Set the flag to prefer tensor order when interpolating a cell 2087 2088 Logically Collective 2089 2090 Input Parameters: 2091 + dm - The `DMPLEX` object 2092 - preferTensor - Flag to prefer tensor order 2093 2094 Level: intermediate 2095 2096 Notes: 2097 Currently, this just affects triangular prisms, 2098 2099 .seealso: `DMPlexGetInterpolatePreferTensor()`, `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()` 2100 @*/ 2101 PetscErrorCode DMPlexSetInterpolatePreferTensor(DM dm, PetscBool preferTensor) 2102 { 2103 DM_Plex *plex = (DM_Plex *)dm->data; 2104 2105 PetscFunctionBegin; 2106 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2107 PetscValidLogicalCollectiveBool(dm, preferTensor, 2); 2108 plex->interpolatePreferTensor = preferTensor; 2109 PetscFunctionReturn(PETSC_SUCCESS); 2110 } 2111