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