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