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