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