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