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