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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 943 owned: 944 remotePoint->rank = rank; 945 remotePoint->index = localPoint; 946 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 960 PetscCall(PetscFindInt(p, Nl, locals, &idx)); 961 if (idx >= 0) { 962 *isShared = PETSC_TRUE; 963 PetscFunctionReturn(PETSC_SUCCESS); 964 } 965 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 966 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 967 if (rootdegree[p] > 0) *isShared = PETSC_TRUE; 968 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 1090 } 1091 1092 /*@ 1093 DMPlexInterpolatePointSF - Insert interpolated points in the overlap into the `PointSF` in parallel, following local interpolation 1094 1095 Collective 1096 1097 Input Parameters: 1098 + dm - The interpolated `DMPLEX` 1099 - pointSF - The initial `PetscSF` without interpolated points 1100 1101 Level: developer 1102 1103 Note: 1104 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` 1105 1106 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexUninterpolate()` 1107 @*/ 1108 PetscErrorCode DMPlexInterpolatePointSF(DM dm, PetscSF pointSF) 1109 { 1110 MPI_Comm comm; 1111 PetscHMapIJ remoteHash; 1112 PetscHMapI claimshash; 1113 PetscSection candidateSection, candidateRemoteSection, claimSection; 1114 PetscSFNode *candidates, *candidatesRemote, *claims; 1115 const PetscInt *localPoints, *rootdegree; 1116 const PetscSFNode *remotePoints; 1117 PetscInt ov, Nr, r, Nl, l; 1118 PetscInt candidatesSize, candidatesRemoteSize, claimsSize; 1119 PetscBool flg, debug = PETSC_FALSE; 1120 PetscMPIInt rank; 1121 1122 PetscFunctionBegin; 1123 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1124 PetscValidHeaderSpecific(pointSF, PETSCSF_CLASSID, 2); 1125 PetscCall(DMPlexIsDistributed(dm, &flg)); 1126 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 1127 /* Set initial SF so that lower level queries work */ 1128 PetscCall(DMSetPointSF(dm, pointSF)); 1129 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1130 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1131 PetscCall(DMPlexGetOverlap(dm, &ov)); 1132 PetscCheck(!ov, comm, PETSC_ERR_SUP, "Interpolation of overlapped DMPlex not implemented yet"); 1133 PetscCall(PetscOptionsHasName(NULL, ((PetscObject)dm)->prefix, "-dmplex_interp_debug", &debug)); 1134 PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_interp_pre_view")); 1135 PetscCall(PetscObjectViewFromOptions((PetscObject)pointSF, NULL, "-petscsf_interp_pre_view")); 1136 PetscCall(PetscLogEventBegin(DMPLEX_InterpolateSF, dm, 0, 0, 0)); 1137 /* Step 0: Precalculations */ 1138 PetscCall(PetscSFGetGraph(pointSF, &Nr, &Nl, &localPoints, &remotePoints)); 1139 PetscCheck(Nr >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "This DMPlex is distributed but input PointSF has no graph set"); 1140 PetscCall(PetscHMapIJCreate(&remoteHash)); 1141 for (l = 0; l < Nl; ++l) { 1142 PetscHashIJKey key; 1143 key.i = remotePoints[l].index; 1144 key.j = remotePoints[l].rank; 1145 PetscCall(PetscHMapIJSet(remoteHash, key, l)); 1146 } 1147 /* Compute root degree to identify shared points */ 1148 PetscCall(PetscSFComputeDegreeBegin(pointSF, &rootdegree)); 1149 PetscCall(PetscSFComputeDegreeEnd(pointSF, &rootdegree)); 1150 PetscCall(IntArrayViewFromOptions(comm, "-interp_root_degree_view", "Root degree", "point", "degree", Nr, rootdegree)); 1151 /* 1152 1) Loop over each leaf point $p$ at depth $d$ in the SF 1153 \item Get set $F(p)$ of faces $f$ in the support of $p$ for which 1154 \begin{itemize} 1155 \item all cone points of $f$ are shared 1156 \item $p$ is the cone point with smallest canonical number 1157 \end{itemize} 1158 \item Send $F(p)$ and the cone of each face to the active root point $r(p)$ 1159 \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 1160 \item Send the root face from the root back to all leaf process 1161 \item Leaf processes add the shared face to the SF 1162 */ 1163 /* Step 1: Construct section+SFNode array 1164 The section has entries for all shared faces for which we have a leaf point in the cone 1165 The array holds candidate shared faces, each face is referred to by the leaf point */ 1166 PetscCall(PetscSectionCreate(comm, &candidateSection)); 1167 PetscCall(PetscSectionSetChart(candidateSection, 0, Nr)); 1168 { 1169 PetscHMapIJ faceHash; 1170 1171 PetscCall(PetscHMapIJCreate(&faceHash)); 1172 for (l = 0; l < Nl; ++l) { 1173 const PetscInt p = localPoints[l]; 1174 1175 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] First pass leaf point %" PetscInt_FMT "\n", rank, p)); 1176 PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, NULL, faceHash, p, debug)); 1177 } 1178 PetscCall(PetscHMapIJClear(faceHash)); 1179 PetscCall(PetscSectionSetUp(candidateSection)); 1180 PetscCall(PetscSectionGetStorageSize(candidateSection, &candidatesSize)); 1181 PetscCall(PetscMalloc1(candidatesSize, &candidates)); 1182 for (l = 0; l < Nl; ++l) { 1183 const PetscInt p = localPoints[l]; 1184 1185 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Second pass leaf point %" PetscInt_FMT "\n", rank, p)); 1186 PetscCall(DMPlexAddSharedFace_Private(dm, candidateSection, candidates, faceHash, p, debug)); 1187 } 1188 PetscCall(PetscHMapIJDestroy(&faceHash)); 1189 if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL)); 1190 } 1191 PetscCall(PetscObjectSetName((PetscObject)candidateSection, "Candidate Section")); 1192 PetscCall(PetscObjectViewFromOptions((PetscObject)candidateSection, NULL, "-petscsection_interp_candidate_view")); 1193 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_view", "Candidates", NULL, candidatesSize, candidates)); 1194 /* Step 2: Gather candidate section / array pair into the root partition via inverse(multi(pointSF)). */ 1195 /* Note that this section is indexed by offsets into leaves, not by point number */ 1196 { 1197 PetscSF sfMulti, sfInverse, sfCandidates; 1198 PetscInt *remoteOffsets; 1199 1200 PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti)); 1201 PetscCall(PetscSFCreateInverseSF(sfMulti, &sfInverse)); 1202 PetscCall(PetscSectionCreate(comm, &candidateRemoteSection)); 1203 PetscCall(PetscSFDistributeSection(sfInverse, candidateSection, &remoteOffsets, candidateRemoteSection)); 1204 PetscCall(PetscSFCreateSectionSF(sfInverse, candidateSection, remoteOffsets, candidateRemoteSection, &sfCandidates)); 1205 PetscCall(PetscSectionGetStorageSize(candidateRemoteSection, &candidatesRemoteSize)); 1206 PetscCall(PetscMalloc1(candidatesRemoteSize, &candidatesRemote)); 1207 PetscCall(PetscSFBcastBegin(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE)); 1208 PetscCall(PetscSFBcastEnd(sfCandidates, MPIU_2INT, candidates, candidatesRemote, MPI_REPLACE)); 1209 PetscCall(PetscSFDestroy(&sfInverse)); 1210 PetscCall(PetscSFDestroy(&sfCandidates)); 1211 PetscCall(PetscFree(remoteOffsets)); 1212 1213 PetscCall(PetscObjectSetName((PetscObject)candidateRemoteSection, "Remote Candidate Section")); 1214 PetscCall(PetscObjectViewFromOptions((PetscObject)candidateRemoteSection, NULL, "-petscsection_interp_candidate_remote_view")); 1215 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_candidate_remote_view", "Remote Candidates", NULL, candidatesRemoteSize, candidatesRemote)); 1216 } 1217 /* 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 */ 1218 { 1219 PetscHashIJKLRemote faceTable; 1220 PetscInt idx, idx2; 1221 1222 PetscCall(PetscHashIJKLRemoteCreate(&faceTable)); 1223 /* There is a section point for every leaf attached to a given root point */ 1224 for (r = 0, idx = 0, idx2 = 0; r < Nr; ++r) { 1225 PetscInt deg; 1226 1227 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx) { 1228 PetscInt offset, dof, d; 1229 1230 PetscCall(PetscSectionGetDof(candidateRemoteSection, idx, &dof)); 1231 PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx, &offset)); 1232 /* dof may include many faces from the remote process */ 1233 for (d = 0; d < dof; ++d) { 1234 const PetscInt hidx = offset + d; 1235 const PetscInt Np = candidatesRemote[hidx].index + 1; 1236 const PetscSFNode rface = candidatesRemote[hidx + 1]; 1237 const PetscSFNode *fcone = &candidatesRemote[hidx + 2]; 1238 PetscSFNode fcp0; 1239 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1240 const PetscInt *join = NULL; 1241 PetscHashIJKLRemoteKey key; 1242 PetscHashIter iter; 1243 PetscBool missing, mapToLocalPointFailed = PETSC_FALSE; 1244 PetscInt points[1024], p, joinSize; 1245 1246 if (debug) 1247 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, 1248 rface.index, r, idx, d, Np)); 1249 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); 1250 fcp0.rank = rank; 1251 fcp0.index = r; 1252 d += Np; 1253 /* Put remote face in hash table */ 1254 key.i = fcp0; 1255 key.j = fcone[0]; 1256 key.k = Np > 2 ? fcone[1] : pmax; 1257 key.l = Np > 3 ? fcone[2] : pmax; 1258 PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key)); 1259 PetscCall(PetscHashIJKLRemotePut(faceTable, key, &iter, &missing)); 1260 if (missing) { 1261 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Setting remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank)); 1262 PetscCall(PetscHashIJKLRemoteIterSet(faceTable, iter, rface)); 1263 } else { 1264 PetscSFNode oface; 1265 1266 PetscCall(PetscHashIJKLRemoteIterGet(faceTable, iter, &oface)); 1267 if ((rface.rank < oface.rank) || (rface.rank == oface.rank && rface.index < oface.index)) { 1268 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Replacing with remote face (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, rface.index, rface.rank)); 1269 PetscCall(PetscHashIJKLRemoteIterSet(faceTable, iter, rface)); 1270 } 1271 } 1272 /* Check for local face */ 1273 points[0] = r; 1274 for (p = 1; p < Np; ++p) { 1275 PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, fcone[p - 1], &points[p], &mapToLocalPointFailed)); 1276 if (mapToLocalPointFailed) break; /* We got a point not in our overlap */ 1277 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Checking local candidate %" PetscInt_FMT "\n", rank, points[p])); 1278 } 1279 if (mapToLocalPointFailed) continue; 1280 PetscCall(DMPlexGetJoin(dm, Np, points, &joinSize, &join)); 1281 if (joinSize == 1) { 1282 PetscSFNode lface; 1283 PetscSFNode oface; 1284 1285 /* Always replace with local face */ 1286 lface.rank = rank; 1287 lface.index = join[0]; 1288 PetscCall(PetscHashIJKLRemoteIterGet(faceTable, iter, &oface)); 1289 if (debug) 1290 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)); 1291 PetscCall(PetscHashIJKLRemoteIterSet(faceTable, iter, lface)); 1292 } 1293 PetscCall(DMPlexRestoreJoin(dm, Np, points, &joinSize, &join)); 1294 } 1295 } 1296 /* Put back faces for this root */ 1297 for (deg = 0; deg < rootdegree[r]; ++deg, ++idx2) { 1298 PetscInt offset, dof, d; 1299 1300 PetscCall(PetscSectionGetDof(candidateRemoteSection, idx2, &dof)); 1301 PetscCall(PetscSectionGetOffset(candidateRemoteSection, idx2, &offset)); 1302 /* dof may include many faces from the remote process */ 1303 for (d = 0; d < dof; ++d) { 1304 const PetscInt hidx = offset + d; 1305 const PetscInt Np = candidatesRemote[hidx].index + 1; 1306 const PetscSFNode *fcone = &candidatesRemote[hidx + 2]; 1307 PetscSFNode fcp0; 1308 const PetscSFNode pmax = {PETSC_MAX_INT, PETSC_MAX_INT}; 1309 PetscHashIJKLRemoteKey key; 1310 PetscHashIter iter; 1311 PetscBool missing; 1312 1313 if (debug) PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)dm), "[%d] Entering face at (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, r, idx)); 1314 PetscCheck(Np <= 4, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle faces with %" PetscInt_FMT " cone points", Np); 1315 fcp0.rank = rank; 1316 fcp0.index = r; 1317 d += Np; 1318 /* Find remote face in hash table */ 1319 key.i = fcp0; 1320 key.j = fcone[0]; 1321 key.k = Np > 2 ? fcone[1] : pmax; 1322 key.l = Np > 3 ? fcone[2] : pmax; 1323 PetscCall(PetscSortSFNode(Np, (PetscSFNode *)&key)); 1324 if (debug) 1325 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, 1326 key.i.rank, key.i.index, key.j.rank, key.j.index, key.k.rank, key.k.index, key.l.rank, key.l.index)); 1327 PetscCall(PetscHashIJKLRemotePut(faceTable, key, &iter, &missing)); 1328 PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root %" PetscInt_FMT " Idx %" PetscInt_FMT " ought to have an associated face", r, idx2); 1329 PetscCall(PetscHashIJKLRemoteIterGet(faceTable, iter, &candidatesRemote[hidx])); 1330 } 1331 } 1332 } 1333 if (debug) PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)dm), NULL)); 1334 PetscCall(PetscHashIJKLRemoteDestroy(&faceTable)); 1335 } 1336 /* Step 4: Push back owned faces */ 1337 { 1338 PetscSF sfMulti, sfClaims, sfPointNew; 1339 PetscSFNode *remotePointsNew; 1340 PetscInt *remoteOffsets, *localPointsNew; 1341 PetscInt pStart, pEnd, r, NlNew, p; 1342 1343 /* 4) Push claims back to receiver via the MultiSF and derive new pointSF mapping on receiver */ 1344 PetscCall(PetscSFGetMultiSF(pointSF, &sfMulti)); 1345 PetscCall(PetscSectionCreate(comm, &claimSection)); 1346 PetscCall(PetscSFDistributeSection(sfMulti, candidateRemoteSection, &remoteOffsets, claimSection)); 1347 PetscCall(PetscSFCreateSectionSF(sfMulti, candidateRemoteSection, remoteOffsets, claimSection, &sfClaims)); 1348 PetscCall(PetscSectionGetStorageSize(claimSection, &claimsSize)); 1349 PetscCall(PetscMalloc1(claimsSize, &claims)); 1350 for (p = 0; p < claimsSize; ++p) claims[p].rank = -1; 1351 PetscCall(PetscSFBcastBegin(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE)); 1352 PetscCall(PetscSFBcastEnd(sfClaims, MPIU_2INT, candidatesRemote, claims, MPI_REPLACE)); 1353 PetscCall(PetscSFDestroy(&sfClaims)); 1354 PetscCall(PetscFree(remoteOffsets)); 1355 PetscCall(PetscObjectSetName((PetscObject)claimSection, "Claim Section")); 1356 PetscCall(PetscObjectViewFromOptions((PetscObject)claimSection, NULL, "-petscsection_interp_claim_view")); 1357 PetscCall(SFNodeArrayViewFromOptions(comm, "-petscsection_interp_claim_view", "Claims", NULL, claimsSize, claims)); 1358 /* Step 5) Walk the original section of local supports and add an SF entry for each updated item */ 1359 /* TODO I should not have to do a join here since I already put the face and its cone in the candidate section */ 1360 PetscCall(PetscHMapICreate(&claimshash)); 1361 for (r = 0; r < Nr; ++r) { 1362 PetscInt dof, off, d; 1363 1364 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Checking root for claims %" PetscInt_FMT "\n", rank, r)); 1365 PetscCall(PetscSectionGetDof(candidateSection, r, &dof)); 1366 PetscCall(PetscSectionGetOffset(candidateSection, r, &off)); 1367 for (d = 0; d < dof;) { 1368 if (claims[off + d].rank >= 0) { 1369 const PetscInt faceInd = off + d; 1370 const PetscInt Np = candidates[off + d].index; 1371 const PetscInt *join = NULL; 1372 PetscInt joinSize, points[1024], c; 1373 1374 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found claim for remote point (%" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, claims[faceInd].rank, claims[faceInd].index)); 1375 points[0] = r; 1376 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[0])); 1377 for (c = 0, d += 2; c < Np; ++c, ++d) { 1378 PetscCall(DMPlexMapToLocalPoint(dm, remoteHash, candidates[off + d], &points[c + 1], NULL)); 1379 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT "\n", rank, points[c + 1])); 1380 } 1381 PetscCall(DMPlexGetJoin(dm, Np + 1, points, &joinSize, &join)); 1382 if (joinSize == 1) { 1383 if (claims[faceInd].rank == rank) { 1384 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Ignoring local face %" PetscInt_FMT " for non-remote partner\n", rank, join[0])); 1385 } else { 1386 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found local face %" PetscInt_FMT "\n", rank, join[0])); 1387 PetscCall(PetscHMapISet(claimshash, join[0], faceInd)); 1388 } 1389 } else { 1390 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Failed to find face\n", rank)); 1391 } 1392 PetscCall(DMPlexRestoreJoin(dm, Np + 1, points, &joinSize, &join)); 1393 } else { 1394 if (debug) PetscCall(PetscSynchronizedPrintf(comm, "[%d] No claim for point %" PetscInt_FMT "\n", rank, r)); 1395 d += claims[off + d].index + 1; 1396 } 1397 } 1398 } 1399 if (debug) PetscCall(PetscSynchronizedFlush(comm, NULL)); 1400 /* Step 6) Create new pointSF from hashed claims */ 1401 PetscCall(PetscHMapIGetSize(claimshash, &NlNew)); 1402 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1403 PetscCall(PetscMalloc1(Nl + NlNew, &localPointsNew)); 1404 PetscCall(PetscMalloc1(Nl + NlNew, &remotePointsNew)); 1405 for (l = 0; l < Nl; ++l) { 1406 localPointsNew[l] = localPoints[l]; 1407 remotePointsNew[l].index = remotePoints[l].index; 1408 remotePointsNew[l].rank = remotePoints[l].rank; 1409 } 1410 p = Nl; 1411 PetscCall(PetscHMapIGetKeys(claimshash, &p, localPointsNew)); 1412 /* We sort new points, and assume they are numbered after all existing points */ 1413 PetscCall(PetscSortInt(NlNew, &localPointsNew[Nl])); 1414 for (p = Nl; p < Nl + NlNew; ++p) { 1415 PetscInt off; 1416 PetscCall(PetscHMapIGet(claimshash, localPointsNew[p], &off)); 1417 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); 1418 remotePointsNew[p] = claims[off]; 1419 } 1420 PetscCall(PetscSFCreate(comm, &sfPointNew)); 1421 PetscCall(PetscSFSetGraph(sfPointNew, pEnd - pStart, Nl + NlNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1422 PetscCall(PetscSFSetUp(sfPointNew)); 1423 PetscCall(DMSetPointSF(dm, sfPointNew)); 1424 PetscCall(PetscObjectViewFromOptions((PetscObject)sfPointNew, NULL, "-petscsf_interp_view")); 1425 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPointNew, PETSC_FALSE)); 1426 PetscCall(PetscSFDestroy(&sfPointNew)); 1427 PetscCall(PetscHMapIDestroy(&claimshash)); 1428 } 1429 PetscCall(PetscHMapIJDestroy(&remoteHash)); 1430 PetscCall(PetscSectionDestroy(&candidateSection)); 1431 PetscCall(PetscSectionDestroy(&candidateRemoteSection)); 1432 PetscCall(PetscSectionDestroy(&claimSection)); 1433 PetscCall(PetscFree(candidates)); 1434 PetscCall(PetscFree(candidatesRemote)); 1435 PetscCall(PetscFree(claims)); 1436 PetscCall(PetscLogEventEnd(DMPLEX_InterpolateSF, dm, 0, 0, 0)); 1437 PetscFunctionReturn(PETSC_SUCCESS); 1438 } 1439 1440 /*@ 1441 DMPlexInterpolate - Take in a cell-vertex mesh and return one with all intermediate faces, edges, etc. 1442 1443 Collective 1444 1445 Input Parameter: 1446 . dm - The `DMPLEX` object with only cells and vertices 1447 1448 Output Parameter: 1449 . dmInt - The complete `DMPLEX` object 1450 1451 Level: intermediate 1452 1453 Note: 1454 Labels and coordinates are copied. 1455 1456 Developer Note: 1457 It sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`. 1458 1459 .seealso: `DMPLEX`, `DMPlexUninterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()` 1460 @*/ 1461 PetscErrorCode DMPlexInterpolate(DM dm, DM *dmInt) 1462 { 1463 DMPlexInterpolatedFlag interpolated; 1464 DM idm, odm = dm; 1465 PetscSF sfPoint; 1466 PetscInt depth, dim, d; 1467 const char *name; 1468 PetscBool flg = PETSC_TRUE; 1469 1470 PetscFunctionBegin; 1471 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1472 PetscValidPointer(dmInt, 2); 1473 PetscCall(PetscLogEventBegin(DMPLEX_Interpolate, dm, 0, 0, 0)); 1474 PetscCall(DMPlexGetDepth(dm, &depth)); 1475 PetscCall(DMGetDimension(dm, &dim)); 1476 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1477 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1478 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 1479 PetscCall(PetscObjectReference((PetscObject)dm)); 1480 idm = dm; 1481 } else { 1482 for (d = 1; d < dim; ++d) { 1483 /* Create interpolated mesh */ 1484 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &idm)); 1485 PetscCall(DMSetType(idm, DMPLEX)); 1486 PetscCall(DMSetDimension(idm, dim)); 1487 if (depth > 0) { 1488 PetscCall(DMPlexInterpolateFaces_Internal(odm, 1, idm)); 1489 PetscCall(DMGetPointSF(odm, &sfPoint)); 1490 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(odm, sfPoint, PETSC_FALSE)); 1491 { 1492 /* TODO: We need to systematically fix cases of distributed Plexes with no graph set */ 1493 PetscInt nroots; 1494 PetscCall(PetscSFGetGraph(sfPoint, &nroots, NULL, NULL, NULL)); 1495 if (nroots >= 0) PetscCall(DMPlexInterpolatePointSF(idm, sfPoint)); 1496 } 1497 } 1498 if (odm != dm) PetscCall(DMDestroy(&odm)); 1499 odm = idm; 1500 } 1501 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 1502 PetscCall(PetscObjectSetName((PetscObject)idm, name)); 1503 PetscCall(DMPlexCopyCoordinates(dm, idm)); 1504 PetscCall(DMCopyLabels(dm, idm, PETSC_COPY_VALUES, PETSC_FALSE, DM_COPY_LABELS_FAIL)); 1505 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_interpolate_orient_interfaces", &flg, NULL)); 1506 if (flg) PetscCall(DMPlexOrientInterface_Internal(idm)); 1507 } 1508 /* This function makes the mesh fully interpolated on all ranks */ 1509 { 1510 DM_Plex *plex = (DM_Plex *)idm->data; 1511 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_FULL; 1512 } 1513 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, idm)); 1514 *dmInt = idm; 1515 PetscCall(PetscLogEventEnd(DMPLEX_Interpolate, dm, 0, 0, 0)); 1516 PetscFunctionReturn(PETSC_SUCCESS); 1517 } 1518 1519 /*@ 1520 DMPlexCopyCoordinates - Copy coordinates from one mesh to another with the same vertices 1521 1522 Collective 1523 1524 Input Parameter: 1525 . dmA - The `DMPLEX` object with initial coordinates 1526 1527 Output Parameter: 1528 . dmB - The `DMPLEX` object with copied coordinates 1529 1530 Level: intermediate 1531 1532 Note: 1533 This is typically used when adding pieces other than vertices to a mesh 1534 1535 .seealso: `DMPLEX`, `DMCopyLabels()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMGetCoordinateSection()` 1536 @*/ 1537 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB) 1538 { 1539 Vec coordinatesA, coordinatesB; 1540 VecType vtype; 1541 PetscSection coordSectionA, coordSectionB; 1542 PetscScalar *coordsA, *coordsB; 1543 PetscInt spaceDim, Nf, vStartA, vStartB, vEndA, vEndB, coordSizeB, v, d; 1544 PetscInt cStartA, cEndA, cStartB, cEndB, cS, cE, cdim; 1545 PetscBool lc = PETSC_FALSE; 1546 1547 PetscFunctionBegin; 1548 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 1549 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 1550 if (dmA == dmB) PetscFunctionReturn(PETSC_SUCCESS); 1551 PetscCall(DMGetCoordinateDim(dmA, &cdim)); 1552 PetscCall(DMSetCoordinateDim(dmB, cdim)); 1553 PetscCall(DMPlexGetDepthStratum(dmA, 0, &vStartA, &vEndA)); 1554 PetscCall(DMPlexGetDepthStratum(dmB, 0, &vStartB, &vEndB)); 1555 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); 1556 /* Copy over discretization if it exists */ 1557 { 1558 DM cdmA, cdmB; 1559 PetscDS dsA, dsB; 1560 PetscObject objA, objB; 1561 PetscClassId idA, idB; 1562 const PetscScalar *constants; 1563 PetscInt cdim, Nc; 1564 1565 PetscCall(DMGetCoordinateDM(dmA, &cdmA)); 1566 PetscCall(DMGetCoordinateDM(dmB, &cdmB)); 1567 PetscCall(DMGetField(cdmA, 0, NULL, &objA)); 1568 PetscCall(DMGetField(cdmB, 0, NULL, &objB)); 1569 PetscCall(PetscObjectGetClassId(objA, &idA)); 1570 PetscCall(PetscObjectGetClassId(objB, &idB)); 1571 if ((idA == PETSCFE_CLASSID) && (idA != idB)) { 1572 PetscCall(DMSetField(cdmB, 0, NULL, objA)); 1573 PetscCall(DMCreateDS(cdmB)); 1574 PetscCall(DMGetDS(cdmA, &dsA)); 1575 PetscCall(DMGetDS(cdmB, &dsB)); 1576 PetscCall(PetscDSGetCoordinateDimension(dsA, &cdim)); 1577 PetscCall(PetscDSSetCoordinateDimension(dsB, cdim)); 1578 PetscCall(PetscDSGetConstants(dsA, &Nc, &constants)); 1579 PetscCall(PetscDSSetConstants(dsB, Nc, (PetscScalar *)constants)); 1580 } 1581 } 1582 PetscCall(DMPlexGetHeightStratum(dmA, 0, &cStartA, &cEndA)); 1583 PetscCall(DMPlexGetHeightStratum(dmB, 0, &cStartB, &cEndB)); 1584 PetscCall(DMGetCoordinateSection(dmA, &coordSectionA)); 1585 PetscCall(DMGetCoordinateSection(dmB, &coordSectionB)); 1586 if (coordSectionA == coordSectionB) PetscFunctionReturn(PETSC_SUCCESS); 1587 PetscCall(PetscSectionGetNumFields(coordSectionA, &Nf)); 1588 if (!Nf) PetscFunctionReturn(PETSC_SUCCESS); 1589 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of coordinate fields must be 1, not %" PetscInt_FMT, Nf); 1590 if (!coordSectionB) { 1591 PetscInt dim; 1592 1593 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coordSectionA), &coordSectionB)); 1594 PetscCall(DMGetCoordinateDim(dmA, &dim)); 1595 PetscCall(DMSetCoordinateSection(dmB, dim, coordSectionB)); 1596 PetscCall(PetscObjectDereference((PetscObject)coordSectionB)); 1597 } 1598 PetscCall(PetscSectionSetNumFields(coordSectionB, 1)); 1599 PetscCall(PetscSectionGetFieldComponents(coordSectionA, 0, &spaceDim)); 1600 PetscCall(PetscSectionSetFieldComponents(coordSectionB, 0, spaceDim)); 1601 PetscCall(PetscSectionGetChart(coordSectionA, &cS, &cE)); 1602 if (cStartA <= cS && cS < cEndA) { /* localized coordinates */ 1603 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); 1604 cS = cS - cStartA + cStartB; 1605 cE = vEndB; 1606 lc = PETSC_TRUE; 1607 } else { 1608 cS = vStartB; 1609 cE = vEndB; 1610 } 1611 PetscCall(PetscSectionSetChart(coordSectionB, cS, cE)); 1612 for (v = vStartB; v < vEndB; ++v) { 1613 PetscCall(PetscSectionSetDof(coordSectionB, v, spaceDim)); 1614 PetscCall(PetscSectionSetFieldDof(coordSectionB, v, 0, spaceDim)); 1615 } 1616 if (lc) { /* localized coordinates */ 1617 PetscInt c; 1618 1619 for (c = cS - cStartB; c < cEndB - cStartB; c++) { 1620 PetscInt dof; 1621 1622 PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof)); 1623 PetscCall(PetscSectionSetDof(coordSectionB, c + cStartB, dof)); 1624 PetscCall(PetscSectionSetFieldDof(coordSectionB, c + cStartB, 0, dof)); 1625 } 1626 } 1627 PetscCall(PetscSectionSetUp(coordSectionB)); 1628 PetscCall(PetscSectionGetStorageSize(coordSectionB, &coordSizeB)); 1629 PetscCall(DMGetCoordinatesLocal(dmA, &coordinatesA)); 1630 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinatesB)); 1631 PetscCall(PetscObjectSetName((PetscObject)coordinatesB, "coordinates")); 1632 PetscCall(VecSetSizes(coordinatesB, coordSizeB, PETSC_DETERMINE)); 1633 PetscCall(VecGetBlockSize(coordinatesA, &d)); 1634 PetscCall(VecSetBlockSize(coordinatesB, d)); 1635 PetscCall(VecGetType(coordinatesA, &vtype)); 1636 PetscCall(VecSetType(coordinatesB, vtype)); 1637 PetscCall(VecGetArray(coordinatesA, &coordsA)); 1638 PetscCall(VecGetArray(coordinatesB, &coordsB)); 1639 for (v = 0; v < vEndB - vStartB; ++v) { 1640 PetscInt offA, offB; 1641 1642 PetscCall(PetscSectionGetOffset(coordSectionA, v + vStartA, &offA)); 1643 PetscCall(PetscSectionGetOffset(coordSectionB, v + vStartB, &offB)); 1644 for (d = 0; d < spaceDim; ++d) coordsB[offB + d] = coordsA[offA + d]; 1645 } 1646 if (lc) { /* localized coordinates */ 1647 PetscInt c; 1648 1649 for (c = cS - cStartB; c < cEndB - cStartB; c++) { 1650 PetscInt dof, offA, offB; 1651 1652 PetscCall(PetscSectionGetOffset(coordSectionA, c + cStartA, &offA)); 1653 PetscCall(PetscSectionGetOffset(coordSectionB, c + cStartB, &offB)); 1654 PetscCall(PetscSectionGetDof(coordSectionA, c + cStartA, &dof)); 1655 PetscCall(PetscArraycpy(coordsB + offB, coordsA + offA, dof)); 1656 } 1657 } 1658 PetscCall(VecRestoreArray(coordinatesA, &coordsA)); 1659 PetscCall(VecRestoreArray(coordinatesB, &coordsB)); 1660 PetscCall(DMSetCoordinatesLocal(dmB, coordinatesB)); 1661 PetscCall(VecDestroy(&coordinatesB)); 1662 PetscFunctionReturn(PETSC_SUCCESS); 1663 } 1664 1665 /*@ 1666 DMPlexUninterpolate - Take in a mesh with all intermediate faces, edges, etc. and return a cell-vertex mesh 1667 1668 Collective 1669 1670 Input Parameter: 1671 . dm - The complete `DMPLEX` object 1672 1673 Output Parameter: 1674 . dmUnint - The `DMPLEX` object with only cells and vertices 1675 1676 Level: intermediate 1677 1678 Note: 1679 It does not copy over the coordinates. 1680 1681 Developer Note: 1682 Sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`. 1683 1684 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexCreateFromCellListPetsc()`, `DMPlexCopyCoordinates()` 1685 @*/ 1686 PetscErrorCode DMPlexUninterpolate(DM dm, DM *dmUnint) 1687 { 1688 DMPlexInterpolatedFlag interpolated; 1689 DM udm; 1690 PetscInt dim, vStart, vEnd, cStart, cEnd, c, maxConeSize = 0, *cone; 1691 1692 PetscFunctionBegin; 1693 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1694 PetscValidPointer(dmUnint, 2); 1695 PetscCall(DMGetDimension(dm, &dim)); 1696 PetscCall(DMPlexIsInterpolated(dm, &interpolated)); 1697 PetscCheck(interpolated != DMPLEX_INTERPOLATED_PARTIAL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not for partially interpolated meshes"); 1698 if (interpolated == DMPLEX_INTERPOLATED_NONE || dim <= 1) { 1699 /* in case dim <= 1 just keep the DMPLEX_INTERPOLATED_FULL flag */ 1700 PetscCall(PetscObjectReference((PetscObject)dm)); 1701 *dmUnint = dm; 1702 PetscFunctionReturn(PETSC_SUCCESS); 1703 } 1704 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1705 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1706 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &udm)); 1707 PetscCall(DMSetType(udm, DMPLEX)); 1708 PetscCall(DMSetDimension(udm, dim)); 1709 PetscCall(DMPlexSetChart(udm, cStart, vEnd)); 1710 for (c = cStart; c < cEnd; ++c) { 1711 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1712 1713 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1714 for (cl = 0; cl < closureSize * 2; cl += 2) { 1715 const PetscInt p = closure[cl]; 1716 1717 if ((p >= vStart) && (p < vEnd)) ++coneSize; 1718 } 1719 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1720 PetscCall(DMPlexSetConeSize(udm, c, coneSize)); 1721 maxConeSize = PetscMax(maxConeSize, coneSize); 1722 } 1723 PetscCall(DMSetUp(udm)); 1724 PetscCall(PetscMalloc1(maxConeSize, &cone)); 1725 for (c = cStart; c < cEnd; ++c) { 1726 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 1727 1728 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1729 for (cl = 0; cl < closureSize * 2; cl += 2) { 1730 const PetscInt p = closure[cl]; 1731 1732 if ((p >= vStart) && (p < vEnd)) cone[coneSize++] = p; 1733 } 1734 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 1735 PetscCall(DMPlexSetCone(udm, c, cone)); 1736 } 1737 PetscCall(PetscFree(cone)); 1738 PetscCall(DMPlexSymmetrize(udm)); 1739 PetscCall(DMPlexStratify(udm)); 1740 /* Reduce SF */ 1741 { 1742 PetscSF sfPoint, sfPointUn; 1743 const PetscSFNode *remotePoints; 1744 const PetscInt *localPoints; 1745 PetscSFNode *remotePointsUn; 1746 PetscInt *localPointsUn; 1747 PetscInt numRoots, numLeaves, l; 1748 PetscInt numLeavesUn = 0, n = 0; 1749 1750 /* Get original SF information */ 1751 PetscCall(DMGetPointSF(dm, &sfPoint)); 1752 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(dm, sfPoint, PETSC_FALSE)); 1753 PetscCall(DMGetPointSF(udm, &sfPointUn)); 1754 PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1755 if (numRoots >= 0) { 1756 /* Allocate space for cells and vertices */ 1757 for (l = 0; l < numLeaves; ++l) { 1758 const PetscInt p = localPoints[l]; 1759 1760 if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) numLeavesUn++; 1761 } 1762 /* Fill in leaves */ 1763 PetscCall(PetscMalloc1(numLeavesUn, &remotePointsUn)); 1764 PetscCall(PetscMalloc1(numLeavesUn, &localPointsUn)); 1765 for (l = 0; l < numLeaves; l++) { 1766 const PetscInt p = localPoints[l]; 1767 1768 if ((vStart <= p && p < vEnd) || (cStart <= p && p < cEnd)) { 1769 localPointsUn[n] = p; 1770 remotePointsUn[n].rank = remotePoints[l].rank; 1771 remotePointsUn[n].index = remotePoints[l].index; 1772 ++n; 1773 } 1774 } 1775 PetscCheck(n == numLeavesUn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent number of leaves %" PetscInt_FMT " != %" PetscInt_FMT, n, numLeavesUn); 1776 PetscCall(PetscSFSetGraph(sfPointUn, cEnd - cStart + vEnd - vStart, numLeavesUn, localPointsUn, PETSC_OWN_POINTER, remotePointsUn, PETSC_OWN_POINTER)); 1777 } 1778 } 1779 /* This function makes the mesh fully uninterpolated on all ranks */ 1780 { 1781 DM_Plex *plex = (DM_Plex *)udm->data; 1782 plex->interpolated = plex->interpolatedCollective = DMPLEX_INTERPOLATED_NONE; 1783 } 1784 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, udm)); 1785 if (PetscDefined(USE_DEBUG)) PetscCall(DMPlexCheckPointSF(udm, NULL, PETSC_FALSE)); 1786 *dmUnint = udm; 1787 PetscFunctionReturn(PETSC_SUCCESS); 1788 } 1789 1790 static PetscErrorCode DMPlexIsInterpolated_Internal(DM dm, DMPlexInterpolatedFlag *interpolated) 1791 { 1792 PetscInt coneSize, depth, dim, h, p, pStart, pEnd; 1793 MPI_Comm comm; 1794 1795 PetscFunctionBegin; 1796 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1797 PetscCall(DMPlexGetDepth(dm, &depth)); 1798 PetscCall(DMGetDimension(dm, &dim)); 1799 1800 if (depth == dim) { 1801 *interpolated = DMPLEX_INTERPOLATED_FULL; 1802 if (!dim) goto finish; 1803 1804 /* Check points at height = dim are vertices (have no cones) */ 1805 PetscCall(DMPlexGetHeightStratum(dm, dim, &pStart, &pEnd)); 1806 for (p = pStart; p < pEnd; p++) { 1807 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1808 if (coneSize) { 1809 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1810 goto finish; 1811 } 1812 } 1813 1814 /* Check points at height < dim have cones */ 1815 for (h = 0; h < dim; h++) { 1816 PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 1817 for (p = pStart; p < pEnd; p++) { 1818 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 1819 if (!coneSize) { 1820 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1821 goto finish; 1822 } 1823 } 1824 } 1825 } else if (depth == 1) { 1826 *interpolated = DMPLEX_INTERPOLATED_NONE; 1827 } else { 1828 *interpolated = DMPLEX_INTERPOLATED_PARTIAL; 1829 } 1830 finish: 1831 PetscFunctionReturn(PETSC_SUCCESS); 1832 } 1833 1834 /*@ 1835 DMPlexIsInterpolated - Find out to what extent the `DMPLEX` is topologically interpolated. 1836 1837 Not Collective 1838 1839 Input Parameter: 1840 . dm - The `DMPLEX` object 1841 1842 Output Parameter: 1843 . interpolated - Flag whether the `DM` is interpolated 1844 1845 Level: intermediate 1846 1847 Notes: 1848 Unlike `DMPlexIsInterpolatedCollective()`, this is NOT collective 1849 so the results can be different on different ranks in special cases. 1850 However, `DMPlexInterpolate()` guarantees the result is the same on all. 1851 1852 Unlike `DMPlexIsInterpolatedCollective()`, this cannot return `DMPLEX_INTERPOLATED_MIXED`. 1853 1854 Developer Notes: 1855 Initially, plex->interpolated = `DMPLEX_INTERPOLATED_INVALID`. 1856 1857 If plex->interpolated == `DMPLEX_INTERPOLATED_INVALID`, `DMPlexIsInterpolated_Internal()` is called. 1858 It checks the actual topology and sets plex->interpolated on each rank separately to one of 1859 `DMPLEX_INTERPOLATED_NONE`, `DMPLEX_INTERPOLATED_PARTIAL` or `DMPLEX_INTERPOLATED_FULL`. 1860 1861 If plex->interpolated != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolated. 1862 1863 `DMPlexInterpolate()` sets plex->interpolated = `DMPLEX_INTERPOLATED_FULL`, 1864 and DMPlexUninterpolate() sets plex->interpolated = `DMPLEX_INTERPOLATED_NONE`. 1865 1866 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolatedCollective()` 1867 @*/ 1868 PetscErrorCode DMPlexIsInterpolated(DM dm, DMPlexInterpolatedFlag *interpolated) 1869 { 1870 DM_Plex *plex = (DM_Plex *)dm->data; 1871 1872 PetscFunctionBegin; 1873 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1874 PetscValidPointer(interpolated, 2); 1875 if (plex->interpolated < 0) { 1876 PetscCall(DMPlexIsInterpolated_Internal(dm, &plex->interpolated)); 1877 } else if (PetscDefined(USE_DEBUG)) { 1878 DMPlexInterpolatedFlag flg; 1879 1880 PetscCall(DMPlexIsInterpolated_Internal(dm, &flg)); 1881 PetscCheck(flg == plex->interpolated, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stashed DMPlexInterpolatedFlag %s is inconsistent with current %s", DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[flg]); 1882 } 1883 *interpolated = plex->interpolated; 1884 PetscFunctionReturn(PETSC_SUCCESS); 1885 } 1886 1887 /*@ 1888 DMPlexIsInterpolatedCollective - Find out to what extent the `DMPLEX` is topologically interpolated (in collective manner). 1889 1890 Collective 1891 1892 Input Parameter: 1893 . dm - The `DMPLEX` object 1894 1895 Output Parameter: 1896 . interpolated - Flag whether the `DM` is interpolated 1897 1898 Level: intermediate 1899 1900 Notes: 1901 Unlike `DMPlexIsInterpolated()`, this is collective so the results are guaranteed to be the same on all ranks. 1902 1903 This function will return `DMPLEX_INTERPOLATED_MIXED` if the results of `DMPlexIsInterpolated()` are different on different ranks. 1904 1905 Developer Notes: 1906 Initially, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_INVALID`. 1907 1908 If plex->interpolatedCollective == `DMPLEX_INTERPOLATED_INVALID`, this function calls `DMPlexIsInterpolated()` which sets plex->interpolated. 1909 `MPI_Allreduce()` is then called and collectively consistent flag plex->interpolatedCollective is set and returned; 1910 if plex->interpolated varies on different ranks, plex->interpolatedCollective = `DMPLEX_INTERPOLATED_MIXED`, 1911 otherwise sets plex->interpolatedCollective = plex->interpolated. 1912 1913 If plex->interpolatedCollective != `DMPLEX_INTERPOLATED_INVALID`, this function just returns plex->interpolatedCollective. 1914 1915 .seealso: `DMPLEX`, `DMPlexInterpolate()`, `DMPlexIsInterpolated()` 1916 @*/ 1917 PetscErrorCode DMPlexIsInterpolatedCollective(DM dm, DMPlexInterpolatedFlag *interpolated) 1918 { 1919 DM_Plex *plex = (DM_Plex *)dm->data; 1920 PetscBool debug = PETSC_FALSE; 1921 1922 PetscFunctionBegin; 1923 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1924 PetscValidPointer(interpolated, 2); 1925 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_is_interpolated_collective_debug", &debug, NULL)); 1926 if (plex->interpolatedCollective < 0) { 1927 DMPlexInterpolatedFlag min, max; 1928 MPI_Comm comm; 1929 1930 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1931 PetscCall(DMPlexIsInterpolated(dm, &plex->interpolatedCollective)); 1932 PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &min, 1, MPIU_ENUM, MPI_MIN, comm)); 1933 PetscCall(MPIU_Allreduce(&plex->interpolatedCollective, &max, 1, MPIU_ENUM, MPI_MAX, comm)); 1934 if (min != max) plex->interpolatedCollective = DMPLEX_INTERPOLATED_MIXED; 1935 if (debug) { 1936 PetscMPIInt rank; 1937 1938 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1939 PetscCall(PetscSynchronizedPrintf(comm, "[%d] interpolated=%s interpolatedCollective=%s\n", rank, DMPlexInterpolatedFlags[plex->interpolated], DMPlexInterpolatedFlags[plex->interpolatedCollective])); 1940 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 1941 } 1942 } 1943 *interpolated = plex->interpolatedCollective; 1944 PetscFunctionReturn(PETSC_SUCCESS); 1945 } 1946