1 // Copyright(C) 1999-2021 National Technology & Engineering Solutions 2 // of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with 3 // NTESS, the U.S. Government retains certain rights in this software. 4 // 5 // See packages/seacas/LICENSE for details 6 7 #include <cgns/Iocgns_Defines.h> 8 9 #include <Ioss_CodeTypes.h> 10 #include <Ioss_ParallelUtils.h> 11 #include <Ioss_SmartAssert.h> 12 #include <Ioss_Sort.h> 13 #include <Ioss_StructuredBlock.h> 14 #include <Ioss_Utils.h> 15 #include <cgns/Iocgns_DecompositionData.h> 16 #include <cgns/Iocgns_Utils.h> 17 #include <fmt/color.h> 18 #include <fmt/format.h> 19 #include <fmt/ostream.h> 20 #include <tokenize.h> 21 22 #include <vtk_cgns.h> // xxx(kitware) 23 #ifdef SEACAS_HAVE_MPI 24 #include VTK_CGNS(pcgnslib.h) 25 #else 26 #include VTK_CGNS(cgnslib.h) 27 #endif 28 29 #include <algorithm> 30 #include <cassert> 31 #include <iomanip> 32 #include <numeric> 33 34 namespace { 35 int rank = 0; 36 37 // ZOLTAN Callback functions... 38 39 #if !defined(NO_ZOLTAN_SUPPORT) zoltan_num_dim(void * data,int * ierr)40 int zoltan_num_dim(void *data, int *ierr) 41 { 42 // Return dimensionality of coordinate data. 43 auto *zdata = (Iocgns::DecompositionDataBase *)(data); 44 45 *ierr = ZOLTAN_OK; 46 return zdata->spatial_dimension(); 47 } 48 zoltan_num_obj(void * data,int * ierr)49 int zoltan_num_obj(void *data, int *ierr) 50 { 51 // Return number of objects (element count) on this processor... 52 auto *zdata = (Iocgns::DecompositionDataBase *)(data); 53 54 *ierr = ZOLTAN_OK; 55 return zdata->decomp_elem_count(); 56 } 57 zoltan_obj_list(void * data,int ngid_ent,int,ZOLTAN_ID_PTR gids,ZOLTAN_ID_PTR lids,int wdim,float * wgts,int * ierr)58 void zoltan_obj_list(void *data, int ngid_ent, int /* nlid_ent */, ZOLTAN_ID_PTR gids, 59 ZOLTAN_ID_PTR lids, int wdim, float *wgts, int *ierr) 60 { 61 // Return list of object IDs, both local and global. 62 auto *zdata = (Iocgns::DecompositionDataBase *)(data); 63 64 // At the time this is called, we don't have much information 65 // These routines are the ones that are developing that 66 // information... 67 size_t element_count = zdata->decomp_elem_count(); 68 size_t element_offset = zdata->decomp_elem_offset(); 69 70 *ierr = ZOLTAN_OK; 71 72 if (lids != nullptr) { 73 std::iota(lids, lids + element_count, 0); 74 } 75 76 if (wdim != 0) { 77 std::fill(wgts, wgts + element_count, 1.0); 78 } 79 80 if (ngid_ent == 1) { 81 std::iota(gids, gids + element_count, element_offset); 82 } 83 else if (ngid_ent == 2) { 84 auto *global_ids = (int64_t *)gids; 85 std::iota(global_ids, global_ids + element_count, element_offset); 86 } 87 else { 88 *ierr = ZOLTAN_FATAL; 89 } 90 } 91 zoltan_geom(void * data,int,int,int,ZOLTAN_ID_PTR,ZOLTAN_ID_PTR,int,double * geom,int * ierr)92 void zoltan_geom(void *data, int /* ngid_ent */, int /* nlid_ent */, int /* nobj */, 93 ZOLTAN_ID_PTR /* gids */, ZOLTAN_ID_PTR /* lids */, int /* ndim */, double *geom, 94 int *ierr) 95 { 96 // Return coordinates for objects. 97 auto *zdata = (Iocgns::DecompositionDataBase *)(data); 98 99 std::copy(zdata->centroids().begin(), zdata->centroids().end(), &geom[0]); 100 101 *ierr = ZOLTAN_OK; 102 } 103 #endif 104 105 // These are used for structured parallel decomposition... create_zone_data(int cgns_file_ptr,std::vector<Iocgns::StructuredZoneData * > & zones,MPI_Comm comm)106 void create_zone_data(int cgns_file_ptr, std::vector<Iocgns::StructuredZoneData *> &zones, 107 MPI_Comm comm) 108 { 109 Ioss::ParallelUtils par_util(comm); 110 int myProcessor = par_util.parallel_rank(); // To make error macro work... 111 int base = 1; 112 int num_zones = 0; 113 114 CGCHECK(cg_nzones(cgns_file_ptr, base, &num_zones)); 115 116 std::map<std::string, int> zone_name_map; 117 118 for (cgsize_t zone = 1; zone <= num_zones; zone++) { 119 cgsize_t size[9]; 120 char zone_name[CGNS_MAX_NAME_LENGTH + 1]; 121 CGCHECK(cg_zone_read(cgns_file_ptr, base, zone, zone_name, size)); 122 zone_name_map[zone_name] = zone; 123 124 SMART_ASSERT(size[0] - 1 == size[3])(size[0])(size[3]); 125 SMART_ASSERT(size[1] - 1 == size[4])(size[1])(size[4]); 126 SMART_ASSERT(size[2] - 1 == size[5])(size[2])(size[5]); 127 128 assert(size[6] == 0); 129 assert(size[7] == 0); 130 assert(size[8] == 0); 131 132 auto *zone_data = new Iocgns::StructuredZoneData(zone_name, zone, size[3], size[4], size[5]); 133 zones.push_back(zone_data); 134 135 // Handle zone-grid-connectivity... 136 int nconn = 0; 137 CGCHECK(cg_n1to1(cgns_file_ptr, base, zone, &nconn)); 138 for (int i = 0; i < nconn; i++) { 139 char connectname[CGNS_MAX_NAME_LENGTH + 1]; 140 char donorname[CGNS_MAX_NAME_LENGTH + 1]; 141 std::array<cgsize_t, 6> range; 142 std::array<cgsize_t, 6> donor_range; 143 Ioss::IJK_t transform; 144 145 CGCHECK(cg_1to1_read(cgns_file_ptr, base, zone, i + 1, connectname, donorname, range.data(), 146 donor_range.data(), transform.data())); 147 148 // Get number of nodes shared with other "previous" zones... 149 // A "previous" zone will have a lower zone number this this zone... 150 int donor_zone = -1; 151 auto donor_iter = zone_name_map.find(donorname); 152 if (donor_iter != zone_name_map.end()) { 153 donor_zone = (*donor_iter).second; 154 } 155 Ioss::IJK_t range_beg{{(int)range[0], (int)range[1], (int)range[2]}}; 156 Ioss::IJK_t range_end{{(int)range[3], (int)range[4], (int)range[5]}}; 157 Ioss::IJK_t donor_beg{{(int)donor_range[0], (int)donor_range[1], (int)donor_range[2]}}; 158 Ioss::IJK_t donor_end{{(int)donor_range[3], (int)donor_range[4], (int)donor_range[5]}}; 159 160 #if IOSS_DEBUG_OUTPUT 161 if (rank == 0) { 162 fmt::print(Ioss::DEBUG(), "Adding zgc {} to {} donor: {}\n", connectname, zone_name, 163 donorname); 164 } 165 #endif 166 zone_data->m_zoneConnectivity.emplace_back(connectname, zone, donorname, donor_zone, 167 transform, range_beg, range_end, donor_beg, 168 donor_end); 169 } 170 } 171 172 // If there are any Structured blocks, need to iterate them and their 1-to-1 connections 173 // and update the donor_zone id for zones that had not yet been processed at the time of 174 // definition... 175 for (auto &zone : zones) { 176 for (auto &conn : zone->m_zoneConnectivity) { 177 if (conn.m_donorZone < 0) { 178 auto donor_iter = zone_name_map.find(conn.m_donorName); 179 assert(donor_iter != zone_name_map.end()); 180 conn.m_donorZone = (*donor_iter).second; 181 } 182 } 183 } 184 } 185 } // namespace 186 187 namespace Iocgns { 188 template DecompositionData<int>::DecompositionData(const Ioss::PropertyManager &props, 189 MPI_Comm communicator); 190 template DecompositionData<int64_t>::DecompositionData(const Ioss::PropertyManager &props, 191 MPI_Comm communicator); 192 193 template <typename INT> DecompositionData(const Ioss::PropertyManager & props,MPI_Comm communicator)194 DecompositionData<INT>::DecompositionData(const Ioss::PropertyManager &props, 195 MPI_Comm communicator) 196 : DecompositionDataBase(), m_decomposition(props, communicator) 197 { 198 rank = m_decomposition.m_processor; 199 200 if (props.exists("LOAD_BALANCE_THRESHOLD")) { 201 if (props.get("LOAD_BALANCE_THRESHOLD").get_type() == Ioss::Property::STRING) { 202 std::string lb_thresh = props.get("LOAD_BALANCE_THRESHOLD").get_string(); 203 m_loadBalanceThreshold = std::stod(lb_thresh); 204 } 205 else if (props.get("LOAD_BALANCE_THRESHOLD").get_type() == Ioss::Property::REAL) { 206 m_loadBalanceThreshold = props.get("LOAD_BALANCE_THRESHOLD").get_real(); 207 } 208 } 209 if (props.exists("LINE_DECOMPOSITION")) { 210 m_lineDecomposition = props.get("LINE_DECOMPOSITION").get_string(); 211 } 212 } 213 214 template <typename INT> decompose_model(int filePtr,Ioss::MeshType mesh_type)215 void DecompositionData<INT>::decompose_model(int filePtr, Ioss::MeshType mesh_type) 216 { 217 if (mesh_type == Ioss::MeshType::UNSTRUCTURED) { 218 decompose_unstructured(filePtr); 219 } 220 else if (mesh_type == Ioss::MeshType::STRUCTURED) { 221 decompose_structured(filePtr); 222 } 223 #if IOSS_ENABLE_HYBRID 224 else if (mesh_type == Ioss::MeshType::HYBRID) { 225 std::ostringstream errmsg; 226 fmt::print(errmsg, "ERROR: CGNS: The mesh type is HYBRID which is not supported for parallel " 227 "decomposition yet."); 228 IOSS_ERROR(errmsg); 229 } 230 #endif 231 else { 232 std::ostringstream errmsg; 233 fmt::print(errmsg, "ERROR: CGNS: The mesh type is not Unstructured or Structured " 234 "which are the only types currently supported"); 235 IOSS_ERROR(errmsg); 236 } 237 } 238 decompose_structured(int filePtr)239 template <typename INT> void DecompositionData<INT>::decompose_structured(int filePtr) 240 { 241 m_decomposition.show_progress(__func__); 242 create_zone_data(filePtr, m_structuredZones, m_decomposition.m_comm); 243 if (m_structuredZones.empty()) { 244 return; 245 } 246 247 #if IOSS_DEBUG_OUTPUT 248 bool verbose = true; 249 #else 250 bool verbose = false; 251 #endif 252 253 // Determine whether user has specified "line decompositions" for any of the zones. 254 // The line decomposition is an ordinal which will not be split during the 255 // decomposition. 256 if (!m_lineDecomposition.empty()) { 257 // See if the ordinal is specified as "__ordinal_{ijk}" which is used for testing... 258 if (m_lineDecomposition.find("__ordinal_") == 0) { 259 auto sub = m_lineDecomposition.substr(10); 260 unsigned int ord = 0; 261 for (size_t i = 0; i < sub.size(); i++) { 262 char ordinal = sub[i]; 263 ord |= ordinal == 'i' ? Ordinal::I : ordinal == 'j' ? Ordinal::J : Ordinal::K; 264 } 265 for (auto zone : m_structuredZones) { 266 if (zone->is_active()) { 267 zone->m_lineOrdinal |= ord; 268 } 269 } 270 } 271 else { 272 Utils::set_line_decomposition(filePtr, m_lineDecomposition, m_structuredZones, rank, 273 verbose); 274 } 275 } 276 277 // Do the processor decomposition. 278 Utils::decompose_model(m_structuredZones, m_decomposition.m_processorCount, rank, 279 m_loadBalanceThreshold, verbose); 280 281 Ioss::sort(m_structuredZones.begin(), m_structuredZones.end(), 282 [](Iocgns::StructuredZoneData *a, Iocgns::StructuredZoneData *b) { 283 return a->m_zone < b->m_zone; 284 }); 285 286 for (auto zone : m_structuredZones) { 287 if (zone->is_active()) { 288 zone->resolve_zgc_split_donor(m_structuredZones); 289 } 290 } 291 292 // Update and Output the processor assignments 293 for (auto &zone : m_structuredZones) { 294 if (zone->is_active()) { 295 zone->update_zgc_processor(m_structuredZones); 296 #if IOSS_DEBUG_OUTPUT 297 if (rank == 0) { 298 auto zone_node_count = 299 (zone->m_ordinal[0] + 1) * (zone->m_ordinal[1] + 1) * (zone->m_ordinal[2] + 1); 300 fmt::print( 301 Ioss::DEBUG(), 302 "Zone {}({}) assigned to processor {}, Adam zone = {}, Cells = {}, Nodes = {}\n", 303 zone->m_name, zone->m_zone, zone->m_proc, zone->m_adam->m_zone, zone->work(), 304 zone_node_count); 305 auto zgcs = zone->m_zoneConnectivity; 306 #if 0 307 // This should work, but doesn't... 308 fmt::print(Ioss::DEBUG(), "{}\n", fmt::join(zgcs, "\n")); 309 #else 310 for (auto &zgc : zgcs) { 311 fmt::print(Ioss::DEBUG(), "{}\n", zgc); 312 } 313 #endif 314 } 315 #endif 316 } 317 } 318 319 // Output the processor assignments in form similar to 'split' file 320 if (rank == 0) { 321 int z = 1; 322 fmt::print( 323 Ioss::OUTPUT(), 324 " n proc parent imin imax jmin jmax kmin kmax work\n"); 325 auto tmp_zone(m_structuredZones); 326 Ioss::sort(tmp_zone.begin(), tmp_zone.end(), 327 [](Iocgns::StructuredZoneData *a, Iocgns::StructuredZoneData *b) { 328 return a->m_proc < b->m_proc; 329 }); 330 331 for (auto &zone : tmp_zone) { 332 if (zone->is_active()) { 333 fmt::print(Ioss::OUTPUT(), "{:6d}{:8d}{:8d}{:8d}{:8d}{:8d}{:8d}{:8d}{:8d}{:14L}\n", z++, 334 zone->m_proc, zone->m_adam->m_zone, zone->m_offset[0] + 1, 335 zone->m_ordinal[0] + zone->m_offset[0] + 1, zone->m_offset[1] + 1, 336 zone->m_ordinal[1] + zone->m_offset[1] + 1, zone->m_offset[2] + 1, 337 zone->m_ordinal[2] + zone->m_offset[2] + 1, zone->work()); 338 } 339 } 340 } 341 342 for (auto &zone : m_structuredZones) { 343 if (!zone->is_active()) { 344 zone->m_proc = -1; 345 } 346 } 347 348 #if IOSS_DEBUG_OUTPUT 349 MPI_Barrier(m_decomposition.m_comm); 350 if (rank == 0) { 351 fmt::print(Ioss::DEBUG(), "{}", 352 fmt::format(fg(fmt::color::green), "Returning from decomposition\n")); 353 } 354 #endif 355 } 356 decompose_unstructured(int filePtr)357 template <typename INT> void DecompositionData<INT>::decompose_unstructured(int filePtr) 358 { 359 m_decomposition.show_progress(__func__); 360 361 // Initial decomposition is linear where processor #p contains 362 // elements from (#p * #element/#proc) to (#p+1 * #element/#proc) 363 364 // ======================================================================== 365 // Get the number of zones (element blocks) in the mesh... 366 int num_zones = 0; 367 int base = 1; // Only single base supported so far. 368 369 { 370 int cell_dimension = 0; 371 int phys_dimension = 0; 372 char base_name[CGNS_MAX_NAME_LENGTH + 1]; 373 CGCHECK2(cg_base_read(filePtr, base, base_name, &cell_dimension, &phys_dimension)); 374 m_decomposition.m_spatialDimension = phys_dimension; 375 } 376 377 CGCHECK2(cg_nzones(filePtr, base, &num_zones)); 378 m_zones.resize(num_zones + 1); // Use 1-based zones. 379 380 size_t global_cell_node_count = 0; 381 size_t global_element_count = 0; 382 for (int zone = 1; zone <= num_zones; zone++) { 383 // All zones are "Unstructured" since this was checked prior to 384 // calling this function... 385 cgsize_t size[3]; 386 char zone_name[CGNS_MAX_NAME_LENGTH + 1]; 387 CGCHECK2(cg_zone_read(filePtr, base, zone, zone_name, size)); 388 389 INT total_block_nodes = size[0]; 390 INT total_block_elem = size[1]; 391 392 m_zones[zone].m_nodeCount = total_block_nodes; 393 m_zones[zone].m_nodeOffset = global_cell_node_count; 394 m_zones[zone].m_name = zone_name; 395 m_zones[zone].m_elementOffset = global_element_count; 396 global_cell_node_count += total_block_nodes; 397 global_element_count += total_block_elem; 398 } 399 400 if (global_element_count < (size_t)m_decomposition.m_processorCount) { 401 std::ostringstream errmsg; 402 fmt::print(errmsg, 403 "ERROR: CGNS: Element Count ({}) is less than Processor Count ({}). No " 404 "decomposition possible.", 405 global_element_count, m_decomposition.m_processorCount); 406 IOSS_ERROR(errmsg); 407 } 408 409 // Generate element_dist/node_dist -- size m_decomposition.m_processorCount + 1 410 // processor p contains all elements/nodes from X_dist[p] .. X_dist[p+1] 411 m_decomposition.generate_entity_distributions(global_cell_node_count, global_element_count); 412 413 generate_adjacency_list(filePtr, m_decomposition); 414 415 // Get min and max node used on this processor... 416 auto min_max = 417 std::minmax_element(m_decomposition.m_adjacency.begin(), m_decomposition.m_adjacency.end()); 418 INT min_node = *(min_max.first); 419 INT max_node = *(min_max.second); 420 generate_zone_shared_nodes(filePtr, min_node, max_node); 421 422 // Now iterate adjacency list and update any "zone_shared_node" nodes 423 // with their "sharee" 424 if (!m_zoneSharedMap.empty()) { 425 for (auto &node : m_decomposition.m_adjacency) { 426 auto alias = m_zoneSharedMap.find(node); 427 if (alias != m_zoneSharedMap.end()) { 428 node = (*alias).second; 429 } 430 } 431 } 432 433 #if IOSS_DEBUG_OUTPUT 434 if (rank == 0) { 435 fmt::print(Ioss::DEBUG(), 436 "Processor {0} has {1} elements; offset = {2}\n" 437 "Processor {0} has {3} nodes; offset = {4}.\n", 438 m_decomposition.m_processor, decomp_elem_count(), decomp_elem_offset(), 439 decomp_node_count(), decomp_node_offset()); 440 } 441 #endif 442 443 if (m_decomposition.needs_centroids()) { 444 // Get my coordinate data using direct cgns calls 445 std::vector<double> x(decomp_node_count()); 446 std::vector<double> y; 447 std::vector<double> z; 448 449 get_file_node_coordinates(filePtr, 0, x.data()); 450 if (m_decomposition.m_spatialDimension > 1) { 451 y.resize(decomp_node_count()); 452 get_file_node_coordinates(filePtr, 1, y.data()); 453 } 454 if (m_decomposition.m_spatialDimension > 2) { 455 z.resize(decomp_node_count()); 456 get_file_node_coordinates(filePtr, 2, z.data()); 457 } 458 459 m_decomposition.calculate_element_centroids(x, y, z); 460 } 461 462 #if !defined(NO_ZOLTAN_SUPPORT) 463 float version = 0.0; 464 Zoltan_Initialize(0, nullptr, &version); 465 466 Zoltan zz(m_decomposition.m_comm); 467 468 // Register Zoltan Callback functions... 469 zz.Set_Num_Obj_Fn(zoltan_num_obj, this); 470 zz.Set_Obj_List_Fn(zoltan_obj_list, this); 471 zz.Set_Num_Geom_Fn(zoltan_num_dim, this); 472 zz.Set_Geom_Multi_Fn(zoltan_geom, this); 473 #endif 474 475 m_decomposition.decompose_model( 476 #if !defined(NO_ZOLTAN_SUPPORT) 477 zz, 478 #endif 479 m_elementBlocks); 480 481 if (!m_sideSets.empty()) { 482 // Create elemGTL map which is used for sidesets (also element sets) 483 build_global_to_local_elem_map(); 484 } 485 486 get_sideset_data(filePtr); 487 488 // Have all the decomposition data needed 489 // Can now populate the Ioss metadata... 490 } 491 492 template <typename INT> generate_zone_shared_nodes(int filePtr,INT min_node,INT max_node)493 void DecompositionData<INT>::generate_zone_shared_nodes(int filePtr, INT min_node, INT max_node) 494 { 495 // Begin of Zone-Shared node information 496 497 // Modify adjacency list based on shared nodes between zones... 498 // Need the map from "global" to "global-shared" 499 // * This is not necessarily nodes only on my processor since connectivity can include 500 // nodes other than those I own. 501 // * Potentially large number of shared nodes; practically small(?) 502 503 // * Maintain hash map from old id to new (if any) 504 // * TODO: Make more scalable 505 506 int base = 1; // Only single base supported so far. 507 508 // Donor zone is always lower numbered, so zone 1 has no donor zone. Start at zone 2. 509 for (cgsize_t zone = 2; zone < (cgsize_t)m_zones.size(); zone++) { 510 511 // Determine number of "shared" nodes (shared with other zones) 512 int nconn = 0; 513 CGCHECK2(cg_nconns(filePtr, base, zone, &nconn)); 514 for (int i = 0; i < nconn; i++) { 515 char connectname[CGNS_MAX_NAME_LENGTH + 1]; 516 CGNS_ENUMT(GridLocation_t) location; 517 CGNS_ENUMT(GridConnectivityType_t) connect_type; 518 CGNS_ENUMT(PointSetType_t) ptset_type; 519 cgsize_t npnts = 0; 520 char donorname[CGNS_MAX_NAME_LENGTH + 1]; 521 CGNS_ENUMT(ZoneType_t) donor_zonetype; 522 CGNS_ENUMT(PointSetType_t) donor_ptset_type; 523 CGNS_ENUMT(DataType_t) donor_datatype; 524 cgsize_t ndata_donor; 525 526 CGCHECK2(cg_conn_info(filePtr, base, zone, i + 1, connectname, &location, &connect_type, 527 &ptset_type, &npnts, donorname, &donor_zonetype, &donor_ptset_type, 528 &donor_datatype, &ndata_donor)); 529 530 if (connect_type != CGNS_ENUMV(Abutting1to1) || ptset_type != CGNS_ENUMV(PointList) || 531 donor_ptset_type != CGNS_ENUMV(PointListDonor)) { 532 std::ostringstream errmsg; 533 fmt::print(errmsg, 534 "ERROR: CGNS: Zone {} adjacency data is not correct type. Require " 535 "Abutting1to1 and PointList. {}\t{}\t{}", 536 zone, connect_type, ptset_type, donor_ptset_type); 537 IOSS_ERROR(errmsg); 538 } 539 540 // Verify data consistency... 541 if (npnts != ndata_donor) { 542 std::ostringstream errmsg; 543 fmt::print(errmsg, 544 "ERROR: CGNS: Zone {} point count ({}) does not match donor point count ({}).", 545 zone, npnts, ndata_donor); 546 IOSS_ERROR(errmsg); 547 } 548 549 // Get number of nodes shared with other "previous" zones... 550 // A "previous" zone will have a lower zone number this this zone... 551 std::string dz_name(donorname); 552 int dz = 1; 553 for (; dz < zone; dz++) { 554 if (m_zones[dz].m_name == dz_name) 555 break; 556 } 557 558 if (dz != zone) { 559 #if IOSS_DEBUG_OUTPUT 560 if (m_decomposition.m_processor == 0) { 561 fmt::print(Ioss::DEBUG(), "Zone {} shares {} nodes with {}\n", zone, npnts, donorname); 562 } 563 #endif 564 // The 'ids' in 'points' and 'donors' will be zone-local 1-based. 565 CGNSIntVector points(npnts); 566 CGNSIntVector donors(npnts); 567 568 CGCHECK2(cg_conn_read(filePtr, base, zone, i + 1, points.data(), donor_datatype, 569 donors.data())); 570 571 for (int j = 0; j < npnts; j++) { 572 // Convert to 0-based global id by subtracting 1 and adding zone.m_nodeOffset 573 cgsize_t point = points[j] - 1 + m_zones[zone].m_nodeOffset; 574 cgsize_t donor = donors[j] - 1 + m_zones[dz].m_nodeOffset; 575 576 // See if 'donor' is mapped to a different node already 577 auto donor_map = m_zoneSharedMap.find(donor); 578 if (donor_map != m_zoneSharedMap.end()) { 579 donor = (*donor_map).second; 580 } 581 m_zoneSharedMap.insert({point, donor}); 582 #if IOSS_DEBUG_OUTPUT 583 if (m_decomposition.m_processor == 0) { 584 fmt::print(Ioss::DEBUG(), "Inserted {} to {}\n", point, donor); 585 } 586 #endif 587 } 588 } 589 } 590 } 591 // Filter m_zoneSharedMap down to nodes on this processor... 592 // This processor contains global zone ids from `min_node` to `max_node` 593 // global zone ids are the first entry in m_zoneShardedMap. 594 for (auto it = m_zoneSharedMap.cbegin(); it != m_zoneSharedMap.cend(); /* no increment */) { 595 if ((*it).first < min_node || (*it).first > max_node) { 596 it = m_zoneSharedMap.erase(it); 597 } 598 else { 599 ++it; 600 } 601 } 602 } 603 604 template <typename INT> generate_adjacency_list(int filePtr,Ioss::Decomposition<INT> & decomposition)605 void DecompositionData<INT>::generate_adjacency_list(int filePtr, 606 Ioss::Decomposition<INT> &decomposition) 607 { 608 int base = 1; // Only single base supported so far. 609 610 // Range of elements currently handled by this processor [) 611 size_t p_start = decomp_elem_offset(); 612 size_t p_end = p_start + decomp_elem_count(); 613 614 size_t sum = 0; // Size of adjacency vector. 615 size_t offset = 0; 616 617 int num_zones = 0; 618 INT zone_node_offset = 0; 619 620 CGCHECK2(cg_nzones(filePtr, base, &num_zones)); 621 for (int zone = 1; zone <= num_zones; zone++) { 622 // ======================================================================== 623 // Read the ZoneBC_t node to get list of SideBlocks to define on this zone 624 // The BC_t nodes in the ZoneBC_t give the element range for each SideBlock 625 // which can be matched up below with the Elements_t nodes to get contents 626 // of the SideBlocks. 627 auto zonebc = 628 Utils::parse_zonebc_sideblocks(filePtr, base, zone, m_decomposition.m_processor); 629 630 cgsize_t size[3]; 631 char zone_name[CGNS_MAX_NAME_LENGTH + 1]; 632 CGCHECK2(cg_zone_read(filePtr, base, zone, zone_name, size)); 633 634 INT total_elements = size[1]; 635 // NOTE: A Zone will have a single set of nodes, but can have 636 // multiple sections each with their own element type... 637 // Keep treating sections as element blocks until we 638 // have handled 'size[1]' number of elements; the remaining 639 // sections are then the boundary faces (?) 640 int num_sections = 0; 641 CGCHECK2(cg_nsections(filePtr, base, zone, &num_sections)); 642 643 size_t last_blk_location = 0; 644 for (int is = 1; is <= num_sections; is++) { 645 char section_name[CGNS_MAX_NAME_LENGTH + 1]; 646 CGNS_ENUMT(ElementType_t) e_type; 647 cgsize_t el_start = 0; 648 cgsize_t el_end = 0; 649 int num_bndry = 0; 650 int parent_flag = 0; 651 652 // Get the type of elements in this section... 653 CGCHECK2(cg_section_read(filePtr, base, zone, is, section_name, &e_type, &el_start, &el_end, 654 &num_bndry, &parent_flag)); 655 656 INT num_entity = el_end - el_start + 1; 657 658 if (parent_flag == 0 && total_elements > 0) { 659 total_elements -= num_entity; 660 661 // Range of elements in element block b [) 662 size_t b_start = offset; // offset is index of first element in this block... 663 offset += num_entity; 664 size_t b_end = offset; 665 666 int element_nodes; 667 CGCHECK2(cg_npe(e_type, &element_nodes)); 668 669 if (b_start < p_end && p_start < b_end) { 670 // Some of this blocks elements are on this processor... 671 size_t overlap = std::min(b_end, p_end) - std::max(b_start, p_start); 672 sum += overlap * element_nodes; 673 } 674 675 Ioss::BlockDecompositionData block; 676 block.zone_ = zone; 677 block.section_ = is; 678 block.name_ = zone_name; 679 block.topologyType = Utils::map_cgns_to_topology_type(e_type); 680 block.nodesPerEntity = element_nodes; 681 block.fileCount = num_entity; 682 block.zoneNodeOffset = zone_node_offset; 683 684 last_blk_location = m_elementBlocks.size(); 685 m_elementBlocks.push_back(block); 686 } 687 else { 688 // This is a boundary-condition -- sideset (?) 689 std::string bc_name(section_name); 690 std::string ss_name; 691 // Search zonebc (if it exists) for an entry such that the element ranges overlap. 692 if (!zonebc.empty()) { 693 size_t i = 0; 694 for (; i < zonebc.size(); i++) { 695 if (zonebc[i].range_beg >= el_start && zonebc[i].range_end <= el_end) { 696 break; 697 } 698 } 699 if (i < zonebc.size()) { 700 ss_name = zonebc[i].name; 701 } 702 } 703 else { 704 ss_name = section_name; 705 } 706 707 Ioss::SetDecompositionData sset; 708 sset.zone_ = zone; 709 sset.section_ = is; 710 sset.name_ = bc_name; 711 sset.ss_name_ = ss_name; 712 sset.fileCount = num_entity; 713 sset.topologyType = Utils::map_cgns_to_topology_type(e_type); 714 sset.parentBlockIndex = last_blk_location; 715 m_sideSets.emplace_back(std::move(sset)); 716 } 717 } 718 zone_node_offset += size[0]; 719 } 720 int block_count = (int)m_elementBlocks.size(); 721 722 // Get the global element block index list at this time also. 723 // The global element at index 'I' (0-based) is on block B 724 // if global_block_index[B] <= I && global_block_index[B+1] < I 725 // allocate and TODO: Fill 726 m_decomposition.m_fileBlockIndex.reserve(block_count + 1); 727 for (auto block : m_elementBlocks) { 728 m_decomposition.m_fileBlockIndex.push_back(block.file_count()); 729 } 730 m_decomposition.m_fileBlockIndex.push_back(0); 731 Ioss::Utils::generate_index(m_decomposition.m_fileBlockIndex); 732 733 // Make sure 'sum' can fit in INT... 734 INT tmp_sum = (INT)sum; 735 if ((size_t)tmp_sum != sum) { 736 std::ostringstream errmsg; 737 fmt::print( 738 errmsg, 739 "ERROR: The decomposition of this mesh requires 64-bit integers, but is being\n" 740 " run with 32-bit integer code. Please rerun with the property INTEGER_SIZE_API\n" 741 " set to 8. The details of how to do this vary with the code that is being run.\n" 742 " Contact gdsjaar@sandia.gov for more details.\n"); 743 IOSS_ERROR(errmsg); 744 } 745 746 // Now, populate the vectors... 747 decomposition.m_pointer.reserve(decomp_elem_count() + 1); 748 decomposition.m_adjacency.reserve(sum); 749 offset = 0; 750 sum = 0; // Size of adjacency vector. 751 752 for (auto &block : m_elementBlocks) { 753 // Range of elements in element block b [) 754 size_t b_start = offset; // offset is index of first element in this block... 755 offset += block.file_count(); 756 size_t b_end = b_start + block.file_count(); 757 758 ssize_t overlap = std::min(b_end, p_end) - std::max(b_start, p_start); 759 overlap = std::max(overlap, (ssize_t)0); 760 block.fileCount = overlap; 761 size_t element_nodes = block.nodesPerEntity; 762 int zone = block.zone_; 763 int section = block.section_; 764 765 // Get the connectivity (raw) for this portion of elements... 766 CGNSIntVector connectivity(overlap * element_nodes); 767 INT blk_start = std::max(b_start, p_start) - b_start + 1; 768 INT blk_end = blk_start + overlap - 1; 769 blk_start = blk_start < 0 ? 0 : blk_start; 770 blk_end = blk_end < 0 ? 0 : blk_end; 771 #if IOSS_DEBUG_OUTPUT 772 if (rank == 0) { 773 fmt::print(Ioss::DEBUG(), "Processor {} has {} elements on element block {}\t({} to {})\n", 774 m_decomposition.m_processor, overlap, block.name(), blk_start, blk_end); 775 } 776 #endif 777 block.fileSectionOffset = blk_start; 778 CGCHECK2(cgp_elements_read_data(filePtr, base, zone, section, blk_start, blk_end, 779 connectivity.data())); 780 size_t el = 0; 781 INT zone_offset = block.zoneNodeOffset; 782 783 for (ssize_t elem = 0; elem < overlap; elem++) { 784 decomposition.m_pointer.push_back(decomposition.m_adjacency.size()); 785 for (size_t k = 0; k < element_nodes; k++) { 786 INT node = connectivity[el++] - 1 + zone_offset; // 0-based node 787 decomposition.m_adjacency.push_back(node); 788 } 789 } 790 sum += overlap * element_nodes; 791 } 792 decomposition.m_pointer.push_back(decomposition.m_adjacency.size()); 793 } 794 get_sideset_data(int filePtr)795 template <typename INT> void DecompositionData<INT>::get_sideset_data(int filePtr) 796 { 797 int base = 1; // Only single base supported so far. 798 799 // NOTE: Not currently used; assume can read all on single processor... 800 // Calculate the max "buffer" size usable for storing sideset 801 // elemlists. This is basically the space used to store the file 802 // decomposition nodal coordinates. The "decomp_node_count()/2*2" is to 803 // equalize the decomp_node_count() among processors since some procs have 1 804 // more node than others. For small models, assume we can handle 805 // at least 10000 nodes. 806 // size_t max_size = std::max(10000, (decomp_node_count() / 2) * 2 * 3 *sizeof(double) / 807 // sizeof(cgsize_t)); 808 809 bool subsetting = false; 810 811 if (subsetting) { 812 assert(1 == 0); 813 } 814 else { 815 for (auto &sset : m_sideSets) { 816 817 auto topology = Ioss::ElementTopology::factory(sset.topologyType, true); 818 int nodes_per_face = topology->number_nodes(); 819 CGNSIntVector nodes(nodes_per_face * sset.file_count()); 820 821 // We get: 822 // * num_to_get parent elements, 823 // * num_to_get zeros (other parent element for face, but on boundary so 0) 824 // * num_to_get face_on_element 825 // * num_to_get zeros (face on other parent element) 826 CGNSIntVector parent(4 * sset.file_count()); 827 828 CGCHECK2(cg_elements_read(filePtr, base, sset.zone(), sset.section(), nodes.data(), 829 parent.data())); 830 831 if (parent[0] == 0) { 832 // Get rid of 'parent' list -- not used. 833 Ioss::Utils::clear(parent); 834 835 // This model does not contain parent/face data; it only contains the 836 // face connectivity (nodes) data. Need to construct parent/face data 837 // The faces in the list should all be boundaries of the block and should 838 // therefore, not be shared with another block or shared across processors. 839 // If we check the 'corner_node_cnt' nodes of the face and they are all 840 // on this processor, then assume the face is on this processor... 841 if (m_boundaryFaces[sset.zone()].empty()) { 842 // Need map of proc-zone-implicit element ids to global-zone-implicit ids 843 // so we can assign the correct element id to the faces. 844 auto blk = m_elementBlocks[sset.zone() - 1]; 845 std::vector<INT> file_data(blk.fileCount); 846 std::iota(file_data.begin(), file_data.end(), blk.fileSectionOffset); 847 std::vector<INT> zone_local_zone_global(blk.iossCount); 848 communicate_element_data(file_data.data(), zone_local_zone_global.data(), 1); 849 Ioss::Utils::clear(file_data); 850 851 std::vector<INT> connectivity(blk.ioss_count() * blk.nodesPerEntity); 852 get_block_connectivity(filePtr, connectivity.data(), sset.zone() - 1, true); 853 854 auto topo = Ioss::ElementTopology::factory(blk.topologyType, true); 855 // Should map the connectivity from cgns to ioss, but only use the lower order which is 856 // same. 857 Iocgns::Utils::generate_block_faces(topo, blk.ioss_count(), connectivity, 858 m_boundaryFaces[sset.zone()], 859 zone_local_zone_global); 860 } 861 862 // TODO: Should we filter down to just corner nodes? 863 // Now, iterate the face connectivity vector and see if 864 // there is a match in `m_boundaryFaces` 865 size_t offset = 0; 866 auto & boundary = m_boundaryFaces[sset.zone()]; 867 int num_corner_nodes = topology->number_corner_nodes(); 868 SMART_ASSERT(num_corner_nodes == 3 || num_corner_nodes == 4)(num_corner_nodes); 869 870 for (size_t iface = 0; iface < sset.file_count(); iface++) { 871 std::array<size_t, 4> conn = {{0, 0, 0, 0}}; 872 873 for (int i = 0; i < num_corner_nodes; i++) { 874 conn[i] = nodes[offset + i]; 875 } 876 offset += nodes_per_face; 877 878 Ioss::Face face(conn); 879 // See if face is in m_boundaryFaces 880 // If not, then owned by another rank 881 // If so, then get parent element and element side. 882 auto it = boundary.find(face); 883 if (it != boundary.end()) { 884 sset.entitylist_map.push_back(iface); 885 } 886 } 887 } 888 else { 889 size_t zone_element_id_offset = m_zones[sset.zone()].m_elementOffset; 890 for (size_t i = 0; i < sset.file_count(); i++) { 891 cgsize_t elem = parent[i] + zone_element_id_offset; 892 // See if elem owned by this processor... 893 if (i_own_elem(elem)) { 894 // Save elem in this processors element list for this set. 895 // The saved data is this elements location in the global 896 // element list (parent) for this set. 897 sset.entitylist_map.push_back(i); 898 } 899 } 900 } 901 } 902 903 // Each processor knows how many of the sideset elems it owns; 904 // broadcast that information (the count) to the other 905 // processors. The first processor with non-zero elem count is 906 // the "root" for this sideset. 907 { 908 std::vector<int> has_elems_local(m_sideSets.size()); 909 for (size_t i = 0; i < m_sideSets.size(); i++) { 910 has_elems_local[i] = m_sideSets[i].entitylist_map.empty() ? 0 : 1; 911 } 912 913 std::vector<int> has_elems(m_sideSets.size() * m_decomposition.m_processorCount); 914 MPI_Allgather(has_elems_local.data(), has_elems_local.size(), MPI_INT, has_elems.data(), 915 has_elems_local.size(), MPI_INT, m_decomposition.m_comm); 916 917 for (size_t i = 0; i < m_sideSets.size(); i++) { 918 m_sideSets[i].hasEntities.resize(m_decomposition.m_processorCount); 919 m_sideSets[i].root_ = m_decomposition.m_processorCount; 920 for (int p = 0; p < m_decomposition.m_processorCount; p++) { 921 if (p < m_sideSets[i].root_ && has_elems[p * m_sideSets.size() + i] != 0) { 922 m_sideSets[i].root_ = p; 923 } 924 m_sideSets[i].hasEntities[p] = has_elems[p * m_sideSets.size() + i]; 925 } 926 int color = m_sideSets[i].hasEntities[m_decomposition.m_processor] ? 1 : MPI_UNDEFINED; 927 MPI_Comm_split(m_decomposition.m_comm, color, m_decomposition.m_processor, 928 &m_sideSets[i].setComm_); 929 } 930 } 931 } 932 } 933 934 template <typename INT> get_file_node_coordinates(int filePtr,int direction,double * data)935 void DecompositionData<INT>::get_file_node_coordinates(int filePtr, int direction, 936 double *data) const 937 { 938 int base = 1; // Only single base supported so far. 939 cgsize_t beg = 0; 940 cgsize_t end = 0; 941 cgsize_t offset = 0; 942 cgsize_t node_count = decomp_node_count(); 943 cgsize_t node_offset = decomp_node_offset(); 944 945 int num_zones = (int)m_zones.size() - 1; 946 for (int zone = 1; zone <= num_zones; zone++) { 947 end += m_zones[zone].m_nodeCount; 948 949 cgsize_t start = std::max(node_offset, beg); 950 cgsize_t finish = std::min(end, node_offset + node_count); 951 cgsize_t count = (finish > start) ? finish - start : 0; 952 953 // Now adjust start for 1-based node numbering and the start of this zone... 954 start = start - beg + 1; 955 finish = finish - beg; 956 if (count == 0) { 957 start = 0; 958 finish = 0; 959 } 960 #if IOSS_DEBUG_OUTPUT 961 if (rank == 0) { 962 fmt::print( 963 Ioss::DEBUG(), 964 "{}: reading {} nodes from zone {} starting at {} with an offset of {} ending at {}\n", 965 m_decomposition.m_processor, count, zone, start, offset, finish); 966 } 967 #endif 968 double *coords = nullptr; 969 if (count > 0) { 970 coords = &data[offset]; 971 } 972 CGCHECK2(cgp_coord_read_data(filePtr, base, zone, direction + 1, &start, &finish, coords)); 973 offset += count; 974 beg = end; 975 } 976 } 977 978 template <typename INT> get_node_coordinates(int filePtr,double * ioss_data,const Ioss::Field & field)979 void DecompositionData<INT>::get_node_coordinates(int filePtr, double *ioss_data, 980 const Ioss::Field &field) const 981 { 982 std::vector<double> tmp(decomp_node_count()); 983 if (field.get_name() == "mesh_model_coordinates_x") { 984 get_file_node_coordinates(filePtr, 0, tmp.data()); 985 communicate_node_data(tmp.data(), ioss_data, 1); 986 } 987 988 else if (field.get_name() == "mesh_model_coordinates_y") { 989 get_file_node_coordinates(filePtr, 1, tmp.data()); 990 communicate_node_data(tmp.data(), ioss_data, 1); 991 } 992 993 else if (field.get_name() == "mesh_model_coordinates_z") { 994 get_file_node_coordinates(filePtr, 2, tmp.data()); 995 communicate_node_data(tmp.data(), ioss_data, 1); 996 } 997 998 else if (field.get_name() == "mesh_model_coordinates") { 999 // Data required by upper classes store x0, y0, z0, ... xn, 1000 // yn, zn. Data stored in cgns file is x0, ..., xn, y0, 1001 // ..., yn, z0, ..., zn so we have to allocate some scratch 1002 // memory to read in the data and then map into supplied 1003 // 'data' 1004 1005 std::vector<double> ioss_tmp(ioss_node_count()); 1006 1007 // This implementation trades off extra communication for 1008 // reduced memory overhead. 1009 // * This method uses 'ioss_node_count' extra memory; 3 1010 // reads; and 3 communicate_node_data calls. 1011 // 1012 // * Other method uses 6*ioss_node_count extra memory; 3 reads; 1013 // and 1 communicate_node_data call. 1014 // 1015 for (int d = 0; d < m_decomposition.m_spatialDimension; d++) { 1016 get_file_node_coordinates(filePtr, d, tmp.data()); 1017 communicate_node_data(tmp.data(), ioss_tmp.data(), 1); 1018 1019 size_t index = d; 1020 for (size_t i = 0; i < ioss_node_count(); i++) { 1021 ioss_data[index] = ioss_tmp[i]; 1022 index += m_decomposition.m_spatialDimension; 1023 } 1024 } 1025 } 1026 } 1027 1028 template <typename INT> get_node_field(int filePtr,int step,int field_offset,double * ioss_data)1029 void DecompositionData<INT>::get_node_field(int filePtr, int step, int field_offset, 1030 double *ioss_data) const 1031 { 1032 std::vector<double> tmp(decomp_node_count()); 1033 1034 int base = 1; // Only single base supported so far. 1035 cgsize_t beg = 0; 1036 cgsize_t end = 0; 1037 cgsize_t offset = 0; 1038 cgsize_t node_count = decomp_node_count(); 1039 cgsize_t node_offset = decomp_node_offset(); 1040 1041 int num_zones = (int)m_zones.size() - 1; 1042 for (int zone = 1; zone <= num_zones; zone++) { 1043 end += m_zones[zone].m_nodeCount; 1044 1045 int solution_index = Utils::find_solution_index(filePtr, base, zone, step, CGNS_ENUMV(Vertex)); 1046 1047 cgsize_t start = std::max(node_offset, beg); 1048 cgsize_t finish = std::min(end, node_offset + node_count); 1049 cgsize_t count = (finish > start) ? finish - start : 0; 1050 1051 // Now adjust start for 1-based node numbering and the start of this zone... 1052 start = (count == 0) ? 0 : start - beg + 1; 1053 finish = (count == 0) ? 0 : finish - beg; 1054 1055 double * data = (count > 0) ? &tmp[offset] : nullptr; 1056 cgsize_t range_min[1] = {start}; 1057 cgsize_t range_max[1] = {finish}; 1058 1059 CGCHECK2(cgp_field_read_data(filePtr, base, zone, solution_index, field_offset, range_min, 1060 range_max, data)); 1061 1062 offset += count; 1063 beg = end; 1064 } 1065 communicate_node_data(tmp.data(), ioss_data, 1); 1066 } 1067 1068 template void DecompositionData<int>::get_sideset_element_side( 1069 int filePtr, const Ioss::SetDecompositionData &sset, int *data) const; 1070 template void DecompositionData<int64_t>::get_sideset_element_side( 1071 int filePtr, const Ioss::SetDecompositionData &sset, int64_t *data) const; 1072 template <typename INT> get_sideset_element_side(int filePtr,const Ioss::SetDecompositionData & sset,INT * ioss_data)1073 void DecompositionData<INT>::get_sideset_element_side(int filePtr, 1074 const Ioss::SetDecompositionData &sset, 1075 INT *ioss_data) const 1076 { 1077 std::vector<INT> element_side; 1078 int base = 1; 1079 1080 auto topology = Ioss::ElementTopology::factory(sset.topologyType, true); 1081 int nodes_per_face = topology->number_nodes(); 1082 CGNSIntVector nodes(nodes_per_face * sset.file_count()); 1083 1084 CGNSIntVector parent(4 * sset.file_count()); 1085 1086 CGCHECK2( 1087 cg_elements_read(filePtr, base, sset.zone(), sset.section(), nodes.data(), parent.data())); 1088 1089 if (parent[0] == 0) { 1090 // Get rid of 'parent' list -- not used. 1091 Ioss::Utils::clear(parent); 1092 1093 // This model does not contain parent/face data; it only contains the 1094 // face connectivity (nodes) data. Need to construct parent/face data 1095 // The faces in the list should all be boundaries of the block and should 1096 // therefore, not be shared with another block or shared across processors. 1097 // If we check the 'corner_node_cnt' nodes of the face and they are all 1098 // on this processor, then assume the face is on this processor... 1099 1100 // TODO: Should we filter down to just corner nodes? 1101 CGNSIntVector face_nodes(sset.entitylist_map.size() * nodes_per_face); 1102 communicate_set_data(nodes.data(), face_nodes.data(), sset, nodes_per_face); 1103 1104 // Now, iterate the face connectivity vector and find a match in `m_boundaryFaces` 1105 size_t offset = 0; 1106 size_t j = 0; 1107 1108 // NOTE: The boundary face generation doesn't filter proc-boundary faces, 1109 // so all zones will have boundary faces generated in 'get_sideset_data` 1110 assert(!m_boundaryFaces[sset.zone()].empty()); 1111 auto &boundary = m_boundaryFaces[sset.zone()]; 1112 1113 int num_corner_nodes = topology->number_corner_nodes(); 1114 SMART_ASSERT(num_corner_nodes == 3 || num_corner_nodes == 4)(num_corner_nodes); 1115 1116 for (size_t iface = 0; iface < sset.ioss_count(); iface++) { 1117 std::array<size_t, 4> conn = {{0, 0, 0, 0}}; 1118 1119 for (int i = 0; i < num_corner_nodes; i++) { 1120 conn[i] = face_nodes[offset + i]; 1121 } 1122 offset += nodes_per_face; 1123 1124 size_t zone_element_id_offset = m_zones[sset.zone()].m_elementOffset; 1125 Ioss::Face face(conn); 1126 // See if face is in m_boundaryFaces 1127 // If not, error 1128 // If so, then get parent element and element side. 1129 auto it = boundary.find(face); 1130 if (it != boundary.end()) { 1131 cgsize_t fid = (*it).element[0]; 1132 #if IOSS_DEBUG_OUTPUT 1133 fmt::print(Ioss::DEBUG(), "Connectivity: {} {} {} {} maps to element {}, face {}\n", 1134 conn[0], conn[1], conn[2], conn[3], fid / 10, fid % 10 + 1); 1135 #endif 1136 ioss_data[j++] = fid / 10 + zone_element_id_offset; 1137 ioss_data[j++] = fid % 10 + 1; 1138 } 1139 else { 1140 std::ostringstream errmsg; 1141 fmt::print(errmsg, 1142 "ERROR: CGNS: Could not find face with connectivity {} {} {} {} on " 1143 "sideblock {}.", 1144 conn[0], conn[1], conn[2], conn[3], sset.name()); 1145 IOSS_ERROR(errmsg); 1146 } 1147 } 1148 } 1149 else { 1150 // Get rid of 'nodes' list -- not used. 1151 Ioss::Utils::clear(nodes); 1152 1153 // `parent` contains: 1154 // `num_to_get` parent elements, followed by -- 1155 // `num_to_get` zeros (other parent element for face, but on boundary so 0) 1156 // `num_to_get` face_on_element 1157 // `num_to_get` zeros (face on other parent element) 1158 1159 // Move from 'parent' to 'element_side' and interleave. element, side, element, side, ... 1160 element_side.reserve(sset.file_count() * 2); 1161 size_t zone_element_id_offset = m_zones[sset.zone()].m_elementOffset; 1162 for (size_t i = 0; i < sset.file_count(); i++) { 1163 element_side.push_back(parent[0 * sset.file_count() + i] + zone_element_id_offset); 1164 element_side.push_back(parent[2 * sset.file_count() + i]); 1165 } 1166 auto blk = m_elementBlocks[sset.zone() - 1]; 1167 auto topo = Ioss::ElementTopology::factory(blk.topologyType, true); 1168 Utils::map_cgns_face_to_ioss(topo, sset.file_count(), element_side.data()); 1169 // The above was all on root processor for this side set, now need to send data to other 1170 // processors that own any of the elements in the sideset. 1171 communicate_set_data(element_side.data(), ioss_data, sset, 2); 1172 } 1173 } 1174 1175 #ifndef DOXYGEN_SKIP_THIS 1176 template void DecompositionData<int>::get_block_connectivity(int filePtr, int *data, int blk_seq, 1177 bool raw_ids) const; 1178 template void DecompositionData<int64_t>::get_block_connectivity(int filePtr, int64_t *data, 1179 int blk_seq, bool raw_ids) const; 1180 #endif 1181 1182 template <typename INT> get_block_connectivity(int filePtr,INT * data,int blk_seq,bool raw_ids)1183 void DecompositionData<INT>::get_block_connectivity(int filePtr, INT *data, int blk_seq, 1184 bool raw_ids) const 1185 { 1186 auto blk = m_elementBlocks[blk_seq]; 1187 CGNSIntVector file_conn(blk.file_count() * blk.nodesPerEntity); 1188 int base = 1; 1189 CGCHECK2(cgp_elements_read_data(filePtr, base, blk.zone(), blk.section(), blk.fileSectionOffset, 1190 blk.fileSectionOffset + blk.file_count() - 1, 1191 file_conn.data())); 1192 1193 if (!raw_ids) { 1194 // Map from zone-local node numbers to global implicit 1195 if (blk.zoneNodeOffset != 0) { 1196 for (auto &node : file_conn) { 1197 node += blk.zoneNodeOffset; 1198 } 1199 } 1200 1201 if (!m_zoneSharedMap.empty()) { 1202 for (auto &node : file_conn) { 1203 ZoneSharedMap::const_iterator alias = m_zoneSharedMap.find(node - 1); 1204 if (alias != m_zoneSharedMap.end()) { 1205 node = (*alias).second + 1; 1206 } 1207 } 1208 } 1209 } 1210 1211 communicate_block_data(file_conn.data(), data, blk, (size_t)blk.nodesPerEntity); 1212 } 1213 1214 #ifndef DOXYGEN_SKIP_THIS 1215 template void DecompositionData<int>::get_element_field(int filePtr, int solution_index, 1216 int blk_seq, int field_index, 1217 double *data) const; 1218 template void DecompositionData<int64_t>::get_element_field(int filePtr, int solution_index, 1219 int blk_seq, int field_index, 1220 double *data) const; 1221 #endif 1222 1223 template <typename INT> get_element_field(int filePtr,int solution_index,int blk_seq,int field_index,double * data)1224 void DecompositionData<INT>::get_element_field(int filePtr, int solution_index, int blk_seq, 1225 int field_index, double *data) const 1226 { 1227 const auto blk = m_elementBlocks[blk_seq]; 1228 std::vector<double> cgns_data(blk.file_count()); 1229 int base = 1; 1230 cgsize_t range_min[1] = {(cgsize_t)blk.fileSectionOffset}; 1231 cgsize_t range_max[1] = {(cgsize_t)(blk.fileSectionOffset + blk.file_count() - 1)}; 1232 1233 CGCHECK2(cgp_field_read_data(filePtr, base, blk.zone(), solution_index, field_index, range_min, 1234 range_max, cgns_data.data())); 1235 1236 communicate_block_data(cgns_data.data(), data, blk, (size_t)1); 1237 } 1238 ~DecompositionDataBase()1239 DecompositionDataBase::~DecompositionDataBase() 1240 { 1241 for (auto &zone : m_structuredZones) { 1242 delete zone; 1243 } 1244 } 1245 1246 #ifndef DOXYGEN_SKIP_THIS 1247 template void DecompositionDataBase::communicate_node_data(int *file_data, int *ioss_data, 1248 size_t comp_count) const; 1249 template void DecompositionDataBase::communicate_node_data(int64_t *file_data, int64_t *ioss_data, 1250 size_t comp_count) const; 1251 template void DecompositionDataBase::communicate_node_data(double *file_data, double *ioss_data, 1252 size_t comp_count) const; 1253 #endif 1254 1255 template <typename T> communicate_node_data(T * file_data,T * ioss_data,size_t comp_count)1256 void DecompositionDataBase::communicate_node_data(T *file_data, T *ioss_data, 1257 size_t comp_count) const 1258 { 1259 if (int_size() == sizeof(int)) { 1260 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1261 Ioss::Utils::check_dynamic_cast(this32); 1262 this32->communicate_node_data(file_data, ioss_data, comp_count); 1263 } 1264 else { 1265 const DecompositionData<int64_t> *this64 = 1266 dynamic_cast<const DecompositionData<int64_t> *>(this); 1267 Ioss::Utils::check_dynamic_cast(this64); 1268 this64->communicate_node_data(file_data, ioss_data, comp_count); 1269 } 1270 } 1271 1272 #ifndef DOXYGEN_SKIP_THIS 1273 template void DecompositionDataBase::communicate_element_data(int *file_data, int *ioss_data, 1274 size_t comp_count) const; 1275 template void DecompositionDataBase::communicate_element_data(int64_t *file_data, 1276 int64_t *ioss_data, 1277 size_t comp_count) const; 1278 template void DecompositionDataBase::communicate_element_data(double *file_data, 1279 double *ioss_data, 1280 size_t comp_count) const; 1281 #endif 1282 1283 template <typename T> communicate_element_data(T * file_data,T * ioss_data,size_t comp_count)1284 void DecompositionDataBase::communicate_element_data(T *file_data, T *ioss_data, 1285 size_t comp_count) const 1286 { 1287 if (int_size() == sizeof(int)) { 1288 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1289 Ioss::Utils::check_dynamic_cast(this32); 1290 this32->communicate_element_data(file_data, ioss_data, comp_count); 1291 } 1292 else { 1293 const DecompositionData<int64_t> *this64 = 1294 dynamic_cast<const DecompositionData<int64_t> *>(this); 1295 Ioss::Utils::check_dynamic_cast(this64); 1296 this64->communicate_element_data(file_data, ioss_data, comp_count); 1297 } 1298 } 1299 get_node_entity_proc_data(void * entity_proc,const Ioss::MapContainer & node_map,bool do_map)1300 void DecompositionDataBase::get_node_entity_proc_data(void * entity_proc, 1301 const Ioss::MapContainer &node_map, 1302 bool do_map) const 1303 { 1304 if (int_size() == sizeof(int)) { 1305 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1306 Ioss::Utils::check_dynamic_cast(this32); 1307 this32->m_decomposition.get_node_entity_proc_data((int *)entity_proc, node_map, do_map); 1308 } 1309 else { 1310 const DecompositionData<int64_t> *this64 = 1311 dynamic_cast<const DecompositionData<int64_t> *>(this); 1312 Ioss::Utils::check_dynamic_cast(this64); 1313 this64->m_decomposition.get_node_entity_proc_data((int64_t *)entity_proc, node_map, do_map); 1314 } 1315 } 1316 get_block_connectivity(int filePtr,void * data,int blk_seq,bool raw_ids)1317 void DecompositionDataBase::get_block_connectivity(int filePtr, void *data, int blk_seq, 1318 bool raw_ids) const 1319 { 1320 if (int_size() == sizeof(int)) { 1321 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1322 Ioss::Utils::check_dynamic_cast(this32); 1323 this32->get_block_connectivity(filePtr, (int *)data, blk_seq, raw_ids); 1324 } 1325 else { 1326 const DecompositionData<int64_t> *this64 = 1327 dynamic_cast<const DecompositionData<int64_t> *>(this); 1328 Ioss::Utils::check_dynamic_cast(this64); 1329 this64->get_block_connectivity(filePtr, (int64_t *)data, blk_seq, raw_ids); 1330 } 1331 } 1332 get_element_field(int filePtr,int solution_index,int blk_seq,int field_index,double * data)1333 void DecompositionDataBase::get_element_field(int filePtr, int solution_index, int blk_seq, 1334 int field_index, double *data) const 1335 { 1336 if (int_size() == sizeof(int)) { 1337 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1338 Ioss::Utils::check_dynamic_cast(this32); 1339 this32->get_element_field(filePtr, solution_index, blk_seq, field_index, data); 1340 } 1341 else { 1342 const DecompositionData<int64_t> *this64 = 1343 dynamic_cast<const DecompositionData<int64_t> *>(this); 1344 Ioss::Utils::check_dynamic_cast(this64); 1345 this64->get_element_field(filePtr, solution_index, blk_seq, field_index, data); 1346 } 1347 } 1348 get_node_field(int filePtr,int step,int field_index,double * data)1349 void DecompositionDataBase::get_node_field(int filePtr, int step, int field_index, 1350 double *data) const 1351 { 1352 if (int_size() == sizeof(int)) { 1353 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1354 Ioss::Utils::check_dynamic_cast(this32); 1355 this32->get_node_field(filePtr, step, field_index, data); 1356 } 1357 else { 1358 const DecompositionData<int64_t> *this64 = 1359 dynamic_cast<const DecompositionData<int64_t> *>(this); 1360 Ioss::Utils::check_dynamic_cast(this64); 1361 this64->get_node_field(filePtr, step, field_index, data); 1362 } 1363 } 1364 get_sideset_element_side(int filePtr,const Ioss::SetDecompositionData & sset,void * data)1365 void DecompositionDataBase::get_sideset_element_side(int filePtr, 1366 const Ioss::SetDecompositionData &sset, 1367 void * data) const 1368 { 1369 if (int_size() == sizeof(int)) { 1370 const DecompositionData<int> *this32 = dynamic_cast<const DecompositionData<int> *>(this); 1371 Ioss::Utils::check_dynamic_cast(this32); 1372 this32->get_sideset_element_side(filePtr, sset, (int *)data); 1373 } 1374 else { 1375 const DecompositionData<int64_t> *this64 = 1376 dynamic_cast<const DecompositionData<int64_t> *>(this); 1377 Ioss::Utils::check_dynamic_cast(this64); 1378 this64->get_sideset_element_side(filePtr, sset, (int64_t *)data); 1379 } 1380 } 1381 } // namespace Iocgns 1382