1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/viewercgnsimpl.h> 3 4 #include <pcgnslib.h> 5 #include <cgns_io.h> 6 7 #if !defined(CGNS_ENUMT) 8 #define CGNS_ENUMT(a) a 9 #endif 10 #if !defined(CGNS_ENUMV) 11 #define CGNS_ENUMV(a) a 12 #endif 13 // Permute plex closure ordering to CGNS 14 static PetscErrorCode DMPlexCGNSGetPermutation_Internal(DMPolytopeType cell_type, PetscInt closure_size, CGNS_ENUMT(ElementType_t) * element_type, const int **perm) 15 { 16 CGNS_ENUMT(ElementType_t) element_type_tmp; 17 18 // https://cgns.github.io/CGNS_docs_current/sids/conv.html#unstructgrid 19 static const int bar_2[2] = {0, 1}; 20 static const int bar_3[3] = {1, 2, 0}; 21 static const int bar_4[4] = {2, 3, 0, 1}; 22 static const int bar_5[5] = {3, 4, 0, 1, 2}; 23 static const int tri_3[3] = {0, 1, 2}; 24 static const int tri_6[6] = {3, 4, 5, 0, 1, 2}; 25 static const int tri_10[10] = {7, 8, 9, 1, 2, 3, 4, 5, 6, 0}; 26 static const int quad_4[4] = {0, 1, 2, 3}; 27 static const int quad_9[9] = { 28 5, 6, 7, 8, // vertices 29 1, 2, 3, 4, // edges 30 0, // center 31 }; 32 static const int quad_16[] = { 33 12, 13, 14, 15, // vertices 34 4, 5, 6, 7, 8, 9, 10, 11, // edges 35 0, 1, 3, 2, // centers 36 }; 37 static const int quad_25[] = { 38 21, 22, 23, 24, // vertices 39 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, // edges 40 0, 1, 2, 5, 8, 7, 6, 3, 4, // centers 41 }; 42 static const int tetra_4[4] = {0, 2, 1, 3}; 43 static const int tetra_10[10] = {6, 8, 7, 9, 2, 1, 0, 3, 5, 4}; 44 static const int tetra_20[20] = { 45 16, 18, 17, 19, // vertices 46 9, 8, 7, 6, 5, 4, // bottom edges 47 10, 11, 14, 15, 13, 12, // side edges 48 0, 2, 3, 1, // faces 49 }; 50 static const int hexa_8[8] = {0, 3, 2, 1, 4, 5, 6, 7}; 51 static const int hexa_27[27] = { 52 19, 22, 21, 20, 23, 24, 25, 26, // vertices 53 10, 9, 8, 7, // bottom edges 54 16, 15, 18, 17, // mid edges 55 11, 12, 13, 14, // top edges 56 1, 3, 5, 4, 6, 2, // faces 57 0, // center 58 }; 59 static const int hexa_64[64] = { 60 // 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 61 56, 59, 58, 57, 60, 61, 62, 63, // vertices 62 39, 38, 37, 36, 35, 34, 33, 32, // bottom edges 63 51, 50, 48, 49, 52, 53, 55, 54, // mid edges; Paraview needs edge 21-22 swapped with 23-24 64 40, 41, 42, 43, 44, 45, 46, 47, // top edges 65 8, 10, 11, 9, // z-minus face 66 16, 17, 19, 18, // y-minus face 67 24, 25, 27, 26, // x-plus face 68 20, 21, 23, 22, // y-plus face 69 30, 28, 29, 31, // x-minus face 70 12, 13, 15, 14, // z-plus face 71 0, 1, 3, 2, 4, 5, 7, 6, // center 72 }; 73 74 PetscFunctionBegin; 75 element_type_tmp = CGNS_ENUMV(ElementTypeNull); 76 *perm = NULL; 77 switch (cell_type) { 78 case DM_POLYTOPE_SEGMENT: 79 switch (closure_size) { 80 case 2: 81 element_type_tmp = CGNS_ENUMV(BAR_2); 82 *perm = bar_2; 83 break; 84 case 3: 85 element_type_tmp = CGNS_ENUMV(BAR_3); 86 *perm = bar_3; 87 break; 88 case 4: 89 element_type_tmp = CGNS_ENUMV(BAR_4); 90 *perm = bar_4; 91 break; 92 case 5: 93 element_type_tmp = CGNS_ENUMV(BAR_5); 94 *perm = bar_5; 95 break; 96 default: 97 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 98 } 99 break; 100 case DM_POLYTOPE_TRIANGLE: 101 switch (closure_size) { 102 case 3: 103 element_type_tmp = CGNS_ENUMV(TRI_3); 104 *perm = tri_3; 105 break; 106 case 6: 107 element_type_tmp = CGNS_ENUMV(TRI_6); 108 *perm = tri_6; 109 break; 110 case 10: 111 element_type_tmp = CGNS_ENUMV(TRI_10); 112 *perm = tri_10; 113 break; 114 default: 115 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 116 } 117 break; 118 case DM_POLYTOPE_QUADRILATERAL: 119 switch (closure_size) { 120 case 4: 121 element_type_tmp = CGNS_ENUMV(QUAD_4); 122 *perm = quad_4; 123 break; 124 case 9: 125 element_type_tmp = CGNS_ENUMV(QUAD_9); 126 *perm = quad_9; 127 break; 128 case 16: 129 element_type_tmp = CGNS_ENUMV(QUAD_16); 130 *perm = quad_16; 131 break; 132 case 25: 133 element_type_tmp = CGNS_ENUMV(QUAD_25); 134 *perm = quad_25; 135 break; 136 default: 137 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 138 } 139 break; 140 case DM_POLYTOPE_TETRAHEDRON: 141 switch (closure_size) { 142 case 4: 143 element_type_tmp = CGNS_ENUMV(TETRA_4); 144 *perm = tetra_4; 145 break; 146 case 10: 147 element_type_tmp = CGNS_ENUMV(TETRA_10); 148 *perm = tetra_10; 149 break; 150 case 20: 151 element_type_tmp = CGNS_ENUMV(TETRA_20); 152 *perm = tetra_20; 153 break; 154 default: 155 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 156 } 157 break; 158 case DM_POLYTOPE_HEXAHEDRON: 159 switch (closure_size) { 160 case 8: 161 element_type_tmp = CGNS_ENUMV(HEXA_8); 162 *perm = hexa_8; 163 break; 164 case 27: 165 element_type_tmp = CGNS_ENUMV(HEXA_27); 166 *perm = hexa_27; 167 break; 168 case 64: 169 element_type_tmp = CGNS_ENUMV(HEXA_64); 170 *perm = hexa_64; 171 break; 172 default: 173 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 174 } 175 break; 176 default: 177 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size); 178 } 179 if (element_type) *element_type = element_type_tmp; 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 /* 184 Input Parameters: 185 + cellType - The CGNS-defined element type 186 187 Output Parameters: 188 + dmcelltype - The equivalent DMPolytopeType for the cellType 189 . numCorners - Number of corners of the polytope 190 - dim - The topological dimension of the polytope 191 192 CGNS elements defined in: https://cgns.github.io/CGNS_docs_current/sids/conv.html#unstructgrid 193 */ 194 static inline PetscErrorCode CGNSElementTypeGetTopologyInfo(CGNS_ENUMT(ElementType_t) cellType, DMPolytopeType *dmcelltype, PetscInt *numCorners, PetscInt *dim) 195 { 196 DMPolytopeType _dmcelltype; 197 198 PetscFunctionBegin; 199 switch (cellType) { 200 case CGNS_ENUMV(BAR_2): 201 case CGNS_ENUMV(BAR_3): 202 case CGNS_ENUMV(BAR_4): 203 case CGNS_ENUMV(BAR_5): 204 _dmcelltype = DM_POLYTOPE_SEGMENT; 205 break; 206 case CGNS_ENUMV(TRI_3): 207 case CGNS_ENUMV(TRI_6): 208 case CGNS_ENUMV(TRI_9): 209 case CGNS_ENUMV(TRI_10): 210 case CGNS_ENUMV(TRI_12): 211 case CGNS_ENUMV(TRI_15): 212 _dmcelltype = DM_POLYTOPE_TRIANGLE; 213 break; 214 case CGNS_ENUMV(QUAD_4): 215 case CGNS_ENUMV(QUAD_8): 216 case CGNS_ENUMV(QUAD_9): 217 case CGNS_ENUMV(QUAD_12): 218 case CGNS_ENUMV(QUAD_16): 219 case CGNS_ENUMV(QUAD_P4_16): 220 case CGNS_ENUMV(QUAD_25): 221 _dmcelltype = DM_POLYTOPE_QUADRILATERAL; 222 break; 223 case CGNS_ENUMV(TETRA_4): 224 case CGNS_ENUMV(TETRA_10): 225 case CGNS_ENUMV(TETRA_16): 226 case CGNS_ENUMV(TETRA_20): 227 case CGNS_ENUMV(TETRA_22): 228 case CGNS_ENUMV(TETRA_34): 229 case CGNS_ENUMV(TETRA_35): 230 _dmcelltype = DM_POLYTOPE_TETRAHEDRON; 231 break; 232 case CGNS_ENUMV(PYRA_5): 233 case CGNS_ENUMV(PYRA_13): 234 case CGNS_ENUMV(PYRA_14): 235 case CGNS_ENUMV(PYRA_21): 236 case CGNS_ENUMV(PYRA_29): 237 case CGNS_ENUMV(PYRA_P4_29): 238 case CGNS_ENUMV(PYRA_30): 239 case CGNS_ENUMV(PYRA_50): 240 case CGNS_ENUMV(PYRA_55): 241 _dmcelltype = DM_POLYTOPE_PYRAMID; 242 break; 243 case CGNS_ENUMV(PENTA_6): 244 case CGNS_ENUMV(PENTA_15): 245 case CGNS_ENUMV(PENTA_18): 246 case CGNS_ENUMV(PENTA_24): 247 case CGNS_ENUMV(PENTA_33): 248 case CGNS_ENUMV(PENTA_38): 249 case CGNS_ENUMV(PENTA_40): 250 case CGNS_ENUMV(PENTA_66): 251 case CGNS_ENUMV(PENTA_75): 252 _dmcelltype = DM_POLYTOPE_TRI_PRISM; 253 break; 254 case CGNS_ENUMV(HEXA_8): 255 case CGNS_ENUMV(HEXA_20): 256 case CGNS_ENUMV(HEXA_27): 257 case CGNS_ENUMV(HEXA_32): 258 case CGNS_ENUMV(HEXA_44): 259 case CGNS_ENUMV(HEXA_56): 260 case CGNS_ENUMV(HEXA_64): 261 case CGNS_ENUMV(HEXA_98): 262 case CGNS_ENUMV(HEXA_125): 263 _dmcelltype = DM_POLYTOPE_HEXAHEDRON; 264 break; 265 case CGNS_ENUMV(MIXED): 266 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid CGNS ElementType_t: MIXED"); 267 default: 268 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid CGNS ElementType_t: %d", (int)cellType); 269 } 270 271 if (dmcelltype) *dmcelltype = _dmcelltype; 272 if (numCorners) *numCorners = DMPolytopeTypeGetNumVertices(_dmcelltype); 273 if (dim) *dim = DMPolytopeTypeGetDim(_dmcelltype); 274 PetscFunctionReturn(PETSC_SUCCESS); 275 } 276 277 /* 278 Input Parameters: 279 + cellType - The CGNS-defined cell type 280 281 Output Parameters: 282 + numClosure - Number of nodes that define the function space on the cell 283 - pOrder - The polynomial order of the cell 284 285 CGNS elements defined in: https://cgns.github.io/CGNS_docs_current/sids/conv.html#unstructgrid 286 287 Note: we only support "full" elements, ie. not seredipity elements 288 */ 289 static inline PetscErrorCode CGNSElementTypeGetDiscretizationInfo(CGNS_ENUMT(ElementType_t) cellType, PetscInt *numClosure, PetscInt *pOrder) 290 { 291 PetscInt _numClosure, _pOrder; 292 293 PetscFunctionBegin; 294 switch (cellType) { 295 case CGNS_ENUMV(BAR_2): 296 _numClosure = 2; 297 _pOrder = 1; 298 break; 299 case CGNS_ENUMV(BAR_3): 300 _numClosure = 3; 301 _pOrder = 2; 302 break; 303 case CGNS_ENUMV(BAR_4): 304 _numClosure = 4; 305 _pOrder = 3; 306 break; 307 case CGNS_ENUMV(BAR_5): 308 _numClosure = 5; 309 _pOrder = 4; 310 break; 311 case CGNS_ENUMV(TRI_3): 312 _numClosure = 3; 313 _pOrder = 1; 314 break; 315 case CGNS_ENUMV(TRI_6): 316 _numClosure = 6; 317 _pOrder = 2; 318 break; 319 case CGNS_ENUMV(TRI_10): 320 _numClosure = 10; 321 _pOrder = 3; 322 break; 323 case CGNS_ENUMV(TRI_15): 324 _numClosure = 15; 325 _pOrder = 4; 326 break; 327 case CGNS_ENUMV(QUAD_4): 328 _numClosure = 4; 329 _pOrder = 1; 330 break; 331 case CGNS_ENUMV(QUAD_9): 332 _numClosure = 9; 333 _pOrder = 2; 334 break; 335 case CGNS_ENUMV(QUAD_16): 336 _numClosure = 16; 337 _pOrder = 3; 338 break; 339 case CGNS_ENUMV(QUAD_25): 340 _numClosure = 25; 341 _pOrder = 4; 342 break; 343 case CGNS_ENUMV(TETRA_4): 344 _numClosure = 4; 345 _pOrder = 1; 346 break; 347 case CGNS_ENUMV(TETRA_10): 348 _numClosure = 10; 349 _pOrder = 2; 350 break; 351 case CGNS_ENUMV(TETRA_20): 352 _numClosure = 20; 353 _pOrder = 3; 354 break; 355 case CGNS_ENUMV(TETRA_35): 356 _numClosure = 35; 357 _pOrder = 4; 358 break; 359 case CGNS_ENUMV(PYRA_5): 360 _numClosure = 5; 361 _pOrder = 1; 362 break; 363 case CGNS_ENUMV(PYRA_14): 364 _numClosure = 14; 365 _pOrder = 2; 366 break; 367 case CGNS_ENUMV(PYRA_30): 368 _numClosure = 30; 369 _pOrder = 3; 370 break; 371 case CGNS_ENUMV(PYRA_55): 372 _numClosure = 55; 373 _pOrder = 4; 374 break; 375 case CGNS_ENUMV(PENTA_6): 376 _numClosure = 6; 377 _pOrder = 1; 378 break; 379 case CGNS_ENUMV(PENTA_18): 380 _numClosure = 18; 381 _pOrder = 2; 382 break; 383 case CGNS_ENUMV(PENTA_40): 384 _numClosure = 40; 385 _pOrder = 3; 386 break; 387 case CGNS_ENUMV(PENTA_75): 388 _numClosure = 75; 389 _pOrder = 4; 390 break; 391 case CGNS_ENUMV(HEXA_8): 392 _numClosure = 8; 393 _pOrder = 1; 394 break; 395 case CGNS_ENUMV(HEXA_27): 396 _numClosure = 27; 397 _pOrder = 2; 398 break; 399 case CGNS_ENUMV(HEXA_64): 400 _numClosure = 64; 401 _pOrder = 3; 402 break; 403 case CGNS_ENUMV(HEXA_125): 404 _numClosure = 125; 405 _pOrder = 4; 406 break; 407 case CGNS_ENUMV(MIXED): 408 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid CGNS ElementType_t: MIXED"); 409 default: 410 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported or Invalid cell type %d", (int)cellType); 411 } 412 if (numClosure) *numClosure = _numClosure; 413 if (pOrder) *pOrder = _pOrder; 414 PetscFunctionReturn(PETSC_SUCCESS); 415 } 416 417 static PetscErrorCode PetscCGNSDataType(PetscDataType pd, CGNS_ENUMT(DataType_t) * cd) 418 { 419 PetscFunctionBegin; 420 switch (pd) { 421 case PETSC_FLOAT: 422 *cd = CGNS_ENUMV(RealSingle); 423 break; 424 case PETSC_DOUBLE: 425 *cd = CGNS_ENUMV(RealDouble); 426 break; 427 case PETSC_COMPLEX: 428 *cd = CGNS_ENUMV(ComplexDouble); 429 break; 430 default: 431 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Data type %s", PetscDataTypes[pd]); 432 } 433 PetscFunctionReturn(PETSC_SUCCESS); 434 } 435 436 PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 437 { 438 int cgid = -1; 439 PetscBool use_parallel_viewer = PETSC_FALSE; 440 441 PetscFunctionBegin; 442 PetscAssertPointer(filename, 2); 443 PetscCall(PetscViewerCGNSRegisterLogEvents_Internal()); 444 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_cgns_parallel", &use_parallel_viewer, NULL)); 445 446 if (use_parallel_viewer) { 447 PetscCallCGNS(cgp_mpi_comm(comm)); 448 PetscCallCGNSOpen(cgp_open(filename, CG_MODE_READ, &cgid), 0, 0); 449 PetscCheck(cgid > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "cgp_open(\"%s\",...) did not return a valid file ID", filename); 450 PetscCall(DMPlexCreateCGNS(comm, cgid, interpolate, dm)); 451 PetscCallCGNSClose(cgp_close(cgid), 0, 0); 452 } else { 453 PetscCallCGNSOpen(cg_open(filename, CG_MODE_READ, &cgid), 0, 0); 454 PetscCheck(cgid > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename); 455 PetscCall(DMPlexCreateCGNS(comm, cgid, interpolate, dm)); 456 PetscCallCGNSClose(cg_close(cgid), 0, 0); 457 } 458 PetscFunctionReturn(PETSC_SUCCESS); 459 } 460 461 PetscErrorCode DMPlexCreateCGNS_Internal_Serial(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm) 462 { 463 PetscMPIInt num_proc, rank; 464 DM cdm; 465 DMLabel label; 466 PetscSection coordSection; 467 Vec coordinates; 468 PetscScalar *coords; 469 PetscInt *cellStart, *vertStart, v; 470 PetscInt labelIdRange[2], labelId; 471 /* Read from file */ 472 char basename[CGIO_MAX_NAME_LENGTH + 1]; 473 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 474 int dim = 0, physDim = 0, coordDim = 0, numVertices = 0, numCells = 0; 475 int nzones = 0; 476 const int B = 1; // Only support single base 477 478 PetscFunctionBegin; 479 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 480 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 481 PetscCall(DMCreate(comm, dm)); 482 PetscCall(DMSetType(*dm, DMPLEX)); 483 484 /* Open CGNS II file and read basic information on rank 0, then broadcast to all processors */ 485 if (rank == 0) { 486 int nbases, z; 487 488 PetscCallCGNSRead(cg_nbases(cgid, &nbases), *dm, 0); 489 PetscCheck(nbases <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single base, not %d", nbases); 490 PetscCallCGNSRead(cg_base_read(cgid, B, basename, &dim, &physDim), *dm, 0); 491 PetscCallCGNSRead(cg_nzones(cgid, B, &nzones), *dm, 0); 492 PetscCall(PetscCalloc2(nzones + 1, &cellStart, nzones + 1, &vertStart)); 493 for (z = 1; z <= nzones; ++z) { 494 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 495 496 PetscCallCGNSRead(cg_zone_read(cgid, B, z, buffer, sizes), *dm, 0); 497 numVertices += sizes[0]; 498 numCells += sizes[1]; 499 cellStart[z] += sizes[1] + cellStart[z - 1]; 500 vertStart[z] += sizes[0] + vertStart[z - 1]; 501 } 502 for (z = 1; z <= nzones; ++z) vertStart[z] += numCells; 503 coordDim = dim; 504 } 505 PetscCallMPI(MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH + 1, MPI_CHAR, 0, comm)); 506 PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm)); 507 PetscCallMPI(MPI_Bcast(&coordDim, 1, MPI_INT, 0, comm)); 508 PetscCallMPI(MPI_Bcast(&nzones, 1, MPI_INT, 0, comm)); 509 510 PetscCall(PetscObjectSetName((PetscObject)*dm, basename)); 511 PetscCall(DMSetDimension(*dm, dim)); 512 PetscCall(DMCreateLabel(*dm, "celltype")); 513 PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices)); 514 515 /* Read zone information */ 516 if (rank == 0) { 517 int z, c, c_loc; 518 519 /* Read the cell set connectivity table and build mesh topology 520 CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */ 521 /* First set sizes */ 522 for (z = 1, c = 0; z <= nzones; ++z) { 523 CGNS_ENUMT(ZoneType_t) zonetype; 524 int nsections; 525 CGNS_ENUMT(ElementType_t) cellType; 526 cgsize_t start, end; 527 int nbndry, parentFlag; 528 PetscInt numCorners, pOrder; 529 DMPolytopeType ctype; 530 const int S = 1; // Only support single section 531 532 PetscCallCGNSRead(cg_zone_type(cgid, B, z, &zonetype), *dm, 0); 533 PetscCheck(zonetype != CGNS_ENUMV(Structured), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can only handle Unstructured zones for CGNS"); 534 PetscCallCGNSRead(cg_nsections(cgid, B, z, &nsections), *dm, 0); 535 PetscCheck(nsections <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single section, not %d", nsections); 536 PetscCallCGNSRead(cg_section_read(cgid, B, z, S, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 537 if (cellType == CGNS_ENUMV(MIXED)) { 538 cgsize_t elementDataSize, *elements; 539 PetscInt off; 540 541 PetscCallCGNSRead(cg_ElementDataSize(cgid, B, z, S, &elementDataSize), *dm, 0); 542 PetscCall(PetscMalloc1(elementDataSize, &elements)); 543 PetscCallCGNSReadData(cg_poly_elements_read(cgid, B, z, S, elements, NULL, NULL), *dm, 0); 544 for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) { 545 PetscCall(CGNSElementTypeGetTopologyInfo((CGNS_ENUMT(ElementType_t))elements[off], &ctype, &numCorners, NULL)); 546 PetscCall(CGNSElementTypeGetDiscretizationInfo((CGNS_ENUMT(ElementType_t))elements[off], NULL, &pOrder)); 547 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 548 PetscCall(DMPlexSetConeSize(*dm, c, numCorners)); 549 PetscCall(DMPlexSetCellType(*dm, c, ctype)); 550 off += numCorners + 1; 551 } 552 PetscCall(PetscFree(elements)); 553 } else { 554 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, &ctype, &numCorners, NULL)); 555 PetscCall(CGNSElementTypeGetDiscretizationInfo(cellType, NULL, &pOrder)); 556 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 557 for (c_loc = start; c_loc <= end; ++c_loc, ++c) { 558 PetscCall(DMPlexSetConeSize(*dm, c, numCorners)); 559 PetscCall(DMPlexSetCellType(*dm, c, ctype)); 560 } 561 } 562 } 563 for (v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 564 } 565 566 PetscCall(DMSetUp(*dm)); 567 568 PetscCall(DMCreateLabel(*dm, "zone")); 569 if (rank == 0) { 570 int z, c, c_loc, v_loc; 571 572 PetscCall(DMGetLabel(*dm, "zone", &label)); 573 for (z = 1, c = 0; z <= nzones; ++z) { 574 CGNS_ENUMT(ElementType_t) cellType; 575 cgsize_t elementDataSize, *elements, start, end; 576 int nbndry, parentFlag; 577 PetscInt *cone, numc, numCorners, maxCorners = 27, pOrder; 578 const int S = 1; // Only support single section 579 580 PetscCallCGNSRead(cg_section_read(cgid, B, z, S, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 581 numc = end - start; 582 PetscCallCGNSRead(cg_ElementDataSize(cgid, B, z, S, &elementDataSize), *dm, 0); 583 PetscCall(PetscMalloc2(elementDataSize, &elements, maxCorners, &cone)); 584 PetscCallCGNSReadData(cg_poly_elements_read(cgid, B, z, S, elements, NULL, NULL), *dm, 0); 585 if (cellType == CGNS_ENUMV(MIXED)) { 586 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 587 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 588 PetscCall(CGNSElementTypeGetTopologyInfo((CGNS_ENUMT(ElementType_t))elements[v], NULL, &numCorners, NULL)); 589 PetscCall(CGNSElementTypeGetDiscretizationInfo((CGNS_ENUMT(ElementType_t))elements[v], NULL, &pOrder)); 590 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 591 ++v; 592 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) cone[v_loc] = elements[v] + numCells - 1; 593 PetscCall(DMPlexReorderCell(*dm, c, cone)); 594 PetscCall(DMPlexSetCone(*dm, c, cone)); 595 PetscCall(DMLabelSetValue(label, c, z)); 596 } 597 } else { 598 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, NULL, &numCorners, NULL)); 599 PetscCall(CGNSElementTypeGetDiscretizationInfo(cellType, NULL, &pOrder)); 600 PetscCheck(pOrder == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Serial CGNS reader only supports first order elements, not %" PetscInt_FMT " order", pOrder); 601 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 602 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 603 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) cone[v_loc] = elements[v] + numCells - 1; 604 PetscCall(DMPlexReorderCell(*dm, c, cone)); 605 PetscCall(DMPlexSetCone(*dm, c, cone)); 606 PetscCall(DMLabelSetValue(label, c, z)); 607 } 608 } 609 PetscCall(PetscFree2(elements, cone)); 610 } 611 } 612 613 PetscCall(DMPlexSymmetrize(*dm)); 614 PetscCall(DMPlexStratify(*dm)); 615 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(*dm)); 616 617 /* Read coordinates */ 618 PetscCall(DMSetCoordinateDim(*dm, coordDim)); 619 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 620 PetscCall(DMGetLocalSection(cdm, &coordSection)); 621 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 622 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, coordDim)); 623 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 624 for (v = numCells; v < numCells + numVertices; ++v) { 625 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 626 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, coordDim)); 627 } 628 PetscCall(PetscSectionSetUp(coordSection)); 629 630 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 631 PetscCall(VecGetArray(coordinates, &coords)); 632 if (rank == 0) { 633 PetscInt off = 0; 634 float *x[3]; 635 int z, d; 636 637 PetscCall(PetscMalloc3(numVertices, &x[0], numVertices, &x[1], numVertices, &x[2])); 638 for (z = 1; z <= nzones; ++z) { 639 CGNS_ENUMT(DataType_t) datatype; 640 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 641 cgsize_t range_min[3] = {1, 1, 1}; 642 cgsize_t range_max[3] = {1, 1, 1}; 643 int ngrids, ncoords; 644 645 PetscCallCGNSRead(cg_zone_read(cgid, B, z, buffer, sizes), *dm, 0); 646 range_max[0] = sizes[0]; 647 PetscCallCGNSRead(cg_ngrids(cgid, B, z, &ngrids), *dm, 0); 648 PetscCheck(ngrids <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single grid, not %d", ngrids); 649 PetscCallCGNSRead(cg_ncoords(cgid, B, z, &ncoords), *dm, 0); 650 PetscCheck(ncoords == coordDim, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a coordinate array for each dimension, not %d", ncoords); 651 for (d = 0; d < coordDim; ++d) { 652 PetscCallCGNSRead(cg_coord_info(cgid, B, z, 1 + d, &datatype, buffer), *dm, 0); 653 PetscCallCGNSReadData(cg_coord_read(cgid, B, z, buffer, CGNS_ENUMV(RealSingle), range_min, range_max, x[d]), *dm, 0); 654 } 655 if (coordDim >= 1) { 656 for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 0] = x[0][v]; 657 } 658 if (coordDim >= 2) { 659 for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 1] = x[1][v]; 660 } 661 if (coordDim >= 3) { 662 for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 2] = x[2][v]; 663 } 664 off += sizes[0]; 665 } 666 PetscCall(PetscFree3(x[0], x[1], x[2])); 667 } 668 PetscCall(VecRestoreArray(coordinates, &coords)); 669 670 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 671 PetscCall(VecSetBlockSize(coordinates, coordDim)); 672 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 673 PetscCall(VecDestroy(&coordinates)); 674 675 /* Read boundary conditions */ 676 PetscCall(DMGetNumLabels(*dm, &labelIdRange[0])); 677 if (rank == 0) { 678 CGNS_ENUMT(BCType_t) bctype; 679 CGNS_ENUMT(DataType_t) datatype; 680 CGNS_ENUMT(PointSetType_t) pointtype; 681 cgsize_t *points; 682 PetscReal *normals; 683 int normal[3]; 684 char *bcname = buffer; 685 cgsize_t npoints, nnormals; 686 int z, nbc, bc, c, ndatasets; 687 688 for (z = 1; z <= nzones; ++z) { 689 PetscCallCGNSRead(cg_nbocos(cgid, B, z, &nbc), *dm, 0); 690 for (bc = 1; bc <= nbc; ++bc) { 691 PetscCallCGNSRead(cg_boco_info(cgid, B, z, bc, bcname, &bctype, &pointtype, &npoints, normal, &nnormals, &datatype, &ndatasets), *dm, 0); 692 PetscCall(DMCreateLabel(*dm, bcname)); 693 PetscCall(DMGetLabel(*dm, bcname, &label)); 694 PetscCall(PetscMalloc2(npoints, &points, nnormals, &normals)); 695 PetscCallCGNSReadData(cg_boco_read(cgid, B, z, bc, points, (void *)normals), *dm, 0); 696 if (pointtype == CGNS_ENUMV(ElementRange)) { 697 // Range of cells: assuming half-open interval 698 for (c = points[0]; c < points[1]; ++c) PetscCall(DMLabelSetValue(label, c - cellStart[z - 1], 1)); 699 } else if (pointtype == CGNS_ENUMV(ElementList)) { 700 // List of cells 701 for (c = 0; c < npoints; ++c) PetscCall(DMLabelSetValue(label, points[c] - cellStart[z - 1], 1)); 702 } else if (pointtype == CGNS_ENUMV(PointRange)) { 703 CGNS_ENUMT(GridLocation_t) gridloc; 704 705 // List of points: 706 PetscCallCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 707 PetscCallCGNSRead(cg_gridlocation_read(&gridloc), *dm, 0); 708 // Range of points: assuming half-open interval 709 for (c = points[0]; c < points[1]; ++c) { 710 if (gridloc == CGNS_ENUMV(Vertex)) PetscCall(DMLabelSetValue(label, c - vertStart[z - 1], 1)); 711 else PetscCall(DMLabelSetValue(label, c - cellStart[z - 1], 1)); 712 } 713 } else if (pointtype == CGNS_ENUMV(PointList)) { 714 CGNS_ENUMT(GridLocation_t) gridloc; 715 716 // List of points: 717 PetscCallCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 718 PetscCallCGNSRead(cg_gridlocation_read(&gridloc), *dm, 0); 719 for (c = 0; c < npoints; ++c) { 720 if (gridloc == CGNS_ENUMV(Vertex)) PetscCall(DMLabelSetValue(label, points[c] - vertStart[z - 1], 1)); 721 else PetscCall(DMLabelSetValue(label, points[c] - cellStart[z - 1], 1)); 722 } 723 } else SETERRQ(comm, PETSC_ERR_SUP, "Unsupported point set type %d", (int)pointtype); 724 PetscCall(PetscFree2(points, normals)); 725 } 726 } 727 PetscCall(PetscFree2(cellStart, vertStart)); 728 } 729 PetscCall(DMGetNumLabels(*dm, &labelIdRange[1])); 730 PetscCallMPI(MPI_Bcast(labelIdRange, 2, MPIU_INT, 0, comm)); 731 732 /* Create BC labels at all processes */ 733 for (labelId = labelIdRange[0]; labelId < labelIdRange[1]; ++labelId) { 734 char *labelName = buffer; 735 size_t len = sizeof(buffer); 736 const char *locName; 737 738 if (rank == 0) { 739 PetscCall(DMGetLabelByNum(*dm, labelId, &label)); 740 PetscCall(PetscObjectGetName((PetscObject)label, &locName)); 741 PetscCall(PetscStrncpy(labelName, locName, len)); 742 } 743 PetscCallMPI(MPI_Bcast(labelName, (PetscMPIInt)len, MPIU_INT, 0, comm)); 744 PetscCallMPI(DMCreateLabel(*dm, labelName)); 745 } 746 PetscFunctionReturn(PETSC_SUCCESS); 747 } 748 749 PetscErrorCode DMPlexCreateCGNS_Internal_Parallel(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm) 750 { 751 PetscMPIInt num_proc, rank; 752 /* Read from file */ 753 char basename[CGIO_MAX_NAME_LENGTH + 1]; 754 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 755 int dim = 0, physDim = 0, coordDim = 0; 756 PetscInt NVertices = 0, NCells = 0; 757 int nzones = 0, nbases; 758 int zone = 1; // Only supports single zone files 759 int base = 1; // Only supports single base 760 761 PetscFunctionBegin; 762 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 763 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 764 PetscCall(DMCreate(comm, dm)); 765 PetscCall(DMSetType(*dm, DMPLEX)); 766 767 PetscCallCGNSRead(cg_nbases(cgid, &nbases), *dm, 0); 768 PetscCheck(nbases <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single base, not %d", nbases); 769 // From the CGNS web page cell_dim phys_dim (embedding space in PETSc) CGNS defines as length of spatial vectors/components) 770 PetscCallCGNSRead(cg_base_read(cgid, base, basename, &dim, &physDim), *dm, 0); 771 PetscCallCGNSRead(cg_nzones(cgid, base, &nzones), *dm, 0); 772 PetscCheck(nzones == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Parallel reader limited to one zone, not %d", nzones); 773 { 774 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 775 776 PetscCallCGNSRead(cg_zone_read(cgid, base, zone, buffer, sizes), *dm, 0); 777 NVertices = sizes[0]; 778 NCells = sizes[1]; 779 } 780 781 PetscCall(PetscObjectSetName((PetscObject)*dm, basename)); 782 PetscCall(DMSetDimension(*dm, dim)); 783 coordDim = dim; 784 785 // 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 786 // call to get this right but continuing for now with single section, single topology, one zone. 787 // establish element ranges for my rank 788 PetscInt mystarte, myende, mystartv, myendv, myownede, myownedv; 789 PetscLayout elem_map, vtx_map; 790 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NCells, 1, &elem_map)); 791 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NVertices, 1, &vtx_map)); 792 PetscCall(PetscLayoutGetRange(elem_map, &mystarte, &myende)); 793 PetscCall(PetscLayoutGetLocalSize(elem_map, &myownede)); 794 PetscCall(PetscLayoutGetRange(vtx_map, &mystartv, &myendv)); 795 PetscCall(PetscLayoutGetLocalSize(vtx_map, &myownedv)); 796 797 // -- Build Plex in parallel 798 DMPolytopeType dm_cell_type = DM_POLYTOPE_UNKNOWN; 799 PetscInt pOrder = 1, numClosure = -1; 800 cgsize_t *elements; 801 { 802 int nsections; 803 PetscInt *elementsQ1, numCorners = -1; 804 const int *perm; 805 cgsize_t start, end; // Throwaway 806 807 cg_nsections(cgid, base, zone, &nsections); 808 // Read element connectivity 809 for (int index_sect = 1; index_sect <= nsections; index_sect++) { 810 int nbndry, parentFlag; 811 PetscInt cell_dim; 812 CGNS_ENUMT(ElementType_t) cellType; 813 814 PetscCallCGNSRead(cg_section_read(cgid, base, zone, index_sect, buffer, &cellType, &start, &end, &nbndry, &parentFlag), *dm, 0); 815 816 PetscCall(CGNSElementTypeGetTopologyInfo(cellType, &dm_cell_type, &numCorners, &cell_dim)); 817 // Skip over element that are not max dimension (ie. boundary elements) 818 if (cell_dim != dim) continue; 819 PetscCall(CGNSElementTypeGetDiscretizationInfo(cellType, &numClosure, &pOrder)); 820 PetscCall(PetscMalloc1(myownede * numClosure, &elements)); 821 PetscCallCGNSReadData(cgp_elements_read_data(cgid, base, zone, index_sect, mystarte + 1, myende, elements), *dm, 0); 822 for (PetscInt v = 0; v < myownede * numClosure; ++v) elements[v] -= 1; // 0 based 823 break; 824 } 825 826 // Create corners-only connectivity 827 PetscCall(PetscMalloc1(myownede * numCorners, &elementsQ1)); 828 PetscCall(DMPlexCGNSGetPermutation_Internal(dm_cell_type, numCorners, NULL, &perm)); 829 for (PetscInt e = 0; e < myownede; ++e) { 830 for (PetscInt v = 0; v < numCorners; ++v) elementsQ1[e * numCorners + perm[v]] = elements[e * numClosure + v]; 831 } 832 833 // Build cell-vertex Plex 834 PetscCall(DMPlexBuildFromCellListParallel(*dm, myownede, myownedv, NVertices, numCorners, elementsQ1, NULL, NULL)); 835 PetscCall(DMViewFromOptions(*dm, NULL, "-corner_dm_view")); 836 PetscCall(PetscFree(elementsQ1)); 837 } 838 839 if (interpolate) PetscCall(DMPlexInterpolateInPlace_Internal(*dm)); 840 841 // -- Create SF for naive nodal-data read to elements 842 PetscSF plex_to_cgns_sf; 843 { 844 PetscInt nleaves, num_comp; 845 PetscInt *leaf, num_leaves = 0; 846 PetscInt cStart, cEnd; 847 const int *perm; 848 PetscSF cgns_to_local_sf; 849 PetscSection local_section; 850 PetscFE fe; 851 852 // sfNatural requires PetscSection to handle DMDistribute, so we use PetscFE to define the section 853 // Use number of components = 1 to work with just the nodes themselves 854 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, dm_cell_type, pOrder, PETSC_DETERMINE, &fe)); 855 PetscCall(PetscObjectSetName((PetscObject)fe, "FE for sfNatural")); 856 PetscCall(DMAddField(*dm, NULL, (PetscObject)fe)); 857 PetscCall(DMCreateDS(*dm)); 858 PetscCall(PetscFEDestroy(&fe)); 859 860 PetscCall(DMGetLocalSection(*dm, &local_section)); 861 PetscCall(PetscSectionViewFromOptions(local_section, NULL, "-fe_natural_section_view")); 862 PetscCall(PetscSectionGetFieldComponents(local_section, 0, &num_comp)); 863 PetscCall(PetscSectionGetStorageSize(local_section, &nleaves)); 864 nleaves /= num_comp; 865 PetscCall(PetscMalloc1(nleaves, &leaf)); 866 for (PetscInt i = 0; i < nleaves; i++) leaf[i] = -1; 867 868 // Get the permutation from CGNS closure numbering to PLEX closure numbering 869 PetscCall(DMPlexCGNSGetPermutation_Internal(dm_cell_type, numClosure, NULL, &perm)); 870 PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd)); 871 for (PetscInt cell = cStart; cell < cEnd; ++cell) { 872 PetscInt num_closure_dof, *closure_idx = NULL; 873 874 PetscCall(DMPlexGetClosureIndices(*dm, local_section, local_section, cell, PETSC_FALSE, &num_closure_dof, &closure_idx, NULL, NULL)); 875 PetscAssert(numClosure * num_comp == num_closure_dof, comm, PETSC_ERR_PLIB, "Closure dof size does not match polytope"); 876 for (PetscInt i = 0; i < numClosure; i++) { 877 PetscInt li = closure_idx[perm[i] * num_comp] / num_comp; 878 if (li < 0) continue; 879 880 PetscInt cgns_idx = elements[cell * numClosure + i]; 881 if (leaf[li] == -1) { 882 leaf[li] = cgns_idx; 883 num_leaves++; 884 } else PetscAssert(leaf[li] == cgns_idx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "leaf does not match previously set"); 885 } 886 PetscCall(DMPlexRestoreClosureIndices(*dm, local_section, local_section, cell, PETSC_FALSE, &num_closure_dof, &closure_idx, NULL, NULL)); 887 } 888 PetscAssert(num_leaves == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "leaf count in closure does not match nleaves"); 889 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)*dm), &cgns_to_local_sf)); 890 PetscCall(PetscSFSetGraphLayout(cgns_to_local_sf, vtx_map, nleaves, NULL, PETSC_USE_POINTER, leaf)); 891 PetscCall(PetscObjectSetName((PetscObject)cgns_to_local_sf, "CGNS to Plex SF")); 892 PetscCall(PetscSFViewFromOptions(cgns_to_local_sf, NULL, "-CGNStoPlex_sf_view")); 893 PetscCall(PetscFree(leaf)); 894 PetscCall(PetscFree(elements)); 895 896 { // Convert cgns_to_local to global_to_cgns 897 PetscSF sectionsf, cgns_to_global_sf; 898 899 PetscCall(DMGetSectionSF(*dm, §ionsf)); 900 PetscCall(PetscSFComposeInverse(cgns_to_local_sf, sectionsf, &cgns_to_global_sf)); 901 PetscCall(PetscSFDestroy(&cgns_to_local_sf)); 902 PetscCall(PetscSFCreateInverseSF(cgns_to_global_sf, &plex_to_cgns_sf)); 903 PetscCall(PetscObjectSetName((PetscObject)plex_to_cgns_sf, "Global Plex to CGNS")); 904 PetscCall(PetscSFDestroy(&cgns_to_global_sf)); 905 } 906 } 907 908 { // -- Set coordinates for DM 909 PetscScalar *coords; 910 float *x[3]; 911 double *xd[3]; 912 PetscBool read_with_double; 913 PetscFE cfe; 914 915 // Setup coordinate space first. Use pOrder here for isoparametric; revisit with CPEX-0045 High Order. 916 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, coordDim, dm_cell_type, pOrder, PETSC_DETERMINE, &cfe)); 917 PetscCall(DMSetCoordinateDisc(*dm, cfe, PETSC_FALSE, PETSC_FALSE)); 918 PetscCall(PetscFEDestroy(&cfe)); 919 920 { // Determine if coords are written in single or double precision 921 CGNS_ENUMT(DataType_t) datatype; 922 923 PetscCallCGNSRead(cg_coord_info(cgid, base, zone, 1, &datatype, buffer), *dm, 0); 924 read_with_double = datatype == CGNS_ENUMV(RealDouble) ? PETSC_TRUE : PETSC_FALSE; 925 } 926 927 // Read coords from file and set into component-major ordering 928 if (read_with_double) PetscCall(PetscMalloc3(myownedv, &xd[0], myownedv, &xd[1], myownedv, &xd[2])); 929 else PetscCall(PetscMalloc3(myownedv, &x[0], myownedv, &x[1], myownedv, &x[2])); 930 PetscCall(PetscMalloc1(myownedv * coordDim, &coords)); 931 { 932 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 933 cgsize_t range_min[3] = {mystartv + 1, 1, 1}; 934 cgsize_t range_max[3] = {myendv, 1, 1}; 935 int ngrids, ncoords; 936 937 PetscCallCGNSRead(cg_zone_read(cgid, base, zone, buffer, sizes), *dm, 0); 938 PetscCallCGNSRead(cg_ngrids(cgid, base, zone, &ngrids), *dm, 0); 939 PetscCheck(ngrids <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single grid, not %d", ngrids); 940 PetscCallCGNSRead(cg_ncoords(cgid, base, zone, &ncoords), *dm, 0); 941 PetscCheck(ncoords == coordDim, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a coordinate array for each dimension, not %d", ncoords); 942 if (read_with_double) { 943 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); 944 if (coordDim >= 1) { 945 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 0] = xd[0][v]; 946 } 947 if (coordDim >= 2) { 948 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 1] = xd[1][v]; 949 } 950 if (coordDim >= 3) { 951 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 2] = xd[2][v]; 952 } 953 } else { 954 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); 955 if (coordDim >= 1) { 956 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 0] = x[0][v]; 957 } 958 if (coordDim >= 2) { 959 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 1] = x[1][v]; 960 } 961 if (coordDim >= 3) { 962 for (PetscInt v = 0; v < myownedv; ++v) coords[v * coordDim + 2] = x[2][v]; 963 } 964 } 965 } 966 if (read_with_double) PetscCall(PetscFree3(xd[0], xd[1], xd[2])); 967 else PetscCall(PetscFree3(x[0], x[1], x[2])); 968 969 { // Reduce CGNS-ordered coordinate nodes to Plex ordering, and set DM's coordinates 970 Vec coord_global; 971 MPI_Datatype unit; 972 PetscScalar *coord_global_array; 973 DM cdm; 974 975 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 976 PetscCall(DMCreateGlobalVector(cdm, &coord_global)); 977 PetscCall(VecGetArrayWrite(coord_global, &coord_global_array)); 978 PetscCallMPI(MPI_Type_contiguous(coordDim, MPIU_SCALAR, &unit)); 979 PetscCallMPI(MPI_Type_commit(&unit)); 980 PetscCall(PetscSFReduceBegin(plex_to_cgns_sf, unit, coords, coord_global_array, MPI_REPLACE)); 981 PetscCall(PetscSFReduceEnd(plex_to_cgns_sf, unit, coords, coord_global_array, MPI_REPLACE)); 982 PetscCall(VecRestoreArrayWrite(coord_global, &coord_global_array)); 983 PetscCallMPI(MPI_Type_free(&unit)); 984 PetscCall(DMSetCoordinates(*dm, coord_global)); 985 PetscCall(VecDestroy(&coord_global)); 986 } 987 PetscCall(PetscFree(coords)); 988 } 989 990 // -- Set sfNatural for solution vectors in CGNS file 991 // NOTE: We set sfNatural to be the map between the original CGNS ordering of nodes and the Plex ordering of nodes. 992 PetscCall(PetscSFViewFromOptions(plex_to_cgns_sf, NULL, "-sfNatural_init_view")); 993 PetscCall(DMSetNaturalSF(*dm, plex_to_cgns_sf)); 994 PetscCall(DMSetUseNatural(*dm, PETSC_TRUE)); 995 PetscCall(PetscSFDestroy(&plex_to_cgns_sf)); 996 997 PetscCall(PetscLayoutDestroy(&elem_map)); 998 PetscCall(PetscLayoutDestroy(&vtx_map)); 999 PetscFunctionReturn(PETSC_SUCCESS); 1000 } 1001 1002 // node_l2g must be freed 1003 static PetscErrorCode DMPlexCreateNodeNumbering(DM dm, PetscInt *num_local_nodes, PetscInt *num_global_nodes, PetscInt *nStart, PetscInt *nEnd, const PetscInt **node_l2g) 1004 { 1005 PetscSection local_section; 1006 PetscSF point_sf; 1007 PetscInt pStart, pEnd, spStart, spEnd, *points, nleaves, ncomp, *nodes; 1008 PetscMPIInt comm_size; 1009 const PetscInt *ilocal, field = 0; 1010 1011 PetscFunctionBegin; 1012 *num_local_nodes = -1; 1013 *num_global_nodes = -1; 1014 *nStart = -1; 1015 *nEnd = -1; 1016 *node_l2g = NULL; 1017 1018 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &comm_size)); 1019 PetscCall(DMGetLocalSection(dm, &local_section)); 1020 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1021 PetscCall(PetscSectionGetChart(local_section, &spStart, &spEnd)); 1022 PetscCall(DMGetPointSF(dm, &point_sf)); 1023 if (comm_size == 1) nleaves = 0; 1024 else PetscCall(PetscSFGetGraph(point_sf, NULL, &nleaves, &ilocal, NULL)); 1025 PetscCall(PetscSectionGetFieldComponents(local_section, field, &ncomp)); 1026 1027 PetscInt local_node = 0, owned_node = 0, owned_start = 0; 1028 for (PetscInt p = spStart, leaf = 0; p < spEnd; p++) { 1029 PetscInt dof; 1030 PetscCall(PetscSectionGetFieldDof(local_section, p, field, &dof)); 1031 PetscAssert(dof % ncomp == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field dof %" PetscInt_FMT " must be divisible by components %" PetscInt_FMT, dof, ncomp); 1032 local_node += dof / ncomp; 1033 if (leaf < nleaves && p == ilocal[leaf]) { // skip points owned by a different process 1034 leaf++; 1035 } else { 1036 owned_node += dof / ncomp; 1037 } 1038 } 1039 PetscCallMPI(MPI_Exscan(&owned_node, &owned_start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm))); 1040 PetscCall(PetscMalloc1(pEnd - pStart, &points)); 1041 owned_node = 0; 1042 for (PetscInt p = spStart, leaf = 0; p < spEnd; p++) { 1043 if (leaf < nleaves && p == ilocal[leaf]) { // skip points owned by a different process 1044 points[p - pStart] = -1; 1045 leaf++; 1046 continue; 1047 } 1048 PetscInt dof, offset; 1049 PetscCall(PetscSectionGetFieldDof(local_section, p, field, &dof)); 1050 PetscCall(PetscSectionGetFieldOffset(local_section, p, field, &offset)); 1051 PetscAssert(offset % ncomp == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field offset %" PetscInt_FMT " must be divisible by components %" PetscInt_FMT, offset, ncomp); 1052 points[p - pStart] = owned_start + owned_node; 1053 owned_node += dof / ncomp; 1054 } 1055 if (comm_size > 1) { 1056 PetscCall(PetscSFBcastBegin(point_sf, MPIU_INT, points, points, MPI_REPLACE)); 1057 PetscCall(PetscSFBcastEnd(point_sf, MPIU_INT, points, points, MPI_REPLACE)); 1058 } 1059 1060 // Set up global indices for each local node 1061 PetscCall(PetscMalloc1(local_node, &nodes)); 1062 for (PetscInt p = spStart; p < spEnd; p++) { 1063 PetscInt dof, offset; 1064 PetscCall(PetscSectionGetFieldDof(local_section, p, field, &dof)); 1065 PetscCall(PetscSectionGetFieldOffset(local_section, p, field, &offset)); 1066 for (PetscInt n = 0; n < dof / ncomp; n++) nodes[offset / ncomp + n] = points[p - pStart] + n; 1067 } 1068 PetscCall(PetscFree(points)); 1069 *num_local_nodes = local_node; 1070 *nStart = owned_start; 1071 *nEnd = owned_start + owned_node; 1072 PetscCallMPI(MPIU_Allreduce(&owned_node, num_global_nodes, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm))); 1073 *node_l2g = nodes; 1074 PetscFunctionReturn(PETSC_SUCCESS); 1075 } 1076 1077 PetscErrorCode DMView_PlexCGNS(DM dm, PetscViewer viewer) 1078 { 1079 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 1080 PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data; 1081 PetscInt fvGhostStart; 1082 PetscInt topo_dim, coord_dim, num_global_elems; 1083 PetscInt cStart, cEnd, num_local_nodes, num_global_nodes, nStart, nEnd; 1084 const PetscInt *node_l2g; 1085 Vec coord; 1086 DM colloc_dm, cdm; 1087 PetscMPIInt size; 1088 const char *dm_name; 1089 int base, zone; 1090 cgsize_t isize[3]; 1091 1092 PetscFunctionBegin; 1093 if (!cgv->file_num) { 1094 PetscInt time_step; 1095 PetscCall(DMGetOutputSequenceNumber(dm, &time_step, NULL)); 1096 PetscCall(PetscViewerCGNSFileOpen_Internal(viewer, time_step)); 1097 } 1098 PetscCallMPI(MPI_Comm_size(comm, &size)); 1099 PetscCall(DMGetDimension(dm, &topo_dim)); 1100 PetscCall(DMGetCoordinateDim(dm, &coord_dim)); 1101 PetscCall(PetscObjectGetName((PetscObject)dm, &dm_name)); 1102 PetscCallCGNSWrite(cg_base_write(cgv->file_num, dm_name, topo_dim, coord_dim, &base), dm, viewer); 1103 PetscCallCGNS(cg_goto(cgv->file_num, base, NULL)); 1104 PetscCallCGNSWrite(cg_dataclass_write(CGNS_ENUMV(NormalizedByDimensional)), dm, viewer); 1105 1106 { 1107 PetscFE fe, fe_coord; 1108 PetscClassId ds_id; 1109 PetscDualSpace dual_space, dual_space_coord; 1110 PetscInt num_fields, field_order = -1, field_order_coord; 1111 PetscBool is_simplex; 1112 PetscCall(DMGetNumFields(dm, &num_fields)); 1113 if (num_fields > 0) { 1114 PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe)); 1115 PetscCall(PetscObjectGetClassId((PetscObject)fe, &ds_id)); 1116 if (ds_id != PETSCFE_CLASSID) { 1117 fe = NULL; 1118 if (ds_id == PETSCFV_CLASSID) field_order = -1; // use whatever is present for coords; field will be CellCenter 1119 else field_order = 1; // assume vertex-based linear elements 1120 } 1121 } else fe = NULL; 1122 if (fe) { 1123 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 1124 PetscCall(PetscDualSpaceGetOrder(dual_space, &field_order)); 1125 } 1126 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1127 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe_coord)); 1128 { 1129 PetscClassId id; 1130 PetscCall(PetscObjectGetClassId((PetscObject)fe_coord, &id)); 1131 if (id != PETSCFE_CLASSID) fe_coord = NULL; 1132 } 1133 if (fe_coord) { 1134 PetscCall(PetscFEGetDualSpace(fe_coord, &dual_space_coord)); 1135 PetscCall(PetscDualSpaceGetOrder(dual_space_coord, &field_order_coord)); 1136 } else field_order_coord = 1; 1137 if (field_order > 0 && field_order != field_order_coord) { 1138 PetscInt quadrature_order = field_order; 1139 PetscCall(DMClone(dm, &colloc_dm)); 1140 { // Inform the new colloc_dm that it is a coordinate DM so isoperiodic affine corrections can be applied 1141 const PetscSF *face_sfs; 1142 PetscInt num_face_sfs; 1143 PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &num_face_sfs, &face_sfs)); 1144 PetscCall(DMPlexSetIsoperiodicFaceSF(colloc_dm, num_face_sfs, (PetscSF *)face_sfs)); 1145 if (face_sfs) colloc_dm->periodic.setup = DMPeriodicCoordinateSetUp_Internal; 1146 } 1147 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 1148 PetscCall(PetscFECreateLagrange(comm, topo_dim, coord_dim, is_simplex, field_order, quadrature_order, &fe)); 1149 PetscCall(DMSetCoordinateDisc(colloc_dm, fe, PETSC_FALSE, PETSC_TRUE)); 1150 PetscCall(PetscFEDestroy(&fe)); 1151 } else { 1152 PetscCall(PetscObjectReference((PetscObject)dm)); 1153 colloc_dm = dm; 1154 } 1155 } 1156 PetscCall(DMGetCoordinateDM(colloc_dm, &cdm)); 1157 PetscCall(DMPlexCreateNodeNumbering(cdm, &num_local_nodes, &num_global_nodes, &nStart, &nEnd, &node_l2g)); 1158 PetscCall(DMGetCoordinatesLocal(colloc_dm, &coord)); 1159 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1160 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &fvGhostStart, NULL)); 1161 if (fvGhostStart >= 0) cEnd = fvGhostStart; 1162 num_global_elems = cEnd - cStart; 1163 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &num_global_elems, 1, MPIU_INT, MPI_SUM, comm)); 1164 isize[0] = num_global_nodes; 1165 isize[1] = num_global_elems; 1166 isize[2] = 0; 1167 PetscCallCGNSWrite(cg_zone_write(cgv->file_num, base, "Zone", isize, CGNS_ENUMV(Unstructured), &zone), dm, viewer); 1168 1169 { 1170 const PetscScalar *X; 1171 PetscScalar *x; 1172 int coord_ids[3]; 1173 1174 PetscCall(VecGetArrayRead(coord, &X)); 1175 for (PetscInt d = 0; d < coord_dim; d++) { 1176 const double exponents[] = {0, 1, 0, 0, 0}; 1177 char coord_name[64]; 1178 PetscCall(PetscSNPrintf(coord_name, sizeof coord_name, "Coordinate%c", 'X' + (int)d)); 1179 PetscCallCGNSWrite(cgp_coord_write(cgv->file_num, base, zone, CGNS_ENUMV(RealDouble), coord_name, &coord_ids[d]), dm, viewer); 1180 PetscCallCGNS(cg_goto(cgv->file_num, base, "Zone_t", zone, "GridCoordinates", 0, coord_name, 0, NULL)); 1181 PetscCallCGNSWrite(cg_exponents_write(CGNS_ENUMV(RealDouble), exponents), dm, viewer); 1182 } 1183 1184 int section; 1185 cgsize_t e_owned, e_global, e_start, *conn = NULL; 1186 const int *perm; 1187 CGNS_ENUMT(ElementType_t) element_type = CGNS_ENUMV(ElementTypeNull); 1188 { 1189 PetscCall(PetscMalloc1(nEnd - nStart, &x)); 1190 for (PetscInt d = 0; d < coord_dim; d++) { 1191 for (PetscInt n = 0; n < num_local_nodes; n++) { 1192 PetscInt gn = node_l2g[n]; 1193 if (gn < nStart || nEnd <= gn) continue; 1194 x[gn - nStart] = X[n * coord_dim + d]; 1195 } 1196 // CGNS nodes use 1-based indexing 1197 cgsize_t start = nStart + 1, end = nEnd; 1198 PetscCallCGNSWriteData(cgp_coord_write_data(cgv->file_num, base, zone, coord_ids[d], &start, &end, x), dm, viewer); 1199 } 1200 PetscCall(PetscFree(x)); 1201 PetscCall(VecRestoreArrayRead(coord, &X)); 1202 } 1203 1204 e_owned = cEnd - cStart; 1205 if (e_owned > 0) { 1206 DMPolytopeType cell_type; 1207 1208 PetscCall(DMPlexGetCellType(dm, cStart, &cell_type)); 1209 for (PetscInt i = cStart, c = 0; i < cEnd; i++) { 1210 PetscInt closure_dof, *closure_indices, elem_size; 1211 1212 PetscCall(DMPlexGetClosureIndices(cdm, cdm->localSection, cdm->localSection, i, PETSC_FALSE, &closure_dof, &closure_indices, NULL, NULL)); 1213 elem_size = closure_dof / coord_dim; 1214 if (!conn) PetscCall(PetscMalloc1(e_owned * elem_size, &conn)); 1215 PetscCall(DMPlexCGNSGetPermutation_Internal(cell_type, closure_dof / coord_dim, &element_type, &perm)); 1216 for (PetscInt j = 0; j < elem_size; j++) conn[c++] = node_l2g[closure_indices[perm[j] * coord_dim] / coord_dim] + 1; 1217 PetscCall(DMPlexRestoreClosureIndices(cdm, cdm->localSection, cdm->localSection, i, PETSC_FALSE, &closure_dof, &closure_indices, NULL, NULL)); 1218 } 1219 } 1220 1221 { // Get global element_type (for ranks that do not have owned elements) 1222 PetscInt local_element_type, global_element_type; 1223 1224 local_element_type = e_owned > 0 ? element_type : -1; 1225 PetscCallMPI(MPIU_Allreduce(&local_element_type, &global_element_type, 1, MPIU_INT, MPI_MAX, comm)); 1226 if (local_element_type != -1) PetscCheck(local_element_type == global_element_type, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ranks with different element types not supported"); 1227 element_type = (CGNS_ENUMT(ElementType_t))global_element_type; 1228 } 1229 PetscCallMPI(MPIU_Allreduce(&e_owned, &e_global, 1, MPIU_CGSIZE, MPI_SUM, comm)); 1230 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); 1231 e_start = 0; 1232 PetscCallMPI(MPI_Exscan(&e_owned, &e_start, 1, MPIU_CGSIZE, MPI_SUM, comm)); 1233 PetscCallCGNSWrite(cgp_section_write(cgv->file_num, base, zone, "Elem", element_type, 1, e_global, 0, §ion), dm, viewer); 1234 PetscCallCGNSWriteData(cgp_elements_write_data(cgv->file_num, base, zone, section, e_start + 1, e_start + e_owned, conn), dm, viewer); 1235 PetscCall(PetscFree(conn)); 1236 1237 cgv->base = base; 1238 cgv->zone = zone; 1239 cgv->node_l2g = node_l2g; 1240 cgv->num_local_nodes = num_local_nodes; 1241 cgv->nStart = nStart; 1242 cgv->nEnd = nEnd; 1243 cgv->eStart = e_start; 1244 cgv->eEnd = e_start + e_owned; 1245 if (1) { 1246 PetscMPIInt rank; 1247 int *efield; 1248 int sol, field; 1249 DMLabel label; 1250 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1251 PetscCall(PetscMalloc1(e_owned, &efield)); 1252 for (PetscInt i = 0; i < e_owned; i++) efield[i] = rank; 1253 PetscCallCGNSWrite(cg_sol_write(cgv->file_num, base, zone, "CellInfo", CGNS_ENUMV(CellCenter), &sol), dm, viewer); 1254 PetscCallCGNSWrite(cgp_field_write(cgv->file_num, base, zone, sol, CGNS_ENUMV(Integer), "Rank", &field), dm, viewer); 1255 cgsize_t start = e_start + 1, end = e_start + e_owned; 1256 PetscCallCGNSWriteData(cgp_field_write_data(cgv->file_num, base, zone, sol, field, &start, &end, efield), dm, viewer); 1257 PetscCall(DMGetLabel(dm, "Cell Sets", &label)); 1258 if (label) { 1259 for (PetscInt c = cStart; c < cEnd; c++) { 1260 PetscInt value; 1261 PetscCall(DMLabelGetValue(label, c, &value)); 1262 efield[c - cStart] = value; 1263 } 1264 PetscCallCGNSWrite(cgp_field_write(cgv->file_num, base, zone, sol, CGNS_ENUMV(Integer), "CellSet", &field), dm, viewer); 1265 PetscCallCGNSWriteData(cgp_field_write_data(cgv->file_num, base, zone, sol, field, &start, &end, efield), dm, viewer); 1266 } 1267 PetscCall(PetscFree(efield)); 1268 } 1269 } 1270 PetscCall(DMDestroy(&colloc_dm)); 1271 PetscFunctionReturn(PETSC_SUCCESS); 1272 } 1273 1274 PetscErrorCode VecView_Plex_Local_CGNS(Vec V, PetscViewer viewer) 1275 { 1276 PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data; 1277 DM dm; 1278 PetscSection section; 1279 PetscInt time_step, num_fields, pStart, pEnd, fvGhostStart; 1280 PetscReal time, *time_slot; 1281 size_t *step_slot; 1282 const PetscScalar *v; 1283 char solution_name[PETSC_MAX_PATH_LEN]; 1284 int sol; 1285 1286 PetscFunctionBegin; 1287 PetscCall(VecGetDM(V, &dm)); 1288 PetscCall(DMGetLocalSection(dm, §ion)); 1289 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 1290 PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &fvGhostStart, NULL)); 1291 if (fvGhostStart >= 0) pEnd = fvGhostStart; 1292 1293 if (!cgv->node_l2g) PetscCall(DMView(dm, viewer)); 1294 if (!cgv->grid_loc) { // Determine if writing to cell-centers or to nodes 1295 PetscInt cStart, cEnd; 1296 PetscInt local_grid_loc, global_grid_loc; 1297 1298 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1299 if (fvGhostStart >= 0) cEnd = fvGhostStart; 1300 if (cgv->num_local_nodes == 0) local_grid_loc = -1; 1301 else if (cStart == pStart && cEnd == pEnd) local_grid_loc = CGNS_ENUMV(CellCenter); 1302 else local_grid_loc = CGNS_ENUMV(Vertex); 1303 1304 PetscCallMPI(MPIU_Allreduce(&local_grid_loc, &global_grid_loc, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)viewer))); 1305 if (local_grid_loc != -1) 1306 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); 1307 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); 1308 cgv->grid_loc = (CGNS_ENUMT(GridLocation_t))global_grid_loc; 1309 } 1310 if (!cgv->nodal_field) { 1311 switch (cgv->grid_loc) { 1312 case CGNS_ENUMV(Vertex): { 1313 PetscCall(PetscMalloc1(cgv->nEnd - cgv->nStart, &cgv->nodal_field)); 1314 } break; 1315 case CGNS_ENUMV(CellCenter): { 1316 PetscCall(PetscMalloc1(cgv->eEnd - cgv->eStart, &cgv->nodal_field)); 1317 } break; 1318 default: 1319 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only write for Vertex and CellCenter grid locations"); 1320 } 1321 } 1322 if (!cgv->output_times) PetscCall(PetscSegBufferCreate(sizeof(PetscReal), 20, &cgv->output_times)); 1323 if (!cgv->output_steps) PetscCall(PetscSegBufferCreate(sizeof(size_t), 20, &cgv->output_steps)); 1324 1325 PetscCall(DMGetOutputSequenceNumber(dm, &time_step, &time)); 1326 if (time_step < 0) { 1327 time_step = 0; 1328 time = 0.; 1329 } 1330 // Avoid "Duplicate child name found" error by not writing an already-written solution. 1331 // This usually occurs when a solution is written and then diverges on the very next timestep. 1332 if (time_step == cgv->previous_output_step) PetscFunctionReturn(PETSC_SUCCESS); 1333 1334 PetscCall(PetscSegBufferGet(cgv->output_times, 1, &time_slot)); 1335 *time_slot = time; 1336 PetscCall(PetscSegBufferGet(cgv->output_steps, 1, &step_slot)); 1337 *step_slot = cgv->previous_output_step = time_step; 1338 PetscCall(PetscSNPrintf(solution_name, sizeof solution_name, "FlowSolution%" PetscInt_FMT, time_step)); 1339 PetscCallCGNSWrite(cg_sol_write(cgv->file_num, cgv->base, cgv->zone, solution_name, cgv->grid_loc, &sol), V, viewer); 1340 PetscCall(VecGetArrayRead(V, &v)); 1341 PetscCall(PetscSectionGetNumFields(section, &num_fields)); 1342 for (PetscInt field = 0; field < num_fields; field++) { 1343 PetscInt ncomp; 1344 const char *field_name; 1345 PetscCall(PetscSectionGetFieldName(section, field, &field_name)); 1346 PetscCall(PetscSectionGetFieldComponents(section, field, &ncomp)); 1347 for (PetscInt comp = 0; comp < ncomp; comp++) { 1348 int cgfield; 1349 const char *comp_name; 1350 char cgns_field_name[32]; // CGNS max field name is 32 1351 CGNS_ENUMT(DataType_t) datatype; 1352 PetscCall(PetscSectionGetComponentName(section, field, comp, &comp_name)); 1353 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)); 1354 else if (field_name[0] == '\0') PetscCall(PetscStrncpy(cgns_field_name, comp_name, sizeof cgns_field_name)); 1355 else PetscCall(PetscSNPrintf(cgns_field_name, sizeof cgns_field_name, "%s.%s", field_name, comp_name)); 1356 PetscCall(PetscCGNSDataType(PETSC_SCALAR, &datatype)); 1357 PetscCallCGNSWrite(cgp_field_write(cgv->file_num, cgv->base, cgv->zone, sol, datatype, cgns_field_name, &cgfield), V, viewer); 1358 for (PetscInt p = pStart, n = 0; p < pEnd; p++) { 1359 PetscInt off, dof; 1360 PetscCall(PetscSectionGetFieldDof(section, p, field, &dof)); 1361 if (dof == 0) continue; 1362 PetscCall(PetscSectionGetFieldOffset(section, p, field, &off)); 1363 for (PetscInt c = comp; c < dof; c += ncomp, n++) { 1364 switch (cgv->grid_loc) { 1365 case CGNS_ENUMV(Vertex): { 1366 PetscInt gn = cgv->node_l2g[n]; 1367 if (gn < cgv->nStart || cgv->nEnd <= gn) continue; 1368 cgv->nodal_field[gn - cgv->nStart] = v[off + c]; 1369 } break; 1370 case CGNS_ENUMV(CellCenter): { 1371 cgv->nodal_field[n] = v[off + c]; 1372 } break; 1373 default: 1374 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only pack for Vertex and CellCenter grid locations"); 1375 } 1376 } 1377 } 1378 // CGNS nodes use 1-based indexing 1379 cgsize_t start, end; 1380 switch (cgv->grid_loc) { 1381 case CGNS_ENUMV(Vertex): { 1382 start = cgv->nStart + 1; 1383 end = cgv->nEnd; 1384 } break; 1385 case CGNS_ENUMV(CellCenter): { 1386 start = cgv->eStart + 1; 1387 end = cgv->eEnd; 1388 } break; 1389 default: 1390 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Can only write for Vertex and CellCenter grid locations"); 1391 } 1392 PetscCallCGNSWriteData(cgp_field_write_data(cgv->file_num, cgv->base, cgv->zone, sol, cgfield, &start, &end, cgv->nodal_field), V, viewer); 1393 } 1394 } 1395 PetscCall(VecRestoreArrayRead(V, &v)); 1396 PetscCall(PetscViewerCGNSCheckBatch_Internal(viewer)); 1397 PetscFunctionReturn(PETSC_SUCCESS); 1398 } 1399 1400 PetscErrorCode VecLoad_Plex_CGNS_Internal(Vec V, PetscViewer viewer) 1401 { 1402 MPI_Comm comm; 1403 char buffer[CGIO_MAX_NAME_LENGTH + 1]; 1404 PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data; 1405 int cgid = cgv->file_num; 1406 PetscBool use_parallel_viewer = PETSC_FALSE; 1407 int z = 1; // Only support one zone 1408 int B = 1; // Only support one base 1409 int numComp; 1410 PetscInt V_numComps, mystartv, myendv, myownedv; 1411 1412 PetscFunctionBegin; 1413 PetscCall(PetscObjectGetComm((PetscObject)V, &comm)); 1414 1415 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_cgns_parallel", &use_parallel_viewer, NULL)); 1416 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"); 1417 1418 { // Get CGNS node ownership information 1419 int nbases, nzones; 1420 PetscInt NVertices; 1421 PetscLayout vtx_map; 1422 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 1423 1424 PetscCallCGNSRead(cg_nbases(cgid, &nbases), V, viewer); 1425 PetscCheck(nbases <= 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "CGNS file must have a single base, not %d", nbases); 1426 PetscCallCGNSRead(cg_nzones(cgid, B, &nzones), V, viewer); 1427 PetscCheck(nzones == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "limited to one zone %d", (int)nzones); 1428 1429 PetscCallCGNSRead(cg_zone_read(cgid, B, z, buffer, sizes), V, viewer); 1430 NVertices = sizes[0]; 1431 1432 PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, NVertices, 1, &vtx_map)); 1433 PetscCall(PetscLayoutGetRange(vtx_map, &mystartv, &myendv)); 1434 PetscCall(PetscLayoutGetLocalSize(vtx_map, &myownedv)); 1435 PetscCall(PetscLayoutDestroy(&vtx_map)); 1436 } 1437 1438 { // -- Read data from file into Vec 1439 PetscScalar *fields = NULL; 1440 PetscSF sfNatural; 1441 1442 { // Check compatibility between sfNatural and the data source and sink 1443 DM dm; 1444 PetscInt nleaves, nroots, V_local_size; 1445 1446 PetscCall(VecGetDM(V, &dm)); 1447 PetscCall(DMGetNaturalSF(dm, &sfNatural)); 1448 PetscCheck(sfNatural, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM of Vec must have sfNatural"); 1449 PetscCall(PetscSFGetGraph(sfNatural, &nroots, &nleaves, NULL, NULL)); 1450 PetscCall(VecGetLocalSize(V, &V_local_size)); 1451 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); 1452 if (nroots == 0) { 1453 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); 1454 V_numComps = 0; 1455 } else { 1456 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); 1457 V_numComps = V_local_size / nroots; 1458 } 1459 } 1460 1461 { // Read data into component-major ordering 1462 int isol, numSols; 1463 CGNS_ENUMT(DataType_t) datatype; 1464 double *fields_CGNS; 1465 1466 PetscCallCGNSRead(cg_nsols(cgid, B, z, &numSols), V, viewer); 1467 PetscCall(PetscViewerCGNSGetSolutionFileIndex_Internal(viewer, &isol)); 1468 PetscCallCGNSRead(cg_nfields(cgid, B, z, isol, &numComp), V, viewer); 1469 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); 1470 1471 cgsize_t range_min[3] = {mystartv + 1, 1, 1}; 1472 cgsize_t range_max[3] = {myendv, 1, 1}; 1473 PetscCall(PetscMalloc1(myownedv * numComp, &fields_CGNS)); 1474 PetscCall(PetscMalloc1(myownedv * numComp, &fields)); 1475 for (int d = 0; d < numComp; ++d) { 1476 PetscCallCGNSRead(cg_field_info(cgid, B, z, isol, (d + 1), &datatype, buffer), V, viewer); 1477 PetscCheck(datatype == CGNS_ENUMV(RealDouble), PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Field %s in file is not of type double", buffer); 1478 PetscCallCGNSReadData(cgp_field_read_data(cgid, B, z, isol, (d + 1), range_min, range_max, &fields_CGNS[d * myownedv]), V, viewer); 1479 } 1480 for (int d = 0; d < numComp; ++d) { 1481 for (PetscInt v = 0; v < myownedv; ++v) fields[v * numComp + d] = fields_CGNS[d * myownedv + v]; 1482 } 1483 PetscCall(PetscFree(fields_CGNS)); 1484 } 1485 1486 { // Reduce fields into Vec array 1487 PetscScalar *V_array; 1488 MPI_Datatype fieldtype; 1489 1490 PetscCall(VecGetArrayWrite(V, &V_array)); 1491 PetscCallMPI(MPI_Type_contiguous(numComp, MPIU_SCALAR, &fieldtype)); 1492 PetscCallMPI(MPI_Type_commit(&fieldtype)); 1493 PetscCall(PetscSFReduceBegin(sfNatural, fieldtype, fields, V_array, MPI_REPLACE)); 1494 PetscCall(PetscSFReduceEnd(sfNatural, fieldtype, fields, V_array, MPI_REPLACE)); 1495 PetscCallMPI(MPI_Type_free(&fieldtype)); 1496 PetscCall(VecRestoreArrayWrite(V, &V_array)); 1497 } 1498 PetscCall(PetscFree(fields)); 1499 } 1500 PetscFunctionReturn(PETSC_SUCCESS); 1501 } 1502