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 static 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, d, pStart, Np, cStart, cEnd, c, fStart, fEnd; 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, cellDepth, &cStart, &cEnd)); 506 /* Number new faces and save face vertices in hash table */ 507 PetscCall(DMPlexGetDepthStratum(dm, depth > cellDepth ? cellDepth : 0, NULL, &fStart)); 508 { 509 PetscInt opEnd; 510 511 // We need to account for existing faces in non-manifold meshes 512 PetscCall(DMPlexGetChart(dm, NULL, &opEnd)); 513 fStart = PetscMax(opEnd, fStart); 514 } 515 fEnd = fStart; 516 517 minCone = PETSC_MAX_INT; 518 for (c = cStart, cntFaces = 0; c < cEnd; ++c) { 519 const PetscInt *cone; 520 DMPolytopeType ct; 521 PetscInt numFaces = 0, coneSize; 522 523 PetscCall(DMPlexGetCellType(dm, c, &ct)); 524 PetscCall(DMPlexGetCone(dm, c, &cone)); 525 PetscCall(DMPlexGetConeSize(dm, c, &coneSize)); 526 for (PetscInt j = 0; j < coneSize; j++) minCone = PetscMin(cone[j], minCone); 527 // Ignore faces since they are interpolated 528 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, NULL, NULL, NULL)); 529 cntFaces += numFaces; 530 } 531 // Encode so that we can use 0 as an excluded value, instead of PETSC_MAX_INT 532 minCone = -(minCone - 1); 533 534 PetscCall(PetscMalloc1(cntFaces, &facesId)); 535 536 for (c = cStart, cntFaces = 0; c < cEnd; ++c) { 537 const PetscInt *cone, *faceSizes, *faces; 538 const DMPolytopeType *faceTypes; 539 DMPolytopeType ct; 540 PetscInt numFaces, cf, 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 (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 for (c = cStart, cntFaces = 0; c < cEnd; ++c) { 587 const PetscInt *cone, *faceSizes, *faces; 588 const DMPolytopeType *faceTypes; 589 DMPolytopeType ct; 590 PetscInt numFaces, cf, foff = 0; 591 592 PetscCall(DMPlexGetCellType(dm, c, &ct)); 593 PetscCall(DMPlexGetCone(dm, c, &cone)); 594 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) { 595 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 596 } else { 597 numFaces = 0; 598 } 599 for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) { 600 const PetscInt faceSize = faceSizes[cf]; 601 const DMPolytopeType faceType = faceTypes[cf]; 602 const PetscInt *face = &faces[foff]; 603 PetscHashIJKLKey key; 604 PetscHashIter iter; 605 PetscBool missing; 606 607 key.i = face[0] + minCone; 608 key.j = faceSize > 1 ? face[1] + minCone : 0; 609 key.k = faceSize > 2 ? face[2] + minCone : 0; 610 key.l = faceSize > 3 ? face[3] + minCone : 0; 611 PetscCall(PetscSortInt(faceSize, (PetscInt *)&key)); 612 PetscCall(PetscHMapIJKLPut(faceTable, key, &iter, &missing)); 613 if (missing) { 614 facesId[cntFaces] = faceTypeStart[faceType]; 615 PetscCall(PetscHMapIJKLIterSet(faceTable, iter, faceTypeStart[faceType]++)); 616 } else PetscCall(PetscHMapIJKLIterGet(faceTable, iter, &facesId[cntFaces])); 617 cntFaces++; 618 } 619 if (ct != DM_POLYTOPE_SEGMENT && ct != DM_POLYTOPE_POINT_PRISM_TENSOR) PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 620 } 621 for (ct = 1; ct < DM_NUM_POLYTOPES; ++ct) { 622 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]); 623 } 624 } 625 } 626 PetscCall(PetscHMapIJKLDestroy(&faceTable)); 627 628 /* Add new points, always at the end of the numbering */ 629 PetscCall(DMPlexGetChart(dm, &pStart, &Np)); 630 PetscCall(DMPlexSetChart(idm, pStart, Np + (fEnd - fStart))); 631 /* Set cone sizes */ 632 /* Must create the celltype label here so that we do not automatically try to compute the types */ 633 PetscCall(DMCreateLabel(idm, "celltype")); 634 PetscCall(DMPlexGetCellTypeLabel(idm, &ctLabel)); 635 for (d = 0; d <= depth; ++d) { 636 DMPolytopeType ct; 637 PetscInt coneSize, pStart, pEnd, p; 638 639 if (d == cellDepth) continue; 640 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 641 for (p = pStart; p < pEnd; ++p) { 642 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 643 PetscCall(DMPlexSetConeSize(idm, p, coneSize)); 644 PetscCall(DMPlexGetCellType(dm, p, &ct)); 645 PetscCall(DMPlexSetCellType(idm, p, ct)); 646 } 647 } 648 for (c = cStart, cntFaces = 0; c < cEnd; ++c) { 649 const PetscInt *cone, *faceSizes; 650 const DMPolytopeType *faceTypes; 651 DMPolytopeType ct; 652 PetscInt numFaces, cf; 653 654 PetscCall(DMPlexGetCellType(dm, c, &ct)); 655 PetscCall(DMPlexGetCone(dm, c, &cone)); 656 if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) { 657 PetscCall(DMPlexSetCellType(idm, c, ct)); 658 PetscCall(DMPlexSetConeSize(idm, c, 2)); 659 continue; 660 } 661 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL)); 662 PetscCall(DMPlexSetCellType(idm, c, ct)); 663 PetscCall(DMPlexSetConeSize(idm, c, numFaces)); 664 for (cf = 0; cf < numFaces; ++cf) { 665 const PetscInt f = facesId[cntFaces]; 666 DMPolytopeType faceType = faceTypes[cf]; 667 const PetscInt faceSize = faceSizes[cf]; 668 PetscCall(DMPlexSetConeSize(idm, f, faceSize)); 669 PetscCall(DMPlexSetCellType(idm, f, faceType)); 670 cntFaces++; 671 } 672 PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, NULL)); 673 } 674 PetscCall(DMSetUp(idm)); 675 /* Initialize cones so we do not need the bash table to tell us that a cone has been set */ 676 { 677 PetscSection cs; 678 PetscInt *cones, csize; 679 680 PetscCall(DMPlexGetConeSection(idm, &cs)); 681 PetscCall(DMPlexGetCones(idm, &cones)); 682 PetscCall(PetscSectionGetStorageSize(cs, &csize)); 683 for (c = 0; c < csize; ++c) cones[c] = -1; 684 } 685 /* Set cones */ 686 for (d = 0; d <= depth; ++d) { 687 const PetscInt *cone; 688 PetscInt pStart, pEnd, p; 689 690 if (d == cellDepth) continue; 691 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 692 for (p = pStart; p < pEnd; ++p) { 693 PetscCall(DMPlexGetCone(dm, p, &cone)); 694 PetscCall(DMPlexSetCone(idm, p, cone)); 695 PetscCall(DMPlexGetConeOrientation(dm, p, &cone)); 696 PetscCall(DMPlexSetConeOrientation(idm, p, cone)); 697 } 698 } 699 for (c = cStart, cntFaces = 0; c < cEnd; ++c) { 700 const PetscInt *cone, *faceSizes, *faces; 701 const DMPolytopeType *faceTypes; 702 DMPolytopeType ct; 703 PetscInt numFaces, cf, foff = 0; 704 705 PetscCall(DMPlexGetCellType(dm, c, &ct)); 706 PetscCall(DMPlexGetCone(dm, c, &cone)); 707 if (ct == DM_POLYTOPE_SEGMENT || ct == DM_POLYTOPE_POINT_PRISM_TENSOR) { 708 PetscCall(DMPlexSetCone(idm, c, cone)); 709 PetscCall(DMPlexGetConeOrientation(dm, c, &cone)); 710 PetscCall(DMPlexSetConeOrientation(idm, c, cone)); 711 continue; 712 } 713 PetscCall(DMPlexGetRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 714 for (cf = 0; cf < numFaces; foff += faceSizes[cf], ++cf) { 715 DMPolytopeType faceType = faceTypes[cf]; 716 const PetscInt faceSize = faceSizes[cf]; 717 const PetscInt f = facesId[cntFaces]; 718 const PetscInt *face = &faces[foff]; 719 const PetscInt *fcone; 720 721 PetscCall(DMPlexInsertCone(idm, c, cf, f)); 722 PetscCall(DMPlexGetCone(idm, f, &fcone)); 723 if (fcone[0] < 0) PetscCall(DMPlexSetCone(idm, f, face)); 724 { 725 const PetscInt *cone; 726 PetscInt coneSize, ornt; 727 728 PetscCall(DMPlexGetConeSize(idm, f, &coneSize)); 729 PetscCall(DMPlexGetCone(idm, f, &cone)); 730 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); 731 /* Notice that we have to use vertices here because the lower dimensional faces have not been created yet */ 732 PetscCall(DMPolytopeGetVertexOrientation(faceType, cone, face, &ornt)); 733 PetscCall(DMPlexInsertConeOrientation(idm, c, cf, ornt)); 734 } 735 cntFaces++; 736 } 737 PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, cone, &numFaces, &faceTypes, &faceSizes, &faces)); 738 } 739 PetscCall(PetscFree(facesId)); 740 PetscCall(DMPlexSymmetrize(idm)); 741 PetscCall(DMPlexStratify(idm)); 742 PetscFunctionReturn(PETSC_SUCCESS); 743 } 744 745 static PetscErrorCode SortRmineRremoteByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 746 { 747 PetscInt nleaves; 748 PetscInt nranks; 749 const PetscMPIInt *ranks = NULL; 750 const PetscInt *roffset = NULL, *rmine = NULL, *rremote = NULL; 751 PetscInt n, o, r; 752 753 PetscFunctionBegin; 754 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote)); 755 nleaves = roffset[nranks]; 756 PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1)); 757 for (r = 0; r < nranks; r++) { 758 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 759 - to unify order with the other side */ 760 o = roffset[r]; 761 n = roffset[r + 1] - o; 762 PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n)); 763 PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n)); 764 PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o])); 765 } 766 PetscFunctionReturn(PETSC_SUCCESS); 767 } 768 769 PetscErrorCode DMPlexOrientInterface_Internal(DM dm) 770 { 771 PetscSF sf; 772 const PetscInt *locals; 773 const PetscSFNode *remotes; 774 const PetscMPIInt *ranks; 775 const PetscInt *roffset; 776 PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 777 PetscInt nroots, p, nleaves, nranks, r, maxConeSize = 0; 778 PetscInt(*roots)[4], (*leaves)[4], mainCone[4]; 779 PetscMPIInt(*rootsRanks)[4], (*leavesRanks)[4]; 780 MPI_Comm comm; 781 PetscMPIInt rank, size; 782 PetscInt debug = 0; 783 784 PetscFunctionBegin; 785 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 786 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 787 PetscCallMPI(MPI_Comm_size(comm, &size)); 788 PetscCall(DMGetPointSF(dm, &sf)); 789 PetscCall(DMViewFromOptions(dm, NULL, "-before_orient_interface_dm_view")); 790 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sf, PETSC_FALSE)); 791 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &locals, &remotes)); 792 if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS); 793 PetscCall(PetscSFSetUp(sf)); 794 PetscCall(SortRmineRremoteByRemote_Private(sf, &rmine1, &rremote1)); 795 for (p = 0; p < nleaves; ++p) { 796 PetscInt coneSize; 797 PetscCall(DMPlexGetConeSize(dm, locals[p], &coneSize)); 798 maxConeSize = PetscMax(maxConeSize, coneSize); 799 } 800 PetscCheck(maxConeSize <= 4, comm, PETSC_ERR_SUP, "This method does not support cones of size %" PetscInt_FMT, maxConeSize); 801 PetscCall(PetscMalloc4(nroots, &roots, nroots, &leaves, nroots, &rootsRanks, nroots, &leavesRanks)); 802 for (p = 0; p < nroots; ++p) { 803 const PetscInt *cone; 804 PetscInt coneSize, c, ind0; 805 806 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 807 PetscCall(DMPlexGetCone(dm, p, &cone)); 808 /* Ignore vertices */ 809 if (coneSize < 2) { 810 for (c = 0; c < 4; c++) { 811 roots[p][c] = -1; 812 rootsRanks[p][c] = -1; 813 } 814 continue; 815 } 816 /* Translate all points to root numbering */ 817 for (c = 0; c < PetscMin(coneSize, 4); c++) { 818 PetscCall(PetscFindInt(cone[c], nleaves, locals, &ind0)); 819 if (ind0 < 0) { 820 roots[p][c] = cone[c]; 821 rootsRanks[p][c] = rank; 822 } else { 823 roots[p][c] = remotes[ind0].index; 824 rootsRanks[p][c] = remotes[ind0].rank; 825 } 826 } 827 for (c = coneSize; c < 4; c++) { 828 roots[p][c] = -1; 829 rootsRanks[p][c] = -1; 830 } 831 } 832 for (p = 0; p < nroots; ++p) { 833 PetscInt c; 834 for (c = 0; c < 4; c++) { 835 leaves[p][c] = -2; 836 leavesRanks[p][c] = -2; 837 } 838 } 839 PetscCall(PetscSFBcastBegin(sf, MPIU_4INT, roots, leaves, MPI_REPLACE)); 840 PetscCall(PetscSFBcastBegin(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE)); 841 PetscCall(PetscSFBcastEnd(sf, MPIU_4INT, roots, leaves, MPI_REPLACE)); 842 PetscCall(PetscSFBcastEnd(sf, MPI_4INT, rootsRanks, leavesRanks, MPI_REPLACE)); 843 if (debug) { 844 PetscCall(PetscSynchronizedFlush(comm, NULL)); 845 if (rank == 0) PetscCall(PetscSynchronizedPrintf(comm, "Referenced roots\n")); 846 } 847 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL)); 848 for (p = 0; p < nroots; ++p) { 849 DMPolytopeType ct; 850 const PetscInt *cone; 851 PetscInt coneSize, c, ind0, o; 852 853 if (leaves[p][0] < 0) continue; /* Ignore vertices */ 854 PetscCall(DMPlexGetCellType(dm, p, &ct)); 855 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 856 PetscCall(DMPlexGetCone(dm, p, &cone)); 857 if (debug) { 858 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])); 859 } 860 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]) { 861 /* Translate these leaves to my cone points; mainCone means desired order p's cone points */ 862 for (c = 0; c < PetscMin(coneSize, 4); ++c) { 863 PetscInt rS, rN; 864 865 if (leavesRanks[p][c] == rank) { 866 /* A local leaf is just taken as it is */ 867 mainCone[c] = leaves[p][c]; 868 continue; 869 } 870 /* Find index of rank leavesRanks[p][c] among remote ranks */ 871 /* No need for PetscMPIIntCast because these integers were originally cast from PetscMPIInt. */ 872 PetscCall(PetscFindMPIInt((PetscMPIInt)leavesRanks[p][c], nranks, ranks, &r)); 873 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]); 874 PetscCheck(ranks[r] >= 0 && ranks[r] < size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "p=%" PetscInt_FMT " c=%" PetscInt_FMT " commsize=%d: ranks[%" PetscInt_FMT "] = %d makes no sense", p, c, size, r, ranks[r]); 875 /* Find point leaves[p][c] among remote points aimed at rank leavesRanks[p][c] */ 876 rS = roffset[r]; 877 rN = roffset[r + 1] - rS; 878 PetscCall(PetscFindInt(leaves[p][c], rN, &rremote1[rS], &ind0)); 879 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]); 880 /* Get the corresponding local point */ 881 mainCone[c] = rmine1[rS + ind0]; 882 } 883 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])); 884 /* Set the desired order of p's cone points and fix orientations accordingly */ 885 PetscCall(DMPolytopeGetOrientation(ct, cone, mainCone, &o)); 886 PetscCall(DMPlexOrientPoint(dm, p, o)); 887 } else if (debug) PetscCall(PetscSynchronizedPrintf(comm, " ==\n")); 888 } 889 if (debug) { 890 PetscCall(PetscSynchronizedFlush(comm, NULL)); 891 PetscCallMPI(MPI_Barrier(comm)); 892 } 893 PetscCall(DMViewFromOptions(dm, NULL, "-after_orient_interface_dm_view")); 894 PetscCall(PetscFree4(roots, leaves, rootsRanks, leavesRanks)); 895 PetscCall(PetscFree2(rmine1, rremote1)); 896 PetscFunctionReturn(PETSC_SUCCESS); 897 } 898 899 static PetscErrorCode IntArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], const char valname[], PetscInt n, const PetscInt a[]) 900 { 901 PetscInt idx; 902 PetscMPIInt rank; 903 PetscBool flg; 904 905 PetscFunctionBegin; 906 PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg)); 907 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 908 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 909 PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name)); 910 for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s %" PetscInt_FMT " %s %" PetscInt_FMT "\n", rank, idxname, idx, valname, a[idx])); 911 PetscCall(PetscSynchronizedFlush(comm, NULL)); 912 PetscFunctionReturn(PETSC_SUCCESS); 913 } 914 915 static PetscErrorCode SFNodeArrayViewFromOptions(MPI_Comm comm, const char opt[], const char name[], const char idxname[], PetscInt n, const PetscSFNode a[]) 916 { 917 PetscInt idx; 918 PetscMPIInt rank; 919 PetscBool flg; 920 921 PetscFunctionBegin; 922 PetscCall(PetscOptionsHasName(NULL, NULL, opt, &flg)); 923 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 924 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 925 PetscCall(PetscSynchronizedPrintf(comm, "[%d]%s:\n", rank, name)); 926 if (idxname) { 927 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)); 928 } else { 929 for (idx = 0; idx < n; ++idx) PetscCall(PetscSynchronizedPrintf(comm, "[%d]rank %" PetscInt_FMT " index %" PetscInt_FMT "\n", rank, a[idx].rank, a[idx].index)); 930 } 931 PetscCall(PetscSynchronizedFlush(comm, NULL)); 932 PetscFunctionReturn(PETSC_SUCCESS); 933 } 934 935 static PetscErrorCode DMPlexMapToLocalPoint(DM dm, PetscHMapIJ remotehash, PetscSFNode remotePoint, PetscInt *localPoint, PetscBool *mapFailed) 936 { 937 PetscSF sf; 938 const PetscInt *locals; 939 PetscMPIInt rank; 940 941 PetscFunctionBegin; 942 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 943 PetscCall(DMGetPointSF(dm, &sf)); 944 PetscCall(PetscSFGetGraph(sf, NULL, NULL, &locals, NULL)); 945 if (mapFailed) *mapFailed = PETSC_FALSE; 946 if (remotePoint.rank == rank) { 947 *localPoint = remotePoint.index; 948 } else { 949 PetscHashIJKey key; 950 PetscInt l; 951 952 key.i = remotePoint.index; 953 key.j = remotePoint.rank; 954 PetscCall(PetscHMapIJGet(remotehash, key, &l)); 955 if (l >= 0) { 956 *localPoint = locals[l]; 957 } else if (mapFailed) *mapFailed = PETSC_TRUE; 958 } 959 PetscFunctionReturn(PETSC_SUCCESS); 960 } 961 962 static PetscErrorCode DMPlexMapToGlobalPoint(DM dm, PetscInt localPoint, PetscSFNode *remotePoint, PetscBool *mapFailed) 963 { 964 PetscSF sf; 965 const PetscInt *locals, *rootdegree; 966 const PetscSFNode *remotes; 967 PetscInt Nl, l; 968 PetscMPIInt rank; 969 970 PetscFunctionBegin; 971 if (mapFailed) *mapFailed = PETSC_FALSE; 972 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 973 PetscCall(DMGetPointSF(dm, &sf)); 974 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, &remotes)); 975 if (Nl < 0) goto owned; 976 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 977 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 978 if (rootdegree[localPoint]) goto owned; 979 PetscCall(PetscFindInt(localPoint, Nl, locals, &l)); 980 if (l < 0) { 981 if (mapFailed) *mapFailed = PETSC_TRUE; 982 } else *remotePoint = remotes[l]; 983 PetscFunctionReturn(PETSC_SUCCESS); 984 owned: 985 remotePoint->rank = rank; 986 remotePoint->index = localPoint; 987 PetscFunctionReturn(PETSC_SUCCESS); 988 } 989 990 static PetscErrorCode DMPlexPointIsShared(DM dm, PetscInt p, PetscBool *isShared) 991 { 992 PetscSF sf; 993 const PetscInt *locals, *rootdegree; 994 PetscInt Nl, idx; 995 996 PetscFunctionBegin; 997 *isShared = PETSC_FALSE; 998 PetscCall(DMGetPointSF(dm, &sf)); 999 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &locals, NULL)); 1000 if (Nl < 0) PetscFunctionReturn(PETSC_SUCCESS); 1001 PetscCall(PetscFindInt(p, Nl, locals, &idx)); 1002 if (idx >= 0) { 1003 *isShared = PETSC_TRUE; 1004 PetscFunctionReturn(PETSC_SUCCESS); 1005 } 1006 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 1007 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 1008 if (rootdegree[p] > 0) *isShared = PETSC_TRUE; 1009 PetscFunctionReturn(PETSC_SUCCESS); 1010 } 1011 1012 static PetscErrorCode DMPlexConeIsShared(DM dm, PetscInt p, PetscBool *isShared) 1013 { 1014 const PetscInt *cone; 1015 PetscInt coneSize, c; 1016 PetscBool cShared = PETSC_TRUE; 1017 1018 PetscFunctionBegin; 1019 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1020 PetscCall(DMPlexGetCone(dm, p, &cone)); 1021 for (c = 0; c < coneSize; ++c) { 1022 PetscBool pointShared; 1023 1024 PetscCall(DMPlexPointIsShared(dm, cone[c], &pointShared)); 1025 cShared = (PetscBool)(cShared && pointShared); 1026 } 1027 *isShared = coneSize ? cShared : PETSC_FALSE; 1028 PetscFunctionReturn(PETSC_SUCCESS); 1029 } 1030 1031 static PetscErrorCode DMPlexGetConeMinimum(DM dm, PetscInt p, PetscSFNode *cpmin) 1032 { 1033 const PetscInt *cone; 1034 PetscInt coneSize, c; 1035 PetscSFNode cmin = {PETSC_MAX_INT, PETSC_MAX_INT}, missing = {-1, -1}; 1036 1037 PetscFunctionBegin; 1038 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1039 PetscCall(DMPlexGetCone(dm, p, &cone)); 1040 for (c = 0; c < coneSize; ++c) { 1041 PetscSFNode rcp; 1042 PetscBool mapFailed; 1043 1044 PetscCall(DMPlexMapToGlobalPoint(dm, cone[c], &rcp, &mapFailed)); 1045 if (mapFailed) { 1046 cmin = missing; 1047 } else { 1048 cmin = (rcp.rank < cmin.rank) || (rcp.rank == cmin.rank && rcp.index < cmin.index) ? rcp : cmin; 1049 } 1050 } 1051 *cpmin = coneSize ? cmin : missing; 1052 PetscFunctionReturn(PETSC_SUCCESS); 1053 } 1054 1055 /* 1056 Each shared face has an entry in the candidates array: 1057 (-1, coneSize-1), {(global cone point)} 1058 where the set is missing the point p which we use as the key for the face 1059 */ 1060 static PetscErrorCode DMPlexAddSharedFace_Private(DM dm, PetscSection candidateSection, PetscSFNode candidates[], PetscHMapIJ faceHash, PetscInt p, PetscBool debug) 1061 { 1062 MPI_Comm comm; 1063 const PetscInt *support; 1064 PetscInt supportSize, s, off = 0, idx = 0, overlap, cellHeight, height; 1065 PetscMPIInt rank; 1066 1067 PetscFunctionBegin; 1068 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1069 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1070 PetscCall(DMPlexGetOverlap(dm, &overlap)); 1071 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 1072 PetscCall(DMPlexGetPointHeight(dm, p, &height)); 1073 if (!overlap && height <= cellHeight + 1) { 1074 /* cells can't be shared for non-overlapping meshes */ 1075 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Skipping face %" PetscInt_FMT " to avoid adding cell to hashmap since this is nonoverlapping mesh\n", rank, p)); 1076 PetscFunctionReturn(PETSC_SUCCESS); 1077 } 1078 PetscCall(DMPlexGetSupportSize(dm, p, &supportSize)); 1079 PetscCall(DMPlexGetSupport(dm, p, &support)); 1080 if (candidates) PetscCall(PetscSectionGetOffset(candidateSection, p, &off)); 1081 for (s = 0; s < supportSize; ++s) { 1082 const PetscInt face = support[s]; 1083 const PetscInt *cone; 1084 PetscSFNode cpmin = {-1, -1}, rp = {-1, -1}; 1085 PetscInt coneSize, c, f; 1086 PetscBool isShared = PETSC_FALSE; 1087 PetscHashIJKey key; 1088 1089 /* Only add point once */ 1090 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Support face %" PetscInt_FMT "\n", rank, face)); 1091 key.i = p; 1092 key.j = face; 1093 PetscCall(PetscHMapIJGet(faceHash, key, &f)); 1094 if (f >= 0) continue; 1095 PetscCall(DMPlexConeIsShared(dm, face, &isShared)); 1096 PetscCall(DMPlexGetConeMinimum(dm, face, &cpmin)); 1097 PetscCall(DMPlexMapToGlobalPoint(dm, p, &rp, NULL)); 1098 if (debug) { 1099 PetscCall(PetscSynchronizedPrintf(comm, "[%d] Face point %" PetscInt_FMT " is shared: %d\n", rank, face, (int)isShared)); 1100 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)); 1101 } 1102 if (isShared && (rp.rank == cpmin.rank && rp.index == cpmin.index)) { 1103 PetscCall(PetscHMapIJSet(faceHash, key, p)); 1104 if (candidates) { 1105 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Adding shared face %" PetscInt_FMT " at idx %" PetscInt_FMT "\n[%d] ", rank, face, idx, rank)); 1106 PetscCall(DMPlexGetConeSize(dm, face, &coneSize)); 1107 PetscCall(DMPlexGetCone(dm, face, &cone)); 1108 candidates[off + idx].rank = -1; 1109 candidates[off + idx++].index = coneSize - 1; 1110 candidates[off + idx].rank = rank; 1111 candidates[off + idx++].index = face; 1112 for (c = 0; c < coneSize; ++c) { 1113 const PetscInt cp = cone[c]; 1114 1115 if (cp == p) continue; 1116 PetscCall(DMPlexMapToGlobalPoint(dm, cp, &candidates[off + idx], NULL)); 1117 if (debug) PetscCall(PetscSynchronizedPrintf(comm, " (%" PetscInt_FMT ",%" PetscInt_FMT ")", candidates[off + idx].rank, candidates[off + idx].index)); 1118 ++idx; 1119 } 1120 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "\n")); 1121 } else { 1122 /* Add cone size to section */ 1123 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Scheduling shared face %" PetscInt_FMT "\n", rank, face)); 1124 PetscCall(DMPlexGetConeSize(dm, face, &coneSize)); 1125 PetscCall(PetscHMapIJSet(faceHash, key, p)); 1126 PetscCall(PetscSectionAddDof(candidateSection, p, coneSize + 1)); 1127 } 1128 } 1129 } 1130 PetscFunctionReturn(PETSC_SUCCESS); 1131 } 1132 1133 /*@ 1134 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation 1135 1136 Collective 1137 1138 Input Parameters: 1139 + dm - The interpolated `DMPLEX` 1140 - pointSF - The initial `PetscSF` without interpolated points 1141 1142 Level: developer 1143 1144 Note: 1145 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` 1146 1147 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()` 1148 @*/ 1149 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 1150 { 1151 MPI_Comm comm; 1152 PetscHMapIJ remoteHash; 1153 PetscHMapI claimshash; 1154 PetscSection candidateSection, candidateRemoteSection, claimSection; 1155 PetscSFNode *candidates, *candidatesRemote, *claims; 1156 const PetscInt *localPoints, *rootdegree; 1157 const PetscSFNode *remotePoints; 1158 PetscInt ov, Nr, r, Nl, l; 1159 PetscInt candidatesSize, candidatesRemoteSize, claimsSize; 1160 PetscBool flg, debug = PETSC_FALSE; 1161 PetscMPIInt rank; 1162 1163 PetscFunctionBegin; 1164 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1165 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 2); 1166 PetscCall(DMPlexIsDistributed(dm, &flg)); 1167 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 1168 /* Set initial SF so that lower level queries work */ 1169 PetscCall(DMSetPointSF(dm, pointSF)); 1170 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1171 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1172 PetscCall(DMPlexGetOverlap(dm, &ov)); 1173 PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1174 PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug)); 1175 PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view")); 1176 PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view")); 1177 PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0)); 1178 /* Step 0: Precalculations */ 1179 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints)); 1180 PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 1181 PetscCall(PetscHMapIJCreate(&remoteHash)); 1182 for (l = 0; l < Nl; ++l) { 1183 PetscHashIJKey key; 1184 key.i = remotePoints[l].index; 1185 key.j = remotePoints[l].rank; 1186 PetscCall(PetscHMapIJSet(remoteHash, key, l)); 1187 } 1188 /* Compute root degree to identify shared points */ 1189 PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree)); 1190 PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree)); 1191 PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree)); 1192 /* 1193 1) Loop over each leaf point $p$ at depth $d$ in the SF 1194 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which 1195 \begin{itemize} 1196 \item all cone points of $f$ are shared 1197 \item $p$ is the cone point with smallest canonical number 1198 \end{itemize} 1199 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$ 1200 \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 1201 \item Send the root face from the root back to all leaf process 1202 \item Leaf processes add the shared face to the SF 1203 */ 1204 /* Step 1: Construct section+SFNode array 1205 The section has entries for all shared faces for which we have a leaf point in the cone 1206 The array holds candidate shared faces, each face is referred to by the leaf point */ 1207 PetscCall(PetscSectionCreate(comm, &candidateSection)); 1208 PetscCall(PetscSectionSetChart(candidateSection, 0, Nr)); 1209 { 1210 PetscHMapIJ faceHash; 1211 1212 PetscCall(PetscHMapIJCreate(&faceHash)); 1213 for (l = 0; l < Nl; ++l) { 1214 const PetscInt p = localPoints[l]; 1215 1216 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %" PetscInt_FMT "\n", rank, p)); 1217 PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug)); 1218 } 1219 PetscCall(PetscHMapIJClear(faceHash)); 1220 PetscCall(PetscSectionSetUp(candidateSection)); 1221 PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize)); 1222 PetscCall(PetscMalloc1(candidatesSize, &candidates)); 1223 for (l = 0; l < Nl; ++l) { 1224 const PetscInt p = localPoints[l]; 1225 1226 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %" PetscInt_FMT "\n", rank, p)); 1227 PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug)); 1228 } 1229 PetscCall(PetscHMapIJDestroy(&faceHash)); 1230 if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL)); 1231 } 1232 PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section")); 1233 PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view")); 1234 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates)); 1235 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1236 /* Note that this section is indexed by offsets into leaves, not by point number */ 1237 { 1238 PetscSF sfMulti, sfInverse, sfCandidates; 1239 PetscInt *remoteOffsets; 1240 1241 PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti)); 1242 PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse)); 1243 PetscCall(PetscSectionCreate(comm, &candidateRemoteSection)); 1244 PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection)); 1245 PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates)); 1246 PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize)); 1247 PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote)); 1248 PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE)); 1249 PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE)); 1250 PetscCall(PetscSFDestroy(&sfInverse)); 1251 PetscCall(PetscSFDestroy(&sfCandidates)); 1252 PetscCall(PetscFree(remoteOffsets)); 1253 1254 PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section")); 1255 PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view")); 1256 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote)); 1257 } 1258 /* 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 */ 1259 { 1260 PetscHMapIJKLRemote faceTable; 1261 PetscInt idx, idx2; 1262 1263 PetscCall(PetscHMapIJKLRemoteCreate(&faceTable)); 1264 /* There is a section point for every leaf attached to a given root point */ 1265 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) { 1266 PetscInt deg; 1267 1268 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1269 PetscInt offset, dof, d; 1270 1271 PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof)); 1272 PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset)); 1273 /* dof may include many faces from the remote process */ 1274 for (d = 0; d < dof; ++d) { 1275 const PetscInt hidx = offset + d; 1276 const PetscInt Np = candidatesRemote[hidx].index + 1; 1277 const PetscSFNode rface = candidatesRemote[hidx + 1]; 1278 const PetscSFNode *fcone = &candidatesRemote[hidx + 2]; 1279 PetscSFNode fcp0; 1280 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1281 const PetscInt *join = NULL; 1282 PetscHMapIJKLRemoteKey key; 1283 PetscHashIter iter; 1284 PetscBool missing, mapToLocalPointFailed = PETSC_FALSE; 1285 PetscInt points[1024], p, joinSize; 1286 1287 if (debug) 1288 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, 1289 rface.index, r, idx, d, Np)); 1290 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); 1291 fcp0.rank = rank; 1292 fcp0.index = r; 1293 d += Np; 1294 /* Put remote face in hash table */ 1295 key.i = fcp0; 1296 key.j = fcone[0]; 1297 key.k = Np > 2 ? fcone[1] : pmax; 1298 key.l = Np > 3 ? fcone[2] : pmax; 1299 PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key)); 1300 PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing)); 1301 if (missing) { 1302 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank)); 1303 PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface)); 1304 } else { 1305 PetscSFNode oface; 1306 1307 PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface)); 1308 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1309 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank)); 1310 PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, rface)); 1311 } 1312 } 1313 /* Check for local face */ 1314 points[0] = r; 1315 for (p = 1; p < Np; ++p) { 1316 PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed)); 1317 if (mapToLocalPointFailed) break; /* We got a point not in our overlap */ 1318 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking local candidate %" PetscInt_FMT "\n", rank, points[p])); 1319 } 1320 if (mapToLocalPointFailed) continue; 1321 PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join)); 1322 if (joinSize == 1) { 1323 PetscSFNode lface; 1324 PetscSFNode oface; 1325 1326 /* Always replace with local face */ 1327 lface.rank = rank; 1328 lface.index = join[0]; 1329 PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &oface)); 1330 if (debug) 1331 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)); 1332 PetscCall(PetscHMapIJKLRemoteIterSet(faceTable, iter, lface)); 1333 } 1334 PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join)); 1335 } 1336 } 1337 /* Put back faces for this root */ 1338 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) { 1339 PetscInt offset, dof, d; 1340 1341 PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof)); 1342 PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset)); 1343 /* dof may include many faces from the remote process */ 1344 for (d = 0; d < dof; ++d) { 1345 const PetscInt hidx = offset + d; 1346 const PetscInt Np = candidatesRemote[hidx].index + 1; 1347 const PetscSFNode *fcone = &candidatesRemote[hidx + 2]; 1348 PetscSFNode fcp0; 1349 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1350 PetscHMapIJKLRemoteKey key; 1351 PetscHashIter iter; 1352 PetscBool missing; 1353 1354 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx)); 1355 PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np); 1356 fcp0.rank = rank; 1357 fcp0.index = r; 1358 d += Np; 1359 /* Find remote face in hash table */ 1360 key.i = fcp0; 1361 key.j = fcone[0]; 1362 key.k = Np > 2 ? fcone[1] : pmax; 1363 key.l = Np > 3 ? fcone[2] : pmax; 1364 PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key)); 1365 if (debug) 1366 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, 1367 key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index)); 1368 PetscCall(PetscHMapIJKLRemotePut(faceTable, key, &iter, &missing)); 1369 PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2); 1370 PetscCall(PetscHMapIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx])); 1371 } 1372 } 1373 } 1374 if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL)); 1375 PetscCall(PetscHMapIJKLRemoteDestroy(&faceTable)); 1376 } 1377 /* Step 4: Push back owned faces */ 1378 { 1379 PetscSF sfMulti, sfClaims, sfPointNew; 1380 PetscSFNode *remotePointsNew; 1381 PetscInt *remoteOffsets, *localPointsNew; 1382 PetscInt pStart, pEnd, r, NlNew, p; 1383 1384 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1385 PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti)); 1386 PetscCall(PetscSectionCreate(comm, &claimSection)); 1387 PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection)); 1388 PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims)); 1389 PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize)); 1390 PetscCall(PetscMalloc1(claimsSize, &claims)); 1391 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1392 PetscCall(PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE)); 1393 PetscCall(PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE)); 1394 PetscCall(PetscSFDestroy(&sfClaims)); 1395 PetscCall(PetscFree(remoteOffsets)); 1396 PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section")); 1397 PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view")); 1398 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims)); 1399 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */ 1400 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */ 1401 PetscCall(PetscHMapICreate(&claimshash)); 1402 for (r = 0; r < Nr; ++r) { 1403 PetscInt dof, off, d; 1404 1405 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %" PetscInt_FMT "\n", rank, r)); 1406 PetscCall(PetscSectionGetDof(candidateSection, r, &dof)); 1407 PetscCall(PetscSectionGetOffset(candidateSection, r, &off)); 1408 for (d = 0; d < dof;) { 1409 if (claims[off + d].rank >= 0) { 1410 const PetscInt faceInd = off + d; 1411 const PetscInt Np = candidates[off + d].index; 1412 const PetscInt *join = NULL; 1413 PetscInt joinSize, points[1024], c; 1414 1415 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index)); 1416 points[0] = r; 1417 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[0])); 1418 for (c = 0, d += 2; c < Np; ++c, ++d) { 1419 PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL)); 1420 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[c + 1])); 1421 } 1422 PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join)); 1423 if (joinSize == 1) { 1424 if (claims[faceInd].rank == rank) { 1425 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0])); 1426 } else { 1427 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found local face %" PetscInt_FMT "\n", rank, join[0])); 1428 PetscCall(PetscHMapISet(claimshash, join[0], faceInd)); 1429 } 1430 } else { 1431 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank)); 1432 } 1433 PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join)); 1434 } else { 1435 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] No claim for point %" PetscInt_FMT "\n", rank, r)); 1436 d += claims[off + d].index + 1; 1437 } 1438 } 1439 } 1440 if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL)); 1441 /* Step 6) Create new pointSF from hashed claims */ 1442 PetscCall(PetscHMapIGetSize(claimshash, &NlNew)); 1443 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1444 PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew)); 1445 PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew)); 1446 for (l = 0; l < Nl; ++l) { 1447 localPointsNew[l] = localPoints[l]; 1448 remotePointsNew[l].index = remotePoints[l].index; 1449 remotePointsNew[l].rank = remotePoints[l].rank; 1450 } 1451 p = Nl; 1452 PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew)); 1453 /* We sort new points, and assume they are numbered after all existing points */ 1454 PetscCall(PetscSortInt(NlNew, PetscSafePointerPlusOffset(localPointsNew, Nl))); 1455 for (p = Nl; p < Nl + NlNew; ++p) { 1456 PetscInt off; 1457 PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off)); 1458 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); 1459 remotePointsNew[p] = claims[off]; 1460 } 1461 PetscCall(PetscSFCreate(comm, &sfPointNew)); 1462 PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1463 PetscCall(PetscSFSetUp(sfPointNew)); 1464 PetscCall(DMSetPointSF(dm, sfPointNew)); 1465 PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view")); 1466 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE)); 1467 PetscCall(PetscSFDestroy(&sfPointNew)); 1468 PetscCall(PetscHMapIDestroy(&claimshash)); 1469 } 1470 PetscCall(PetscHMapIJDestroy(&remoteHash)); 1471 PetscCall(PetscSectionDestroy(&candidateSection)); 1472 PetscCall(PetscSectionDestroy(&candidateRemoteSection)); 1473 PetscCall(PetscSectionDestroy(&claimSection)); 1474 PetscCall(PetscFree(candidates)); 1475 PetscCall(PetscFree(candidatesRemote)); 1476 PetscCall(PetscFree(claims)); 1477 PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0)); 1478 PetscFunctionReturn(PETSC_SUCCESS); 1479 } 1480 1481 /*@ 1482 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1483 1484 Collective 1485 1486 Input Parameter: 1487 . dm - The `DMPLEX` object with only cells and vertices 1488 1489 Output Parameter: 1490 . dmInt - The complete `DMPLEX` object 1491 1492 Level: intermediate 1493 1494 Note: 1495 Labels and coordinates are copied. 1496 1497 Developer Notes: 1498 It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`. 1499 1500 .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()` 1501 @*/ 1502 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1503 { 1504 DMPlexInterpolatedFlag interpolated; 1505 DM idm, odm = dm; 1506 PetscSF sfPoint; 1507 PetscInt depth, dim, d; 1508 const char *name; 1509 PetscBool flg = PETSC_TRUE; 1510 1511 PetscFunctionBegin; 1512 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1513 PetscAssertPointer(dmInt, 2); 1514 PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0)); 1515 PetscCall(DMPlexGetDepth(dm, &depth)); 1516 PetscCall(DMGetDimension(dm, &dim)); 1517 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1518 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1519 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1520 PetscCall(PetscObjectReference((PetscObject)dm)); 1521 idm = dm; 1522 } else { 1523 for (d = 1; d < dim; ++d) { 1524 const char *prefix; 1525 1526 /* Create interpolated mesh */ 1527 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm)); 1528 PetscCall(DMSetType(idm, DMPLEX)); 1529 PetscCall(DMSetDimension(idm, dim)); 1530 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 1531 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)idm, prefix)); 1532 if (depth > 0) { 1533 PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm)); 1534 PetscCall(DMGetPointSF(odm, &sfPoint)); 1535 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE)); 1536 { 1537 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1538 PetscInt nroots; 1539 PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1540 if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint)); 1541 } 1542 } 1543 if (odm != dm) PetscCall(DMDestroy(&odm)); 1544 odm = idm; 1545 } 1546 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 1547 PetscCall(PetscObjectSetName((PetscObject)idm, name)); 1548 PetscCall(DMPlexCopyCoordinates(dm, idm)); 1549 PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1550 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL)); 1551 if (flg) PetscCall(DMPlexOrientInterface_Internal(idm)); 1552 } 1553 /* This function makes the mesh fully interpolated on all ranks */ 1554 { 1555 DM_Plex *plex = (DM_Plex *)idm->data; 1556 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1557 } 1558 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm)); 1559 *dmInt = idm; 1560 PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0)); 1561 PetscFunctionReturn(PETSC_SUCCESS); 1562 } 1563 1564 /*@ 1565 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1566 1567 Collective 1568 1569 Input Parameter: 1570 . dmA - The `DMPLEX` object with initial coordinates 1571 1572 Output Parameter: 1573 . dmB - The `DMPLEX` object with copied coordinates 1574 1575 Level: intermediate 1576 1577 Notes: 1578 This is typically used when adding pieces other than vertices to a mesh 1579 1580 This function does not copy localized coordinates. 1581 1582 .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()` 1583 @*/ 1584 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1585 { 1586 Vec coordinatesA, coordinatesB; 1587 VecType vtype; 1588 PetscSection coordSectionA, coordSectionB; 1589 PetscScalar *coordsA, *coordsB; 1590 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1591 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim; 1592 1593 PetscFunctionBegin; 1594 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1595 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1596 if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS); 1597 PetscCall(DMGetCoordinateDim(dmA, &cdim)); 1598 PetscCall(DMSetCoordinateDim(dmB, cdim)); 1599 PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA)); 1600 PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB)); 1601 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); 1602 /* Copy over discretization if it exists */ 1603 { 1604 DM cdmA, cdmB; 1605 PetscDS dsA, dsB; 1606 PetscObject objA, objB; 1607 PetscClassId idA, idB; 1608 const PetscScalar *constants; 1609 PetscInt cdim, Nc; 1610 1611 PetscCall(DMGetCoordinateDM(dmA, &cdmA)); 1612 PetscCall(DMGetCoordinateDM(dmB, &cdmB)); 1613 PetscCall(DMGetField(cdmA, 0, NULL, &objA)); 1614 PetscCall(DMGetField(cdmB, 0, NULL, &objB)); 1615 PetscCall(PetscObjectGetClassId(objA, &idA)); 1616 PetscCall(PetscObjectGetClassId(objB, &idB)); 1617 if ((idA == PETSCFE_CLASSID) && (idA != idB)) { 1618 PetscCall(DMSetField(cdmB, 0, NULL, objA)); 1619 PetscCall(DMCreateDS(cdmB)); 1620 PetscCall(DMGetDS(cdmA, &dsA)); 1621 PetscCall(DMGetDS(cdmB, &dsB)); 1622 PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim)); 1623 PetscCall(PetscDSSetCoordinateDimension(dsB, cdim)); 1624 PetscCall(PetscDSGetConstants(dsA, &Nc, &constants)); 1625 PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants)); 1626 } 1627 } 1628 PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA)); 1629 PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB)); 1630 PetscCall(DMGetCoordinateSection(dmA, &coordSectionA)); 1631 PetscCall(DMGetCoordinateSection(dmB, &coordSectionB)); 1632 if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS); 1633 PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf)); 1634 if (!Nf) PetscFunctionReturn(PETSC_SUCCESS); 1635 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf); 1636 if (!coordSectionB) { 1637 PetscInt dim; 1638 1639 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB)); 1640 PetscCall(DMGetCoordinateDim(dmA, &dim)); 1641 PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB)); 1642 PetscCall(PetscObjectDereference((PetscObject)coordSectionB)); 1643 } 1644 PetscCall(PetscSectionSetNumFields(coordSectionB, 1)); 1645 PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim)); 1646 PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim)); 1647 PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE)); 1648 cS = vStartB; 1649 cE = vEndB; 1650 PetscCall(PetscSectionSetChart(coordSectionB, cS, cE)); 1651 for (v = vStartB; v < vEndB; ++v) { 1652 PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim)); 1653 PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim)); 1654 } 1655 PetscCall(PetscSectionSetUp(coordSectionB)); 1656 PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB)); 1657 PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA)); 1658 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB)); 1659 PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates")); 1660 PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE)); 1661 PetscCall(VecGetBlockSize(coordinatesA, &d)); 1662 PetscCall(VecSetBlockSize(coordinatesB, d)); 1663 PetscCall(VecGetType(coordinatesA, &vtype)); 1664 PetscCall(VecSetType(coordinatesB, vtype)); 1665 PetscCall(VecGetArray(coordinatesA, &coordsA)); 1666 PetscCall(VecGetArray(coordinatesB, &coordsB)); 1667 for (v = 0; v < vEndB - vStartB; ++v) { 1668 PetscInt offA, offB; 1669 1670 PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA)); 1671 PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB)); 1672 for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d]; 1673 } 1674 PetscCall(VecRestoreArray(coordinatesA, &coordsA)); 1675 PetscCall(VecRestoreArray(coordinatesB, &coordsB)); 1676 PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB)); 1677 PetscCall(VecDestroy(&coordinatesB)); 1678 PetscFunctionReturn(PETSC_SUCCESS); 1679 } 1680 1681 /*@ 1682 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1683 1684 Collective 1685 1686 Input Parameter: 1687 . dm - The complete `DMPLEX` object 1688 1689 Output Parameter: 1690 . dmUnint - The `DMPLEX` object with only cells and vertices 1691 1692 Level: intermediate 1693 1694 Note: 1695 It does not copy over the coordinates. 1696 1697 Developer Notes: 1698 Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`. 1699 1700 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()` 1701 @*/ 1702 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1703 { 1704 DMPlexInterpolatedFlag interpolated; 1705 DM udm; 1706 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 1707 1708 PetscFunctionBegin; 1709 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1710 PetscAssertPointer(dmUnint, 2); 1711 PetscCall(PetscLogEventBegin(DMPLEX_Uninterpolate, dm, 0, 0, 0)); 1712 PetscCall(DMGetDimension(dm, &dim)); 1713 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1714 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1715 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1716 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1717 PetscCall(PetscObjectReference((PetscObject)dm)); 1718 *dmUnint = dm; 1719 PetscFunctionReturn(PETSC_SUCCESS); 1720 } 1721 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1722 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1723 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm)); 1724 PetscCall(DMSetType(udm, DMPLEX)); 1725 PetscCall(DMSetDimension(udm, dim)); 1726 PetscCall(DMPlexSetChart(udm, cStart, vEnd)); 1727 for (c = cStart; c < cEnd; ++c) { 1728 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1729 1730 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1731 for (cl = 0; cl < closureSize * 2; cl += 2) { 1732 const PetscInt p = closure[cl]; 1733 1734 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1735 } 1736 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1737 PetscCall(DMPlexSetConeSize(udm, c, coneSize)); 1738 maxConeSize = PetscMax(maxConeSize, coneSize); 1739 } 1740 PetscCall(DMSetUp(udm)); 1741 PetscCall(PetscMalloc1(maxConeSize, &cone)); 1742 for (c = cStart; c < cEnd; ++c) { 1743 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1744 1745 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1746 for (cl = 0; cl < closureSize * 2; cl += 2) { 1747 const PetscInt p = closure[cl]; 1748 1749 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1750 } 1751 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1752 PetscCall(DMPlexSetCone(udm, c, cone)); 1753 } 1754 PetscCall(PetscFree(cone)); 1755 PetscCall(DMPlexSymmetrize(udm)); 1756 PetscCall(DMPlexStratify(udm)); 1757 /* Reduce SF */ 1758 { 1759 PetscSF sfPoint, sfPointUn; 1760 const PetscSFNode *remotePoints; 1761 const PetscInt *localPoints; 1762 PetscSFNode *remotePointsUn; 1763 PetscInt *localPointsUn; 1764 PetscInt numRoots, numLeaves, l; 1765 PetscInt numLeavesUn = 0, n = 0; 1766 1767 /* Get original SF information */ 1768 PetscCall(DMGetPointSF(dm, &sfPoint)); 1769 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE)); 1770 PetscCall(DMGetPointSF(udm, &sfPointUn)); 1771 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1772 if (numRoots >= 0) { 1773 /* Allocate space for cells and vertices */ 1774 for (l = 0; l < numLeaves; ++l) { 1775 const PetscInt p = localPoints[l]; 1776 1777 if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++; 1778 } 1779 /* Fill in leaves */ 1780 PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn)); 1781 PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn)); 1782 for (l = 0; l < numLeaves; l++) { 1783 const PetscInt p = localPoints[l]; 1784 1785 if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) { 1786 localPointsUn[n] = p; 1787 remotePointsUn[n].rank = remotePoints[l].rank; 1788 remotePointsUn[n].index = remotePoints[l].index; 1789 ++n; 1790 } 1791 } 1792 PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn); 1793 PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER)); 1794 } 1795 } 1796 /* This function makes the mesh fully uninterpolated on all ranks */ 1797 { 1798 DM_Plex *plex = (DM_Plex *)udm->data; 1799 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1800 } 1801 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm)); 1802 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE)); 1803 *dmUnint = udm; 1804 PetscCall(PetscLogEventEnd(DMPLEX_Uninterpolate, dm, 0, 0, 0)); 1805 PetscFunctionReturn(PETSC_SUCCESS); 1806 } 1807 1808 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1809 { 1810 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1811 MPI_Comm comm; 1812 1813 PetscFunctionBegin; 1814 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1815 PetscCall(DMPlexGetDepth(dm, &depth)); 1816 PetscCall(DMGetDimension(dm, &dim)); 1817 1818 if (depth == dim) { 1819 *interpolated = DMPLEX_INTERPOLATED_FULL; 1820 if (!dim) goto finish; 1821 1822 /* Check points at height = dim are vertices (have no cones) */ 1823 PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd)); 1824 for (p = pStart; p < pEnd; p++) { 1825 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1826 if (coneSize) { 1827 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1828 goto finish; 1829 } 1830 } 1831 1832 /* Check points at height < dim have cones */ 1833 for (h = 0; h < dim; h++) { 1834 PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 1835 for (p = pStart; p < pEnd; p++) { 1836 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1837 if (!coneSize) { 1838 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1839 goto finish; 1840 } 1841 } 1842 } 1843 } else if (depth == 1) { 1844 *interpolated = DMPLEX_INTERPOLATED_NONE; 1845 } else { 1846 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1847 } 1848 finish: 1849 PetscFunctionReturn(PETSC_SUCCESS); 1850 } 1851 1852 /*@ 1853 DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated. 1854 1855 Not Collective 1856 1857 Input Parameter: 1858 . dm - The `DMPLEX` object 1859 1860 Output Parameter: 1861 . interpolated - Flag whether the `DM` is interpolated 1862 1863 Level: intermediate 1864 1865 Notes: 1866 Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective 1867 so the results can be different on different ranks in special cases. 1868 However, `DMPlexInterpolate()` guarantees the result is the same on all. 1869 1870 Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`. 1871 1872 Developer Notes: 1873 Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`. 1874 1875 If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called. 1876 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1877 `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`. 1878 1879 If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated. 1880 1881 `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`, 1882 and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`. 1883 1884 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()` 1885 @*/ 1886 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1887 { 1888 DM_Plex *plex = (DM_Plex *)dm->data; 1889 1890 PetscFunctionBegin; 1891 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1892 PetscAssertPointer(interpolated, 2); 1893 if (plex->interpolated < 0) { 1894 PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated)); 1895 } else if (PetscDefined(USE_DEBUG)) { 1896 DMPlexInterpolatedFlag flg; 1897 1898 PetscCall(DMPlexIsInterpolated_Internal(dm, &flg)); 1899 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]); 1900 } 1901 *interpolated = plex->interpolated; 1902 PetscFunctionReturn(PETSC_SUCCESS); 1903 } 1904 1905 /*@ 1906 DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner). 1907 1908 Collective 1909 1910 Input Parameter: 1911 . dm - The `DMPLEX` object 1912 1913 Output Parameter: 1914 . interpolated - Flag whether the `DM` is interpolated 1915 1916 Level: intermediate 1917 1918 Notes: 1919 Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks. 1920 1921 This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks. 1922 1923 Developer Notes: 1924 Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`. 1925 1926 If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated. 1927 `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 1928 if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`, 1929 otherwise sets plex->interpolatedCollective = plex->interpolated. 1930 1931 If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective. 1932 1933 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()` 1934 @*/ 1935 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1936 { 1937 DM_Plex *plex = (DM_Plex *)dm->data; 1938 PetscBool debug = PETSC_FALSE; 1939 1940 PetscFunctionBegin; 1941 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1942 PetscAssertPointer(interpolated, 2); 1943 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL)); 1944 if (plex->interpolatedCollective < 0) { 1945 DMPlexInterpolatedFlag min, max; 1946 MPI_Comm comm; 1947 1948 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1949 PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective)); 1950 PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm)); 1951 PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm)); 1952 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1953 if (debug) { 1954 PetscMPIInt rank; 1955 1956 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1957 PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective])); 1958 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 1959 } 1960 } 1961 *interpolated = plex->interpolatedCollective; 1962 PetscFunctionReturn(PETSC_SUCCESS); 1963 } 1964