1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/sfimpl.h> /*I "petscsf.h" I*/ 3 #include <petsc/private/viewercgnsimpl.h> 4 #include <petsc/private/hashseti.h> 5 6 #include <pcgnslib.h> 7 #include <cgns_io.h> 8 9 #if !defined(CGNS_ENUMT) 10 #define CGNS_ENUMT(a) a 11 #endif 12 #if !defined(CGNS_ENUMV) 13 #define CGNS_ENUMV(a) a 14 #endif 15 // Permute plex closure ordering to CGNS 16 static PetscErrorCode DMPlexCGNSGetPermutation_Internal(DMPolytopeType cell_type, PetscInt closure_size, CGNS_ENUMT(ElementType_t) * element_type, const int **perm) 17 { 18 CGNS_ENUMT(ElementType_t) element_type_tmp; 19 20 // https://cgns.github.io/CGNS_docs_current/sids/conv.html#unstructgrid 21 static const int bar_2[2] = {0, 1}; 22 static const int bar_3[3] = {1, 2, 0}; 23 static const int bar_4[4] = {2, 3, 0, 1}; 24 static const int bar_5[5] = {3, 4, 0, 1, 2}; 25 static const int tri_3[3] = {0, 1, 2}; 26 static const int tri_6[6] = {3, 4, 5, 0, 1, 2}; 27 static const int tri_10[10] = {7, 8, 9, 1, 2, 3, 4, 5, 6, 0}; 28 static const int quad_4[4] = {0, 1, 2, 3}; 29 static const int quad_9[9] = { 30 5, 6, 7, 8, // vertices 31 1, 2, 3, 4, // edges 32 0, // center 33 }; 34 static const int quad_16[] = { 35 12, 13, 14, 15, // vertices 36 4, 5, 6, 7, 8, 9, 10, 11, // edges 37 0, 1, 3, 2, // centers 38 }; 39 static const int quad_25[] = { 40 21, 22, 23, 24, // vertices 41 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, // edges 42 0, 1, 2, 5, 8, 7, 6, 3, 4, // centers 43 }; 44 static const int tetra_4[4] = {0, 2, 1, 3}; 45 static const int tetra_10[10] = {6, 8, 7, 9, 2, 1, 0, 3, 5, 4}; 46 static const int tetra_20[20] = { 47 16, 18, 17, 19, // vertices 48 9, 8, 7, 6, 5, 4, // bottom edges 49 10, 11, 14, 15, 13, 12, // side edges 50 0, 2, 3, 1, // faces 51 }; 52 static const int hexa_8[8] = {0, 3, 2, 1, 4, 5, 6, 7}; 53 static const int hexa_27[27] = { 54 19, 22, 21, 20, 23, 24, 25, 26, // vertices 55 10, 9, 8, 7, // bottom edges 56 16, 15, 18, 17, // mid edges 57 11, 12, 13, 14, // top edges 58 1, 3, 5, 4, 6, 2, // faces 59 0, // center 60 }; 61 static const int hexa_64[64] = { 62 // debug with $PETSC_ARCH/tests/dm/impls/plex/tests/ex49 -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -dm_coord_petscspace_degree 3 63 56, 59, 58, 57, 60, 61, 62, 63, // vertices 64 39, 38, 37, 36, 35, 34, 33, 32, // bottom edges 65 51, 50, 48, 49, 52, 53, 55, 54, // mid edges; Paraview needs edge 21-22 swapped with 23-24 66 40, 41, 42, 43, 44, 45, 46, 47, // top edges 67 8, 10, 11, 9, // z-minus face 68 16, 17, 19, 18, // y-minus face 69 24, 25, 27, 26, // x-plus face 70 20, 21, 23, 22, // y-plus face 71 30, 28, 29, 31, // x-minus face 72 12, 13, 15, 14, // z-plus face 73 0, 1, 3, 2, 4, 5, 7, 6, // center 74 }; 75 76 PetscFunctionBegin; 77 element_type_tmp = CGNS_ENUMV(ElementTypeNull); 78 *perm = NULL; 79 switch (cell_type) { 80 case DM_POLYTOPE_SEGMENT: 81 switch (closure_size) { 82 case 2: 83 element_type_tmp = CGNS_ENUMV(BAR_2); 84 *perm = bar_2; 85 break; 86 case 3: 87 element_type_tmp = CGNS_ENUMV(BAR_3); 88 *perm = bar_3; 89 break; 90 case 4: 91 element_type_tmp = CGNS_ENUMV(BAR_4); 92 *perm = bar_4; 93 break; 94 case 5: 95 element_type_tmp = CGNS_ENUMV(BAR_5); 96 *perm = bar_5; 97 break; 98 default: 99 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 100 } 101 break; 102 case DM_POLYTOPE_TRIANGLE: 103 switch (closure_size) { 104 case 3: 105 element_type_tmp = CGNS_ENUMV(TRI_3); 106 *perm = tri_3; 107 break; 108 case 6: 109 element_type_tmp = CGNS_ENUMV(TRI_6); 110 *perm = tri_6; 111 break; 112 case 10: 113 element_type_tmp = CGNS_ENUMV(TRI_10); 114 *perm = tri_10; 115 break; 116 default: 117 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 118 } 119 break; 120 case DM_POLYTOPE_QUADRILATERAL: 121 switch (closure_size) { 122 case 4: 123 element_type_tmp = CGNS_ENUMV(QUAD_4); 124 *perm = quad_4; 125 break; 126 case 9: 127 element_type_tmp = CGNS_ENUMV(QUAD_9); 128 *perm = quad_9; 129 break; 130 case 16: 131 element_type_tmp = CGNS_ENUMV(QUAD_16); 132 *perm = quad_16; 133 break; 134 case 25: 135 element_type_tmp = CGNS_ENUMV(QUAD_25); 136 *perm = quad_25; 137 break; 138 default: 139 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 140 } 141 break; 142 case DM_POLYTOPE_TETRAHEDRON: 143 switch (closure_size) { 144 case 4: 145 element_type_tmp = CGNS_ENUMV(TETRA_4); 146 *perm = tetra_4; 147 break; 148 case 10: 149 element_type_tmp = CGNS_ENUMV(TETRA_10); 150 *perm = tetra_10; 151 break; 152 case 20: 153 element_type_tmp = CGNS_ENUMV(TETRA_20); 154 *perm = tetra_20; 155 break; 156 default: 157 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 158 } 159 break; 160 case DM_POLYTOPE_HEXAHEDRON: 161 switch (closure_size) { 162 case 8: 163 element_type_tmp = CGNS_ENUMV(HEXA_8); 164 *perm = hexa_8; 165 break; 166 case 27: 167 element_type_tmp = CGNS_ENUMV(HEXA_27); 168 *perm = hexa_27; 169 break; 170 case 64: 171 element_type_tmp = CGNS_ENUMV(HEXA_64); 172 *perm = hexa_64; 173 break; 174 default: 175 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 176 } 177 break; 178 default: 179 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 180 } 181 if (element_type) *element_type = element_type_tmp; 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 185 /* 186 Input Parameters: 187 + cellType - The CGNS-defined element type 188 189 Output Parameters: 190 + dmcelltype - The equivalent DMPolytopeType for the cellType 191 . numCorners - Number of corners of the polytope 192 - dim - The topological dimension of the polytope 193 194 CGNS elements defined in: https://cgns.github.io/CGNS_docs_current/sids/conv.html#unstructgrid 195 */ 196 static inline PetscErrorCode CGNSElementTypeGetTopologyInfo(CGNS_ENUMT(ElementType_t) cellType, DMPolytopeType *dmcelltype, PetscInt *numCorners, PetscInt *dim) 197 { 198 DMPolytopeType _dmcelltype; 199 200 PetscFunctionBegin; 201 switch (cellType) { 202 case CGNS_ENUMV(BAR_2): 203 case CGNS_ENUMV(BAR_3): 204 case CGNS_ENUMV(BAR_4): 205 case CGNS_ENUMV(BAR_5): 206 _dmcelltype = DM_POLYTOPE_SEGMENT; 207 break; 208 case CGNS_ENUMV(TRI_3): 209 case CGNS_ENUMV(TRI_6): 210 case CGNS_ENUMV(TRI_9): 211 case CGNS_ENUMV(TRI_10): 212 case CGNS_ENUMV(TRI_12): 213 case CGNS_ENUMV(TRI_15): 214 _dmcelltype = DM_POLYTOPE_TRIANGLE; 215 break; 216 case CGNS_ENUMV(QUAD_4): 217 case CGNS_ENUMV(QUAD_8): 218 case CGNS_ENUMV(QUAD_9): 219 case CGNS_ENUMV(QUAD_12): 220 case CGNS_ENUMV(QUAD_16): 221 case CGNS_ENUMV(QUAD_P4_16): 222 case CGNS_ENUMV(QUAD_25): 223 _dmcelltype = DM_POLYTOPE_QUADRILATERAL; 224 break; 225 case CGNS_ENUMV(TETRA_4): 226 case CGNS_ENUMV(TETRA_10): 227 case CGNS_ENUMV(TETRA_16): 228 case CGNS_ENUMV(TETRA_20): 229 case CGNS_ENUMV(TETRA_22): 230 case CGNS_ENUMV(TETRA_34): 231 case CGNS_ENUMV(TETRA_35): 232 _dmcelltype = DM_POLYTOPE_TETRAHEDRON; 233 break; 234 case CGNS_ENUMV(PYRA_5): 235 case CGNS_ENUMV(PYRA_13): 236 case CGNS_ENUMV(PYRA_14): 237 case CGNS_ENUMV(PYRA_21): 238 case CGNS_ENUMV(PYRA_29): 239 case CGNS_ENUMV(PYRA_P4_29): 240 case CGNS_ENUMV(PYRA_30): 241 case CGNS_ENUMV(PYRA_50): 242 case CGNS_ENUMV(PYRA_55): 243 _dmcelltype = DM_POLYTOPE_PYRAMID; 244 break; 245 case CGNS_ENUMV(PENTA_6): 246 case CGNS_ENUMV(PENTA_15): 247 case CGNS_ENUMV(PENTA_18): 248 case CGNS_ENUMV(PENTA_24): 249 case CGNS_ENUMV(PENTA_33): 250 case CGNS_ENUMV(PENTA_38): 251 case CGNS_ENUMV(PENTA_40): 252 case CGNS_ENUMV(PENTA_66): 253 case CGNS_ENUMV(PENTA_75): 254 _dmcelltype = DM_POLYTOPE_TRI_PRISM; 255 break; 256 case CGNS_ENUMV(HEXA_8): 257 case CGNS_ENUMV(HEXA_20): 258 case CGNS_ENUMV(HEXA_27): 259 case CGNS_ENUMV(HEXA_32): 260 case CGNS_ENUMV(HEXA_44): 261 case CGNS_ENUMV(HEXA_56): 262 case CGNS_ENUMV(HEXA_64): 263 case CGNS_ENUMV(HEXA_98): 264 case CGNS_ENUMV(HEXA_125): 265 _dmcelltype = DM_POLYTOPE_HEXAHEDRON; 266 break; 267 case CGNS_ENUMV(MIXED): 268 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid CGNS ElementType_t: MIXED"); 269 default: 270 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid CGNS ElementType_t: %d", (int)cellType); 271 } 272 273 if (dmcelltype) *dmcelltype = _dmcelltype; 274 if (numCorners) *numCorners = DMPolytopeTypeGetNumVertices(_dmcelltype); 275 if (dim) *dim = DMPolytopeTypeGetDim(_dmcelltype); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 /* 280 Input Parameters: 281 + cellType - The CGNS-defined cell type 282 283 Output Parameters: 284 + numClosure - Number of nodes that define the function space on the cell 285 - pOrder - The polynomial order of the cell 286 287 CGNS elements defined in: https://cgns.github.io/CGNS_docs_current/sids/conv.html#unstructgrid 288 289 Note: we only support "full" elements, ie. not seredipity elements 290 */ 291 static inline PetscErrorCode CGNSElementTypeGetDiscretizationInfo(CGNS_ENUMT(ElementType_t) cellType, PetscInt *numClosure, PetscInt *pOrder) 292 { 293 PetscInt _numClosure, _pOrder; 294 295 PetscFunctionBegin; 296 switch (cellType) { 297 case CGNS_ENUMV(BAR_2): 298 _numClosure = 2; 299 _pOrder = 1; 300 break; 301 case CGNS_ENUMV(BAR_3): 302 _numClosure = 3; 303 _pOrder = 2; 304 break; 305 case CGNS_ENUMV(BAR_4): 306 _numClosure = 4; 307 _pOrder = 3; 308 break; 309 case CGNS_ENUMV(BAR_5): 310 _numClosure = 5; 311 _pOrder = 4; 312 break; 313 case CGNS_ENUMV(TRI_3): 314 _numClosure = 3; 315 _pOrder = 1; 316 break; 317 case CGNS_ENUMV(TRI_6): 318 _numClosure = 6; 319 _pOrder = 2; 320 break; 321 case CGNS_ENUMV(TRI_10): 322 _numClosure = 10; 323 _pOrder = 3; 324 break; 325 case CGNS_ENUMV(TRI_15): 326 _numClosure = 15; 327 _pOrder = 4; 328 break; 329 case CGNS_ENUMV(QUAD_4): 330 _numClosure = 4; 331 _pOrder = 1; 332 break; 333 case CGNS_ENUMV(QUAD_9): 334 _numClosure = 9; 335 _pOrder = 2; 336 break; 337 case CGNS_ENUMV(QUAD_16): 338 _numClosure = 16; 339 _pOrder = 3; 340 break; 341 case CGNS_ENUMV(QUAD_25): 342 _numClosure = 25; 343 _pOrder = 4; 344 break; 345 case CGNS_ENUMV(TETRA_4): 346 _numClosure = 4; 347 _pOrder = 1; 348 break; 349 case CGNS_ENUMV(TETRA_10): 350 _numClosure = 10; 351 _pOrder = 2; 352 break; 353 case CGNS_ENUMV(TETRA_20): 354 _numClosure = 20; 355 _pOrder = 3; 356 break; 357 case CGNS_ENUMV(TETRA_35): 358 _numClosure = 35; 359 _pOrder = 4; 360 break; 361 case CGNS_ENUMV(PYRA_5): 362 _numClosure = 5; 363 _pOrder = 1; 364 break; 365 case CGNS_ENUMV(PYRA_14): 366 _numClosure = 14; 367 _pOrder = 2; 368 break; 369 case CGNS_ENUMV(PYRA_30): 370 _numClosure = 30; 371 _pOrder = 3; 372 break; 373 case CGNS_ENUMV(PYRA_55): 374 _numClosure = 55; 375 _pOrder = 4; 376 break; 377 case CGNS_ENUMV(PENTA_6): 378 _numClosure = 6; 379 _pOrder = 1; 380 break; 381 case CGNS_ENUMV(PENTA_18): 382 _numClosure = 18; 383 _pOrder = 2; 384 break; 385 case CGNS_ENUMV(PENTA_40): 386 _numClosure = 40; 387 _pOrder = 3; 388 break; 389 case CGNS_ENUMV(PENTA_75): 390 _numClosure = 75; 391 _pOrder = 4; 392 break; 393 case CGNS_ENUMV(HEXA_8): 394 _numClosure = 8; 395 _pOrder = 1; 396 break; 397 case CGNS_ENUMV(HEXA_27): 398 _numClosure = 27; 399 _pOrder = 2; 400 break; 401 case CGNS_ENUMV(HEXA_64): 402 _numClosure = 64; 403 _pOrder = 3; 404 break; 405 case CGNS_ENUMV(HEXA_125): 406 _numClosure = 125; 407 _pOrder = 4; 408 break; 409 case CGNS_ENUMV(MIXED): 410 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid CGNS ElementType_t: MIXED"); 411 default: 412 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported or Invalid cell type %d", (int)cellType); 413 } 414 if (numClosure) *numClosure = _numClosure; 415 if (pOrder) *pOrder = _pOrder; 416 PetscFunctionReturn(PETSC_SUCCESS); 417 } 418 419 static PetscErrorCode PetscCGNSDataType(PetscDataType pd, CGNS_ENUMT(DataType_t) * cd) 420 { 421 PetscFunctionBegin; 422 switch (pd) { 423 case PETSC_FLOAT: 424 *cd = CGNS_ENUMV(RealSingle); 425 break; 426 case PETSC_DOUBLE: 427 *cd = CGNS_ENUMV(RealDouble); 428 break; 429 case PETSC_COMPLEX: 430 *cd = CGNS_ENUMV(ComplexDouble); 431 break; 432 default: 433 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Data type %s", PetscDataTypes[pd]); 434 } 435 PetscFunctionReturn(PETSC_SUCCESS); 436 } 437 438 PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 439 { 440 int cgid = -1; 441 PetscBool use_parallel_viewer = PETSC_FALSE; 442 443 PetscFunctionBegin; 444 PetscAssertPointer(filename, 2); 445 PetscCall(PetscViewerCGNSRegisterLogEvents_Internal()); 446 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_cgns_parallel", &use_parallel_viewer, NULL)); 447 448 if (use_parallel_viewer) { 449 PetscCallCGNS(cgp_mpi_comm(comm)); 450 PetscCallCGNSOpen(cgp_open(filename, CG_MODE_READ, &cgid), 0, 0); 451 PetscCheck(cgid > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "cgp_open(\"%s\",...) did not return a valid file ID", filename); 452 PetscCall(DMPlexCreateCGNS(comm, cgid, interpolate, dm)); 453 PetscCallCGNSClose(cgp_close(cgid), 0, 0); 454 } else { 455 PetscCallCGNSOpen(cg_open(filename, CG_MODE_READ, &cgid), 0, 0); 456 PetscCheck(cgid > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename); 457 PetscCall(DMPlexCreateCGNS(comm, cgid, interpolate, dm)); 458 PetscCallCGNSClose(cg_close(cgid), 0, 0); 459 } 460 PetscFunctionReturn(PETSC_SUCCESS); 461 } 462 463 PetscErrorCode DMPlexCreateCGNS_Internal_Serial(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm) 464 { 465 PetscMPIInt num_proc, rank; 466 DM cdm; 467 DMLabel label; 468 PetscSection coordSection; 469 Vec coordinates; 470 PetscScalar *coords; 471 PetscInt *cellStart, *vertStart, v; 472 PetscInt labelIdRange[2], labelId; 473 /* Read from file */ 474 char basename[CGIO_MAX_NAME_LENGTH + 1]; 475 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 476 int dim = 0, physDim = 0, coordDim = 0, numVertices = 0, numCells = 0; 477 int nzones = 0; 478 const int B = 1; // Only support single base 479 480 PetscFunctionBegin; 481 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 482 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 483 PetscCall(DMCreate(comm, dm)); 484 PetscCall(DMSetType(*dm, DMPLEX)); 485 486 /* Open CGNS II file and read basic information on rank 0, then broadcast to all processors */ 487 if (rank == 0) { 488 int nbases, z; 489 490 PetscCallCGNSRead(cg_nbases(cgid, &nbases), *dm, 0); 491 PetscCheck(nbases <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single base, not %d", nbases); 492 PetscCallCGNSRead(cg_base_read(cgid, B, basename, &dim, &physDim), *dm, 0); 493 PetscCallCGNSRead(cg_nzones(cgid, B, &nzones), *dm, 0); 494 PetscCall(PetscCalloc2(nzones + 1, &cellStart, nzones + 1, &vertStart)); 495 for (z = 1; z <= nzones; ++z) { 496 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 497 498 PetscCallCGNSRead(cg_zone_read(cgid, B, z, buffer, sizes), *dm, 0); 499 numVertices += sizes[0]; 500 numCells += sizes[1]; 501 cellStart[z] += sizes[1] + cellStart[z - 1]; 502 vertStart[z] += sizes[0] + vertStart[z - 1]; 503 } 504 for (z = 1; z <= nzones; ++z) vertStart[z] += numCells; 505 coordDim = dim; 506 } 507 PetscCallMPI(MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH + 1, MPI_CHAR, 0, comm)); 508 PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm)); 509 PetscCallMPI(MPI_Bcast(&coordDim, 1, MPI_INT, 0, comm)); 510 PetscCallMPI(MPI_Bcast(&nzones, 1, MPI_INT, 0, comm)); 511 512 PetscCall(PetscObjectSetName((PetscObject)*dm, basename)); 513 PetscCall(DMSetDimension(*dm, dim)); 514 PetscCall(DMCreateLabel(*dm, "celltype")); 515 PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices)); 516 517 /* Read zone information */ 518 if (rank == 0) { 519 int z, c, c_loc; 520 521 /* Read the cell set connectivity table and build mesh topology 522 CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */ 523 /* First set sizes */ 524 for (z = 1, c = 0; z <= nzones; ++z) { 525 CGNS_ENUMT(ZoneType_t) zonetype; 526 int nsections; 527 CGNS_ENUMT(ElementType_t) cellType; 528 cgsize_t start, end; 529 int nbndry, parentFlag; 530 PetscInt numCorners, pOrder; 531 DMPolytopeType ctype; 532 const int S = 1; // Only support single section 533 534 PetscCallCGNSRead(cg_zone_type(cgid, B, z, &zonetype), *dm, 0); 535 PetscCheck(zonetype != CGNS_ENUMV(Structured), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can only handle Unstructured zones for CGNS"); 536 PetscCallCGNSRead(cg_nsections(cgid, B, z, &nsections), *dm, 0); 537 PetscCheck(nsections <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single section, not %d", nsections); 538 PetscCallCGNSRead(cg_section_read(cgid, B, z, S, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 539 if (cellType == CGNS_ENUMV(MIXED)) { 540 cgsize_t elementDataSize, *elements; 541 PetscInt off; 542 543 PetscCallCGNSRead(cg_ElementDataSize(cgid, B, z, S, &elementDataSize), *dm, 0); 544 PetscCall(PetscMalloc1(elementDataSize, &elements)); 545 PetscCallCGNSReadData(cg_poly_elements_read(cgid, B, z, S, elements, NULL, NULL), *dm, 0); 546 for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) { 547 PetscCall(CGNSElementTypeGetTopologyInfo((CGNS_ENUMT(ElementType_t))elements[off], &ctype, &numCorners, NULL)); 548 PetscCall(CGNSElementTypeGetDiscretizationInfo((CGNS_ENUMT(ElementType_t))elements[off], NULL, &pOrder)); 549 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 550 PetscCall(DMPlexSetConeSize(*dm, c, numCorners)); 551 PetscCall(DMPlexSetCellType(*dm, c, ctype)); 552 off += numCorners + 1; 553 } 554 PetscCall(PetscFree(elements)); 555 } else { 556 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, &ctype, &numCorners, NULL)); 557 PetscCall(CGNSElementTypeGetDiscretizationInfo(cellType, NULL, &pOrder)); 558 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 559 for (c_loc = start; c_loc <= end; ++c_loc, ++c) { 560 PetscCall(DMPlexSetConeSize(*dm, c, numCorners)); 561 PetscCall(DMPlexSetCellType(*dm, c, ctype)); 562 } 563 } 564 } 565 for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 566 } 567 568 PetscCall(DMSetUp(*dm)); 569 570 PetscCall(DMCreateLabel(*dm, "zone")); 571 if (rank == 0) { 572 int z, c, c_loc, v_loc; 573 574 PetscCall(DMGetLabel(*dm, "zone", &label)); 575 for (z = 1, c = 0; z <= nzones; ++z) { 576 CGNS_ENUMT(ElementType_t) cellType; 577 cgsize_t elementDataSize, *elements, start, end; 578 int nbndry, parentFlag; 579 PetscInt *cone, numc, numCorners, maxCorners = 27, pOrder; 580 const int S = 1; // Only support single section 581 582 PetscCallCGNSRead(cg_section_read(cgid, B, z, S, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 583 numc = end - start; 584 PetscCallCGNSRead(cg_ElementDataSize(cgid, B, z, S, &elementDataSize), *dm, 0); 585 PetscCall(PetscMalloc2(elementDataSize, &elements, maxCorners, &cone)); 586 PetscCallCGNSReadData(cg_poly_elements_read(cgid, B, z, S, elements, NULL, NULL), *dm, 0); 587 if (cellType == CGNS_ENUMV(MIXED)) { 588 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 589 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 590 PetscCall(CGNSElementTypeGetTopologyInfo((CGNS_ENUMT(ElementType_t))elements[v], NULL, &numCorners, NULL)); 591 PetscCall(CGNSElementTypeGetDiscretizationInfo((CGNS_ENUMT(ElementType_t))elements[v], NULL, &pOrder)); 592 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 593 ++v; 594 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) cone[v_loc] = elements[v] + numCells - 1; 595 PetscCall(DMPlexReorderCell(*dm, c, cone)); 596 PetscCall(DMPlexSetCone(*dm, c, cone)); 597 PetscCall(DMLabelSetValue(label, c, z)); 598 } 599 } else { 600 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, NULL, &numCorners, NULL)); 601 PetscCall(CGNSElementTypeGetDiscretizationInfo(cellType, NULL, &pOrder)); 602 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 603 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 604 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 605 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) cone[v_loc] = elements[v] + numCells - 1; 606 PetscCall(DMPlexReorderCell(*dm, c, cone)); 607 PetscCall(DMPlexSetCone(*dm, c, cone)); 608 PetscCall(DMLabelSetValue(label, c, z)); 609 } 610 } 611 PetscCall(PetscFree2(elements, cone)); 612 } 613 } 614 615 PetscCall(DMPlexSymmetrize(*dm)); 616 PetscCall(DMPlexStratify(*dm)); 617 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(*dm)); 618 619 /* Read coordinates */ 620 PetscCall(DMSetCoordinateDim(*dm, coordDim)); 621 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 622 PetscCall(DMGetLocalSection(cdm, &coordSection)); 623 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 624 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, coordDim)); 625 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 626 for (v = numCells; v < numCells + numVertices; ++v) { 627 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 628 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, coordDim)); 629 } 630 PetscCall(PetscSectionSetUp(coordSection)); 631 632 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 633 PetscCall(VecGetArray(coordinates, &coords)); 634 if (rank == 0) { 635 PetscInt off = 0; 636 float *x[3]; 637 int z, d; 638 639 PetscCall(PetscMalloc3(numVertices, &x[0], numVertices, &x[1], numVertices, &x[2])); 640 for (z = 1; z <= nzones; ++z) { 641 CGNS_ENUMT(DataType_t) datatype; 642 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 643 cgsize_t range_min[3] = {1, 1, 1}; 644 cgsize_t range_max[3] = {1, 1, 1}; 645 int ngrids, ncoords; 646 647 PetscCallCGNSRead(cg_zone_read(cgid, B, z, buffer, sizes), *dm, 0); 648 range_max[0] = sizes[0]; 649 PetscCallCGNSRead(cg_ngrids(cgid, B, z, &ngrids), *dm, 0); 650 PetscCheck(ngrids <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single grid, not %d", ngrids); 651 PetscCallCGNSRead(cg_ncoords(cgid, B, z, &ncoords), *dm, 0); 652 PetscCheck(ncoords == coordDim, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a coordinate array for each dimension, not %d", ncoords); 653 for (d = 0; d < coordDim; ++d) { 654 PetscCallCGNSRead(cg_coord_info(cgid, B, z, 1 + d, &datatype, buffer), *dm, 0); 655 PetscCallCGNSReadData(cg_coord_read(cgid, B, z, buffer, CGNS_ENUMV(RealSingle), range_min, range_max, x[d]), *dm, 0); 656 } 657 if (coordDim >= 1) { 658 for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 0] = x[0][v]; 659 } 660 if (coordDim >= 2) { 661 for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 1] = x[1][v]; 662 } 663 if (coordDim >= 3) { 664 for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 2] = x[2][v]; 665 } 666 off += sizes[0]; 667 } 668 PetscCall(PetscFree3(x[0], x[1], x[2])); 669 } 670 PetscCall(VecRestoreArray(coordinates, &coords)); 671 672 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 673 PetscCall(VecSetBlockSize(coordinates, coordDim)); 674 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 675 PetscCall(VecDestroy(&coordinates)); 676 677 /* Read boundary conditions */ 678 PetscCall(DMGetNumLabels(*dm, &labelIdRange[0])); 679 if (rank == 0) { 680 CGNS_ENUMT(BCType_t) bctype; 681 CGNS_ENUMT(DataType_t) datatype; 682 CGNS_ENUMT(PointSetType_t) pointtype; 683 cgsize_t *points; 684 PetscReal *normals; 685 int normal[3]; 686 char *bcname = buffer; 687 cgsize_t npoints, nnormals; 688 int z, nbc, bc, c, ndatasets; 689 690 for (z = 1; z <= nzones; ++z) { 691 PetscCallCGNSRead(cg_nbocos(cgid, B, z, &nbc), *dm, 0); 692 for (bc = 1; bc <= nbc; ++bc) { 693 PetscCallCGNSRead(cg_boco_info(cgid, B, z, bc, bcname, &bctype, &pointtype, &npoints, normal, &nnormals, &datatype, &ndatasets), *dm, 0); 694 PetscCall(DMCreateLabel(*dm, bcname)); 695 PetscCall(DMGetLabel(*dm, bcname, &label)); 696 PetscCall(PetscMalloc2(npoints, &points, nnormals, &normals)); 697 PetscCallCGNSReadData(cg_boco_read(cgid, B, z, bc, points, (void *)normals), *dm, 0); 698 if (pointtype == CGNS_ENUMV(ElementRange)) { 699 // Range of cells: assuming half-open interval 700 for (c = points[0]; c < points[1]; ++c) PetscCall(DMLabelSetValue(label, c - cellStart[z - 1], 1)); 701 } else if (pointtype == CGNS_ENUMV(ElementList)) { 702 // List of cells 703 for (c = 0; c < npoints; ++c) PetscCall(DMLabelSetValue(label, points[c] - cellStart[z - 1], 1)); 704 } else if (pointtype == CGNS_ENUMV(PointRange)) { 705 CGNS_ENUMT(GridLocation_t) gridloc; 706 707 // List of points: 708 PetscCallCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 709 PetscCallCGNSRead(cg_gridlocation_read(&gridloc), *dm, 0); 710 // Range of points: assuming half-open interval 711 for (c = points[0]; c < points[1]; ++c) { 712 if (gridloc == CGNS_ENUMV(Vertex)) PetscCall(DMLabelSetValue(label, c - vertStart[z - 1], 1)); 713 else PetscCall(DMLabelSetValue(label, c - cellStart[z - 1], 1)); 714 } 715 } else if (pointtype == CGNS_ENUMV(PointList)) { 716 CGNS_ENUMT(GridLocation_t) gridloc; 717 718 // List of points: 719 PetscCallCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 720 PetscCallCGNSRead(cg_gridlocation_read(&gridloc), *dm, 0); 721 for (c = 0; c < npoints; ++c) { 722 if (gridloc == CGNS_ENUMV(Vertex)) PetscCall(DMLabelSetValue(label, points[c] - vertStart[z - 1], 1)); 723 else PetscCall(DMLabelSetValue(label, points[c] - cellStart[z - 1], 1)); 724 } 725 } else SETERRQ(comm, PETSC_ERR_SUP, "Unsupported point set type %d", (int)pointtype); 726 PetscCall(PetscFree2(points, normals)); 727 } 728 } 729 PetscCall(PetscFree2(cellStart, vertStart)); 730 } 731 PetscCall(DMGetNumLabels(*dm, &labelIdRange[1])); 732 PetscCallMPI(MPI_Bcast(labelIdRange, 2, MPIU_INT, 0, comm)); 733 734 /* Create BC labels at all processes */ 735 for (labelId = labelIdRange[0]; labelId < labelIdRange[1]; ++labelId) { 736 char *labelName = buffer; 737 size_t len = sizeof(buffer); 738 const char *locName; 739 740 if (rank == 0) { 741 PetscCall(DMGetLabelByNum(*dm, labelId, &label)); 742 PetscCall(PetscObjectGetName((PetscObject)label, &locName)); 743 PetscCall(PetscStrncpy(labelName, locName, len)); 744 } 745 PetscCallMPI(MPI_Bcast(labelName, (PetscMPIInt)len, MPIU_INT, 0, comm)); 746 PetscCallMPI(DMCreateLabel(*dm, labelName)); 747 } 748 PetscFunctionReturn(PETSC_SUCCESS); 749 } 750 751 typedef struct { 752 cgsize_t start, end; 753 } CGRange; 754 755 // Creates a PetscLayout from the given sizes, but adjusts the ranges by the offset. So the first rank ranges will be [offset, offset + local_size) rather than [0, local_size) 756 static PetscErrorCode PetscLayoutCreateFromSizesAndOffset(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt bs, PetscInt offset, PetscLayout *map) 757 { 758 PetscLayout init; 759 const PetscInt *ranges; 760 PetscInt *new_ranges; 761 PetscMPIInt num_ranks; 762 763 PetscFunctionBegin; 764 PetscCall(PetscLayoutCreateFromSizes(comm, n, N, bs, &init)); 765 PetscCall(PetscLayoutGetRanges(init, &ranges)); 766 PetscCallMPI(MPI_Comm_size(comm, &num_ranks)); 767 PetscCall(PetscMalloc1(num_ranks + 1, &new_ranges)); 768 PetscCall(PetscArraycpy(new_ranges, ranges, num_ranks + 1)); 769 for (PetscInt r = 0; r < num_ranks + 1; r++) new_ranges[r] += offset; 770 PetscCall(PetscLayoutCreateFromRanges(comm, new_ranges, PETSC_OWN_POINTER, bs, map)); 771 PetscCall(PetscLayoutDestroy(&init)); 772 PetscFunctionReturn(PETSC_SUCCESS); 773 } 774 775 /** 776 @brief Creates connectivity array of CGNS elements' corners in Plex ordering, with a `PetscSection` to describe the data layout 777 778 @param[in] dm `DM` 779 @param[in] cgid CGNS file ID 780 @param[in] base CGNS base ID 781 @param[in] zone CGNS zone ID 782 @param[in] num_sections Number of sections to put in connectivity 783 @param[in] section_ids CGNS section IDs to obtain connectivity from 784 @param[out] section Section describing the connectivity for each element 785 @param[out] cellTypes Array specifying the CGNS ElementType_t for each element 786 @param[out] cell_ids CGNS IDs of the cells 787 @param[out] layouts Array of `PetscLayout` that describes the distributed ownership of the cells in each `section_ids` 788 @param[out] connectivity Array of the cell connectivity, described by `section`. The vertices are the local Plex IDs 789 790 @description 791 792 Each point in `section` corresponds to the `cellTypes` and `cell_ids` arrays. 793 The dof and offset of `section` maps into the `connectivity` array. 794 795 The `layouts` array is intended to be used with `PetscLayoutFindOwnerIndex_CGNSSectionLayouts()` 796 **/ 797 static PetscErrorCode DMPlexCGNS_CreateCornersConnectivitySection(DM dm, PetscInt cgid, int base, int zone, PetscInt num_sections, const int section_ids[], PetscSection *section, CGNS_ENUMT(ElementType_t) * cellTypes[], PetscInt *cell_ids[], PetscLayout *layouts[], PetscInt *connectivity[]) 798 { 799 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 800 PetscSection section_; 801 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 802 CGNS_ENUMT(ElementType_t) * sectionCellTypes, *cellTypes_; 803 CGRange *ranges; 804 PetscInt nlocal_cells = 0, global_cell_dim = -1; 805 PetscSegBuffer conn_sb; 806 PetscLayout *layouts_; 807 808 PetscFunctionBegin; 809 PetscCall(PetscMalloc2(num_sections, &ranges, num_sections, §ionCellTypes)); 810 PetscCall(PetscMalloc1(num_sections, &layouts_)); 811 for (PetscInt s = 0; s < num_sections; s++) { 812 int nbndry, parentFlag; 813 PetscInt local_size; 814 815 PetscCallCGNSRead(cg_section_read(cgid, base, zone, section_ids[s], buffer, §ionCellTypes[s], &ranges[s].start, &ranges[s].end, &nbndry, &parentFlag), dm, 0); 816 PetscCheck(sectionCellTypes[s] != CGNS_ENUMV(NGON_n) && sectionCellTypes[s] != CGNS_ENUMV(NFACE_n), comm, PETSC_ERR_SUP, "CGNS reader does not support elements of type NGON_n or NFACE_n"); 817 PetscInt num_section_cells = ranges[s].end - ranges[s].start + 1; 818 PetscCall(PetscLayoutCreateFromSizesAndOffset(comm, PETSC_DECIDE, num_section_cells, 1, ranges[s].start, &layouts_[s])); 819 PetscCall(PetscLayoutGetLocalSize(layouts_[s], &local_size)); 820 nlocal_cells += local_size; 821 } 822 PetscCall(PetscSectionCreate(comm, §ion_)); 823 PetscCall(PetscSectionSetChart(section_, 0, nlocal_cells)); 824 825 PetscCall(PetscMalloc1(nlocal_cells, cell_ids)); 826 PetscCall(PetscMalloc1(nlocal_cells, &cellTypes_)); 827 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), nlocal_cells * 2, &conn_sb)); 828 for (PetscInt s = 0, c = 0; s < num_sections; s++) { 829 PetscInt mystart, myend, myowned; 830 831 PetscCall(PetscLayoutGetRange(layouts_[s], &mystart, &myend)); 832 PetscCall(PetscLayoutGetLocalSize(layouts_[s], &myowned)); 833 if (sectionCellTypes[s] == CGNS_ENUMV(MIXED)) { 834 cgsize_t *offsets, *conn_cg; 835 836 PetscCall(PetscMalloc1(myowned + 1, &offsets)); // The last element in the array is the total size of the connectivity for the given [start,end] range 837 PetscCallCGNSRead(cgp_poly_elements_read_data_offsets(cgid, base, zone, section_ids[s], mystart, myend - 1, offsets), dm, 0); 838 PetscCall(PetscMalloc1(offsets[myowned + 1], &conn_cg)); 839 PetscCallCGNSRead(cgp_poly_elements_read_data_elements(cgid, base, zone, section_ids[s], mystart, myend - 1, offsets, conn_cg), dm, 0); 840 for (PetscInt i = 0; i < myowned; i++) { 841 DMPolytopeType dm_cell_type = DM_POLYTOPE_UNKNOWN; 842 PetscInt numCorners, cell_dim, *conn_sb_seg; 843 const int *perm; 844 845 (*cell_ids)[c] = mystart + i; 846 847 cellTypes_[c] = (CGNS_ENUMT(ElementType_t))conn_cg[offsets[i]]; 848 PetscCall(CGNSElementTypeGetTopologyInfo(cellTypes_[c], &dm_cell_type, &numCorners, &cell_dim)); 849 if (global_cell_dim == -1) global_cell_dim = cell_dim; 850 else 851 PetscCheck(cell_dim == global_cell_dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Can only combine cells of the same dimension. Global cell dimension detected as %" PetscInt_FMT ", but CGNS element %" PetscInt_FMT " is dimension %" PetscInt_FMT, global_cell_dim, (*cell_ids)[c], cell_dim); 852 PetscCall(PetscSegBufferGetInts(conn_sb, numCorners, &conn_sb_seg)); 853 854 PetscCall(DMPlexCGNSGetPermutation_Internal(dm_cell_type, numCorners, NULL, &perm)); 855 for (PetscInt v = 0; v < numCorners; ++v) conn_sb_seg[perm[v]] = conn_cg[offsets[i] + 1 + v]; 856 PetscCall(PetscSectionSetDof(section_, c, numCorners)); 857 c++; 858 } 859 PetscCall(PetscFree(offsets)); 860 PetscCall(PetscFree(conn_cg)); 861 } else { 862 PetscInt numCorners, cell_dim; 863 PetscInt *conn_sb_seg; 864 DMPolytopeType dm_cell_type = DM_POLYTOPE_UNKNOWN; 865 int npe; // nodes per element 866 const int *perm; 867 cgsize_t *conn_cg; 868 869 PetscCall(CGNSElementTypeGetTopologyInfo(sectionCellTypes[s], &dm_cell_type, &numCorners, &cell_dim)); 870 if (global_cell_dim == -1) global_cell_dim = cell_dim; 871 else 872 PetscCheck(cell_dim == global_cell_dim, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Can only combine cells of the same dimension. Global cell dimension detected as %" PetscInt_FMT ", but CGNS element %" PetscInt_FMT " is dimension %" PetscInt_FMT, global_cell_dim, (*cell_ids)[c], cell_dim); 873 874 PetscCallCGNSRead(cg_npe(sectionCellTypes[s], &npe), dm, 0); 875 PetscCall(PetscMalloc1(myowned * npe, &conn_cg)); 876 PetscCallCGNSRead(cgp_elements_read_data(cgid, base, zone, section_ids[s], mystart, myend - 1, conn_cg), dm, 0); 877 PetscCall(DMPlexCGNSGetPermutation_Internal(dm_cell_type, numCorners, NULL, &perm)); 878 PetscCall(PetscSegBufferGetInts(conn_sb, numCorners * myowned, &conn_sb_seg)); 879 for (PetscInt i = 0; i < myowned; i++) { 880 (*cell_ids)[c] = mystart + i; 881 cellTypes_[c] = sectionCellTypes[s]; 882 for (PetscInt v = 0; v < numCorners; ++v) conn_sb_seg[i * numCorners + perm[v]] = conn_cg[i * npe + v]; 883 PetscCall(PetscSectionSetDof(section_, c, numCorners)); 884 c++; 885 } 886 PetscCall(PetscFree(conn_cg)); 887 } 888 } 889 890 PetscCall(PetscSectionSetUp(section_)); 891 *section = section_; 892 PetscCall(PetscSegBufferExtractAlloc(conn_sb, connectivity)); 893 PetscInt connSize; 894 PetscCall(PetscSectionGetStorageSize(section_, &connSize)); 895 for (PetscInt i = 0; i < connSize; i++) (*connectivity)[i] -= 1; // vertices should be 0-based indexing for consistency with DMPlexBuildFromCellListParallel() 896 *layouts = layouts_; 897 if (cellTypes) *cellTypes = cellTypes_; 898 else PetscCall(PetscFree(cellTypes_)); 899 900 PetscCall(PetscSegBufferDestroy(&conn_sb)); 901 PetscCall(PetscFree2(ranges, sectionCellTypes)); 902 PetscFunctionReturn(PETSC_SUCCESS); 903 } 904 905 static inline PetscErrorCode PetscFindIntUnsorted(PetscInt key, PetscInt size, const PetscInt array[], PetscInt *loc) 906 { 907 PetscFunctionBegin; 908 *loc = -1; 909 for (PetscInt i = 0; i < size; i++) { 910 if (array[i] == key) { 911 *loc = i; 912 break; 913 } 914 } 915 PetscFunctionReturn(PETSC_SUCCESS); 916 } 917 918 // Same as `PetscLayoutFindOwnerIndex()`, but does not fail if owner not in PetscLayout 919 static PetscErrorCode PetscLayoutFindOwnerIndex_Internal(PetscLayout map, PetscInt idx, PetscMPIInt *owner, PetscInt *lidx, PetscBool *found_owner) 920 { 921 PetscMPIInt lo = 0, hi, t; 922 923 PetscFunctionBegin; 924 PetscAssert((map->n >= 0) && (map->N >= 0) && (map->range), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PetscLayoutSetUp() must be called first"); 925 if (owner) *owner = -1; 926 if (lidx) *lidx = -1; 927 if (idx < map->range[0] && idx >= map->range[map->size + 1]) { 928 *found_owner = PETSC_FALSE; 929 PetscFunctionReturn(PETSC_SUCCESS); 930 } 931 hi = map->size; 932 while (hi - lo > 1) { 933 t = lo + (hi - lo) / 2; 934 if (idx < map->range[t]) hi = t; 935 else lo = t; 936 } 937 if (owner) *owner = lo; 938 if (lidx) *lidx = idx - map->range[lo]; 939 if (found_owner) *found_owner = PETSC_TRUE; 940 PetscFunctionReturn(PETSC_SUCCESS); 941 } 942 943 // This function assumes that there is an array that `maps` describes the layout too. Locally, the range of each map is concatenated onto each other. 944 // So [maps[0].start, ..., maps[0].end - 1, maps[1].start, ..., maps[1].end - 1, ...] 945 // The returned index is the index into this array 946 static PetscErrorCode PetscLayoutFindOwnerIndex_CGNSSectionLayouts(PetscLayout maps[], PetscInt nmaps, PetscInt idx, PetscMPIInt *owner, PetscInt *lidx, PetscInt *mapidx) 947 { 948 PetscFunctionBegin; 949 for (PetscInt m = 0; m < nmaps; m++) { 950 PetscBool found_owner = PETSC_FALSE; 951 PetscCall(PetscLayoutFindOwnerIndex_Internal(maps[m], idx, owner, lidx, &found_owner)); 952 if (found_owner) { 953 // Now loop back through the previous maps to get the local offset for the containing index 954 for (PetscInt mm = m - 1; mm >= 0; mm--) { 955 PetscInt size = maps[mm]->range[*owner + 1] - maps[mm]->range[*owner]; 956 *lidx += size; 957 } 958 if (mapidx) *mapidx = m; 959 PetscFunctionReturn(PETSC_SUCCESS); 960 } 961 } 962 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "CGNS id %" PetscInt_FMT " not found in layouts", idx); 963 } 964 965 /* 966 @brief Matching root and leaf indices across processes, allowing multiple roots per leaf 967 968 Collective 969 970 Input Parameters: 971 + layout - `PetscLayout` defining the global index space and the MPI rank that brokers each index 972 . numRootIndices - size of `rootIndices` 973 . rootIndices - array of global indices of which this process requests ownership 974 . rootLocalIndices - root local index permutation (`NULL` if no permutation) 975 . rootLocalOffset - offset to be added to `rootLocalIndices` 976 . numLeafIndices - size of `leafIndices` 977 . leafIndices - array of global indices with which this process requires data associated 978 . leafLocalIndices - leaf local index permutation (`NULL` if no permutation) 979 - leafLocalOffset - offset to be added to `leafLocalIndices` 980 981 Output Parameters: 982 + matchSection - The `PetscSection` describing the layout of `matches` with respect to the `leafIndices` (the points of the section indexes into `leafIndices`) 983 - matches - Array of `PetscSFNode` denoting the location (rank and index) of the root matching the leaf. 984 985 Example 1: 986 .vb 987 rank : 0 1 2 988 rootIndices : [1 0 2] [3] [3] 989 rootLocalOffset : 100 200 300 990 layout : [0 1] [2] [3] 991 leafIndices : [0] [2] [0 3] 992 leafLocalOffset : 400 500 600 993 994 would build the following result 995 996 rank 0: [0] 400 <- [0] (0,101) 997 ------------------------------ 998 rank 1: [0] 500 <- [0] (0,102) 999 ------------------------------ 1000 rank 2: [0] 600 <- [0] (0,101) 1001 [1] 601 <- [1] (1,200) 1002 [1] 601 <- [2] (2,300) 1003 | | | | 1004 | | | + `matches`, the rank and index of the respective match with the leaf index 1005 | | + index into `matches` array 1006 | + the leaves for the respective root 1007 + The point in `matchSection` (indexes into `leafIndices`) 1008 1009 For rank 2, the `matchSection` would be: 1010 1011 [0]: (1, 0) 1012 [1]: (2, 1) 1013 | | | 1014 | | + offset 1015 | + ndof 1016 + point 1017 .ve 1018 1019 Notes: 1020 This function is identical to `PetscSFCreateByMatchingIndices()` except it includes *all* matching indices instead of designating a single rank as the "owner". 1021 Attempting to create an SF with all matching indices would create an invalid SF, thus we give an array of `matches` and `matchSection` to describe the layout 1022 Compare the examples in this document to those in `PetscSFCreateByMatchingIndices()`. 1023 1024 .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCreate()`, `PetscSFCreateByMatchingIndices()` 1025 */ 1026 static PetscErrorCode PetscSFFindMatchingIndices(PetscLayout layout, PetscInt numRootIndices, const PetscInt rootIndices[], const PetscInt rootLocalIndices[], PetscInt rootLocalOffset, PetscInt numLeafIndices, const PetscInt leafIndices[], const PetscInt leafLocalIndices[], PetscInt leafLocalOffset, PetscSection *matchSection, PetscSFNode *matches[]) 1027 { 1028 MPI_Comm comm = layout->comm; 1029 PetscMPIInt rank; 1030 PetscSF sf1; 1031 PetscSection sectionBuffer, matchSection_; 1032 PetscInt numMatches; 1033 PetscSFNode *roots, *buffer, *matches_; 1034 PetscInt N, n, pStart, pEnd; 1035 PetscBool areIndicesSame; 1036 1037 PetscFunctionBegin; 1038 if (rootIndices) PetscAssertPointer(rootIndices, 3); 1039 if (rootLocalIndices) PetscAssertPointer(rootLocalIndices, 4); 1040 if (leafIndices) PetscAssertPointer(leafIndices, 7); 1041 if (leafLocalIndices) PetscAssertPointer(leafLocalIndices, 8); 1042 PetscCheck(numRootIndices >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numRootIndices (%" PetscInt_FMT ") must be non-negative", numRootIndices); 1043 PetscCheck(numLeafIndices >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numLeafIndices (%" PetscInt_FMT ") must be non-negative", numLeafIndices); 1044 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1045 PetscCall(PetscLayoutSetUp(layout)); 1046 PetscCall(PetscLayoutGetSize(layout, &N)); 1047 PetscCall(PetscLayoutGetLocalSize(layout, &n)); 1048 areIndicesSame = (PetscBool)(leafIndices == rootIndices); 1049 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &areIndicesSame, 1, MPI_C_BOOL, MPI_LAND, comm)); 1050 PetscCheck(!areIndicesSame || numLeafIndices == numRootIndices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "leafIndices == rootIndices, but numLeafIndices (%" PetscInt_FMT ") != numRootIndices(%" PetscInt_FMT ")", numLeafIndices, numRootIndices); 1051 if (PetscDefined(USE_DEBUG)) { 1052 PetscInt N1 = PETSC_INT_MIN; 1053 for (PetscInt i = 0; i < numRootIndices; i++) 1054 if (rootIndices[i] > N1) N1 = rootIndices[i]; 1055 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm)); 1056 PetscCheck(N1 < N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. root index (%" PetscInt_FMT ") out of layout range [0,%" PetscInt_FMT ")", N1, N); 1057 if (!areIndicesSame) { 1058 N1 = PETSC_INT_MIN; 1059 for (PetscInt i = 0; i < numLeafIndices; i++) 1060 if (leafIndices[i] > N1) N1 = leafIndices[i]; 1061 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm)); 1062 PetscCheck(N1 < N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. leaf index (%" PetscInt_FMT ") out of layout range [0,%" PetscInt_FMT ")", N1, N); 1063 } 1064 } 1065 1066 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1067 { /* Reduce: roots -> buffer */ 1068 // Data in buffer described by section_buffer. The chart of `section_buffer` maps onto the local portion of `layout`, with dofs denoting how many matches there are. 1069 PetscInt bufsize; 1070 const PetscInt *root_degree; 1071 1072 PetscCall(PetscSFCreate(comm, &sf1)); 1073 PetscCall(PetscSFSetFromOptions(sf1)); 1074 PetscCall(PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices)); 1075 1076 PetscCall(PetscSFComputeDegreeBegin(sf1, &root_degree)); 1077 PetscCall(PetscSFComputeDegreeEnd(sf1, &root_degree)); 1078 PetscCall(PetscSectionCreate(comm, §ionBuffer)); 1079 PetscCall(PetscSectionSetChart(sectionBuffer, 0, n)); 1080 PetscCall(PetscMalloc1(numRootIndices, &roots)); 1081 for (PetscInt i = 0; i < numRootIndices; ++i) { 1082 roots[i].rank = rank; 1083 roots[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i); 1084 } 1085 for (PetscInt i = 0; i < n; i++) PetscCall(PetscSectionSetDof(sectionBuffer, i, root_degree[i])); 1086 PetscCall(PetscSectionSetUp(sectionBuffer)); 1087 PetscCall(PetscSectionGetStorageSize(sectionBuffer, &bufsize)); 1088 PetscCall(PetscMalloc1(bufsize, &buffer)); 1089 for (PetscInt i = 0; i < bufsize; ++i) { 1090 buffer[i].index = -1; 1091 buffer[i].rank = -1; 1092 } 1093 PetscCall(PetscSFGatherBegin(sf1, MPIU_SF_NODE, roots, buffer)); 1094 PetscCall(PetscSFGatherEnd(sf1, MPIU_SF_NODE, roots, buffer)); 1095 PetscCall(PetscFree(roots)); 1096 } 1097 1098 // Distribute data in buffers to the leaf locations. The chart of `sectionMatches` maps to `leafIndices`, with dofs denoting how many matches there are for each leaf. 1099 if (!areIndicesSame) PetscCall(PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices)); 1100 PetscCall(PetscSectionCreate(comm, &matchSection_)); 1101 PetscCall(PetscSectionMigrateData(sf1, MPIU_SF_NODE, sectionBuffer, buffer, matchSection_, (void **)&matches_, NULL)); 1102 PetscCall(PetscSectionGetChart(matchSection_, &pStart, &pEnd)); 1103 PetscCheck(pEnd - pStart == numLeafIndices, comm, PETSC_ERR_PLIB, "Section of matches has different chart size (%" PetscInt_FMT ") than number of leaf indices %" PetscInt_FMT ". Section chart is [%" PetscInt_FMT ", %" PetscInt_FMT ")", pEnd - pStart, numLeafIndices, pStart, pEnd); 1104 PetscCall(PetscSectionGetStorageSize(matchSection_, &numMatches)); 1105 for (PetscInt p = pStart; p < pEnd; p++) { 1106 PetscInt ndofs; 1107 PetscCall(PetscSectionGetDof(matchSection_, p, &ndofs)); 1108 PetscCheck(ndofs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No match found for index %" PetscInt_FMT, leafIndices[p]); 1109 } 1110 1111 *matchSection = matchSection_; 1112 *matches = matches_; 1113 1114 PetscCall(PetscFree(buffer)); 1115 PetscCall(PetscSectionDestroy(§ionBuffer)); 1116 PetscCall(PetscSFDestroy(&sf1)); 1117 PetscFunctionReturn(PETSC_SUCCESS); 1118 } 1119 1120 /** 1121 @brief Match CGNS faces to their Plex equivalents 1122 1123 @param[in] dm DM that holds the Plex to match against 1124 @param[in] nuniq_verts Number of unique CGNS vertices on this rank 1125 @param[in] uniq_verts Unique CGNS vertices on this rank 1126 @param[in] plex_vertex_offset Index offset to match `uniq_vertices` to their respective Plex vertices 1127 @param[in] NVertices Number of vertices for Layout for `PetscSFFindMatchingIndices()` 1128 @param[in] connSection PetscSection describing the CGNS face connectivity 1129 @param[in] face_ids Array of the CGNS face IDs 1130 @param[in] conn Array of the CGNS face connectivity 1131 @param[out] cg2plexSF PetscSF describing the mapping from owned CGNS faces to remote `plexFaces` 1132 @param[out] plexFaces Matching Plex face IDs 1133 1134 @description 1135 1136 `cg2plexSF` is a mapping from the owned CGNS faces to the rank whose local Plex has that face. 1137 `plexFaces` holds the actual mesh point in the local Plex that corresponds to the owned CGNS face (which is the root) 1138 1139 cg2plexSF 1140 __________|__________ 1141 | | 1142 1143 [F0_11] -----> [P0_0] [38] 1144 [F0_12] -- Rank 0 1145 ~~~~~~~~~ \ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1146 --> [P1_0] [59] Rank 1 1147 [F1_13] -----> [P1_1] [98] 1148 | | | 1149 | | + plexFaces, maps the leaves of cg2plexSF to local Plex face mesh points 1150 | + Leaves of cg2plexSF. P(rank)_(root_index) 1151 + Roots of cg2plexSF, F(rank)_(CGNS face ID) 1152 1153 Note that, unlike a pointSF, the leaves of `cg2plexSF` do not map onto chart of the local Plex, but just onto an array. 1154 The plexFaces array is then what maps the leaves to the actual local Plex mesh points. 1155 1156 `plex_vertex_offset` is used to map the CGNS vertices in `uniq_vertices` to their respective Plex vertices. 1157 From `DMPlexBuildFromCellListParallel()`, the mapping of CGNS vertices to Plex vertices is uniq_vert[i] -> i + plex_vertex_offset, where the right side is the Plex point ID. 1158 So with `plex_vertex_offset = 5`, 1159 uniq_vertices: [19, 52, 1, 89] 1160 plex_point_ids: [5, 6, 7, 8] 1161 **/ 1162 static PetscErrorCode DMPlexCGNS_MatchCGNSFacesToPlexFaces(DM dm, PetscInt nuniq_verts, const PetscInt uniq_verts[], PetscInt plex_vertex_offset, PetscInt NVertices, PetscSection connSection, const PetscInt face_ids[], const PetscInt conn[], PetscSF *cg2plexSF, PetscInt *plexFaces[]) 1163 { 1164 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 1165 PetscMPIInt myrank, nranks; 1166 PetscInt fownedStart, fownedEnd, fownedSize; 1167 1168 PetscFunctionBegin; 1169 PetscCallMPI(MPI_Comm_rank(comm, &myrank)); 1170 PetscCallMPI(MPI_Comm_size(comm, &nranks)); 1171 1172 { // -- Create cg2plexSF 1173 PetscInt nuniq_face_verts, *uniq_face_verts; 1174 PetscSection fvert2mvertSection; 1175 PetscSFNode *fvert2mvert = NULL; 1176 1177 { // -- Create fvert2mvert, which map CGNS vertices in the owned-face connectivity to the CGNS vertices in the global mesh 1178 PetscLayout layout; 1179 1180 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NVertices, 1, &layout)); 1181 { // Count locally unique vertices in the face connectivity 1182 PetscHSetI vhash; 1183 PetscInt off = 0, conn_size; 1184 1185 PetscCall(PetscHSetICreate(&vhash)); 1186 PetscCall(PetscSectionGetStorageSize(connSection, &conn_size)); 1187 for (PetscInt v = 0; v < conn_size; ++v) PetscCall(PetscHSetIAdd(vhash, conn[v])); 1188 PetscCall(PetscHSetIGetSize(vhash, &nuniq_face_verts)); 1189 PetscCall(PetscMalloc1(nuniq_face_verts, &uniq_face_verts)); 1190 PetscCall(PetscHSetIGetElems(vhash, &off, uniq_face_verts)); 1191 PetscCall(PetscHSetIDestroy(&vhash)); 1192 PetscCheck(off == nuniq_face_verts, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %" PetscInt_FMT " should be %" PetscInt_FMT, off, nuniq_face_verts); 1193 } 1194 PetscCall(PetscSortInt(nuniq_face_verts, uniq_face_verts)); 1195 PetscCall(PetscSFFindMatchingIndices(layout, nuniq_verts, uniq_verts, NULL, 0, nuniq_face_verts, uniq_face_verts, NULL, 0, &fvert2mvertSection, &fvert2mvert)); 1196 1197 PetscCall(PetscLayoutDestroy(&layout)); 1198 } 1199 1200 PetscSFNode *plexFaceRemotes, *ownedFaceRemotes; 1201 PetscCount nPlexFaceRemotes; 1202 PetscInt *local_rank_count; 1203 1204 PetscCall(PetscSectionGetChart(connSection, &fownedStart, &fownedEnd)); 1205 fownedSize = fownedEnd - fownedStart; 1206 PetscCall(PetscCalloc1(nranks, &local_rank_count)); 1207 1208 { // Find the rank(s) whose local Plex has the owned CGNS face 1209 // We determine ownership by determining which ranks contain all the vertices in a face's connectivity 1210 PetscInt maxRanksPerVert; 1211 PetscInt *face_ranks; 1212 PetscSegBuffer plexFaceRemotes_SB, ownedFaceRemotes_SB; 1213 1214 PetscCall(PetscSegBufferCreate(sizeof(PetscSFNode), fownedSize, &plexFaceRemotes_SB)); 1215 PetscCall(PetscSegBufferCreate(sizeof(PetscSFNode), fownedSize, &ownedFaceRemotes_SB)); 1216 PetscCall(PetscSectionGetMaxDof(fvert2mvertSection, &maxRanksPerVert)); 1217 PetscCall(PetscMalloc1(maxRanksPerVert, &face_ranks)); 1218 for (PetscInt f = fownedStart, f_i = 0; f < fownedEnd; f++, f_i++) { 1219 PetscInt fndof, foffset, lndof, loffset, idx, nface_ranks = 0; 1220 1221 PetscCall(PetscSectionGetDof(connSection, f, &fndof)); 1222 PetscCall(PetscSectionGetOffset(connSection, f, &foffset)); 1223 PetscCall(PetscFindInt(conn[foffset + 0], nuniq_face_verts, uniq_face_verts, &idx)); 1224 PetscCall(PetscSectionGetDof(fvert2mvertSection, idx, &lndof)); 1225 PetscCall(PetscSectionGetOffset(fvert2mvertSection, idx, &loffset)); 1226 // Loop over ranks of the first vertex in the face connectivity 1227 for (PetscInt l = 0; l < lndof; l++) { 1228 PetscInt rank = fvert2mvert[loffset + l].rank; 1229 // Loop over vertices of face (except for the first) to see if those vertices have the same candidate rank 1230 for (PetscInt v = 1; v < fndof; v++) { 1231 PetscInt ldndof, ldoffset, idx; 1232 PetscBool vert_has_rank = PETSC_FALSE; 1233 1234 PetscCall(PetscFindInt(conn[foffset + v], nuniq_face_verts, uniq_face_verts, &idx)); 1235 PetscCall(PetscSectionGetDof(fvert2mvertSection, idx, &ldndof)); 1236 PetscCall(PetscSectionGetOffset(fvert2mvertSection, idx, &ldoffset)); 1237 // Loop over ranks of the vth vertex to see if it has the candidate rank 1238 for (PetscInt ld = 0; ld < ldndof; ld++) vert_has_rank = (fvert2mvert[ldoffset + ld].rank == rank) || vert_has_rank; 1239 if (vert_has_rank) continue; // This vertex has the candidate rank, proceed to the next vertex 1240 else goto next_candidate_rank; 1241 } 1242 face_ranks[nface_ranks++] = rank; 1243 next_candidate_rank: 1244 continue; 1245 } 1246 PetscCheck(nface_ranks > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Rank containing CGNS face %" PetscInt_FMT " could not be found", face_ids[f_i]); 1247 1248 // Once we have found the rank(s) whose Plex has the vertices of CGNS face `f`, we can begin to build the information for the SF. 1249 // We want the roots to be the CGNS faces and the leaves to be the corresponding Plex faces, but we have the opposite; for the owned face on myrank, we know the rank that has the corresponding Plex face. 1250 // To get the inverse, we assign each SF edge with a tuple of PetscSFNodes; one in `plexFaceRemotes` and the other in `ownedFaceRemotes`. 1251 // `plexFaceRemotes` has the rank and index of the Plex face. 1252 // `ownedFaceRemotes` has the rank and index of the owned CGNS face. 1253 // Note that `ownedFaceRemotes` is all on my rank (e.g. rank == myrank). 1254 // 1255 // Then, to build the `cg2plexSF`, we communicate the `ownedFaceRemotes` to the `plexFaceRemotes` (via SFReduce). 1256 // Those `ownedFaceRemotes` then act as the leaves to the roots on this process. 1257 // 1258 // Conceptually, this is the same as calling `PetscSFCreateInverseSF` on an SF with `iremotes = plexFaceRemotes` and the `ilocal = ownedFaceRemotes[:].index`. 1259 // However, we cannot use this much simpler way because when there are multiple matching Plex faces, `PetscSFCreateInverseSF()` will be invalid due to `ownedFaceRemotes[:].index` having repeated values (only root vertices of the SF graph may have degree > 1) 1260 PetscSFNode *plexFaceRemotes_buffer, *ownedFaceRemotes_buffer; 1261 PetscCall(PetscSegBufferGet(plexFaceRemotes_SB, nface_ranks, &plexFaceRemotes_buffer)); 1262 PetscCall(PetscSegBufferGet(ownedFaceRemotes_SB, nface_ranks, &ownedFaceRemotes_buffer)); 1263 for (PetscInt n = 0; n < nface_ranks; n++) { 1264 plexFaceRemotes_buffer[n].rank = face_ranks[n]; 1265 local_rank_count[face_ranks[n]]++; 1266 ownedFaceRemotes_buffer[n] = (PetscSFNode){.rank = myrank, .index = f_i}; 1267 } 1268 } 1269 PetscCall(PetscFree(face_ranks)); 1270 1271 PetscCall(PetscSegBufferGetSize(plexFaceRemotes_SB, &nPlexFaceRemotes)); 1272 PetscCall(PetscSegBufferExtractAlloc(plexFaceRemotes_SB, &plexFaceRemotes)); 1273 PetscCall(PetscSegBufferDestroy(&plexFaceRemotes_SB)); 1274 PetscCall(PetscSegBufferExtractAlloc(ownedFaceRemotes_SB, &ownedFaceRemotes)); 1275 PetscCall(PetscSegBufferDestroy(&ownedFaceRemotes_SB)); 1276 1277 // To get the index for plexFaceRemotes, we partition the leaves on each rank (e.g. the array that will hold the local Plex face mesh points) by each rank that has the CGNS owned rank. 1278 // For r in [0,numranks), local_rank_count[r] holds the number plexFaces that myrank holds. 1279 // This determines how large a partition the leaves on rank r need to create for myrank. 1280 // To get the offset into the leaves, we use Exscan to get rank_start. 1281 // For r in [0, numranks), rank_start[r] holds the offset into rank r's leaves that myrank will index into. 1282 1283 // Below is an example: 1284 // 1285 // myrank: | 0 1 2 1286 // local_rank_count: | [3, 2, 0] [1, 0, 2] [2, 2, 1] 1287 // myrank_total_count: | 6 4 3 1288 // rank_start: | [0, 0, 0] [3, 2, 0] [4, 2, 2] 1289 // | 1290 // plexFaceRemotes: 0 | (0, 0) (2, 0) (2, 2) <-- (rank, index) tuples (e.g. PetscSFNode) 1291 // 1 | (1, 0) (2, 1) (0, 4) 1292 // 2 | (0, 1) (0, 3) (1, 2) 1293 // 3 | (0, 2) (0, 5) 1294 // 4 | (1, 1) (1, 3) 1295 // 1296 // leaves: 0 | (0, 0) (0, 1) (1, 0) (rank and index into plexFaceRemotes) 1297 // 1 | (0, 2) (0, 4) (1, 1) 1298 // 2 | (0, 3) (2, 2) (2, 0) 1299 // 3 | (1, 2) (2, 4) 1300 // 4 | (2, 1) 1301 // 5 | (2, 3) 1302 // 1303 // Note how at the leaves, the ranks are contiguous and in order 1304 PetscInt myrank_total_count; 1305 { 1306 PetscInt *rank_start, *rank_offset; 1307 1308 PetscCall(PetscCalloc2(nranks, &rank_start, nranks, &rank_offset)); 1309 PetscCallMPI(MPIU_Allreduce(local_rank_count, rank_start, nranks, MPIU_INT, MPI_SUM, comm)); 1310 myrank_total_count = rank_start[myrank]; 1311 PetscCall(PetscArrayzero(rank_start, nranks)); 1312 PetscCallMPI(MPI_Exscan(local_rank_count, rank_start, nranks, MPIU_INT, MPI_SUM, comm)); 1313 1314 for (PetscInt r = 0; r < nPlexFaceRemotes; r++) { 1315 PetscInt rank = plexFaceRemotes[r].rank; 1316 plexFaceRemotes[r].index = rank_start[rank] + rank_offset[rank]; 1317 rank_offset[rank]++; 1318 } 1319 PetscCall(PetscFree2(rank_start, rank_offset)); 1320 } 1321 1322 { // Communicate the leaves to roots and build cg2plexSF 1323 PetscSF plexRemotes2ownedRemotesSF; 1324 PetscSFNode *iremote_cg2plexSF; 1325 1326 PetscCall(PetscSFCreate(comm, &plexRemotes2ownedRemotesSF)); 1327 PetscCall(PetscSFSetGraph(plexRemotes2ownedRemotesSF, myrank_total_count, nPlexFaceRemotes, NULL, PETSC_COPY_VALUES, plexFaceRemotes, PETSC_OWN_POINTER)); 1328 PetscCall(PetscMalloc1(myrank_total_count, &iremote_cg2plexSF)); 1329 PetscCall(PetscSFViewFromOptions(plexRemotes2ownedRemotesSF, NULL, "-plex2ownedremotes_sf_view")); 1330 for (PetscInt i = 0; i < myrank_total_count; i++) iremote_cg2plexSF[i] = (PetscSFNode){.rank = -1, .index = -1}; 1331 PetscCall(PetscSFReduceBegin(plexRemotes2ownedRemotesSF, MPIU_SF_NODE, ownedFaceRemotes, iremote_cg2plexSF, MPI_REPLACE)); 1332 PetscCall(PetscSFReduceEnd(plexRemotes2ownedRemotesSF, MPIU_SF_NODE, ownedFaceRemotes, iremote_cg2plexSF, MPI_REPLACE)); 1333 PetscCall(PetscSFDestroy(&plexRemotes2ownedRemotesSF)); 1334 for (PetscInt i = 0; i < myrank_total_count; i++) PetscCheck(iremote_cg2plexSF[i].rank >= 0 && iremote_cg2plexSF[i].index != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Owned face SFNode was not reduced properly"); 1335 1336 PetscCall(PetscSFCreate(comm, cg2plexSF)); 1337 PetscCall(PetscSFSetGraph(*cg2plexSF, fownedSize, myrank_total_count, NULL, PETSC_COPY_VALUES, iremote_cg2plexSF, PETSC_OWN_POINTER)); 1338 1339 PetscCall(PetscFree(ownedFaceRemotes)); 1340 } 1341 1342 PetscCall(PetscFree(local_rank_count)); 1343 } 1344 1345 PetscCall(PetscSectionDestroy(&fvert2mvertSection)); 1346 PetscCall(PetscFree(uniq_face_verts)); 1347 PetscCall(PetscFree(fvert2mvert)); 1348 } 1349 1350 { // -- Find plexFaces 1351 // Distribute owned-CGNS-face connectivity to the ranks which have corresponding Plex faces, and then find the corresponding Plex faces 1352 PetscSection connDistSection; 1353 PetscInt *connDist; 1354 1355 // Distribute the face connectivity to the rank that has that face 1356 PetscCall(PetscSectionCreate(comm, &connDistSection)); 1357 PetscCall(PetscSectionMigrateData(*cg2plexSF, MPIU_INT, connSection, conn, connDistSection, (void **)&connDist, NULL)); 1358 1359 { // Translate CGNS vertex numbering to local Plex numbering 1360 PetscInt *dmplex_verts, *uniq_verts_sorted; 1361 PetscInt connDistSize; 1362 1363 PetscCall(PetscMalloc2(nuniq_verts, &dmplex_verts, nuniq_verts, &uniq_verts_sorted)); 1364 PetscCall(PetscArraycpy(uniq_verts_sorted, uniq_verts, nuniq_verts)); 1365 // uniq_verts are one-to-one with the DMPlex vertices with an offset, see DMPlexBuildFromCellListParallel() 1366 for (PetscInt v = 0; v < nuniq_verts; v++) dmplex_verts[v] = v + plex_vertex_offset; 1367 PetscCall(PetscSortIntWithArray(nuniq_verts, uniq_verts_sorted, dmplex_verts)); 1368 1369 PetscCall(PetscSectionGetStorageSize(connDistSection, &connDistSize)); 1370 for (PetscInt v = 0; v < connDistSize; v++) { 1371 PetscInt idx; 1372 PetscCall(PetscFindInt(connDist[v], nuniq_verts, uniq_verts_sorted, &idx)); 1373 PetscCheck(idx >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find CGNS vertex id (from face connectivity) in local plex"); 1374 connDist[v] = dmplex_verts[idx]; 1375 } 1376 PetscCall(PetscFree2(dmplex_verts, uniq_verts_sorted)); 1377 } 1378 1379 // Debugging info 1380 PetscBool view_connectivity = PETSC_FALSE; 1381 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_cgns_view_face_connectivity", &view_connectivity, NULL)); 1382 if (view_connectivity) { 1383 PetscSection conesSection; 1384 PetscInt *cones; 1385 1386 PetscCall(PetscPrintf(comm, "Distributed CGNS Face Connectivity (in Plex vertex numbering):\n")); 1387 PetscCall(PetscSectionArrayView(connDistSection, connDist, PETSC_INT, NULL)); 1388 PetscCall(DMPlexGetCones(dm, &cones)); 1389 PetscCall(DMPlexGetConeSection(dm, &conesSection)); 1390 PetscCall(PetscPrintf(comm, "Plex Cones:\n")); 1391 PetscCall(PetscSectionArrayView(conesSection, cones, PETSC_INT, NULL)); 1392 } 1393 1394 // For every face in connDistSection, find the transitive support of a vertex in that face connectivity. 1395 // Loop through the faces of the transitive support and find the matching face 1396 PetscBT plex_face_found; 1397 PetscInt fplexStart, fplexEnd, vplexStart, vplexEnd; 1398 PetscInt fdistStart, fdistEnd, numfdist; 1399 PetscCall(DMPlexGetHeightStratum(dm, 1, &fplexStart, &fplexEnd)); 1400 PetscCall(DMPlexGetDepthStratum(dm, 0, &vplexStart, &vplexEnd)); 1401 PetscCall(PetscSectionGetChart(connDistSection, &fdistStart, &fdistEnd)); 1402 numfdist = fdistEnd - fdistStart; 1403 PetscCall(PetscMalloc1(numfdist, plexFaces)); 1404 PetscCall(PetscBTCreate(numfdist, &plex_face_found)); 1405 for (PetscInt i = 0; i < numfdist; i++) (*plexFaces)[i] = -1; 1406 1407 for (PetscInt f = fdistStart, f_i = 0; f < fdistEnd; f++, f_i++) { 1408 PetscInt ndof, offset, support_size; 1409 PetscInt *support = NULL; 1410 1411 PetscCall(PetscSectionGetDof(connDistSection, f, &ndof)); 1412 PetscCall(PetscSectionGetOffset(connDistSection, f, &offset)); 1413 1414 // Loop through transitive support of a vertex in the CGNS face connectivity 1415 PetscCall(DMPlexGetTransitiveClosure(dm, connDist[offset + 0], PETSC_FALSE, &support_size, &support)); 1416 for (PetscInt s = 0; s < support_size; s++) { 1417 PetscInt face_point = support[s * 2]; // closure stores points and orientations, [p_0, o_0, p_1, o_1, ...] 1418 PetscInt trans_cone_size, *trans_cone = NULL; 1419 1420 if (face_point < fplexStart || face_point >= fplexEnd) continue; // Skip non-face points 1421 // See if face_point has the same vertices 1422 PetscCall(DMPlexGetTransitiveClosure(dm, face_point, PETSC_TRUE, &trans_cone_size, &trans_cone)); 1423 for (PetscInt c = 0; c < trans_cone_size; c++) { 1424 PetscInt vertex_point = trans_cone[c * 2], conn_has_vertex; 1425 if (vertex_point < vplexStart || vertex_point >= vplexEnd) continue; // Skip non-vertex points 1426 PetscCall(PetscFindIntUnsorted(vertex_point, ndof, &connDist[offset], &conn_has_vertex)); 1427 if (conn_has_vertex < 0) goto check_next_face; 1428 } 1429 (*plexFaces)[f_i] = face_point; 1430 PetscCall(DMPlexRestoreTransitiveClosure(dm, face_point, PETSC_TRUE, &trans_cone_size, &trans_cone)); 1431 break; 1432 check_next_face: 1433 PetscCall(DMPlexRestoreTransitiveClosure(dm, face_point, PETSC_TRUE, &trans_cone_size, &trans_cone)); 1434 continue; 1435 } 1436 PetscCall(DMPlexRestoreTransitiveClosure(dm, connDist[offset + 0], PETSC_FALSE, &support_size, &support)); 1437 if ((*plexFaces)[f_i] != -1) PetscCall(PetscBTSet(plex_face_found, f_i)); 1438 } 1439 1440 // Some distributed CGNS faces did not find a matching plex face 1441 // This can happen if a partition has all the faces surrounding a distributed CGNS face, but does not have the face itself (it's parent element is owned by a different partition). 1442 // Thus, the partition has the vertices associated with the CGNS face, but doesn't actually have the face itself. 1443 // For example, take the following quad mesh, where the numbers represent the owning rank and CGNS face ID. 1444 // 1445 // 2 3 4 <-- face ID 1446 // ----- ----- ----- 1447 // | 0 | 1 | 0 | <-- rank 1448 // ----- ----- ----- 1449 // 5 6 7 <-- face ID 1450 // 1451 // In this case, rank 0 will have all the vertices of face 3 and 6 in it's Plex, but does not actually have either face. 1452 // 1453 // To address this, we remove the leaves associated with these missing faces from cg2plexSF and then verify that all owned faces did find a matching plex face (e.g. root degree > 1) 1454 PetscCount num_plex_faces_found = PetscBTCountSet(plex_face_found, numfdist); 1455 PetscBool some_faces_not_found = num_plex_faces_found < numfdist; 1456 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &some_faces_not_found, 1, MPI_C_BOOL, MPI_LOR, comm)); 1457 if (some_faces_not_found) { 1458 PetscSFNode *iremote_cg2plex_new; 1459 const PetscInt *root_degree; 1460 PetscInt num_roots, *plexFacesNew; 1461 1462 PetscCall(PetscMalloc1(num_plex_faces_found, &iremote_cg2plex_new)); 1463 PetscCall(PetscCalloc1(num_plex_faces_found, &plexFacesNew)); 1464 { // Get SFNodes with matching faces 1465 const PetscSFNode *iremote_cg2plex_old; 1466 PetscInt num_leaves_old, n = 0; 1467 PetscCall(PetscSFGetGraph(*cg2plexSF, &num_roots, &num_leaves_old, NULL, &iremote_cg2plex_old)); 1468 PetscAssert(num_roots == fownedSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent roots and owned faces."); 1469 PetscAssert(num_leaves_old == numfdist, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent leaves and distributed faces."); 1470 for (PetscInt o = 0; o < num_leaves_old; o++) { 1471 if (PetscBTLookupSet(plex_face_found, o)) { 1472 iremote_cg2plex_new[n] = iremote_cg2plex_old[o]; 1473 plexFacesNew[n] = (*plexFaces)[o]; 1474 n++; 1475 } 1476 } 1477 PetscAssert(n == num_plex_faces_found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Found %" PetscCount_FMT " matching plex faces, but only set %" PetscInt_FMT " SFNodes", num_plex_faces_found, n); 1478 } 1479 PetscCall(PetscSFSetGraph(*cg2plexSF, num_roots, num_plex_faces_found, NULL, PETSC_COPY_VALUES, iremote_cg2plex_new, PETSC_OWN_POINTER)); 1480 1481 // Verify that all CGNS faces have a matching Plex face on any rank 1482 PetscCall(PetscSFComputeDegreeBegin(*cg2plexSF, &root_degree)); 1483 PetscCall(PetscSFComputeDegreeEnd(*cg2plexSF, &root_degree)); 1484 for (PetscInt r = 0; r < num_roots; r++) { 1485 PetscCheck(root_degree[r] > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find plex face for the CGNS face %" PetscInt_FMT, face_ids[r]); 1486 PetscCheck(root_degree[r] == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Found more than one plex face for the CGNS face %" PetscInt_FMT ". Face may be internal rather than at mesh domain boundary", face_ids[r]); 1487 } 1488 1489 if (PetscDefined(USE_DEBUG)) { 1490 for (PetscInt i = 0; i < num_plex_faces_found; i++) 1491 PetscCheck(plexFacesNew[i] >= fplexStart && plexFacesNew[i] < fplexEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Plex face ID %" PetscInt_FMT "outside of face stratum [%" PetscInt_FMT ", %" PetscInt_FMT ")", plexFacesNew[i], fplexStart, fplexEnd); 1492 } 1493 1494 PetscCall(PetscFree(*plexFaces)); 1495 *plexFaces = plexFacesNew; 1496 } 1497 1498 PetscCall(PetscBTDestroy(&plex_face_found)); 1499 PetscCall(PetscSectionDestroy(&connDistSection)); 1500 PetscCall(PetscFree(connDist)); 1501 } 1502 PetscFunctionReturn(PETSC_SUCCESS); 1503 } 1504 1505 // Copied from PetscOptionsStringToInt 1506 static inline PetscErrorCode PetscStrtoInt(const char name[], PetscInt *a) 1507 { 1508 size_t len; 1509 char *endptr; 1510 long strtolval; 1511 1512 PetscFunctionBegin; 1513 PetscCall(PetscStrlen(name, &len)); 1514 PetscCheck(len, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "character string of length zero has no numerical value"); 1515 1516 strtolval = strtol(name, &endptr, 10); 1517 PetscCheck((size_t)(endptr - name) == len, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input string \"%s\" has no integer value (do not include . in it)", name); 1518 1519 #if defined(PETSC_USE_64BIT_INDICES) && defined(PETSC_HAVE_ATOLL) 1520 (void)strtolval; 1521 *a = atoll(name); 1522 #elif defined(PETSC_USE_64BIT_INDICES) && defined(PETSC_HAVE___INT64) 1523 (void)strtolval; 1524 *a = _atoi64(name); 1525 #else 1526 *a = (PetscInt)strtolval; 1527 #endif 1528 PetscFunctionReturn(PETSC_SUCCESS); 1529 } 1530 1531 typedef struct { 1532 char name[CGIO_MAX_NAME_LENGTH + 1]; 1533 int normal[3], ndatasets; 1534 cgsize_t npoints, nnormals; 1535 CGNS_ENUMT(BCType_t) bctype; 1536 CGNS_ENUMT(DataType_t) normal_datatype; 1537 CGNS_ENUMT(PointSetType_t) pointtype; 1538 } CGBCInfo; 1539 1540 PetscErrorCode DMPlexCreateCGNS_Internal_Parallel(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm) 1541 { 1542 PetscMPIInt num_proc, rank; 1543 /* Read from file */ 1544 char basename[CGIO_MAX_NAME_LENGTH + 1]; 1545 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 1546 int dim = 0, physDim = 0, coordDim = 0; 1547 PetscInt NVertices = 0, NCells = 0; 1548 int nzones = 0, nbases; 1549 int zone = 1; // Only supports single zone files 1550 int base = 1; // Only supports single base 1551 1552 PetscFunctionBegin; 1553 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1554 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 1555 PetscCall(DMCreate(comm, dm)); 1556 PetscCall(DMSetType(*dm, DMPLEX)); 1557 1558 PetscCallCGNSRead(cg_nbases(cgid, &nbases), *dm, 0); 1559 PetscCheck(nbases <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single base, not %d", nbases); 1560 // From the CGNS web page cell_dim phys_dim (embedding space in PETSc) CGNS defines as length of spatial vectors/components) 1561 PetscCallCGNSRead(cg_base_read(cgid, base, basename, &dim, &physDim), *dm, 0); 1562 PetscCallCGNSRead(cg_nzones(cgid, base, &nzones), *dm, 0); 1563 PetscCheck(nzones == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Parallel reader limited to one zone, not %d", nzones); 1564 { 1565 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 1566 1567 PetscCallCGNSRead(cg_zone_read(cgid, base, zone, buffer, sizes), *dm, 0); 1568 NVertices = sizes[0]; 1569 NCells = sizes[1]; 1570 } 1571 1572 PetscCall(PetscObjectSetName((PetscObject)*dm, basename)); 1573 PetscCall(DMSetDimension(*dm, dim)); 1574 coordDim = dim; 1575 1576 // This is going to be a headache for mixed-topology and multiple sections. We may have to restore reading the data twice (once before the SetChart 1577 // call to get this right but continuing for now with single section, single topology, one zone. 1578 // establish element ranges for my rank 1579 PetscInt mystarte, myende, mystartv, myendv, myownede, myownedv; 1580 PetscLayout elem_map, vtx_map; 1581 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NCells, 1, &elem_map)); 1582 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NVertices, 1, &vtx_map)); 1583 PetscCall(PetscLayoutGetRange(elem_map, &mystarte, &myende)); 1584 PetscCall(PetscLayoutGetLocalSize(elem_map, &myownede)); 1585 PetscCall(PetscLayoutGetRange(vtx_map, &mystartv, &myendv)); 1586 PetscCall(PetscLayoutGetLocalSize(vtx_map, &myownedv)); 1587 1588 // -- Build Plex in parallel 1589 DMPolytopeType dm_cell_type = DM_POLYTOPE_UNKNOWN; 1590 PetscInt pOrder = 1, numClosure = -1; 1591 cgsize_t *elements = NULL; 1592 int *face_section_ids, *cell_section_ids, num_face_sections = 0, num_cell_sections = 0; 1593 PetscInt *uniq_verts, nuniq_verts; 1594 { 1595 int nsections; 1596 PetscInt *elementsQ1, numCorners = -1; 1597 const int *perm; 1598 cgsize_t start, end; // Throwaway 1599 1600 cg_nsections(cgid, base, zone, &nsections); 1601 PetscCall(PetscMalloc2(nsections, &face_section_ids, nsections, &cell_section_ids)); 1602 // Read element connectivity 1603 for (int index_sect = 1; index_sect <= nsections; index_sect++) { 1604 int nbndry, parentFlag; 1605 PetscInt cell_dim; 1606 CGNS_ENUMT(ElementType_t) cellType; 1607 1608 PetscCallCGNSRead(cg_section_read(cgid, base, zone, index_sect, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 1609 1610 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, &dm_cell_type, &numCorners, &cell_dim)); 1611 // Skip over element that are not max dimension (ie. boundary elements) 1612 if (cell_dim == dim) cell_section_ids[num_cell_sections++] = index_sect; 1613 else if (cell_dim == dim - 1) face_section_ids[num_face_sections++] = index_sect; 1614 } 1615 PetscCheck(num_cell_sections == 1, comm, PETSC_ERR_SUP, "CGNS Reader does not support more than 1 full-dimension cell section"); 1616 1617 { 1618 int index_sect = cell_section_ids[0], nbndry, parentFlag; 1619 CGNS_ENUMT(ElementType_t) cellType; 1620 1621 PetscCallCGNSRead(cg_section_read(cgid, base, zone, index_sect, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 1622 PetscCall(CGNSElementTypeGetDiscretizationInfo(cellType, &numClosure, &pOrder)); 1623 PetscCall(PetscMalloc1(myownede * numClosure, &elements)); 1624 PetscCallCGNSReadData(cgp_elements_read_data(cgid, base, zone, index_sect, mystarte + 1, myende, elements), *dm, 0); 1625 for (PetscInt v = 0; v < myownede * numClosure; ++v) elements[v] -= 1; // 0 based 1626 1627 // Create corners-only connectivity 1628 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, &dm_cell_type, &numCorners, NULL)); 1629 PetscCall(PetscMalloc1(myownede * numCorners, &elementsQ1)); 1630 PetscCall(DMPlexCGNSGetPermutation_Internal(dm_cell_type, numCorners, NULL, &perm)); 1631 for (PetscInt e = 0; e < myownede; ++e) { 1632 for (PetscInt v = 0; v < numCorners; ++v) elementsQ1[e * numCorners + perm[v]] = elements[e * numClosure + v]; 1633 } 1634 } 1635 1636 // Build cell-vertex Plex 1637 PetscCall(DMPlexBuildFromCellListParallel(*dm, myownede, myownedv, NVertices, numCorners, elementsQ1, NULL, &uniq_verts)); 1638 PetscCall(DMViewFromOptions(*dm, NULL, "-corner_dm_view")); 1639 { 1640 PetscInt pStart, pEnd; 1641 PetscCall(DMPlexGetChart(*dm, &pStart, &pEnd)); 1642 nuniq_verts = (pEnd - pStart) - myownede; 1643 } 1644 PetscCall(PetscFree(elementsQ1)); 1645 } 1646 1647 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(*dm)); 1648 1649 // -- Create SF for naive nodal-data read to elements 1650 PetscSF plex_to_cgns_sf; 1651 { 1652 PetscInt nleaves, num_comp; 1653 PetscInt *leaf, num_leaves = 0; 1654 PetscInt cStart, cEnd; 1655 const int *perm; 1656 PetscSF cgns_to_local_sf; 1657 PetscSection local_section; 1658 PetscFE fe; 1659 1660 // sfNatural requires PetscSection to handle DMDistribute, so we use PetscFE to define the section 1661 // Use number of components = 1 to work with just the nodes themselves 1662 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, dm_cell_type, pOrder, PETSC_DETERMINE, &fe)); 1663 PetscCall(PetscObjectSetName((PetscObject)fe, "FE for sfNatural")); 1664 PetscCall(DMAddField(*dm, NULL, (PetscObject)fe)); 1665 PetscCall(DMCreateDS(*dm)); 1666 PetscCall(PetscFEDestroy(&fe)); 1667 1668 PetscCall(DMGetLocalSection(*dm, &local_section)); 1669 PetscCall(PetscSectionViewFromOptions(local_section, NULL, "-fe_natural_section_view")); 1670 PetscCall(PetscSectionGetFieldComponents(local_section, 0, &num_comp)); 1671 PetscCall(PetscSectionGetStorageSize(local_section, &nleaves)); 1672 nleaves /= num_comp; 1673 PetscCall(PetscMalloc1(nleaves, &leaf)); 1674 for (PetscInt i = 0; i < nleaves; i++) leaf[i] = -1; 1675 1676 // Get the permutation from CGNS closure numbering to PLEX closure numbering 1677 PetscCall(DMPlexCGNSGetPermutation_Internal(dm_cell_type, numClosure, NULL, &perm)); 1678 PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd)); 1679 for (PetscInt cell = cStart; cell < cEnd; ++cell) { 1680 PetscInt num_closure_dof, *closure_idx = NULL; 1681 1682 PetscCall(DMPlexGetClosureIndices(*dm, local_section, local_section, cell, PETSC_FALSE, &num_closure_dof, &closure_idx, NULL, NULL)); 1683 PetscAssert(numClosure * num_comp == num_closure_dof, comm, PETSC_ERR_PLIB, "Closure dof size does not match polytope"); 1684 for (PetscInt i = 0; i < numClosure; i++) { 1685 PetscInt li = closure_idx[perm[i] * num_comp] / num_comp; 1686 if (li < 0) continue; 1687 1688 PetscInt cgns_idx = elements[cell * numClosure + i]; 1689 if (leaf[li] == -1) { 1690 leaf[li] = cgns_idx; 1691 num_leaves++; 1692 } else PetscAssert(leaf[li] == cgns_idx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "leaf does not match previously set"); 1693 } 1694 PetscCall(DMPlexRestoreClosureIndices(*dm, local_section, local_section, cell, PETSC_FALSE, &num_closure_dof, &closure_idx, NULL, NULL)); 1695 } 1696 PetscAssert(num_leaves == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "leaf count in closure does not match nleaves"); 1697 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)*dm), &cgns_to_local_sf)); 1698 PetscCall(PetscSFSetGraphLayout(cgns_to_local_sf, vtx_map, nleaves, NULL, PETSC_USE_POINTER, leaf)); 1699 PetscCall(PetscObjectSetName((PetscObject)cgns_to_local_sf, "CGNS to Plex SF")); 1700 PetscCall(PetscSFViewFromOptions(cgns_to_local_sf, NULL, "-CGNStoPlex_sf_view")); 1701 PetscCall(PetscFree(leaf)); 1702 PetscCall(PetscFree(elements)); 1703 1704 { // Convert cgns_to_local to global_to_cgns 1705 PetscSF sectionsf, cgns_to_global_sf; 1706 1707 PetscCall(DMGetSectionSF(*dm, §ionsf)); 1708 PetscCall(PetscSFComposeInverse(cgns_to_local_sf, sectionsf, &cgns_to_global_sf)); 1709 PetscCall(PetscSFDestroy(&cgns_to_local_sf)); 1710 PetscCall(PetscSFCreateInverseSF(cgns_to_global_sf, &plex_to_cgns_sf)); 1711 PetscCall(PetscObjectSetName((PetscObject)plex_to_cgns_sf, "Global Plex to CGNS")); 1712 PetscCall(PetscSFDestroy(&cgns_to_global_sf)); 1713 } 1714 } 1715 1716 { // -- Set coordinates for DM 1717 PetscScalar *coords; 1718 float *x[3]; 1719 double *xd[3]; 1720 PetscBool read_with_double; 1721 PetscFE cfe; 1722 1723 // Setup coordinate space first. Use pOrder here for isoparametric; revisit with CPEX-0045 High Order. 1724 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, coordDim, dm_cell_type, pOrder, PETSC_DETERMINE, &cfe)); 1725 PetscCall(DMSetCoordinateDisc(*dm, cfe, PETSC_FALSE, PETSC_FALSE)); 1726 PetscCall(PetscFEDestroy(&cfe)); 1727 1728 { // Determine if coords are written in single or double precision 1729 CGNS_ENUMT(DataType_t) datatype; 1730 1731 PetscCallCGNSRead(cg_coord_info(cgid, base, zone, 1, &datatype, buffer), *dm, 0); 1732 read_with_double = datatype == CGNS_ENUMV(RealDouble) ? PETSC_TRUE : PETSC_FALSE; 1733 } 1734 1735 // Read coords from file and set into component-major ordering 1736 if (read_with_double) PetscCall(PetscMalloc3(myownedv, &xd[0], myownedv, &xd[1], myownedv, &xd[2])); 1737 else PetscCall(PetscMalloc3(myownedv, &x[0], myownedv, &x[1], myownedv, &x[2])); 1738 PetscCall(PetscMalloc1(myownedv * coordDim, &coords)); 1739 { 1740 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 1741 cgsize_t range_min[3] = {mystartv + 1, 1, 1}; 1742 cgsize_t range_max[3] = {myendv, 1, 1}; 1743 int ngrids, ncoords; 1744 1745 PetscCallCGNSRead(cg_zone_read(cgid, base, zone, buffer, sizes), *dm, 0); 1746 PetscCallCGNSRead(cg_ngrids(cgid, base, zone, &ngrids), *dm, 0); 1747 PetscCheck(ngrids <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single grid, not %d", ngrids); 1748 PetscCallCGNSRead(cg_ncoords(cgid, base, zone, &ncoords), *dm, 0); 1749 PetscCheck(ncoords == coordDim, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a coordinate array for each dimension, not %d", ncoords); 1750 if (read_with_double) { 1751 for (int d = 0; d < coordDim; ++d) PetscCallCGNSReadData(cgp_coord_read_data(cgid, base, zone, (d + 1), range_min, range_max, xd[d]), *dm, 0); 1752 if (coordDim >= 1) { 1753 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 0] = xd[0][v]; 1754 } 1755 if (coordDim >= 2) { 1756 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 1] = xd[1][v]; 1757 } 1758 if (coordDim >= 3) { 1759 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 2] = xd[2][v]; 1760 } 1761 } else { 1762 for (int d = 0; d < coordDim; ++d) PetscCallCGNSReadData(cgp_coord_read_data(cgid, 1, zone, (d + 1), range_min, range_max, x[d]), *dm, 0); 1763 if (coordDim >= 1) { 1764 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 0] = x[0][v]; 1765 } 1766 if (coordDim >= 2) { 1767 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 1] = x[1][v]; 1768 } 1769 if (coordDim >= 3) { 1770 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 2] = x[2][v]; 1771 } 1772 } 1773 } 1774 if (read_with_double) PetscCall(PetscFree3(xd[0], xd[1], xd[2])); 1775 else PetscCall(PetscFree3(x[0], x[1], x[2])); 1776 1777 { // Reduce CGNS-ordered coordinate nodes to Plex ordering, and set DM's coordinates 1778 Vec coord_global; 1779 MPI_Datatype unit; 1780 PetscScalar *coord_global_array; 1781 DM cdm; 1782 1783 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 1784 PetscCall(DMCreateGlobalVector(cdm, &coord_global)); 1785 PetscCall(VecGetArrayWrite(coord_global, &coord_global_array)); 1786 PetscCallMPI(MPI_Type_contiguous(coordDim, MPIU_SCALAR, &unit)); 1787 PetscCallMPI(MPI_Type_commit(&unit)); 1788 PetscCall(PetscSFReduceBegin(plex_to_cgns_sf, unit, coords, coord_global_array, MPI_REPLACE)); 1789 PetscCall(PetscSFReduceEnd(plex_to_cgns_sf, unit, coords, coord_global_array, MPI_REPLACE)); 1790 PetscCall(VecRestoreArrayWrite(coord_global, &coord_global_array)); 1791 PetscCallMPI(MPI_Type_free(&unit)); 1792 PetscCall(DMSetCoordinates(*dm, coord_global)); 1793 PetscCall(VecDestroy(&coord_global)); 1794 } 1795 PetscCall(PetscFree(coords)); 1796 } 1797 1798 PetscCall(DMViewFromOptions(*dm, NULL, "-corner_interpolated_dm_view")); 1799 1800 int nbocos; 1801 PetscCallCGNSRead(cg_nbocos(cgid, base, zone, &nbocos), *dm, 0); 1802 // In order to extract boundary condition (boco) information into DMLabels, each rank holds: 1803 // - The local Plex 1804 // - Naively read CGNS face connectivity 1805 // - Naively read list of CGNS faces for each boco 1806 // 1807 // First, we need to build a mapping from the CGNS faces to the (probably off-rank) Plex face. 1808 // The CGNS faces that each rank owns is known globally via cgnsLayouts. 1809 // The cg2plexSF maps these CGNS face IDs to their (probably off-rank) Plex face. 1810 // The plexFaces array maps the (contiguous) leaves of cg2plexSF to the local Plex face point. 1811 // 1812 // Next, we read the list of CGNS faces for each boco and find the location of that face's owner using cgnsLayouts. 1813 // Then, we can communicate the label value to the local Plex which corresponds to the CGNS face. 1814 if (interpolate && num_face_sections != 0 && nbocos != 0) { 1815 PetscSection connSection; 1816 PetscInt nCgFaces, nPlexFaces; 1817 PetscInt *face_ids, *conn, *plexFaces; 1818 PetscSF cg2plexSF; 1819 PetscLayout *cgnsLayouts; 1820 1821 PetscCall(DMPlexCGNS_CreateCornersConnectivitySection(*dm, cgid, base, zone, num_face_sections, face_section_ids, &connSection, NULL, &face_ids, &cgnsLayouts, &conn)); 1822 { 1823 PetscBool view_connectivity = PETSC_FALSE; 1824 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_cgns_view_face_connectivity", &view_connectivity, NULL)); 1825 if (view_connectivity) PetscCall(PetscSectionArrayView(connSection, conn, PETSC_INT, NULL)); 1826 } 1827 PetscCall(DMPlexCGNS_MatchCGNSFacesToPlexFaces(*dm, nuniq_verts, uniq_verts, myownede, NVertices, connSection, face_ids, conn, &cg2plexSF, &plexFaces)); 1828 PetscCall(PetscSFGetGraph(cg2plexSF, NULL, &nPlexFaces, NULL, NULL)); 1829 { 1830 PetscInt start, end; 1831 PetscCall(PetscSectionGetChart(connSection, &start, &end)); 1832 nCgFaces = end - start; 1833 } 1834 1835 PetscInt *plexFaceValues, *cgFaceValues; 1836 PetscCall(PetscMalloc2(nPlexFaces, &plexFaceValues, nCgFaces, &cgFaceValues)); 1837 for (PetscInt BC = 1; BC <= nbocos; BC++) { 1838 cgsize_t *points; 1839 CGBCInfo bcinfo; 1840 PetscBool is_faceset = PETSC_FALSE; 1841 PetscInt label_value = 1; 1842 1843 PetscCallCGNSRead(cg_boco_info(cgid, base, zone, BC, bcinfo.name, &bcinfo.bctype, &bcinfo.pointtype, &bcinfo.npoints, bcinfo.normal, &bcinfo.nnormals, &bcinfo.normal_datatype, &bcinfo.ndatasets), *dm, 0); 1844 1845 PetscCall(PetscStrbeginswith(bcinfo.name, "FaceSet", &is_faceset)); 1846 if (is_faceset) { 1847 size_t faceset_len; 1848 PetscCall(PetscStrlen("FaceSet", &faceset_len)); 1849 PetscCall(PetscStrtoInt(bcinfo.name + faceset_len, &label_value)); 1850 } 1851 const char *label_name = is_faceset ? "Face Sets" : bcinfo.name; 1852 1853 if (bcinfo.npoints < 1) continue; 1854 1855 PetscLayout bc_layout; 1856 PetscInt bcStart, bcEnd, bcSize; 1857 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, bcinfo.npoints, 1, &bc_layout)); 1858 PetscCall(PetscLayoutGetRange(bc_layout, &bcStart, &bcEnd)); 1859 PetscCall(PetscLayoutGetLocalSize(bc_layout, &bcSize)); 1860 PetscCall(PetscLayoutDestroy(&bc_layout)); 1861 PetscCall(DMGetWorkArray(*dm, bcSize, MPIU_CGSIZE, &points)); 1862 1863 const char *labels[] = {"Zone_t", "ZoneBC_t", "BC_t", "PointList"}; 1864 PetscCallCGNSRead(cg_golist(cgid, base, 4, (char **)labels, (int[]){zone, 1, BC, 0}), *dm, 0); 1865 PetscCallCGNSReadData(cgp_ptlist_read_data(cgid, bcStart + 1, bcEnd, points), *dm, 0); 1866 1867 PetscInt *label_values; 1868 PetscSFNode *remotes; 1869 PetscCall(PetscMalloc2(bcSize, &remotes, bcSize, &label_values)); 1870 for (PetscInt p = 0; p < bcSize; p++) { 1871 PetscMPIInt bcrank; 1872 PetscInt bcidx; 1873 1874 PetscCall(PetscLayoutFindOwnerIndex_CGNSSectionLayouts(cgnsLayouts, num_face_sections, points[p], &bcrank, &bcidx, NULL)); 1875 remotes[p].rank = bcrank; 1876 remotes[p].index = bcidx; 1877 label_values[p] = label_value; 1878 } 1879 PetscCall(DMRestoreWorkArray(*dm, bcSize, MPIU_CGSIZE, &points)); 1880 1881 { // Communicate the BC values to their Plex-face owners 1882 PetscSF cg2bcSF; 1883 DMLabel label; 1884 1885 for (PetscInt i = 0; i < nCgFaces; i++) cgFaceValues[i] = -1; 1886 for (PetscInt i = 0; i < nPlexFaces; i++) plexFaceValues[i] = -1; 1887 1888 PetscCall(PetscSFCreate(comm, &cg2bcSF)); 1889 PetscCall(PetscSFSetGraph(cg2bcSF, nCgFaces, bcSize, NULL, PETSC_COPY_VALUES, remotes, PETSC_USE_POINTER)); 1890 1891 PetscCall(PetscSFReduceBegin(cg2bcSF, MPIU_INT, label_values, cgFaceValues, MPI_REPLACE)); 1892 PetscCall(PetscSFReduceEnd(cg2bcSF, MPIU_INT, label_values, cgFaceValues, MPI_REPLACE)); 1893 PetscCall(PetscSFBcastBegin(cg2plexSF, MPIU_INT, cgFaceValues, plexFaceValues, MPI_REPLACE)); 1894 PetscCall(PetscSFBcastEnd(cg2plexSF, MPIU_INT, cgFaceValues, plexFaceValues, MPI_REPLACE)); 1895 PetscCall(PetscSFDestroy(&cg2bcSF)); 1896 PetscCall(PetscFree2(remotes, label_values)); 1897 1898 // Set the label values for the communicated faces 1899 PetscCall(DMGetLabel(*dm, label_name, &label)); 1900 if (label == NULL) { 1901 PetscCall(DMCreateLabel(*dm, label_name)); 1902 PetscCall(DMGetLabel(*dm, label_name, &label)); 1903 } 1904 for (PetscInt i = 0; i < nPlexFaces; i++) { 1905 if (plexFaceValues[i] == -1) continue; 1906 PetscCall(DMLabelSetValue(label, plexFaces[i], plexFaceValues[i])); 1907 } 1908 } 1909 } 1910 PetscCall(PetscFree2(plexFaceValues, cgFaceValues)); 1911 PetscCall(PetscFree(plexFaces)); 1912 PetscCall(PetscSFDestroy(&cg2plexSF)); 1913 PetscCall(PetscFree(conn)); 1914 for (PetscInt s = 0; s < num_face_sections; s++) PetscCall(PetscLayoutDestroy(&cgnsLayouts[s])); 1915 PetscCall(PetscSectionDestroy(&connSection)); 1916 PetscCall(PetscFree(cgnsLayouts)); 1917 PetscCall(PetscFree(face_ids)); 1918 } 1919 PetscCall(PetscFree(uniq_verts)); 1920 PetscCall(PetscFree2(face_section_ids, cell_section_ids)); 1921 1922 // -- Set sfNatural for solution vectors in CGNS file 1923 // NOTE: We set sfNatural to be the map between the original CGNS ordering of nodes and the Plex ordering of nodes. 1924 PetscCall(PetscSFViewFromOptions(plex_to_cgns_sf, NULL, "-sfNatural_init_view")); 1925 PetscCall(DMSetNaturalSF(*dm, plex_to_cgns_sf)); 1926 PetscCall(DMSetUseNatural(*dm, PETSC_TRUE)); 1927 PetscCall(PetscSFDestroy(&plex_to_cgns_sf)); 1928 1929 PetscCall(PetscLayoutDestroy(&elem_map)); 1930 PetscCall(PetscLayoutDestroy(&vtx_map)); 1931 PetscFunctionReturn(PETSC_SUCCESS); 1932 } 1933 1934 // node_l2g must be freed 1935 static PetscErrorCode DMPlexCreateNodeNumbering(DM dm, PetscInt *num_local_nodes, PetscInt *num_global_nodes, PetscInt *nStart, PetscInt *nEnd, const PetscInt **node_l2g) 1936 { 1937 PetscSection local_section; 1938 PetscSF point_sf; 1939 PetscInt pStart, pEnd, spStart, spEnd, *points, nleaves, ncomp, *nodes; 1940 PetscMPIInt comm_size; 1941 const PetscInt *ilocal, field = 0; 1942 1943 PetscFunctionBegin; 1944 *num_local_nodes = -1; 1945 *num_global_nodes = -1; 1946 *nStart = -1; 1947 *nEnd = -1; 1948 *node_l2g = NULL; 1949 1950 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &comm_size)); 1951 PetscCall(DMGetLocalSection(dm, &local_section)); 1952 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1953 PetscCall(PetscSectionGetChart(local_section, &spStart, &spEnd)); 1954 PetscCall(DMGetPointSF(dm, &point_sf)); 1955 if (comm_size == 1) nleaves = 0; 1956 else PetscCall(PetscSFGetGraph(point_sf, NULL, &nleaves, &ilocal, NULL)); 1957 PetscCall(PetscSectionGetFieldComponents(local_section, field, &ncomp)); 1958 1959 PetscInt local_node = 0, owned_node = 0, owned_start = 0; 1960 for (PetscInt p = spStart, leaf = 0; p < spEnd; p++) { 1961 PetscInt dof; 1962 PetscCall(PetscSectionGetFieldDof(local_section, p, field, &dof)); 1963 PetscAssert(dof % ncomp == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field dof %" PetscInt_FMT " must be divisible by components %" PetscInt_FMT, dof, ncomp); 1964 local_node += dof / ncomp; 1965 if (leaf < nleaves && p == ilocal[leaf]) { // skip points owned by a different process 1966 leaf++; 1967 } else { 1968 owned_node += dof / ncomp; 1969 } 1970 } 1971 PetscCallMPI(MPI_Exscan(&owned_node, &owned_start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm))); 1972 PetscCall(PetscMalloc1(pEnd - pStart, &points)); 1973 owned_node = 0; 1974 for (PetscInt p = spStart, leaf = 0; p < spEnd; p++) { 1975 if (leaf < nleaves && p == ilocal[leaf]) { // skip points owned by a different process 1976 points[p - pStart] = -1; 1977 leaf++; 1978 continue; 1979 } 1980 PetscInt dof, offset; 1981 PetscCall(PetscSectionGetFieldDof(local_section, p, field, &dof)); 1982 PetscCall(PetscSectionGetFieldOffset(local_section, p, field, &offset)); 1983 PetscAssert(offset % ncomp == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field offset %" PetscInt_FMT " must be divisible by components %" PetscInt_FMT, offset, ncomp); 1984 points[p - pStart] = owned_start + owned_node; 1985 owned_node += dof / ncomp; 1986 } 1987 if (comm_size > 1) { 1988 PetscCall(PetscSFBcastBegin(point_sf, MPIU_INT, points, points, MPI_REPLACE)); 1989 PetscCall(PetscSFBcastEnd(point_sf, MPIU_INT, points, points, MPI_REPLACE)); 1990 } 1991 1992 // Set up global indices for each local node 1993 PetscCall(PetscMalloc1(local_node, &nodes)); 1994 for (PetscInt p = spStart; p < spEnd; p++) { 1995 PetscInt dof, offset; 1996 PetscCall(PetscSectionGetFieldDof(local_section, p, field, &dof)); 1997 PetscCall(PetscSectionGetFieldOffset(local_section, p, field, &offset)); 1998 for (PetscInt n = 0; n < dof / ncomp; n++) nodes[offset / ncomp + n] = points[p - pStart] + n; 1999 } 2000 PetscCall(PetscFree(points)); 2001 *num_local_nodes = local_node; 2002 *nStart = owned_start; 2003 *nEnd = owned_start + owned_node; 2004 PetscCallMPI(MPIU_Allreduce(&owned_node, num_global_nodes, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm))); 2005 *node_l2g = nodes; 2006 PetscFunctionReturn(PETSC_SUCCESS); 2007 } 2008 2009 PetscErrorCode DMView_PlexCGNS(DM dm, PetscViewer viewer) 2010 { 2011 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 2012 PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data; 2013 PetscInt fvGhostStart; 2014 PetscInt topo_dim, coord_dim, num_global_elems; 2015 PetscInt cStart, cEnd, num_local_nodes, num_global_nodes, nStart, nEnd, fStart, fEnd; 2016 const PetscInt *node_l2g; 2017 Vec coord; 2018 DM colloc_dm, cdm; 2019 PetscMPIInt size; 2020 const char *dm_name; 2021 int base, zone; 2022 cgsize_t isize[3], elem_offset = 0; 2023 2024 PetscFunctionBegin; 2025 if (!cgv->file_num) { 2026 PetscInt time_step; 2027 PetscCall(DMGetOutputSequenceNumber(dm, &time_step, NULL)); 2028 PetscCall(PetscViewerCGNSFileOpen_Internal(viewer, time_step)); 2029 } 2030 PetscCallMPI(MPI_Comm_size(comm, &size)); 2031 PetscCall(DMGetDimension(dm, &topo_dim)); 2032 PetscCall(DMGetCoordinateDim(dm, &coord_dim)); 2033 PetscCall(PetscObjectGetName((PetscObject)dm, &dm_name)); 2034 PetscCallCGNSWrite(cg_base_write(cgv->file_num, dm_name, topo_dim, coord_dim, &base), dm, viewer); 2035 PetscCallCGNS(cg_goto(cgv->file_num, base, NULL)); 2036 PetscCallCGNSWrite(cg_dataclass_write(CGNS_ENUMV(NormalizedByDimensional)), dm, viewer); 2037 2038 { 2039 PetscFE fe, fe_coord; 2040 PetscClassId ds_id; 2041 PetscDualSpace dual_space, dual_space_coord; 2042 PetscInt num_fields, field_order = -1, field_order_coord; 2043 PetscBool is_simplex; 2044 PetscCall(DMGetNumFields(dm, &num_fields)); 2045 if (num_fields > 0) { 2046 PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe)); 2047 PetscCall(PetscObjectGetClassId((PetscObject)fe, &ds_id)); 2048 if (ds_id != PETSCFE_CLASSID) { 2049 fe = NULL; 2050 if (ds_id == PETSCFV_CLASSID) field_order = -1; // use whatever is present for coords; field will be CellCenter 2051 else field_order = 1; // assume vertex-based linear elements 2052 } 2053 } else fe = NULL; 2054 if (fe) { 2055 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 2056 PetscCall(PetscDualSpaceGetOrder(dual_space, &field_order)); 2057 } 2058 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2059 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe_coord)); 2060 { 2061 PetscClassId id; 2062 PetscCall(PetscObjectGetClassId((PetscObject)fe_coord, &id)); 2063 if (id != PETSCFE_CLASSID) fe_coord = NULL; 2064 } 2065 if (fe_coord) { 2066 PetscCall(PetscFEGetDualSpace(fe_coord, &dual_space_coord)); 2067 PetscCall(PetscDualSpaceGetOrder(dual_space_coord, &field_order_coord)); 2068 } else field_order_coord = 1; 2069 if (field_order > 0 && field_order != field_order_coord) { 2070 PetscInt quadrature_order = field_order; 2071 PetscCall(DMClone(dm, &colloc_dm)); 2072 { // Inform the new colloc_dm that it is a coordinate DM so isoperiodic affine corrections can be applied 2073 const PetscSF *face_sfs; 2074 PetscInt num_face_sfs; 2075 PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &num_face_sfs, &face_sfs)); 2076 PetscCall(DMPlexSetIsoperiodicFaceSF(colloc_dm, num_face_sfs, (PetscSF *)face_sfs)); 2077 if (face_sfs) colloc_dm->periodic.setup = DMPeriodicCoordinateSetUp_Internal; 2078 } 2079 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 2080 PetscCall(PetscFECreateLagrange(comm, topo_dim, coord_dim, is_simplex, field_order, quadrature_order, &fe)); 2081 PetscCall(DMSetCoordinateDisc(colloc_dm, fe, PETSC_FALSE, PETSC_TRUE)); 2082 PetscCall(PetscFEDestroy(&fe)); 2083 } else { 2084 PetscCall(PetscObjectReference((PetscObject)dm)); 2085 colloc_dm = dm; 2086 } 2087 } 2088 PetscCall(DMGetCoordinateDM(colloc_dm, &cdm)); 2089 PetscCall(DMPlexCreateNodeNumbering(cdm, &num_local_nodes, &num_global_nodes, &nStart, &nEnd, &node_l2g)); 2090 PetscCall(DMGetCoordinatesLocal(colloc_dm, &coord)); 2091 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2092 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 2093 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &fvGhostStart, NULL)); 2094 if (fvGhostStart >= 0) cEnd = fvGhostStart; 2095 num_global_elems = cEnd - cStart; 2096 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &num_global_elems, 1, MPIU_INT, MPI_SUM, comm)); 2097 isize[0] = num_global_nodes; 2098 isize[1] = num_global_elems; 2099 isize[2] = 0; 2100 PetscCallCGNSWrite(cg_zone_write(cgv->file_num, base, "Zone", isize, CGNS_ENUMV(Unstructured), &zone), dm, viewer); 2101 2102 cgsize_t e_owned, e_global, e_start; 2103 { 2104 const PetscScalar *X; 2105 PetscScalar *x; 2106 int coord_ids[3]; 2107 2108 PetscCall(VecGetArrayRead(coord, &X)); 2109 for (PetscInt d = 0; d < coord_dim; d++) { 2110 const double exponents[] = {0, 1, 0, 0, 0}; 2111 char coord_name[64]; 2112 PetscCall(PetscSNPrintf(coord_name, sizeof coord_name, "Coordinate%c", 'X' + (int)d)); 2113 PetscCallCGNSWrite(cgp_coord_write(cgv->file_num, base, zone, CGNS_ENUMV(RealDouble), coord_name, &coord_ids[d]), dm, viewer); 2114 PetscCallCGNS(cg_goto(cgv->file_num, base, "Zone_t", zone, "GridCoordinates", 0, coord_name, 0, NULL)); 2115 PetscCallCGNSWrite(cg_exponents_write(CGNS_ENUMV(RealDouble), exponents), dm, viewer); 2116 } 2117 2118 int section; 2119 cgsize_t *conn = NULL; 2120 const int *perm; 2121 CGNS_ENUMT(ElementType_t) element_type = CGNS_ENUMV(ElementTypeNull); 2122 { 2123 PetscCall(PetscMalloc1(nEnd - nStart, &x)); 2124 for (PetscInt d = 0; d < coord_dim; d++) { 2125 for (PetscInt n = 0; n < num_local_nodes; n++) { 2126 PetscInt gn = node_l2g[n]; 2127 if (gn < nStart || nEnd <= gn) continue; 2128 x[gn - nStart] = X[n * coord_dim + d]; 2129 } 2130 // CGNS nodes use 1-based indexing 2131 cgsize_t start = nStart + 1, end = nEnd; 2132 PetscCallCGNSWriteData(cgp_coord_write_data(cgv->file_num, base, zone, coord_ids[d], &start, &end, x), dm, viewer); 2133 } 2134 PetscCall(PetscFree(x)); 2135 PetscCall(VecRestoreArrayRead(coord, &X)); 2136 } 2137 2138 e_owned = cEnd - cStart; 2139 if (e_owned > 0) { 2140 DMPolytopeType cell_type; 2141 2142 PetscCall(DMPlexGetCellType(dm, cStart, &cell_type)); 2143 for (PetscInt i = cStart, c = 0; i < cEnd; i++) { 2144 PetscInt closure_dof, *closure_indices, elem_size; 2145 2146 PetscCall(DMPlexGetClosureIndices(cdm, cdm->localSection, cdm->localSection, i, PETSC_FALSE, &closure_dof, &closure_indices, NULL, NULL)); 2147 elem_size = closure_dof / coord_dim; 2148 if (!conn) PetscCall(PetscMalloc1(e_owned * elem_size, &conn)); 2149 PetscCall(DMPlexCGNSGetPermutation_Internal(cell_type, closure_dof / coord_dim, &element_type, &perm)); 2150 for (PetscInt j = 0; j < elem_size; j++) conn[c++] = node_l2g[closure_indices[perm[j] * coord_dim] / coord_dim] + 1; 2151 PetscCall(DMPlexRestoreClosureIndices(cdm, cdm->localSection, cdm->localSection, i, PETSC_FALSE, &closure_dof, &closure_indices, NULL, NULL)); 2152 } 2153 } 2154 2155 { // Get global element_type (for ranks that do not have owned elements) 2156 PetscInt local_element_type, global_element_type; 2157 2158 local_element_type = e_owned > 0 ? (PetscInt)element_type : -1; 2159 PetscCallMPI(MPIU_Allreduce(&local_element_type, &global_element_type, 1, MPIU_INT, MPI_MAX, comm)); 2160 if (local_element_type != -1) 2161 PetscCheck(local_element_type == global_element_type, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ranks with different element types not supported. Local element type is %s, but global is %s", cg_ElementTypeName(local_element_type), cg_ElementTypeName(global_element_type)); 2162 element_type = (CGNS_ENUMT(ElementType_t))global_element_type; 2163 } 2164 PetscCallMPI(MPIU_Allreduce(&e_owned, &e_global, 1, MPIU_CGSIZE, MPI_SUM, comm)); 2165 PetscCheck(e_global == num_global_elems, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of elements %" PRIdCGSIZE " vs %" PetscInt_FMT, e_global, num_global_elems); 2166 e_start = 0; 2167 PetscCallMPI(MPI_Exscan(&e_owned, &e_start, 1, MPIU_CGSIZE, MPI_SUM, comm)); 2168 e_start += elem_offset; 2169 PetscCallCGNSWrite(cgp_section_write(cgv->file_num, base, zone, "Elem", element_type, 1, e_global, 0, §ion), dm, viewer); 2170 PetscCallCGNSWriteData(cgp_elements_write_data(cgv->file_num, base, zone, section, e_start + 1, e_start + e_owned, conn), dm, viewer); 2171 elem_offset = e_global; 2172 PetscCall(PetscFree(conn)); 2173 2174 cgv->base = base; 2175 cgv->zone = zone; 2176 cgv->node_l2g = node_l2g; 2177 cgv->num_local_nodes = num_local_nodes; 2178 cgv->nStart = nStart; 2179 cgv->nEnd = nEnd; 2180 cgv->eStart = e_start; 2181 cgv->eEnd = e_start + e_owned; 2182 if (1) { 2183 PetscMPIInt rank; 2184 int *efield; 2185 int sol, field; 2186 DMLabel label; 2187 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2188 PetscCall(PetscMalloc1(e_owned, &efield)); 2189 for (PetscInt i = 0; i < e_owned; i++) efield[i] = rank; 2190 PetscCallCGNSWrite(cg_sol_write(cgv->file_num, base, zone, "CellInfo", CGNS_ENUMV(CellCenter), &sol), dm, viewer); 2191 PetscCallCGNSWrite(cgp_field_write(cgv->file_num, base, zone, sol, CGNS_ENUMV(Integer), "Rank", &field), dm, viewer); 2192 cgsize_t start = e_start + 1, end = e_start + e_owned; 2193 PetscCallCGNSWriteData(cgp_field_write_data(cgv->file_num, base, zone, sol, field, &start, &end, efield), dm, viewer); 2194 PetscCall(DMGetLabel(dm, "Cell Sets", &label)); 2195 if (label) { 2196 for (PetscInt c = cStart; c < cEnd; c++) { 2197 PetscInt value; 2198 PetscCall(DMLabelGetValue(label, c, &value)); 2199 efield[c - cStart] = value; 2200 } 2201 PetscCallCGNSWrite(cgp_field_write(cgv->file_num, base, zone, sol, CGNS_ENUMV(Integer), "CellSet", &field), dm, viewer); 2202 PetscCallCGNSWriteData(cgp_field_write_data(cgv->file_num, base, zone, sol, field, &start, &end, efield), dm, viewer); 2203 } 2204 PetscCall(PetscFree(efield)); 2205 } 2206 } 2207 2208 DMLabel fsLabel; 2209 PetscInt num_fs_global; 2210 IS fsValuesGlobalIS; 2211 PetscCall(DMGetLabel(dm, "Face Sets", &fsLabel)); 2212 PetscCall(DMLabelGetValueISGlobal(comm, fsLabel, PETSC_TRUE, &fsValuesGlobalIS)); 2213 PetscCall(ISGetSize(fsValuesGlobalIS, &num_fs_global)); 2214 2215 if (num_fs_global > 0) { 2216 CGNS_ENUMT(ElementType_t) element_type = CGNS_ENUMV(ElementTypeNull); 2217 const PetscInt *fsValuesLocal; 2218 IS stratumIS, fsFacesAll; 2219 int section; 2220 const int *perm; 2221 cgsize_t f_owned = 0, f_global, f_start; 2222 cgsize_t *parents, *conn = NULL; 2223 PetscInt fStart, fEnd; 2224 2225 PetscInt num_fs_local; 2226 IS fsValuesLocalIS; 2227 2228 if (fsLabel) { 2229 PetscCall(DMLabelGetNonEmptyStratumValuesIS(fsLabel, &fsValuesLocalIS)); 2230 PetscCall(ISGetSize(fsValuesLocalIS, &num_fs_local)); 2231 PetscCall(ISGetIndices(fsValuesLocalIS, &fsValuesLocal)); 2232 } else num_fs_local = 0; 2233 2234 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 2235 { // Get single IS without duplicates of the local face IDs in the FaceSets 2236 IS *fsPoints = NULL; 2237 2238 PetscCall(PetscMalloc1(num_fs_local, &fsPoints)); 2239 for (PetscInt fs = 0; fs < num_fs_local; ++fs) PetscCall(DMLabelGetStratumIS(fsLabel, fsValuesLocal[fs], &fsPoints[fs])); 2240 PetscCall(ISConcatenate(PETSC_COMM_SELF, num_fs_local, fsPoints, &fsFacesAll)); 2241 PetscCall(ISSortRemoveDups(fsFacesAll)); 2242 PetscCall(ISGeneralFilter(fsFacesAll, fStart, fEnd)); // Remove non-face mesh points from the IS 2243 { 2244 PetscInt f_owned_int; 2245 PetscCall(ISGetSize(fsFacesAll, &f_owned_int)); 2246 f_owned = f_owned_int; 2247 } 2248 for (PetscInt fs = 0; fs < num_fs_local; ++fs) PetscCall(ISDestroy(&fsPoints[fs])); 2249 PetscCall(PetscFree(fsPoints)); 2250 } 2251 PetscCall(ISRestoreIndices(fsValuesLocalIS, &fsValuesLocal)); 2252 PetscCall(ISDestroy(&fsValuesLocalIS)); 2253 2254 { 2255 const PetscInt *faces; 2256 DMPolytopeType cell_type, cell_type_f; 2257 PetscInt closure_dof = -1, closure_dof_f; 2258 2259 PetscCall(ISGetIndices(fsFacesAll, &faces)); 2260 if (f_owned) PetscCall(DMPlexGetCellType(dm, faces[0], &cell_type)); 2261 PetscCall(PetscCalloc1(f_owned * 2, &parents)); 2262 for (PetscInt f = 0, c = 0; f < f_owned; f++) { 2263 PetscInt *closure_indices, elem_size; 2264 const PetscInt face = faces[f]; 2265 2266 PetscCall(DMPlexGetCellType(dm, face, &cell_type_f)); 2267 PetscCheck(cell_type_f == cell_type, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Only mono-topology face sets are supported currently. Face %" PetscInt_FMT " is %s, which is different than the previous type %s", face, DMPolytopeTypes[cell_type_f], DMPolytopeTypes[cell_type]); 2268 2269 // Get connectivity of the face 2270 PetscCall(DMPlexGetClosureIndices(cdm, cdm->localSection, cdm->localSection, face, PETSC_FALSE, &closure_dof_f, &closure_indices, NULL, NULL)); 2271 elem_size = closure_dof_f / coord_dim; 2272 if (!conn) { 2273 PetscCall(PetscMalloc1(f_owned * elem_size, &conn)); 2274 closure_dof = closure_dof_f; 2275 } 2276 PetscCheck(closure_dof_f == closure_dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Closure of face %" PetscInt_FMT " has %" PetscInt_FMT " dofs instead of the previously written %" PetscInt_FMT " dofs. Only mono-topology face sets are supported currently", face, closure_dof_f, closure_dof); 2277 PetscCall(DMPlexCGNSGetPermutation_Internal(cell_type, elem_size, &element_type, &perm)); 2278 for (PetscInt j = 0; j < elem_size; j++) conn[c++] = node_l2g[closure_indices[perm[j] * coord_dim] / coord_dim] + 1; 2279 PetscCall(DMPlexRestoreClosureIndices(cdm, cdm->localSection, cdm->localSection, face, PETSC_FALSE, &closure_dof_f, &closure_indices, NULL, NULL)); 2280 } 2281 PetscCall(ISRestoreIndices(fsFacesAll, &faces)); 2282 } 2283 2284 { // Write connectivity for face sets 2285 { // Get global element type 2286 PetscInt local_element_type, global_element_type; 2287 2288 local_element_type = f_owned > 0 ? (PetscInt)element_type : -1; 2289 PetscCallMPI(MPIU_Allreduce(&local_element_type, &global_element_type, 1, MPIU_INT, MPI_MAX, comm)); 2290 if (local_element_type != -1) 2291 PetscCheck(local_element_type == global_element_type, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ranks with different element types not supported. Local element type is %s, but global is %s", cg_ElementTypeName(local_element_type), cg_ElementTypeName(global_element_type)); 2292 element_type = (CGNS_ENUMT(ElementType_t))global_element_type; 2293 } 2294 PetscCallMPI(MPIU_Allreduce(&f_owned, &f_global, 1, MPIU_CGSIZE, MPI_SUM, comm)); 2295 f_start = 0; 2296 PetscCallMPI(MPI_Exscan(&f_owned, &f_start, 1, MPIU_CGSIZE, MPI_SUM, comm)); 2297 f_start += elem_offset; 2298 PetscCallCGNSWrite(cgp_section_write(cgv->file_num, base, zone, "Faces", element_type, elem_offset + 1, elem_offset + f_global, 0, §ion), dm, viewer); 2299 PetscCallCGNSWriteData(cgp_elements_write_data(cgv->file_num, base, zone, section, f_start + 1, f_start + f_owned, conn), dm, viewer); 2300 2301 PetscCall(PetscFree(conn)); 2302 PetscCall(PetscFree(parents)); 2303 } 2304 2305 const PetscInt *fsValuesGlobal = NULL; 2306 PetscCall(ISGetIndices(fsValuesGlobalIS, &fsValuesGlobal)); 2307 for (PetscInt fs = 0; fs < num_fs_global; ++fs) { 2308 int BC; 2309 const PetscInt fsID = fsValuesGlobal[fs]; 2310 PetscInt *fs_pnts = NULL; 2311 char bc_name[33]; 2312 cgsize_t fs_start, fs_owned, fs_global; 2313 cgsize_t *fs_pnts_cg; 2314 2315 PetscCall(DMLabelGetStratumIS(fsLabel, fsID, &stratumIS)); 2316 if (stratumIS) { // Get list of only face points 2317 PetscSegBuffer fs_pntsSB; 2318 PetscCount fs_owned_count; 2319 PetscInt nstratumPnts; 2320 const PetscInt *stratumPnts; 2321 2322 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 16, &fs_pntsSB)); 2323 PetscCall(ISGetIndices(stratumIS, &stratumPnts)); 2324 PetscCall(ISGetSize(stratumIS, &nstratumPnts)); 2325 for (PetscInt i = 0; i < nstratumPnts; i++) { 2326 PetscInt *fs_pnts_buffer, stratumPnt = stratumPnts[i]; 2327 if (stratumPnt < fStart || stratumPnt >= fEnd) continue; // Skip non-face points 2328 PetscCall(PetscSegBufferGetInts(fs_pntsSB, 1, &fs_pnts_buffer)); 2329 *fs_pnts_buffer = stratumPnt; 2330 } 2331 PetscCall(PetscSegBufferGetSize(fs_pntsSB, &fs_owned_count)); 2332 fs_owned = fs_owned_count; 2333 PetscCall(PetscSegBufferExtractAlloc(fs_pntsSB, &fs_pnts)); 2334 2335 PetscCall(PetscSegBufferDestroy(&fs_pntsSB)); 2336 PetscCall(ISRestoreIndices(stratumIS, &stratumPnts)); 2337 PetscCall(ISDestroy(&stratumIS)); 2338 } else fs_owned = 0; 2339 2340 PetscCallMPI(MPIU_Allreduce(&fs_owned, &fs_global, 1, MPIU_CGSIZE, MPI_SUM, comm)); 2341 fs_start = 0; 2342 PetscCallMPI(MPI_Exscan(&fs_owned, &fs_start, 1, MPIU_CGSIZE, MPI_SUM, comm)); 2343 PetscCheck(fs_start + fs_owned <= fs_global, PETSC_COMM_SELF, PETSC_ERR_PLIB, "End range of point set (%" PRIdCGSIZE ") greater than global point set size (%" PRIdCGSIZE ")", fs_start + fs_owned, fs_global); 2344 2345 PetscCall(PetscSNPrintf(bc_name, sizeof bc_name, "FaceSet%" PetscInt_FMT, fsID)); 2346 PetscCallCGNSWrite(cg_boco_write(cgv->file_num, base, zone, bc_name, CGNS_ENUMV(BCTypeNull), CGNS_ENUMV(PointList), fs_global, NULL, &BC), dm, viewer); 2347 2348 PetscCall(PetscMalloc1(fs_owned, &fs_pnts_cg)); 2349 for (PetscInt i = 0; i < fs_owned; i++) { 2350 PetscInt is_idx; 2351 2352 PetscCall(ISLocate(fsFacesAll, fs_pnts[i], &is_idx)); 2353 PetscCheck(is_idx >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %" PetscInt_FMT " in list of all local face points", fs_pnts[i]); 2354 fs_pnts_cg[i] = is_idx + f_start + 1; 2355 } 2356 2357 const char *labels[] = {"Zone_t", "ZoneBC_t", "BC_t", "PointList"}; 2358 PetscCallCGNSWrite(cg_golist(cgv->file_num, base, 4, (char **)labels, (int[]){zone, 1, BC, 0}), dm, 0); 2359 PetscCallCGNSWriteData(cgp_ptlist_write_data(cgv->file_num, fs_start + 1, fs_start + fs_owned, fs_pnts_cg), dm, viewer); 2360 2361 CGNS_ENUMT(GridLocation_t) grid_loc = CGNS_ENUMV(GridLocationNull); 2362 if (topo_dim == 3) grid_loc = CGNS_ENUMV(FaceCenter); 2363 else if (topo_dim == 2) grid_loc = CGNS_ENUMV(EdgeCenter); 2364 else if (topo_dim == 1) grid_loc = CGNS_ENUMV(Vertex); 2365 PetscCallCGNSWriteData(cg_boco_gridlocation_write(cgv->file_num, base, zone, BC, grid_loc), dm, viewer); 2366 2367 PetscCall(PetscFree(fs_pnts_cg)); 2368 PetscCall(PetscFree(fs_pnts)); 2369 } 2370 PetscCall(ISDestroy(&fsFacesAll)); 2371 PetscCall(ISRestoreIndices(fsValuesGlobalIS, &fsValuesGlobal)); 2372 elem_offset += f_global; 2373 } 2374 PetscCall(ISDestroy(&fsValuesGlobalIS)); 2375 2376 PetscCall(DMDestroy(&colloc_dm)); 2377 PetscFunctionReturn(PETSC_SUCCESS); 2378 } 2379 2380 PetscErrorCode VecView_Plex_Local_CGNS(Vec V, PetscViewer viewer) 2381 { 2382 PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data; 2383 DM dm; 2384 PetscSection section; 2385 PetscInt time_step, num_fields, pStart, pEnd, fvGhostStart; 2386 PetscReal time, *time_slot; 2387 size_t *step_slot; 2388 const PetscScalar *v; 2389 char solution_name[PETSC_MAX_PATH_LEN]; 2390 int sol; 2391 2392 PetscFunctionBegin; 2393 PetscCall(VecGetDM(V, &dm)); 2394 PetscCall(DMGetLocalSection(dm, §ion)); 2395 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 2396 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &fvGhostStart, NULL)); 2397 if (fvGhostStart >= 0) pEnd = fvGhostStart; 2398 2399 if (!cgv->node_l2g) PetscCall(DMView(dm, viewer)); 2400 if (!cgv->grid_loc) { // Determine if writing to cell-centers or to nodes 2401 PetscInt cStart, cEnd; 2402 PetscInt local_grid_loc, global_grid_loc; 2403 2404 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2405 if (fvGhostStart >= 0) cEnd = fvGhostStart; 2406 if (cgv->num_local_nodes == 0) local_grid_loc = -1; 2407 else if (cStart == pStart && cEnd == pEnd) local_grid_loc = CGNS_ENUMV(CellCenter); 2408 else local_grid_loc = CGNS_ENUMV(Vertex); 2409 2410 PetscCallMPI(MPIU_Allreduce(&local_grid_loc, &global_grid_loc, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)viewer))); 2411 if (local_grid_loc != -1) 2412 PetscCheck(local_grid_loc == global_grid_loc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ranks with different grid locations not supported. Local has %" PetscInt_FMT ", allreduce returned %" PetscInt_FMT, local_grid_loc, global_grid_loc); 2413 PetscCheck((global_grid_loc == CGNS_ENUMV(CellCenter)) || (global_grid_loc == CGNS_ENUMV(Vertex)), PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Grid location should only be CellCenter (%d) or Vertex(%d), but have %" PetscInt_FMT, CGNS_ENUMV(CellCenter), CGNS_ENUMV(Vertex), global_grid_loc); 2414 cgv->grid_loc = (CGNS_ENUMT(GridLocation_t))global_grid_loc; 2415 } 2416 if (!cgv->nodal_field) { 2417 switch (cgv->grid_loc) { 2418 case CGNS_ENUMV(Vertex): { 2419 PetscCall(PetscMalloc1(cgv->nEnd - cgv->nStart, &cgv->nodal_field)); 2420 } break; 2421 case CGNS_ENUMV(CellCenter): { 2422 PetscCall(PetscMalloc1(cgv->eEnd - cgv->eStart, &cgv->nodal_field)); 2423 } break; 2424 default: 2425 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only write for Vertex and CellCenter grid locations"); 2426 } 2427 } 2428 if (!cgv->output_times) PetscCall(PetscSegBufferCreate(sizeof(PetscReal), 20, &cgv->output_times)); 2429 if (!cgv->output_steps) PetscCall(PetscSegBufferCreate(sizeof(size_t), 20, &cgv->output_steps)); 2430 2431 PetscCall(DMGetOutputSequenceNumber(dm, &time_step, &time)); 2432 if (time_step < 0) { 2433 time_step = 0; 2434 time = 0.; 2435 } 2436 // Avoid "Duplicate child name found" error by not writing an already-written solution. 2437 // This usually occurs when a solution is written and then diverges on the very next timestep. 2438 if (time_step == cgv->previous_output_step) PetscFunctionReturn(PETSC_SUCCESS); 2439 2440 PetscCall(PetscSegBufferGet(cgv->output_times, 1, &time_slot)); 2441 *time_slot = time; 2442 PetscCall(PetscSegBufferGet(cgv->output_steps, 1, &step_slot)); 2443 *step_slot = cgv->previous_output_step = time_step; 2444 PetscCall(PetscSNPrintf(solution_name, sizeof solution_name, "FlowSolution%" PetscInt_FMT, time_step)); 2445 PetscCallCGNSWrite(cg_sol_write(cgv->file_num, cgv->base, cgv->zone, solution_name, cgv->grid_loc, &sol), V, viewer); 2446 PetscCall(VecGetArrayRead(V, &v)); 2447 PetscCall(PetscSectionGetNumFields(section, &num_fields)); 2448 for (PetscInt field = 0; field < num_fields; field++) { 2449 PetscInt ncomp; 2450 const char *field_name; 2451 PetscCall(PetscSectionGetFieldName(section, field, &field_name)); 2452 PetscCall(PetscSectionGetFieldComponents(section, field, &ncomp)); 2453 for (PetscInt comp = 0; comp < ncomp; comp++) { 2454 int cgfield; 2455 const char *comp_name; 2456 char cgns_field_name[32]; // CGNS max field name is 32 2457 CGNS_ENUMT(DataType_t) datatype; 2458 PetscCall(PetscSectionGetComponentName(section, field, comp, &comp_name)); 2459 if (ncomp == 1 && comp_name[0] == '0' && comp_name[1] == '\0' && field_name[0] != '\0') PetscCall(PetscStrncpy(cgns_field_name, field_name, sizeof cgns_field_name)); 2460 else if (field_name[0] == '\0') PetscCall(PetscStrncpy(cgns_field_name, comp_name, sizeof cgns_field_name)); 2461 else PetscCall(PetscSNPrintf(cgns_field_name, sizeof cgns_field_name, "%s.%s", field_name, comp_name)); 2462 PetscCall(PetscCGNSDataType(PETSC_SCALAR, &datatype)); 2463 PetscCallCGNSWrite(cgp_field_write(cgv->file_num, cgv->base, cgv->zone, sol, datatype, cgns_field_name, &cgfield), V, viewer); 2464 for (PetscInt p = pStart, n = 0; p < pEnd; p++) { 2465 PetscInt off, dof; 2466 PetscCall(PetscSectionGetFieldDof(section, p, field, &dof)); 2467 if (dof == 0) continue; 2468 PetscCall(PetscSectionGetFieldOffset(section, p, field, &off)); 2469 for (PetscInt c = comp; c < dof; c += ncomp, n++) { 2470 switch (cgv->grid_loc) { 2471 case CGNS_ENUMV(Vertex): { 2472 PetscInt gn = cgv->node_l2g[n]; 2473 if (gn < cgv->nStart || cgv->nEnd <= gn) continue; 2474 cgv->nodal_field[gn - cgv->nStart] = v[off + c]; 2475 } break; 2476 case CGNS_ENUMV(CellCenter): { 2477 cgv->nodal_field[n] = v[off + c]; 2478 } break; 2479 default: 2480 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only pack for Vertex and CellCenter grid locations"); 2481 } 2482 } 2483 } 2484 // CGNS nodes use 1-based indexing 2485 cgsize_t start, end; 2486 switch (cgv->grid_loc) { 2487 case CGNS_ENUMV(Vertex): { 2488 start = cgv->nStart + 1; 2489 end = cgv->nEnd; 2490 } break; 2491 case CGNS_ENUMV(CellCenter): { 2492 start = cgv->eStart + 1; 2493 end = cgv->eEnd; 2494 } break; 2495 default: 2496 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only write for Vertex and CellCenter grid locations"); 2497 } 2498 PetscCallCGNSWriteData(cgp_field_write_data(cgv->file_num, cgv->base, cgv->zone, sol, cgfield, &start, &end, cgv->nodal_field), V, viewer); 2499 } 2500 } 2501 PetscCall(VecRestoreArrayRead(V, &v)); 2502 PetscCall(PetscViewerCGNSCheckBatch_Internal(viewer)); 2503 PetscFunctionReturn(PETSC_SUCCESS); 2504 } 2505 2506 PetscErrorCode VecLoad_Plex_CGNS_Internal(Vec V, PetscViewer viewer) 2507 { 2508 MPI_Comm comm; 2509 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 2510 PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data; 2511 int cgid = cgv->file_num; 2512 PetscBool use_parallel_viewer = PETSC_FALSE; 2513 int z = 1; // Only support one zone 2514 int B = 1; // Only support one base 2515 int numComp; 2516 PetscInt V_numComps, mystartv, myendv, myownedv; 2517 2518 PetscFunctionBegin; 2519 PetscCall(PetscObjectGetComm((PetscObject)V, &comm)); 2520 2521 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_cgns_parallel", &use_parallel_viewer, NULL)); 2522 PetscCheck(use_parallel_viewer, comm, PETSC_ERR_USER_INPUT, "Cannot use VecLoad with CGNS file in serial reader; use -dm_plex_cgns_parallel to enable parallel reader"); 2523 2524 { // Get CGNS node ownership information 2525 int nbases, nzones; 2526 PetscInt NVertices; 2527 PetscLayout vtx_map; 2528 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 2529 2530 PetscCallCGNSRead(cg_nbases(cgid, &nbases), V, viewer); 2531 PetscCheck(nbases <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single base, not %d", nbases); 2532 PetscCallCGNSRead(cg_nzones(cgid, B, &nzones), V, viewer); 2533 PetscCheck(nzones == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "limited to one zone %d", (int)nzones); 2534 2535 PetscCallCGNSRead(cg_zone_read(cgid, B, z, buffer, sizes), V, viewer); 2536 NVertices = sizes[0]; 2537 2538 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NVertices, 1, &vtx_map)); 2539 PetscCall(PetscLayoutGetRange(vtx_map, &mystartv, &myendv)); 2540 PetscCall(PetscLayoutGetLocalSize(vtx_map, &myownedv)); 2541 PetscCall(PetscLayoutDestroy(&vtx_map)); 2542 } 2543 2544 { // -- Read data from file into Vec 2545 PetscScalar *fields = NULL; 2546 PetscSF sfNatural; 2547 2548 { // Check compatibility between sfNatural and the data source and sink 2549 DM dm; 2550 PetscInt nleaves, nroots, V_local_size; 2551 2552 PetscCall(VecGetDM(V, &dm)); 2553 PetscCall(DMGetNaturalSF(dm, &sfNatural)); 2554 PetscCheck(sfNatural, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM of Vec must have sfNatural"); 2555 PetscCall(PetscSFGetGraph(sfNatural, &nroots, &nleaves, NULL, NULL)); 2556 PetscCall(VecGetLocalSize(V, &V_local_size)); 2557 PetscCheck(nleaves == myownedv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Number of locally owned vertices (% " PetscInt_FMT ") must match number of leaves in sfNatural (% " PetscInt_FMT ")", myownedv, nleaves); 2558 if (nroots == 0) { 2559 PetscCheck(V_local_size == nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Local Vec size (% " PetscInt_FMT ") must be zero if number of roots in sfNatural is zero", V_local_size); 2560 V_numComps = 0; 2561 } else { 2562 PetscCheck(V_local_size % nroots == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Local Vec size (% " PetscInt_FMT ") not evenly divisible by number of roots in sfNatural (% " PetscInt_FMT ")", V_local_size, nroots); 2563 V_numComps = V_local_size / nroots; 2564 } 2565 } 2566 2567 { // Read data into component-major ordering 2568 int isol, numSols; 2569 CGNS_ENUMT(DataType_t) datatype; 2570 double *fields_CGNS; 2571 2572 PetscCallCGNSRead(cg_nsols(cgid, B, z, &numSols), V, viewer); 2573 PetscCall(PetscViewerCGNSGetSolutionFileIndex_Internal(viewer, &isol)); 2574 PetscCallCGNSRead(cg_nfields(cgid, B, z, isol, &numComp), V, viewer); 2575 PetscCheck(V_numComps == numComp || V_numComps == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec sized for % " PetscInt_FMT " components per node, but file has %d components per node", V_numComps, numComp); 2576 2577 cgsize_t range_min[3] = {mystartv + 1, 1, 1}; 2578 cgsize_t range_max[3] = {myendv, 1, 1}; 2579 PetscCall(PetscMalloc1(myownedv * numComp, &fields_CGNS)); 2580 PetscCall(PetscMalloc1(myownedv * numComp, &fields)); 2581 for (int d = 0; d < numComp; ++d) { 2582 PetscCallCGNSRead(cg_field_info(cgid, B, z, isol, (d + 1), &datatype, buffer), V, viewer); 2583 PetscCheck(datatype == CGNS_ENUMV(RealDouble), PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Field %s in file is not of type double", buffer); 2584 PetscCallCGNSReadData(cgp_field_read_data(cgid, B, z, isol, (d + 1), range_min, range_max, &fields_CGNS[d * myownedv]), V, viewer); 2585 } 2586 for (int d = 0; d < numComp; ++d) { 2587 for (PetscInt v = 0; v < myownedv; ++v) fields[v * numComp + d] = fields_CGNS[d * myownedv + v]; 2588 } 2589 PetscCall(PetscFree(fields_CGNS)); 2590 } 2591 2592 { // Reduce fields into Vec array 2593 PetscScalar *V_array; 2594 MPI_Datatype fieldtype; 2595 2596 PetscCall(VecGetArrayWrite(V, &V_array)); 2597 PetscCallMPI(MPI_Type_contiguous(numComp, MPIU_SCALAR, &fieldtype)); 2598 PetscCallMPI(MPI_Type_commit(&fieldtype)); 2599 PetscCall(PetscSFReduceBegin(sfNatural, fieldtype, fields, V_array, MPI_REPLACE)); 2600 PetscCall(PetscSFReduceEnd(sfNatural, fieldtype, fields, V_array, MPI_REPLACE)); 2601 PetscCallMPI(MPI_Type_free(&fieldtype)); 2602 PetscCall(VecRestoreArrayWrite(V, &V_array)); 2603 } 2604 PetscCall(PetscFree(fields)); 2605 } 2606 PetscFunctionReturn(PETSC_SUCCESS); 2607 } 2608