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 <Ioss_CodeTypes.h> 8 #include <Ioss_ElementTopology.h> 9 #include <Ioss_FileInfo.h> 10 #include <Ioss_IOFactory.h> 11 #include <Ioss_ParallelUtils.h> 12 #include <Ioss_SerializeIO.h> 13 #include <Ioss_SurfaceSplit.h> 14 #include <Ioss_Utils.h> 15 #include <algorithm> 16 #include <cassert> 17 #include <cctype> 18 #include <cfloat> 19 #include <cstddef> 20 #include <cstdio> 21 #include <cstdlib> 22 #include <cstring> 23 #include <ctime> 24 #include <exodus/Ioex_BaseDatabaseIO.h> 25 #include <exodus/Ioex_Internals.h> 26 #include <exodus/Ioex_Utils.h> 27 #include <vtk_exodusII.h> 28 #include <fmt/ostream.h> 29 #include <functional> 30 #include <iostream> 31 #include <map> 32 #include <set> 33 #include <string> 34 #include <tokenize.h> 35 #include <utility> 36 #include <vector> 37 38 #include "Ioss_Assembly.h" 39 #include "Ioss_Blob.h" 40 #include "Ioss_CommSet.h" 41 #include "Ioss_CoordinateFrame.h" 42 #include "Ioss_DBUsage.h" 43 #include "Ioss_DatabaseIO.h" 44 #include "Ioss_EdgeBlock.h" 45 #include "Ioss_EdgeSet.h" 46 #include "Ioss_ElementBlock.h" 47 #include "Ioss_ElementSet.h" 48 #include "Ioss_EntityBlock.h" 49 #include "Ioss_EntitySet.h" 50 #include "Ioss_EntityType.h" 51 #include "Ioss_FaceBlock.h" 52 #include "Ioss_FaceSet.h" 53 #include "Ioss_Field.h" 54 #include "Ioss_GroupingEntity.h" 55 #include "Ioss_Map.h" 56 #include "Ioss_NodeBlock.h" 57 #include "Ioss_NodeSet.h" 58 #include "Ioss_Property.h" 59 #include "Ioss_Region.h" 60 #include "Ioss_SideBlock.h" 61 #include "Ioss_SideSet.h" 62 #include "Ioss_SmartAssert.h" 63 #include "Ioss_State.h" 64 #include "Ioss_VariableType.h" 65 66 // Transitioning from treating global variables as Ioss::Field::TRANSIENT 67 // to Ioss::Field::REDUCTION. To get the old behavior, define the value 68 // below to '1'. 69 #define GLOBALS_ARE_TRANSIENT 0 70 71 // ======================================================================== 72 // Static internal helper functions 73 // ======================================================================== 74 namespace { 75 static bool sixty_four_bit_message_output = false; 76 77 std::vector<ex_entity_type> exodus_types({EX_GLOBAL, EX_BLOB, EX_ASSEMBLY, EX_NODE_BLOCK, 78 EX_EDGE_BLOCK, EX_FACE_BLOCK, EX_ELEM_BLOCK, 79 EX_NODE_SET, EX_EDGE_SET, EX_FACE_SET, EX_ELEM_SET, 80 EX_SIDE_SET}); 81 82 const size_t max_line_length = MAX_LINE_LENGTH; 83 84 const char *complex_suffix[] = {".re", ".im"}; 85 86 void check_variable_consistency(const ex_var_params &exo_params, int my_processor, 87 const std::string &filename, const Ioss::ParallelUtils &util); 88 89 void check_attribute_index_order(Ioss::GroupingEntity *block); 90 91 template <typename T> 92 void write_attribute_names(int exoid, ex_entity_type type, const std::vector<T *> &entities, 93 char suffix_separator); 94 95 template <typename T> 96 void generate_block_truth_table(Ioex::VariableNameMap &variables, Ioss::IntVector &truth_table, 97 std::vector<T *> &blocks, char field_suffix_separator); 98 99 } // namespace 100 101 namespace Ioex { BaseDatabaseIO(Ioss::Region * region,const std::string & filename,Ioss::DatabaseUsage db_usage,MPI_Comm communicator,const Ioss::PropertyManager & props)102 BaseDatabaseIO::BaseDatabaseIO(Ioss::Region *region, const std::string &filename, 103 Ioss::DatabaseUsage db_usage, MPI_Comm communicator, 104 const Ioss::PropertyManager &props) 105 : Ioss::DatabaseIO(region, filename, db_usage, communicator, props) 106 { 107 m_groupCount[EX_GLOBAL] = 1; // To make some common code work more cleanly. 108 m_groupCount[EX_NODE_BLOCK] = 1; // To make some common code work more cleanly. 109 110 // A history file is only written on processor 0... 111 if (db_usage == Ioss::WRITE_HISTORY) { 112 isParallel = false; 113 } 114 115 timeLastFlush = time(nullptr); 116 dbState = Ioss::STATE_UNKNOWN; 117 118 // Set exodusII warning level. 119 if (util().get_environment("EX_DEBUG", isParallel)) { 120 fmt::print( 121 Ioss::DEBUG(), 122 "IOEX: Setting EX_VERBOSE|EX_DEBUG because EX_DEBUG environment variable is set.\n"); 123 ex_opts(EX_VERBOSE | EX_DEBUG); 124 } 125 126 if (!is_input()) { 127 if (util().get_environment("EX_MODE", exodusMode, isParallel)) { 128 fmt::print( 129 Ioss::OUTPUT(), 130 "IOEX: Exodus create mode set to {} from value of EX_MODE environment variable.\n", 131 exodusMode); 132 } 133 134 if (util().get_environment("EX_MINIMIZE_OPEN_FILES", isParallel)) { 135 fmt::print(Ioss::OUTPUT(), 136 "IOEX: Minimizing open files because EX_MINIMIZE_OPEN_FILES environment " 137 "variable is set.\n"); 138 minimizeOpenFiles = true; 139 } 140 else { 141 Ioss::Utils::check_set_bool_property(properties, "MINIMIZE_OPEN_FILES", minimizeOpenFiles); 142 } 143 144 { 145 bool file_per_state = false; 146 Ioss::Utils::check_set_bool_property(properties, "FILE_PER_STATE", file_per_state); 147 if (file_per_state) { 148 set_file_per_state(true); 149 } 150 } 151 } 152 153 // See if there are any properties that need to (or can) be 154 // handled prior to opening/creating database... 155 bool compress = ((properties.exists("COMPRESSION_LEVEL") && 156 properties.get("COMPRESSION_LEVEL").get_int() > 0) || 157 (properties.exists("COMPRESSION_SHUFFLE") && 158 properties.get("COMPRESSION_SHUFFLE").get_int() > 0)); 159 160 if (compress) { 161 exodusMode |= EX_NETCDF4; 162 } 163 164 if (properties.exists("FILE_TYPE")) { 165 std::string type = properties.get("FILE_TYPE").get_string(); 166 if (type == "netcdf4" || type == "netcdf-4" || type == "hdf5") { 167 exodusMode |= EX_NETCDF4; 168 } 169 else if (type == "netcdf5" || type == "netcdf-5" || type == "cdf5") { 170 exodusMode |= EX_64BIT_DATA; 171 } 172 } 173 174 if (properties.exists("ENABLE_FILE_GROUPS")) { 175 exodusMode |= EX_NETCDF4; 176 exodusMode |= EX_NOCLASSIC; 177 } 178 179 if (properties.exists("MAXIMUM_NAME_LENGTH")) { 180 maximumNameLength = properties.get("MAXIMUM_NAME_LENGTH").get_int(); 181 } 182 183 if (properties.exists("REAL_SIZE_DB")) { 184 int rsize = properties.get("REAL_SIZE_DB").get_int(); 185 if (rsize == 4) { 186 dbRealWordSize = 4; // Only used for file create... 187 } 188 } 189 190 if (properties.exists("INTEGER_SIZE_DB")) { 191 int isize = properties.get("INTEGER_SIZE_DB").get_int(); 192 if (isize == 8) { 193 exodusMode |= EX_ALL_INT64_DB; 194 } 195 } 196 197 if (properties.exists("INTEGER_SIZE_API")) { 198 int isize = properties.get("INTEGER_SIZE_API").get_int(); 199 if (isize == 8) { 200 set_int_byte_size_api(Ioss::USE_INT64_API); 201 } 202 } 203 204 if (!is_input()) { 205 if (properties.exists("FLUSH_INTERVAL")) { 206 int interval = properties.get("FLUSH_INTERVAL").get_int(); 207 flushInterval = interval; 208 } 209 } 210 211 // Don't open output files until they are actually going to be 212 // written to. This is needed for proper support of the topology 213 // files and auto restart so we don't overwrite a file with data we 214 // need to save... 215 } 216 set_int_byte_size_api(Ioss::DataSize size)217 void BaseDatabaseIO::set_int_byte_size_api(Ioss::DataSize size) const 218 { 219 if (m_exodusFilePtr > 0) { 220 int old_status = ex_int64_status(get_file_pointer()); 221 if (size == 8) { 222 ex_set_int64_status(get_file_pointer(), EX_ALL_INT64_API | old_status); 223 } 224 else { 225 // Need to clear EX_ALL_INT64_API if set... 226 if ((old_status & EX_ALL_INT64_API) != 0) { 227 old_status &= ~EX_ALL_INT64_API; 228 assert(!(old_status & EX_ALL_INT64_API)); 229 ex_set_int64_status(m_exodusFilePtr, old_status); 230 } 231 } 232 } 233 else { 234 if (size == 8) { 235 exodusMode |= EX_ALL_INT64_API; 236 } 237 else { 238 exodusMode &= ~EX_ALL_INT64_API; 239 } 240 } 241 dbIntSizeAPI = size; // mutable 242 } 243 244 // Returns byte size of integers stored on the database... int_byte_size_db()245 int BaseDatabaseIO::int_byte_size_db() const 246 { 247 int status = ex_int64_status(get_file_pointer()); 248 if (status & EX_MAPS_INT64_DB || status & EX_IDS_INT64_DB || status & EX_BULK_INT64_DB) { 249 return 8; 250 } 251 else { 252 return 4; 253 } 254 } 255 256 // common ~BaseDatabaseIO()257 BaseDatabaseIO::~BaseDatabaseIO() 258 { 259 try { 260 free_file_pointer(); 261 } 262 catch (...) { 263 } 264 } 265 266 // common entity_field_support()267 unsigned BaseDatabaseIO::entity_field_support() const 268 { 269 return Ioss::NODEBLOCK | Ioss::EDGEBLOCK | Ioss::FACEBLOCK | Ioss::ELEMENTBLOCK | 270 Ioss::NODESET | Ioss::EDGESET | Ioss::FACESET | Ioss::ELEMENTSET | Ioss::SIDESET | 271 Ioss::SIDEBLOCK | Ioss::REGION | Ioss::SUPERELEMENT; 272 } 273 274 // common get_file_pointer()275 int BaseDatabaseIO::get_file_pointer() const 276 { 277 // Returns the file_pointer used to access the file on disk. 278 // Checks that the file is open and if not, opens it first. 279 if (m_exodusFilePtr < 0) { 280 bool write_message = true; 281 bool abort_if_error = true; 282 if (is_input()) { 283 open_input_file(write_message, nullptr, nullptr, abort_if_error); 284 } 285 else { 286 bool overwrite = true; 287 handle_output_file(write_message, nullptr, nullptr, overwrite, abort_if_error); 288 } 289 290 if (!m_groupName.empty()) { 291 ex_get_group_id(m_exodusFilePtr, m_groupName.c_str(), &m_exodusFilePtr); 292 } 293 } 294 assert(m_exodusFilePtr >= 0); 295 fileExists = true; 296 return m_exodusFilePtr; 297 } 298 free_file_pointer()299 int BaseDatabaseIO::free_file_pointer() const 300 { 301 if (m_exodusFilePtr != -1) { 302 bool do_timer = false; 303 if (isParallel) { 304 Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer); 305 } 306 double t_begin = (do_timer ? Ioss::Utils::timer() : 0); 307 308 ex_close(m_exodusFilePtr); 309 closeDW(); 310 if (do_timer && isParallel) { 311 double t_end = Ioss::Utils::timer(); 312 double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX); 313 if (myProcessor == 0) { 314 fmt::print(Ioss::DEBUG(), "File Close Time = {}\n", duration); 315 } 316 } 317 } 318 m_exodusFilePtr = -1; 319 320 return m_exodusFilePtr; 321 } 322 ok__(bool write_message,std::string * error_message,int * bad_count)323 bool BaseDatabaseIO::ok__(bool write_message, std::string *error_message, int *bad_count) const 324 { 325 // For input, we try to open the existing file. 326 327 // For output, we do not want to overwrite or clobber the output 328 // file if it already exists since the app might be reading the restart 329 // data from this file and then later clobbering it and then writing 330 // restart data to the same file. So, for output, we first check 331 // whether the file exists and if it it and is writable, assume 332 // that we can later create a new or append to existing file. 333 334 // Returns the number of processors on which this file is *NOT* ok in 'bad_count' if not null. 335 // Will return 'true' only if file ok on all processors. 336 337 if (fileExists) { 338 // File has already been opened at least once... 339 return dbState != Ioss::STATE_INVALID; 340 } 341 342 bool abort_if_error = false; 343 bool is_ok; 344 if (is_input()) { 345 is_ok = open_input_file(write_message, error_message, bad_count, abort_if_error); 346 } 347 else { 348 // See if file exists... Don't overwrite (yet) it it exists. 349 bool overwrite = false; 350 is_ok = 351 handle_output_file(write_message, error_message, bad_count, overwrite, abort_if_error); 352 // Close all open files... 353 if (m_exodusFilePtr >= 0) { 354 ex_close(m_exodusFilePtr); 355 m_exodusFilePtr = -1; 356 } 357 } 358 return is_ok; 359 } 360 finalize_file_open()361 void BaseDatabaseIO::finalize_file_open() const 362 { 363 assert(m_exodusFilePtr >= 0); 364 // Check byte-size of integers stored on the database... 365 if ((ex_int64_status(m_exodusFilePtr) & EX_ALL_INT64_DB) != 0) { 366 if (myProcessor == 0 && !sixty_four_bit_message_output) { 367 fmt::print(Ioss::OUTPUT(), 368 "IOSS: Input database contains 8-byte integers. Setting Ioss to use " 369 "8-byte integers.\n"); 370 sixty_four_bit_message_output = true; 371 } 372 ex_set_int64_status(m_exodusFilePtr, EX_ALL_INT64_API); 373 set_int_byte_size_api(Ioss::USE_INT64_API); 374 } 375 376 // Check for maximum name length used on the input file. 377 int max_name_length = ex_inquire_int(m_exodusFilePtr, EX_INQ_DB_MAX_USED_NAME_LENGTH); 378 if (max_name_length > maximumNameLength) { 379 maximumNameLength = max_name_length; 380 } 381 382 ex_set_max_name_length(m_exodusFilePtr, maximumNameLength); 383 } 384 open_group__(const std::string & group_name)385 bool BaseDatabaseIO::open_group__(const std::string &group_name) 386 { 387 // Get existing file pointer... 388 bool success = false; 389 390 int exoid = get_file_pointer(); 391 392 m_groupName = group_name; 393 ex_get_group_id(exoid, m_groupName.c_str(), &m_exodusFilePtr); 394 395 if (m_exodusFilePtr < 0) { 396 std::ostringstream errmsg; 397 fmt::print(errmsg, "ERROR: Could not open group named '{}' in file '{}'.\n", m_groupName, 398 get_filename()); 399 IOSS_ERROR(errmsg); 400 } 401 success = true; 402 return success; 403 } 404 create_subgroup__(const std::string & group_name)405 bool BaseDatabaseIO::create_subgroup__(const std::string &group_name) 406 { 407 bool success = false; 408 if (!is_input()) { 409 // Get existing file pointer... 410 int exoid = get_file_pointer(); 411 412 // Check name for '/' which is not allowed since it is the 413 // separator character in a full group path 414 if (group_name.find('/') != std::string::npos) { 415 std::ostringstream errmsg; 416 fmt::print(errmsg, "ERROR: Invalid group name '{}' contains a '/' which is not allowed.\n", 417 m_groupName); 418 IOSS_ERROR(errmsg); 419 } 420 421 m_groupName = group_name; 422 exoid = ex_create_group(exoid, m_groupName.c_str()); 423 if (exoid < 0) { 424 std::ostringstream errmsg; 425 fmt::print(errmsg, "ERROR: Could not create group named '{}' in file '{}'.\n", m_groupName, 426 get_filename()); 427 IOSS_ERROR(errmsg); 428 } 429 m_exodusFilePtr = exoid; 430 success = true; 431 } 432 return success; 433 } 434 435 // common put_qa()436 void BaseDatabaseIO::put_qa() 437 { 438 struct qa_element 439 { 440 char *qa_record[1][4]; 441 }; 442 443 size_t num_qa_records = qaRecords.size() / 4; 444 445 auto *qa = new qa_element[num_qa_records + 1]; 446 for (size_t i = 0; i < num_qa_records + 1; i++) { 447 for (int j = 0; j < 4; j++) { 448 qa[i].qa_record[0][j] = new char[MAX_STR_LENGTH + 1]; 449 } 450 } 451 452 { 453 int j = 0; 454 for (size_t i = 0; i < num_qa_records; i++) { 455 Ioss::Utils::copy_string(qa[i].qa_record[0][0], qaRecords[j++], MAX_STR_LENGTH + 1); 456 Ioss::Utils::copy_string(qa[i].qa_record[0][1], qaRecords[j++], MAX_STR_LENGTH + 1); 457 Ioss::Utils::copy_string(qa[i].qa_record[0][2], qaRecords[j++], MAX_STR_LENGTH + 1); 458 Ioss::Utils::copy_string(qa[i].qa_record[0][3], qaRecords[j++], MAX_STR_LENGTH + 1); 459 } 460 } 461 462 Ioss::Utils::time_and_date(qa[num_qa_records].qa_record[0][3], 463 qa[num_qa_records].qa_record[0][2], MAX_STR_LENGTH); 464 465 std::string codename = "unknown"; 466 std::string version = "unknown"; 467 468 if (get_region()->property_exists("code_name")) { 469 codename = get_region()->get_property("code_name").get_string(); 470 } 471 if (get_region()->property_exists("code_version")) { 472 version = get_region()->get_property("code_version").get_string(); 473 } 474 475 Ioss::Utils::copy_string(qa[num_qa_records].qa_record[0][0], codename, MAX_STR_LENGTH + 1); 476 Ioss::Utils::copy_string(qa[num_qa_records].qa_record[0][1], version, MAX_STR_LENGTH + 1); 477 478 int ierr = ex_put_qa(get_file_pointer(), num_qa_records + 1, qa[0].qa_record); 479 if (ierr < 0) { 480 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 481 } 482 483 for (size_t i = 0; i < num_qa_records + 1; i++) { 484 for (int j = 0; j < 4; j++) { 485 delete[] qa[i].qa_record[0][j]; 486 } 487 } 488 delete[] qa; 489 } 490 491 // common put_info()492 void BaseDatabaseIO::put_info() 493 { 494 // dump info records, include the product_registry 495 // See if the input file was specified as a property on the database... 496 std::string filename; 497 std::vector<std::string> input_lines; 498 if (get_region()->property_exists("input_file_name")) { 499 filename = get_region()->get_property("input_file_name").get_string(); 500 // Determine size of input file so can embed it in info records... 501 Ioss::Utils::input_file(filename, &input_lines, max_line_length); 502 } 503 504 // Get configuration information for IOSS library. 505 // Split into strings and remove empty lines... 506 std::string config = Ioss::IOFactory::show_configuration(); 507 std::replace(std::begin(config), std::end(config), '\t', ' '); 508 auto lines = Ioss::tokenize(config, "\n"); 509 lines.erase(std::remove_if(lines.begin(), lines.end(), 510 [](const std::string &line) { return line == ""; }), 511 lines.end()); 512 513 // See if the client added any "information_records" 514 size_t info_rec_size = informationRecords.size(); 515 size_t in_lines = input_lines.size(); 516 size_t qa_lines = 2; // Platform info and Version info... 517 size_t config_lines = lines.size(); 518 519 size_t total_lines = in_lines + qa_lines + info_rec_size + config_lines; 520 521 char **info = Ioss::Utils::get_name_array( 522 total_lines, max_line_length); // 'total_lines' pointers to char buffers 523 524 int i = 0; 525 Ioss::Utils::copy_string(info[i++], Ioss::Utils::platform_information(), max_line_length + 1); 526 527 Ioss::Utils::copy_string(info[i++], Ioex::Version(), max_line_length + 1); 528 529 // Copy input file lines into 'info' array... 530 for (size_t j = 0; j < input_lines.size(); j++, i++) { 531 Ioss::Utils::copy_string(info[i], input_lines[j], max_line_length + 1); 532 } 533 534 // Copy "information_records" property data ... 535 for (size_t j = 0; j < informationRecords.size(); j++, i++) { 536 Ioss::Utils::copy_string(info[i], informationRecords[j], max_line_length + 1); 537 } 538 539 for (size_t j = 0; j < lines.size(); j++, i++) { 540 Ioss::Utils::copy_string(info[i], lines[j], max_line_length + 1); 541 } 542 543 int ierr = ex_put_info(get_file_pointer(), total_lines, info); 544 if (ierr < 0) { 545 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 546 } 547 548 Ioss::Utils::delete_name_array(info, total_lines); 549 } 550 551 // common get_current_state()552 int BaseDatabaseIO::get_current_state() const 553 { 554 int step = get_region()->get_current_state(); 555 556 if (step <= 0) { 557 std::ostringstream errmsg; 558 fmt::print(errmsg, 559 "ERROR: No currently active state. The calling code must call " 560 "Ioss::Region::begin_state(int step)\n" 561 " to set the database timestep from which to read the transient data.\n" 562 " [{}]\n", 563 get_filename()); 564 IOSS_ERROR(errmsg); 565 } 566 return step; 567 } 568 get_assemblies()569 void BaseDatabaseIO::get_assemblies() 570 { 571 Ioss::SerializeIO serializeIO__(this); 572 // Query number of coordinate frames... 573 int nassem = ex_inquire_int(get_file_pointer(), EX_INQ_ASSEMBLY); 574 575 if (nassem > 0) { 576 std::vector<ex_assembly> assemblies(nassem); 577 int max_name_length = ex_inquire_int(m_exodusFilePtr, EX_INQ_DB_MAX_USED_NAME_LENGTH); 578 for (auto &assembly : assemblies) { 579 assembly.name = new char[max_name_length + 1]; 580 } 581 582 int ierr = ex_get_assemblies(get_file_pointer(), assemblies.data()); 583 if (ierr < 0) { 584 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 585 } 586 587 // Now allocate space for member list and get assemblies again... 588 for (auto &assembly : assemblies) { 589 assembly.entity_list = new int64_t[assembly.entity_count]; 590 } 591 592 ierr = ex_get_assemblies(get_file_pointer(), assemblies.data()); 593 if (ierr < 0) { 594 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 595 } 596 597 for (const auto &assembly : assemblies) { 598 auto *assem = new Ioss::Assembly(get_region()->get_database(), assembly.name); 599 assem->property_add(Ioss::Property("id", assembly.id)); 600 get_region()->add(assem); 601 } 602 603 // Now iterate again and populate member lists... 604 for (const auto &assembly : assemblies) { 605 Ioss::Assembly *assem = get_region()->get_assembly(assembly.name); 606 assert(assem != nullptr); 607 auto type = Ioex::map_exodus_type(assembly.type); 608 for (int j = 0; j < assembly.entity_count; j++) { 609 auto *ge = get_region()->get_entity(assembly.entity_list[j], type); 610 if (ge != nullptr) { 611 assem->add(ge); 612 } 613 else { 614 std::ostringstream errmsg; 615 fmt::print(errmsg, 616 "Error: Failed to find entity of type {} with id {} for Assembly {}.\n", 617 type, assembly.entity_list[j], assem->name()); 618 IOSS_ERROR(errmsg); 619 } 620 } 621 SMART_ASSERT(assem->member_count() == (size_t)assembly.entity_count) 622 (assem->member_count())(assembly.entity_count); 623 624 add_mesh_reduction_fields(EX_ASSEMBLY, assembly.id, assem); 625 // Check for additional variables. 626 int attribute_count = assem->get_property("attribute_count").get_int(); 627 add_attribute_fields(EX_ASSEMBLY, assem, attribute_count, "Assembly"); 628 add_reduction_results_fields(EX_ASSEMBLY, assem); 629 } 630 631 // If there are any reduction results fields ("REDUCTION"), then need to 632 // allocate space for the values to be stored on each timestep... 633 if (!m_reductionVariables[EX_ASSEMBLY].empty()) { 634 size_t size = m_reductionVariables[EX_ASSEMBLY].size(); 635 for (const auto &assembly : assemblies) { 636 m_reductionValues[EX_ASSEMBLY][assembly.id].resize(size); 637 } 638 } 639 for (auto &assembly : assemblies) { 640 delete[] assembly.entity_list; 641 delete[] assembly.name; 642 } 643 } 644 } 645 get_blobs()646 void BaseDatabaseIO::get_blobs() 647 { 648 Ioss::SerializeIO serializeIO__(this); 649 // Query number of coordinate frames... 650 int nblob = ex_inquire_int(get_file_pointer(), EX_INQ_BLOB); 651 652 if (nblob > 0) { 653 std::vector<ex_blob> blobs(nblob); 654 int max_name_length = ex_inquire_int(m_exodusFilePtr, EX_INQ_DB_MAX_USED_NAME_LENGTH); 655 for (auto &bl : blobs) { 656 bl.name = new char[max_name_length + 1]; 657 } 658 659 int ierr = ex_get_blobs(get_file_pointer(), blobs.data()); 660 if (ierr < 0) { 661 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 662 } 663 664 for (const auto &bl : blobs) { 665 #ifdef SEACAS_HAVE_MPI 666 // Each blob is spread across all processors (should support a minimum size...) 667 // Determine size of blob on each rank and offset from beginning of blob. 668 size_t per_proc = bl.num_entry / parallel_size(); 669 size_t extra = bl.num_entry % parallel_size(); 670 size_t count = per_proc + (myProcessor < (int)extra ? 1 : 0); 671 672 size_t offset = 0; 673 if (myProcessor < (int)extra) { 674 offset = (per_proc + 1) * myProcessor; 675 } 676 else { 677 offset = (per_proc + 1) * extra + per_proc * (myProcessor - extra); 678 } 679 Ioss::Blob *blob = new Ioss::Blob(get_region()->get_database(), bl.name, count); 680 blob->property_add(Ioss::Property("_processor_offset", (int64_t)offset)); 681 blob->property_add(Ioss::Property("global_size", (int64_t)bl.num_entry)); 682 #else 683 auto *blob = new Ioss::Blob(get_region()->get_database(), bl.name, bl.num_entry); 684 #endif 685 blob->property_add(Ioss::Property("id", bl.id)); 686 get_region()->add(blob); 687 } 688 689 // Now iterate again and populate member lists... 690 int iblk = 0; 691 for (const auto &bl : blobs) { 692 Ioss::Blob *blob = get_region()->get_blob(bl.name); 693 assert(blob != nullptr); 694 695 add_mesh_reduction_fields(EX_BLOB, bl.id, blob); 696 // Check for additional variables. 697 int attribute_count = blob->get_property("attribute_count").get_int(); 698 add_attribute_fields(EX_BLOB, blob, attribute_count, "Blob"); 699 add_reduction_results_fields(EX_BLOB, blob); 700 add_results_fields(EX_BLOB, blob, iblk++); 701 } 702 703 // If there are any reduction results fields ("REDUCTION"), then need to 704 // allocate space for the values to be stored on each timestep... 705 if (!m_reductionVariables[EX_BLOB].empty()) { 706 size_t size = m_reductionVariables[EX_BLOB].size(); 707 for (const auto &blob : blobs) { 708 m_reductionValues[EX_BLOB][blob.id].resize(size); 709 } 710 } 711 712 for (auto &bl : blobs) { 713 delete[] bl.name; 714 } 715 } 716 } 717 718 // common get_nodeblocks()719 void BaseDatabaseIO::get_nodeblocks() 720 { 721 // For exodusII, there is only a single node block which contains 722 // all of the nodes. 723 // The default id assigned is '1' and the name is 'nodeblock_1' 724 725 std::string block_name = "nodeblock_1"; 726 auto block = new Ioss::NodeBlock(this, block_name, nodeCount, spatialDimension); 727 block->property_add(Ioss::Property("id", 1)); 728 block->property_add(Ioss::Property("guid", util().generate_guid(1))); 729 // Check for results variables. 730 731 int num_attr = 0; 732 { 733 Ioss::SerializeIO serializeIO__(this); 734 int ierr = ex_get_attr_param(get_file_pointer(), EX_NODE_BLOCK, 1, &num_attr); 735 if (ierr < 0) { 736 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 737 } 738 } 739 740 add_attribute_fields(EX_NODE_BLOCK, block, num_attr, ""); 741 // Not supported on nodeblocks at this time 742 // add_reduction_results_fields(EX_NODE_BLOCK, block); 743 add_results_fields(EX_NODE_BLOCK, block); 744 745 // If there are any reduction results fields ("REDUCTION"), then need to 746 // allocate space for the values to be stored on each timestep... 747 if (!m_reductionVariables[EX_NODE_BLOCK].empty()) { 748 size_t size = m_reductionVariables[EX_NODE_BLOCK].size(); 749 m_reductionValues[EX_NODE_BLOCK][1].resize(size); 750 } 751 752 bool added = get_region()->add(block); 753 if (!added) { 754 delete block; 755 } 756 } 757 758 // common 759 // common handle_block_ids(const Ioss::EntityBlock * eb,ex_entity_type map_type,Ioss::Map & entity_map,void * ids,size_t num_to_get,size_t offset)760 size_t BaseDatabaseIO::handle_block_ids(const Ioss::EntityBlock *eb, ex_entity_type map_type, 761 Ioss::Map &entity_map, void *ids, size_t num_to_get, 762 size_t offset) const 763 { 764 /*! 765 * NOTE: "element" is generic for "element", "face", or "edge" 766 * 767 * There are two modes we need to support in this routine: 768 * 1. Initial definition of element map (local->global) and 769 * elemMap.reverse (global->local). 770 * 2. Redefinition of element map via 'reordering' of the original 771 * map when the elements on this processor are the same, but their 772 * order is changed. 773 * 774 * So, there will be two maps the 'elemMap.map' map is a 'direct lookup' 775 * map which maps current local position to global id and the 776 * 'elemMap.reverse' is an associative lookup which maps the 777 * global id to 'original local'. There is also a 778 * 'elemMap.reorder' which is direct lookup and maps current local 779 * position to original local. 780 781 * The ids coming in are the global ids; their position is the 782 * local id -1 (That is, data[0] contains the global id of local 783 * element 1 in this element block). The 'model-local' id is 784 * given by eb_offset + 1 + position: 785 * 786 * int local_position = elemMap.reverse[ElementMap[i+1]] 787 * (the elemMap.map and elemMap.reverse are 1-based) 788 * 789 * But, this assumes 1..numel elements are being output at the same 790 * time; we are actually outputting a blocks worth of elements at a 791 * time, so we need to consider the block offsets. 792 * So... local-in-block position 'i' is index 'eb_offset+i' in 793 * 'elemMap.map' and the 'local_position' within the element 794 * blocks data arrays is 'local_position-eb_offset'. With this, the 795 * position within the data array of this element block is: 796 * 797 * int eb_position = 798 * elemMap.reverse[elemMap.map[eb_offset+i+1]]-eb_offset-1 799 * 800 * To determine which map to update on a call to this function, we 801 * use the following hueristics: 802 * -- If the database state is 'Ioss::STATE_MODEL:', then update the 803 * 'elemMap.reverse'. 804 * -- If the database state is not Ioss::STATE_MODEL, then leave 805 * the 'elemMap.reverse' alone since it corresponds to the 806 * information already written to the database. [May want to add 807 * a Ioss::STATE_REDEFINE_MODEL] 808 * -- Always update elemMap.map to match the passed in 'ids' 809 * array. 810 * 811 * NOTE: the maps are built an element block at a time... 812 * NOTE: The mapping is done on TRANSIENT fields only; MODEL fields 813 * should be in the original order... 814 */ 815 816 // Overwrite this portion of the 'elemMap.map', but keep other 817 // parts as they were. We are adding elements starting at position 818 // 'eb_offset+offset' and ending at 819 // 'eb_offset+offset+num_to_get'. If the entire block is being 820 // processed, this reduces to the range 'eb_offset..eb_offset+my_element_count' 821 822 bool in_define = (dbState == Ioss::STATE_MODEL) || (dbState == Ioss::STATE_DEFINE_MODEL); 823 int64_t eb_offset = eb->get_offset(); 824 if (int_byte_size_api() == 4) { 825 entity_map.set_map(static_cast<int *>(ids), num_to_get, eb_offset, in_define); 826 } 827 else { 828 entity_map.set_map(static_cast<int64_t *>(ids), num_to_get, eb_offset, in_define); 829 } 830 831 // Now, if the state is Ioss::STATE_MODEL, output this portion of 832 // the entity number map... 833 if (in_define) { 834 if (ex_put_partial_id_map(get_file_pointer(), map_type, offset + 1, num_to_get, ids) < 0) { 835 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 836 } 837 } 838 return num_to_get; 839 } 840 841 // common compute_block_membership__(Ioss::SideBlock * efblock,std::vector<std::string> & block_membership)842 void BaseDatabaseIO::compute_block_membership__(Ioss::SideBlock * efblock, 843 std::vector<std::string> &block_membership) const 844 { 845 const Ioss::ElementBlockContainer &element_blocks = get_region()->get_element_blocks(); 846 assert(Ioss::Utils::check_block_order(element_blocks)); 847 848 Ioss::Int64Vector block_ids(element_blocks.size()); 849 if (block_ids.size() == 1) { 850 block_ids[0] = 1; 851 } 852 else { 853 Ioss::Int64Vector element_side; 854 if (int_byte_size_api() == 4) { 855 Ioss::IntVector es32; 856 efblock->get_field_data("element_side", es32); 857 element_side.resize(es32.size()); 858 std::copy(es32.begin(), es32.end(), element_side.begin()); 859 } 860 else { 861 efblock->get_field_data("element_side", element_side); 862 } 863 864 size_t number_sides = element_side.size() / 2; 865 Ioss::ElementBlock *block = nullptr; 866 for (size_t iel = 0; iel < number_sides; iel++) { 867 int64_t elem_id = element_side[2 * iel]; // Vector contains both element and side. 868 elem_id = elemMap.global_to_local(elem_id); 869 if (block == nullptr || !block->contains(elem_id)) { 870 block = get_region()->get_element_block(elem_id); 871 assert(block != nullptr); 872 size_t block_order = block->get_property("original_block_order").get_int(); 873 assert(block_order < block_ids.size()); 874 block_ids[block_order] = 1; 875 } 876 } 877 } 878 879 // Synchronize among all processors.... 880 if (isParallel) { 881 util().global_array_minmax(block_ids, Ioss::ParallelUtils::DO_MAX); 882 } 883 884 for (const auto block : element_blocks) { 885 size_t block_order = block->get_property("original_block_order").get_int(); 886 assert(block_order < block_ids.size()); 887 if (block_ids[block_order] == 1) { 888 if (!Ioss::Utils::block_is_omitted(block)) { 889 block_membership.push_back(block->name()); 890 } 891 } 892 } 893 } 894 895 // common get_field_internal(const Ioss::Region *,const Ioss::Field & field,void * data,size_t data_size)896 int64_t BaseDatabaseIO::get_field_internal(const Ioss::Region * /* region */, 897 const Ioss::Field &field, void *data, 898 size_t data_size) const 899 { 900 // For now, assume that all TRANSIENT fields on a region 901 // are REDUCTION fields (1 value). We need to gather these 902 // and output them all at one time. The storage location is a 903 // 'globalVariables' array 904 { 905 size_t num_to_get = field.verify(data_size); 906 Ioss::SerializeIO serializeIO__(this); 907 908 Ioss::Field::RoleType role = field.get_role(); 909 910 if (role == Ioss::Field::TRANSIENT || role == Ioss::Field::REDUCTION) { 911 get_reduction_field(EX_GLOBAL, field, get_region(), data); 912 } 913 else { 914 std::ostringstream errmsg; 915 fmt::print(errmsg, 916 "ERROR: Can not handle non-TRANSIENT or non-REDUCTION fields on regions"); 917 IOSS_ERROR(errmsg); 918 } 919 return num_to_get; 920 } 921 } 922 923 // common put_field_internal(const Ioss::Region *,const Ioss::Field & field,void * data,size_t data_size)924 int64_t BaseDatabaseIO::put_field_internal(const Ioss::Region * /* region */, 925 const Ioss::Field &field, void *data, 926 size_t data_size) const 927 { 928 // For now, assume that all TRANSIENT fields on a region 929 // are REDUCTION fields (1 value). We need to gather these 930 // and output them all at one time. The storage location is a 931 // 'globalVariables' array 932 { 933 Ioss::SerializeIO serializeIO__(this); 934 935 Ioss::Field::RoleType role = field.get_role(); 936 size_t num_to_get = field.verify(data_size); 937 938 if ((role == Ioss::Field::TRANSIENT || role == Ioss::Field::REDUCTION) && num_to_get == 1) { 939 store_reduction_field(EX_GLOBAL, field, get_region(), data); 940 } 941 else if (num_to_get != 1) { 942 // There should have been a warning/error message printed to the 943 // log file earlier for this, so we won't print anything else 944 // here since it would be printed for each and every timestep.... 945 ; 946 } 947 else { 948 std::ostringstream errmsg; 949 fmt::print( 950 errmsg, 951 "ERROR: The variable named '{}' is of the wrong type. A region variable must be of type" 952 " TRANSIENT or REDUCTION.\n" 953 "This is probably an internal error; please notify gdsjaar@sandia.gov", 954 field.get_name()); 955 IOSS_ERROR(errmsg); 956 } 957 return num_to_get; 958 } 959 } 960 961 namespace { 962 // common 963 template <typename T> generate_block_truth_table(const VariableNameMap & variables,Ioss::IntVector & truth_table,std::vector<T * > & blocks,char field_suffix_separator)964 void generate_block_truth_table(const VariableNameMap &variables, Ioss::IntVector &truth_table, 965 std::vector<T *> &blocks, char field_suffix_separator) 966 { 967 size_t block_count = blocks.size(); 968 size_t var_count = variables.size(); 969 970 if (var_count == 0 || block_count == 0) { 971 return; 972 } 973 974 truth_table.resize(block_count * var_count); 975 976 // Fill in the truth table. It is conceptually a two-dimensional array of 977 // the form 'array[num_blocks][num_element_var]'. In C++, 978 // the values for the first block are first, followed by 979 // next block, ... 980 size_t offset = 0; 981 for (const auto &block : blocks) { 982 // Get names of all transient and reduction fields... 983 Ioss::NameList results_fields; 984 block->field_describe(Ioss::Field::TRANSIENT, &results_fields); 985 block->field_describe(Ioss::Field::REDUCTION, &results_fields); 986 987 for (const auto &fn : results_fields) { 988 Ioss::Field field = block->get_field(fn); 989 const Ioss::VariableType *var_type = field.transformed_storage(); 990 Ioss::Field::BasicType ioss_type = field.get_type(); 991 992 int re_im = 1; 993 if (ioss_type == Ioss::Field::COMPLEX) { 994 re_im = 2; 995 } 996 for (int complex_comp = 0; complex_comp < re_im; complex_comp++) { 997 std::string field_name = field.get_name(); 998 if (re_im == 2) { 999 field_name += complex_suffix[complex_comp]; 1000 } 1001 1002 for (int i = 1; i <= var_type->component_count(); i++) { 1003 std::string var_string = var_type->label_name(field_name, i, field_suffix_separator); 1004 // Find position of 'var_string' in 'variables' 1005 auto VN = variables.find(var_string); 1006 if (VN != variables.end()) { 1007 // Index '(*VN).second' is 1-based... 1008 truth_table[offset + (*VN).second - 1] = 1; 1009 } 1010 } 1011 } 1012 } 1013 offset += var_count; 1014 } 1015 assert(offset == var_count * block_count); 1016 } 1017 } // namespace 1018 // common store_reduction_field(ex_entity_type type,const Ioss::Field & field,const Ioss::GroupingEntity * ge,void * variables)1019 void BaseDatabaseIO::store_reduction_field(ex_entity_type type, const Ioss::Field &field, 1020 const Ioss::GroupingEntity *ge, void *variables) const 1021 { 1022 const Ioss::VariableType *var_type = field.transformed_storage(); 1023 1024 Ioss::Field::BasicType ioss_type = field.get_type(); 1025 assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER || 1026 ioss_type == Ioss::Field::INT64 || ioss_type == Ioss::Field::COMPLEX); 1027 auto *rvar = static_cast<double *>(variables); 1028 auto *ivar = static_cast<int *>(variables); 1029 auto *ivar64 = static_cast<int64_t *>(variables); 1030 1031 auto id = ge->get_optional_property("id", 0); 1032 1033 // Note that if the field's basic type is COMPLEX, then each component of 1034 // the VariableType is a complex variable consisting of a real and 1035 // imaginary part. Since exodus cannot handle complex variables, 1036 // we have to output a (real and imaginary) X (number of 1037 // components) fields. For example, if V is a 3d vector of complex 1038 // data, the data in the 'variables' array are v_x, v.im_x, v_y, 1039 // v.im_y, v_z, v.im_z which need to be output in six separate 1040 // exodus fields. These fields were already defined in 1041 // "write_results_metadata". 1042 1043 // get number of components, cycle through each component 1044 // and add suffix to base 'field_name'. Look up index 1045 // of this name in 'm_variables[EX_GLOBAL]' map 1046 int comp_count = var_type->component_count(); 1047 int var_index = 0; 1048 1049 int re_im = 1; 1050 if (field.get_type() == Ioss::Field::COMPLEX) { 1051 re_im = 2; 1052 } 1053 for (int complex_comp = 0; complex_comp < re_im; complex_comp++) { 1054 std::string field_name = field.get_name(); 1055 if (re_im == 2) { 1056 field_name += complex_suffix[complex_comp]; 1057 } 1058 1059 char field_suffix_separator = get_field_separator(); 1060 for (int i = 0; i < comp_count; i++) { 1061 std::string var_name = var_type->label_name(field_name, i + 1, field_suffix_separator); 1062 1063 #if GLOBALS_ARE_TRANSIENT 1064 if (type == EX_GLOBAL) { 1065 SMART_ASSERT(m_variables[type].find(var_name) != m_variables[type].end())(type)(var_name); 1066 var_index = m_variables[type].find(var_name)->second; 1067 } 1068 else { 1069 SMART_ASSERT(m_reductionVariables[type].find(var_name) != 1070 m_reductionVariables[type].end()) 1071 (type)(var_name); 1072 var_index = m_reductionVariables[type].find(var_name)->second; 1073 } 1074 #else 1075 SMART_ASSERT(m_reductionVariables[type].find(var_name) != m_reductionVariables[type].end()) 1076 (type)(var_name); 1077 var_index = m_reductionVariables[type].find(var_name)->second; 1078 #endif 1079 1080 SMART_ASSERT(static_cast<int>(m_reductionValues[type][id].size()) >= var_index) 1081 (id)(m_reductionValues[type][id].size())(var_index); 1082 1083 // Transfer from 'variables' array. 1084 if (ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::COMPLEX) { 1085 m_reductionValues[type][id][var_index - 1] = rvar[i]; 1086 } 1087 else if (ioss_type == Ioss::Field::INTEGER) { 1088 m_reductionValues[type][id][var_index - 1] = ivar[i]; 1089 } 1090 else if (ioss_type == Ioss::Field::INT64) { 1091 m_reductionValues[type][id][var_index - 1] = ivar64[i]; // FIX 64 UNSAFE 1092 } 1093 } 1094 } 1095 } 1096 1097 // common get_reduction_field(ex_entity_type type,const Ioss::Field & field,const Ioss::GroupingEntity * ge,void * variables)1098 void BaseDatabaseIO::get_reduction_field(ex_entity_type type, const Ioss::Field &field, 1099 const Ioss::GroupingEntity *ge, void *variables) const 1100 { 1101 const Ioss::VariableType *var_type = field.raw_storage(); 1102 1103 auto id = ge->get_optional_property("id", 0); 1104 1105 Ioss::Field::BasicType ioss_type = field.get_type(); 1106 assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER || 1107 ioss_type == Ioss::Field::INT64); 1108 auto *rvar = static_cast<double *>(variables); 1109 auto *ivar = static_cast<int *>(variables); 1110 auto *i64var = static_cast<int64_t *>(variables); 1111 1112 // get number of components, cycle through each component 1113 // and add suffix to base 'field_name'. Look up index 1114 // of this name in 'm_variables[type]' map 1115 char field_suffix_separator = get_field_separator(); 1116 1117 int comp_count = var_type->component_count(); 1118 for (int i = 0; i < comp_count; i++) { 1119 int var_index = 0; 1120 const std::string &field_name = field.get_name(); 1121 std::string var_name = var_type->label_name(field_name, i + 1, field_suffix_separator); 1122 1123 #if GLOBALS_ARE_TRANSIENT 1124 if (type == EX_GLOBAL) { 1125 assert(m_variables[type].find(var_name) != m_variables[type].end()); 1126 var_index = m_variables[type].find(var_name)->second; 1127 } 1128 else { 1129 assert(m_reductionVariables[type].find(var_name) != m_reductionVariables[type].end()); 1130 var_index = m_reductionVariables[type].find(var_name)->second; 1131 } 1132 1133 assert(static_cast<int>(m_reductionValues[type][id].size()) >= var_index); 1134 #else 1135 SMART_ASSERT(m_reductionVariables[type].find(var_name) != m_reductionVariables[type].end()) 1136 (type)(var_name); 1137 var_index = m_reductionVariables[type].find(var_name)->second; 1138 SMART_ASSERT(static_cast<int>(m_reductionValues[type][id].size()) >= var_index); 1139 #endif 1140 // Transfer to 'variables' array. 1141 if (ioss_type == Ioss::Field::REAL) { 1142 rvar[i] = m_reductionValues[type][id][var_index - 1]; 1143 } 1144 else if (ioss_type == Ioss::Field::INT64) { 1145 i64var[i] = static_cast<int64_t>(m_reductionValues[type][id][var_index - 1]); 1146 } 1147 else if (ioss_type == Ioss::Field::INTEGER) { 1148 ivar[i] = static_cast<int>(m_reductionValues[type][id][var_index - 1]); 1149 } 1150 } 1151 } 1152 1153 // common write_reduction_fields()1154 void BaseDatabaseIO::write_reduction_fields() const 1155 { 1156 int step = get_current_state(); 1157 step = get_database_step(step); 1158 for (const auto &type : exodus_types) { 1159 auto &id_values = m_reductionValues[type]; 1160 for (const auto &values : id_values) { 1161 int64_t id = values.first; 1162 auto & vals = values.second; 1163 size_t count = vals.size(); 1164 if (count > 0) { 1165 int ierr = ex_put_reduction_vars(get_file_pointer(), step, type, id, count, vals.data()); 1166 if (ierr < 0) { 1167 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1168 } 1169 } 1170 } 1171 } 1172 } 1173 1174 // common read_reduction_fields()1175 void BaseDatabaseIO::read_reduction_fields() const 1176 { 1177 int step = get_current_state(); 1178 1179 for (const auto &type : exodus_types) { 1180 auto &id_values = m_reductionValues[type]; 1181 for (const auto &values : id_values) { 1182 int64_t id = values.first; 1183 auto & vals = values.second; 1184 size_t count = vals.size(); 1185 if (count > 0) { 1186 int ierr = ex_get_reduction_vars(get_file_pointer(), step, type, id, count, 1187 (double *)vals.data()); 1188 if (ierr < 0) { 1189 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1190 } 1191 } 1192 } 1193 } 1194 } 1195 1196 // common begin__(Ioss::State state)1197 bool BaseDatabaseIO::begin__(Ioss::State state) 1198 { 1199 dbState = state; 1200 return true; 1201 } 1202 1203 // common end__(Ioss::State state)1204 bool BaseDatabaseIO::end__(Ioss::State state) 1205 { 1206 // Transitioning out of state 'state' 1207 assert(state == dbState); 1208 switch (state) { 1209 case Ioss::STATE_DEFINE_MODEL: 1210 if (!is_input()) { 1211 write_meta_data(open_create_behavior()); 1212 } 1213 break; 1214 case Ioss::STATE_DEFINE_TRANSIENT: 1215 if (!is_input()) { 1216 write_results_metadata(true, open_create_behavior()); 1217 } 1218 break; 1219 default: // ignore everything else... 1220 break; 1221 } 1222 1223 { 1224 Ioss::SerializeIO serializeIO__(this); 1225 1226 if (!is_input()) { 1227 ex_update(get_file_pointer()); 1228 if (minimizeOpenFiles) { 1229 free_file_pointer(); 1230 } 1231 } 1232 dbState = Ioss::STATE_UNKNOWN; 1233 } 1234 1235 return true; 1236 } 1237 open_state_file(int state)1238 void BaseDatabaseIO::open_state_file(int state) 1239 { 1240 // Close current file... 1241 free_file_pointer(); 1242 1243 // Update filename to append state count... 1244 decodedFilename.clear(); 1245 1246 Ioss::FileInfo db(originalDBFilename); 1247 std::string new_filename; 1248 if (!db.pathname().empty()) { 1249 new_filename += db.pathname() + "/"; 1250 } 1251 1252 if (get_cycle_count() >= 1) { 1253 static const std::string suffix{"ABCDEFGHIJKLMNOPQRSTUVWXYZ"}; 1254 int index = (state - 1) % get_cycle_count(); 1255 new_filename += db.basename() + "-state-" + suffix[index] + "." + db.extension(); 1256 } 1257 else { 1258 new_filename += db.basename() + "-state-" + std::to_string(state) + "." + db.extension(); 1259 } 1260 1261 DBFilename = new_filename; 1262 fileExists = false; 1263 1264 ex_var_params exo_params{}; 1265 #if GLOBALS_ARE_TRANSIENT 1266 exo_params.num_glob = m_variables[EX_GLOBAL].size(); 1267 #else 1268 exo_params.num_glob = m_reductionVariables[EX_GLOBAL].size(); 1269 #endif 1270 exo_params.num_node = m_variables[EX_NODE_BLOCK].size(); 1271 exo_params.num_edge = m_variables[EX_EDGE_BLOCK].size(); 1272 exo_params.num_face = m_variables[EX_FACE_BLOCK].size(); 1273 exo_params.num_elem = m_variables[EX_ELEM_BLOCK].size(); 1274 exo_params.num_nset = m_variables[EX_NODE_SET].size(); 1275 exo_params.num_eset = m_variables[EX_EDGE_SET].size(); 1276 exo_params.num_fset = m_variables[EX_FACE_SET].size(); 1277 exo_params.num_sset = m_variables[EX_SIDE_SET].size(); 1278 exo_params.num_elset = m_variables[EX_ELEM_SET].size(); 1279 1280 char the_title[max_line_length + 1]; 1281 1282 // Title... 1283 if (get_region()->property_exists("title")) { 1284 std::string title_str = get_region()->get_property("title").get_string(); 1285 Ioss::Utils::copy_string(the_title, title_str); 1286 } 1287 else { 1288 Ioss::Utils::copy_string(the_title, "IOSS Default Output Title"); 1289 } 1290 1291 Ioex::Mesh mesh(spatialDimension, the_title, util(), !usingParallelIO); 1292 mesh.populate(get_region()); 1293 1294 // Write the metadata to the exodus file... 1295 Ioex::Internals data(get_file_pointer(), maximumNameLength, util()); 1296 int ierr = data.initialize_state_file(mesh, exo_params, originalDBFilename); 1297 1298 if (ierr < 0) { 1299 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1300 } 1301 } 1302 begin_state__(int state,double time)1303 bool BaseDatabaseIO::begin_state__(int state, double time) 1304 { 1305 Ioss::SerializeIO serializeIO__(this); 1306 1307 time /= timeScaleFactor; 1308 1309 if (!is_input()) { 1310 if (get_file_per_state()) { 1311 // Close current file; create new file and output transient metadata... 1312 open_state_file(state); 1313 write_results_metadata(false, open_create_behavior()); 1314 } 1315 int ierr = ex_put_time(get_file_pointer(), get_database_step(state), &time); 1316 if (ierr < 0) { 1317 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1318 } 1319 1320 // Zero global variable array... 1321 for (auto &type : exodus_types) { 1322 auto &id_values = m_reductionValues[type]; 1323 for (auto &values : id_values) { 1324 auto &vals = values.second; 1325 std::fill(vals.begin(), vals.end(), 0.0); 1326 } 1327 } 1328 } 1329 else { 1330 // Store reduction variables 1331 read_reduction_fields(); 1332 } 1333 return true; 1334 } 1335 1336 // common end_state__(int state,double time)1337 bool BaseDatabaseIO::end_state__(int state, double time) 1338 { 1339 Ioss::SerializeIO serializeIO__(this); 1340 1341 if (!is_input()) { 1342 write_reduction_fields(); 1343 time /= timeScaleFactor; 1344 finalize_write(state, time); 1345 if (minimizeOpenFiles) { 1346 free_file_pointer(); 1347 } 1348 } 1349 return true; 1350 } 1351 1352 // common add_region_fields()1353 void BaseDatabaseIO::add_region_fields() 1354 { 1355 #if GLOBALS_ARE_TRANSIENT 1356 int field_count = add_results_fields(EX_GLOBAL, get_region()); 1357 #else 1358 int field_count = add_reduction_results_fields(EX_GLOBAL, get_region()); 1359 #endif 1360 m_reductionValues[EX_GLOBAL][0].resize(field_count); 1361 add_mesh_reduction_fields(EX_GLOBAL, 0, get_region()); 1362 } 1363 1364 namespace { 1365 // Memory allocated in `ex_get_attributes`, this makes deletion cleaner... 1366 class EX_attribute : public ex_attribute 1367 { 1368 public: EX_attribute()1369 EX_attribute() { values = nullptr; } ~EX_attribute()1370 ~EX_attribute() { free(values); } 1371 }; 1372 } // namespace 1373 add_mesh_reduction_fields(ex_entity_type type,int64_t id,Ioss::GroupingEntity * entity)1374 void BaseDatabaseIO::add_mesh_reduction_fields(ex_entity_type type, int64_t id, 1375 Ioss::GroupingEntity *entity) 1376 { 1377 // Get "global attributes" 1378 // These are single key-value per grouping entity 1379 // Stored as Ioss::Property with origin of ATTRIBUTE 1380 Ioss::SerializeIO serializeIO__(this); 1381 int att_count = ex_get_attribute_count(get_file_pointer(), type, id); 1382 1383 if (att_count > 0) { 1384 std::vector<EX_attribute> attr(att_count); 1385 ex_get_attribute_param(get_file_pointer(), type, id, attr.data()); 1386 ex_get_attributes(get_file_pointer(), att_count, attr.data()); 1387 1388 // Create a property on `entity` for each `attribute` 1389 for (const auto &att : attr) { 1390 if (att.values == nullptr) { 1391 continue; 1392 } 1393 std::string storage = fmt::format("Real[{}]", att.value_count); 1394 switch (att.type) { 1395 case EX_INTEGER: 1396 if (att.value_count == 1) { 1397 entity->property_add( 1398 Ioss::Property(att.name, *(int *)att.values, Ioss::Property::ATTRIBUTE)); 1399 } 1400 else { 1401 const int * idata = static_cast<int *>(att.values); 1402 std::vector<int> tmp(att.value_count); 1403 std::copy(idata, idata + att.value_count, tmp.begin()); 1404 entity->property_add(Ioss::Property(att.name, tmp, Ioss::Property::ATTRIBUTE)); 1405 } 1406 break; 1407 case EX_DOUBLE: 1408 if (att.value_count == 1) { 1409 entity->property_add( 1410 Ioss::Property(att.name, *(double *)att.values, Ioss::Property::ATTRIBUTE)); 1411 } 1412 else { 1413 const double * idata = static_cast<double *>(att.values); 1414 std::vector<double> tmp(att.value_count); 1415 std::copy(idata, idata + att.value_count, tmp.begin()); 1416 entity->property_add(Ioss::Property(att.name, tmp, Ioss::Property::ATTRIBUTE)); 1417 } 1418 break; 1419 case EX_CHAR: 1420 if (att.value_count > 0) { 1421 entity->property_add( 1422 Ioss::Property(att.name, (char *)att.values, Ioss::Property::ATTRIBUTE)); 1423 } 1424 else { 1425 // Just an attribute name. Give it an empty value... 1426 entity->property_add(Ioss::Property(att.name, "", Ioss::Property::ATTRIBUTE)); 1427 } 1428 break; 1429 } 1430 } 1431 } 1432 } 1433 1434 // common add_results_fields(ex_entity_type type,Ioss::GroupingEntity * entity,int64_t position)1435 int64_t BaseDatabaseIO::add_results_fields(ex_entity_type type, Ioss::GroupingEntity *entity, 1436 int64_t position) 1437 { 1438 return internal_add_results_fields(type, entity, position, m_groupCount[type], 1439 m_truthTable[type], m_variables[type]); 1440 } 1441 1442 // common internal_add_results_fields(ex_entity_type type,Ioss::GroupingEntity * entity,int64_t position,int64_t block_count,Ioss::IntVector & truth_table,Ioex::VariableNameMap & variables)1443 int64_t BaseDatabaseIO::internal_add_results_fields(ex_entity_type type, 1444 Ioss::GroupingEntity *entity, 1445 int64_t position, int64_t block_count, 1446 Ioss::IntVector & truth_table, 1447 Ioex::VariableNameMap &variables) 1448 { 1449 int nvar = 0; 1450 { 1451 Ioss::SerializeIO serializeIO__(this); 1452 1453 int ierr = ex_get_variable_param(get_file_pointer(), type, &nvar); 1454 if (ierr < 0) { 1455 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1456 } 1457 } 1458 1459 if (nvar > 0) { 1460 if (truth_table.empty()) { 1461 truth_table.resize(block_count * nvar); 1462 1463 // Read and store the truth table (Should be there since we only 1464 // get to this routine if there are variables...) 1465 1466 if (type == EX_NODE_BLOCK || type == EX_GLOBAL || type == EX_ASSEMBLY) { 1467 // These types don't have a truth table in the exodus api... 1468 // They do in Ioss just for some consistency... 1469 std::fill(truth_table.begin(), truth_table.end(), 1); 1470 } 1471 else { 1472 Ioss::SerializeIO serializeIO__(this); 1473 int ierr = 1474 ex_get_truth_table(get_file_pointer(), type, block_count, nvar, truth_table.data()); 1475 if (ierr < 0) { 1476 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1477 } 1478 } 1479 1480 // If parallel, then synchronize the truth table among all 1481 // processors... Need to know that block_X has variable_Y 1482 // even if block_X is empty on a specific processor... The 1483 // truth table contains 0 if the variable doesn't exist and 1 1484 // if it does, so we just take the maximum at each location... 1485 // This is a collective call... Make sure not in Serialize 1486 if (isParallel) { 1487 util().global_array_minmax(truth_table, Ioss::ParallelUtils::DO_MAX); 1488 } 1489 } 1490 1491 // Get the variable names and add as fields. Need to decode these 1492 // into vector/tensor/... eventually, for now store all as 1493 // scalars. 1494 char **names = Ioss::Utils::get_name_array(nvar, maximumNameLength); 1495 1496 // Read the names... 1497 // (Currently, names are read for every block. We could save them...) 1498 { 1499 Ioss::SerializeIO serializeIO__(this); 1500 1501 int ierr = ex_get_variable_names(get_file_pointer(), type, nvar, names); 1502 if (ierr < 0) { 1503 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1504 } 1505 1506 // Add to VariableNameMap so can determine exodusII index given a 1507 // Sierra field name. exodusII index is just 'i+1' 1508 for (int i = 0; i < nvar; i++) { 1509 if (lowerCaseVariableNames) { 1510 Ioss::Utils::fixup_name(names[i]); 1511 } 1512 variables.insert(VNMValuePair(std::string(names[i]), i + 1)); 1513 } 1514 1515 int offset = position * nvar; 1516 int *local_truth = nullptr; 1517 if (!truth_table.empty()) { 1518 local_truth = &truth_table[offset]; 1519 } 1520 1521 std::vector<Ioss::Field> fields; 1522 int64_t count = entity->entity_count(); 1523 Ioss::Utils::get_fields(count, names, nvar, Ioss::Field::TRANSIENT, get_field_recognition(), 1524 get_field_separator(), local_truth, fields); 1525 1526 for (const auto &field : fields) { 1527 entity->field_add(field); 1528 } 1529 1530 for (int i = 0; i < nvar; i++) { 1531 // Verify that all names were used for a field... 1532 assert(names[i][0] == '\0' || (local_truth && local_truth[i] == 0)); 1533 delete[] names[i]; 1534 } 1535 delete[] names; 1536 } 1537 } 1538 return nvar; 1539 } 1540 1541 // common add_reduction_results_fields(ex_entity_type type,Ioss::GroupingEntity * entity)1542 int64_t BaseDatabaseIO::add_reduction_results_fields(ex_entity_type type, 1543 Ioss::GroupingEntity *entity) 1544 { 1545 int nvar = 0; 1546 { 1547 Ioss::SerializeIO serializeIO__(this); 1548 1549 int ierr = ex_get_reduction_variable_param(get_file_pointer(), type, &nvar); 1550 if (ierr < 0) { 1551 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1552 } 1553 } 1554 1555 if (nvar > 0) { 1556 // Get the variable names and add as fields. Need to decode these 1557 // into vector/tensor/... eventually, for now store all as 1558 // scalars. 1559 char **names = Ioss::Utils::get_name_array(nvar, maximumNameLength); 1560 1561 // Read the names... 1562 // (Currently, names are read for every block. We could save them...) 1563 { 1564 Ioss::SerializeIO serializeIO__(this); 1565 1566 int ierr = ex_get_reduction_variable_names(get_file_pointer(), type, nvar, names); 1567 if (ierr < 0) { 1568 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1569 } 1570 1571 // Add to VariableNameMap so can determine exodusII index given a 1572 // Sierra field name. exodusII index is just 'i+1' 1573 auto &variables = m_reductionVariables[type]; 1574 for (int i = 0; i < nvar; i++) { 1575 if (lowerCaseVariableNames) { 1576 Ioss::Utils::fixup_name(names[i]); 1577 } 1578 variables.insert(VNMValuePair(std::string(names[i]), i + 1)); 1579 } 1580 1581 int *local_truth = nullptr; 1582 1583 std::vector<Ioss::Field> fields; 1584 int64_t count = 1; 1585 Ioss::Utils::get_fields(count, names, nvar, Ioss::Field::REDUCTION, get_field_recognition(), 1586 get_field_separator(), local_truth, fields); 1587 1588 for (const auto &field : fields) { 1589 entity->field_add(field); 1590 } 1591 1592 for (int i = 0; i < nvar; i++) { 1593 // Verify that all names were used for a field... 1594 assert(names[i][0] == '\0' || (local_truth && local_truth[i] == 0)); 1595 delete[] names[i]; 1596 } 1597 delete[] names; 1598 } 1599 } 1600 return nvar; 1601 } 1602 1603 // common write_results_metadata(bool gather_data,Ioss::IfDatabaseExistsBehavior behavior)1604 void BaseDatabaseIO::write_results_metadata(bool gather_data, 1605 Ioss::IfDatabaseExistsBehavior behavior) 1606 { 1607 if (gather_data) { 1608 int glob_index = 0; 1609 #if GLOBALS_ARE_TRANSIENT 1610 glob_index = gather_names(EX_GLOBAL, m_variables[EX_GLOBAL], get_region(), glob_index, true); 1611 #else 1612 glob_index = 1613 gather_names(EX_GLOBAL, m_reductionVariables[EX_GLOBAL], get_region(), glob_index, true); 1614 #endif 1615 m_reductionValues[EX_GLOBAL][0].resize(glob_index); 1616 1617 const Ioss::NodeBlockContainer &node_blocks = get_region()->get_node_blocks(); 1618 assert(node_blocks.size() <= 1); 1619 internal_gather_results_metadata(EX_NODE_BLOCK, node_blocks); 1620 1621 const Ioss::EdgeBlockContainer &edge_blocks = get_region()->get_edge_blocks(); 1622 internal_gather_results_metadata(EX_EDGE_BLOCK, edge_blocks); 1623 1624 const Ioss::FaceBlockContainer &face_blocks = get_region()->get_face_blocks(); 1625 internal_gather_results_metadata(EX_FACE_BLOCK, face_blocks); 1626 1627 const Ioss::ElementBlockContainer &element_blocks = get_region()->get_element_blocks(); 1628 internal_gather_results_metadata(EX_ELEM_BLOCK, element_blocks); 1629 1630 const Ioss::NodeSetContainer &nodesets = get_region()->get_nodesets(); 1631 internal_gather_results_metadata(EX_NODE_SET, nodesets); 1632 1633 const Ioss::EdgeSetContainer &edgesets = get_region()->get_edgesets(); 1634 internal_gather_results_metadata(EX_EDGE_SET, edgesets); 1635 1636 const Ioss::FaceSetContainer &facesets = get_region()->get_facesets(); 1637 internal_gather_results_metadata(EX_FACE_SET, facesets); 1638 1639 const Ioss::ElementSetContainer &elementsets = get_region()->get_elementsets(); 1640 internal_gather_results_metadata(EX_ELEM_SET, elementsets); 1641 1642 const auto &blobs = get_region()->get_blobs(); 1643 internal_gather_results_metadata(EX_BLOB, blobs); 1644 1645 const auto &assemblies = get_region()->get_assemblies(); 1646 internal_gather_results_metadata(EX_ASSEMBLY, assemblies); 1647 1648 { 1649 int index = 0; 1650 const Ioss::SideSetContainer &sidesets = get_region()->get_sidesets(); 1651 for (const auto &sideset : sidesets) { 1652 const Ioss::SideBlockContainer &side_blocks = sideset->get_side_blocks(); 1653 for (const auto &block : side_blocks) { 1654 glob_index = gather_names(EX_SIDE_SET, m_reductionVariables[EX_SIDE_SET], block, 1655 glob_index, true); 1656 index = gather_names(EX_SIDE_SET, m_variables[EX_SIDE_SET], block, index, false); 1657 } 1658 } 1659 generate_sideset_truth_table(); 1660 } 1661 } 1662 1663 if (behavior != Ioss::DB_APPEND && behavior != Ioss::DB_MODIFY) { 1664 ex_var_params exo_params{}; 1665 #if GLOBALS_ARE_TRANSIENT 1666 exo_params.num_glob = m_variables[EX_GLOBAL].size(); 1667 #else 1668 exo_params.num_glob = m_reductionVariables[EX_GLOBAL].size(); 1669 #endif 1670 exo_params.num_node = m_variables[EX_NODE_BLOCK].size(); 1671 exo_params.num_edge = m_variables[EX_EDGE_BLOCK].size(); 1672 exo_params.num_face = m_variables[EX_FACE_BLOCK].size(); 1673 exo_params.num_elem = m_variables[EX_ELEM_BLOCK].size(); 1674 exo_params.num_nset = m_variables[EX_NODE_SET].size(); 1675 exo_params.num_eset = m_variables[EX_EDGE_SET].size(); 1676 exo_params.num_fset = m_variables[EX_FACE_SET].size(); 1677 exo_params.num_sset = m_variables[EX_SIDE_SET].size(); 1678 exo_params.num_elset = m_variables[EX_ELEM_SET].size(); 1679 1680 exo_params.edge_var_tab = m_truthTable[EX_EDGE_BLOCK].data(); 1681 exo_params.face_var_tab = m_truthTable[EX_FACE_BLOCK].data(); 1682 exo_params.elem_var_tab = m_truthTable[EX_ELEM_BLOCK].data(); 1683 exo_params.nset_var_tab = m_truthTable[EX_NODE_SET].data(); 1684 exo_params.eset_var_tab = m_truthTable[EX_EDGE_SET].data(); 1685 exo_params.fset_var_tab = m_truthTable[EX_FACE_SET].data(); 1686 exo_params.sset_var_tab = m_truthTable[EX_SIDE_SET].data(); 1687 exo_params.elset_var_tab = m_truthTable[EX_ELEM_SET].data(); 1688 1689 if (isParallel) { 1690 // Check consistency among all processors. They should all 1691 // have the same number of each variable type... 1692 // The called function will throw an exception if the counts differ. 1693 check_variable_consistency(exo_params, myProcessor, get_filename(), util()); 1694 } 1695 1696 { 1697 Ioss::SerializeIO serializeIO__(this); 1698 1699 int ierr = ex_put_all_var_param_ext(get_file_pointer(), &exo_params); 1700 if (ierr < 0) { 1701 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1702 } 1703 1704 // Blob and Assembly not supported in ex_put_all_var_param_ext... 1705 if (!m_variables[EX_BLOB].empty()) { 1706 ierr = ex_put_variable_param(get_file_pointer(), EX_BLOB, m_variables[EX_BLOB].size()); 1707 if (ierr < 0) { 1708 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1709 } 1710 } 1711 if (!m_variables[EX_ASSEMBLY].empty()) { 1712 ierr = ex_put_variable_param(get_file_pointer(), EX_ASSEMBLY, 1713 m_variables[EX_ASSEMBLY].size()); 1714 if (ierr < 0) { 1715 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1716 } 1717 } 1718 1719 for (const auto &type : exodus_types) { 1720 output_results_names(type, m_variables[type], false); 1721 output_results_names(type, m_reductionVariables[type], true); 1722 } 1723 } 1724 } 1725 } 1726 1727 // common 1728 template <typename T> internal_gather_results_metadata(ex_entity_type type,std::vector<T * > entities)1729 void BaseDatabaseIO::internal_gather_results_metadata(ex_entity_type type, 1730 std::vector<T *> entities) 1731 { 1732 int index = 0; 1733 int red_index = 0; 1734 for (const auto &entity : entities) { 1735 red_index = gather_names(type, m_reductionVariables[type], entity, red_index, true); 1736 index = gather_names(type, m_variables[type], entity, index, false); 1737 } 1738 1739 #if GLOBALS_ARE_TRANSIENT 1740 size_t value_size = 1741 type == EX_GLOBAL ? m_variables[type].size() : m_reductionVariables[type].size(); 1742 #else 1743 size_t value_size = m_reductionVariables[type].size(); 1744 #endif 1745 for (const auto &entity : entities) { 1746 auto id = entity->get_optional_property("id", 0); 1747 m_reductionValues[type][id].resize(value_size); 1748 } 1749 1750 generate_block_truth_table(m_variables[type], m_truthTable[type], entities, 1751 get_field_separator()); 1752 } 1753 1754 // common gather_names(ex_entity_type type,VariableNameMap & variables,const Ioss::GroupingEntity * ge,int index,bool reduction)1755 int BaseDatabaseIO::gather_names(ex_entity_type type, VariableNameMap &variables, 1756 const Ioss::GroupingEntity *ge, int index, bool reduction) 1757 { 1758 int new_index = index; 1759 1760 bool nblock = (type == EX_NODE_BLOCK); 1761 1762 // Get names of all transient and reduction fields... 1763 Ioss::NameList results_fields; 1764 if (reduction) { 1765 ge->field_describe(Ioss::Field::REDUCTION, &results_fields); 1766 } 1767 1768 if (!reduction || type == EX_GLOBAL) { 1769 ge->field_describe(Ioss::Field::TRANSIENT, &results_fields); 1770 } 1771 1772 // NOTE: For exodusII, the convention is that the displacement 1773 // fields are the first 'ndim' fields in the file. 1774 // Try to find a likely displacement field 1775 std::string disp_name; 1776 bool has_disp = false; 1777 if (!reduction && nblock && new_index == 0) { 1778 has_disp = find_displacement_field(results_fields, ge, spatialDimension, &disp_name); 1779 if (has_disp) { 1780 new_index += spatialDimension; 1781 } 1782 } 1783 1784 int save_index = 0; 1785 for (const auto &name : results_fields) { 1786 1787 if (has_disp && name == disp_name && new_index != 0) { 1788 save_index = new_index; 1789 new_index = 0; 1790 } 1791 1792 Ioss::Field field = ge->get_field(name); 1793 const Ioss::VariableType *var_type = field.transformed_storage(); 1794 1795 int re_im = 1; 1796 if (field.get_type() == Ioss::Field::COMPLEX) { 1797 re_im = 2; 1798 } 1799 for (int complex_comp = 0; complex_comp < re_im; complex_comp++) { 1800 std::string field_name = field.get_name(); 1801 if (re_im == 2) { 1802 field_name += complex_suffix[complex_comp]; 1803 } 1804 1805 char field_suffix_separator = get_field_separator(); 1806 for (int i = 1; i <= var_type->component_count(); i++) { 1807 std::string var_string = var_type->label_name(field_name, i, field_suffix_separator); 1808 1809 if (variables.find(var_string) == variables.end()) { 1810 variables.insert(VNMValuePair(var_string, ++new_index)); 1811 } 1812 } 1813 } 1814 if (has_disp && name == disp_name) { 1815 new_index = save_index; 1816 } 1817 } 1818 return new_index; 1819 } 1820 1821 // common generate_sideset_truth_table()1822 void BaseDatabaseIO::generate_sideset_truth_table() 1823 { 1824 size_t var_count = m_variables[EX_SIDE_SET].size(); 1825 1826 if (var_count == 0 || m_groupCount[EX_SIDE_SET] == 0) { 1827 return; 1828 } 1829 1830 // Member variable. Will be deleted in destructor... 1831 m_truthTable[EX_SIDE_SET].resize(m_groupCount[EX_SIDE_SET] * var_count); 1832 1833 // Fill in the truth table. It is conceptually a two-dimensional array of 1834 // the form 'array[num_blocks][num_var]'. In C++, 1835 // the values for the first block are first, followed by 1836 // next block, ... 1837 size_t offset = 0; 1838 1839 char field_suffix_separator = get_field_separator(); 1840 1841 const Ioss::SideSetContainer &sidesets = get_region()->get_sidesets(); 1842 for (const auto &sideset : sidesets) { 1843 1844 const Ioss::SideBlockContainer &side_blocks = sideset->get_side_blocks(); 1845 for (const auto &block : side_blocks) { 1846 // See if this sideblock has a corresponding entry in the sideset list. 1847 if (block->property_exists("invalid")) { 1848 continue; 1849 } 1850 1851 // Get names of all transient and reduction fields... 1852 Ioss::NameList results_fields; 1853 block->field_describe(Ioss::Field::TRANSIENT, &results_fields); 1854 block->field_describe(Ioss::Field::REDUCTION, &results_fields); 1855 1856 for (const auto &fn : results_fields) { 1857 Ioss::Field field = block->get_field(fn); 1858 const Ioss::VariableType *var_type = field.transformed_storage(); 1859 Ioss::Field::BasicType ioss_type = field.get_type(); 1860 1861 int re_im = 1; 1862 if (ioss_type == Ioss::Field::COMPLEX) { 1863 re_im = 2; 1864 } 1865 for (int complex_comp = 0; complex_comp < re_im; complex_comp++) { 1866 std::string field_name = field.get_name(); 1867 if (re_im == 2) { 1868 field_name += complex_suffix[complex_comp]; 1869 } 1870 1871 for (int i = 1; i <= var_type->component_count(); i++) { 1872 std::string var_string = var_type->label_name(field_name, i, field_suffix_separator); 1873 // Find position of 'var_string' in 'm_variables[]' 1874 auto VN = m_variables[EX_SIDE_SET].find(var_string); 1875 if (VN != m_variables[EX_SIDE_SET].end()) { 1876 // Index '(*VN).second' is 1-based... 1877 m_truthTable[EX_SIDE_SET][offset + (*VN).second - 1] = 1; 1878 } 1879 } 1880 } 1881 } 1882 } 1883 offset += var_count; 1884 } 1885 assert(offset == var_count * m_groupCount[EX_SIDE_SET]); 1886 } 1887 1888 // common output_results_names(ex_entity_type type,VariableNameMap & variables,bool reduction)1889 void BaseDatabaseIO::output_results_names(ex_entity_type type, VariableNameMap &variables, 1890 bool reduction) const 1891 { 1892 bool lowercase_names = 1893 (properties.exists("VARIABLE_NAME_CASE") && 1894 Ioss::Utils::lowercase(properties.get("VARIABLE_NAME_CASE").get_string()) == "lower"); 1895 bool uppercase_names = 1896 (properties.exists("VARIABLE_NAME_CASE") && 1897 Ioss::Utils::lowercase(properties.get("VARIABLE_NAME_CASE").get_string()) == "upper"); 1898 1899 size_t var_count = variables.size(); 1900 1901 if (var_count > 0) { 1902 size_t name_length = 0; 1903 // Push into a char** array... 1904 std::vector<char *> var_names(var_count); 1905 std::vector<std::string> variable_names(var_count); 1906 for (const auto &variable : variables) { 1907 size_t index = variable.second; 1908 assert(index > 0 && index <= var_count); 1909 variable_names[index - 1] = variable.first; 1910 if (uppercase_names) { 1911 variable_names[index - 1] = Ioss::Utils::uppercase(variable_names[index - 1]); 1912 } 1913 else if (lowercase_names) { 1914 variable_names[index - 1] = Ioss::Utils::lowercase(variable_names[index - 1]); 1915 } 1916 var_names[index - 1] = const_cast<char *>(variable_names[index - 1].c_str()); 1917 size_t name_len = variable_names[index - 1].length(); 1918 name_length = name_len > name_length ? name_len : name_length; 1919 } 1920 1921 // Should handle this automatically, but by the time we get to defining transient fields, we 1922 // have already created the output database and populated the set/block names. At this point, 1923 // it is too late to change the size of the names stored on the output database... (I think... 1924 // try changing DIM_STR_NAME value and see if works...) 1925 if (name_length > static_cast<size_t>(maximumNameLength)) { 1926 if (myProcessor == 0) { 1927 fmt::print(Ioss::WARNING(), 1928 "There are variables names whose name length ({0}) exceeds the current " 1929 "maximum name length ({1})\n set for this database ({2}).\n" 1930 " You should either reduce the length of the variable name, or " 1931 "set the 'MAXIMUM_NAME_LENGTH' property\n" 1932 " to at least {0}.\n Contact gdsjaar@sandia.gov for more " 1933 "information.\n\n", 1934 name_length, maximumNameLength, get_filename()); 1935 } 1936 } 1937 int ierr = 0; 1938 if (reduction) { 1939 ierr = 1940 ex_put_reduction_variable_names(get_file_pointer(), type, var_count, var_names.data()); 1941 } 1942 else { 1943 ierr = ex_put_variable_names(get_file_pointer(), type, var_count, var_names.data()); 1944 } 1945 if (ierr < 0) { 1946 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1947 } 1948 } 1949 } 1950 1951 // common 1952 // Handle special output time requests -- primarily restart (cycle, overwrite) 1953 // Given the global region step, return the step on the database... get_database_step(int global_step)1954 int BaseDatabaseIO::get_database_step(int global_step) const 1955 { 1956 if (get_file_per_state()) { 1957 return 1; 1958 } 1959 1960 assert(overlayCount >= 0 && cycleCount >= 0); 1961 if (overlayCount == 0 && cycleCount == 0) { 1962 return global_step; 1963 } 1964 1965 int local_step = global_step - 1; 1966 local_step /= (overlayCount + 1); 1967 if (cycleCount > 0) { 1968 local_step %= cycleCount; 1969 } 1970 return local_step + 1; 1971 } 1972 1973 // common flush_database__()1974 void BaseDatabaseIO::flush_database__() const 1975 { 1976 if (!is_input()) { 1977 if (isParallel || myProcessor == 0) { 1978 ex_update(get_file_pointer()); 1979 } 1980 } 1981 } 1982 finalize_write(int state,double sim_time)1983 void BaseDatabaseIO::finalize_write(int state, double sim_time) 1984 { 1985 // Attempt to ensure that all data written up to this point has 1986 // actually made it out to disk. We also write a special attribute 1987 // to the file to indicate that the current timestep should be 1988 // complete on the disk. 1989 // The attribute is a GLOBAL attribute named "last_written_time" 1990 // which is a double value which can be compared to the values in 1991 // the time array to make sure they match. If they don't, then 1992 // hopefully the "last_written_time" is smaller than the time 1993 // array value and indicates that the last step is corrupt. 1994 1995 // Update the attribute. 1996 Ioex::update_last_time_attribute(get_file_pointer(), sim_time); 1997 1998 // Flush the files buffer to disk... 1999 // If a history file, then only flush if there is more 2000 // than 10 seconds since the last flush to avoid 2001 // the flush eating up cpu time for small fast jobs... 2002 // NOTE: If decide to do this on all files, need to sync across 2003 // processors to make sure they all flush at same time. 2004 2005 // GDS: 2011/03/30 -- Use for all non-parallel files, but shorten 2006 // time for non history files. Assume that can afford to lose ~10 2007 // seconds worth of data... (Flush was taking long time on some 2008 // /scratch filesystems at SNL for short regression tests with 2009 // lots of steps) 2010 // GDS: 2011/07/27 -- shorten from 90 to 10. Developers running 2011 // small jobs were not able to view output until job 2012 // finished. Hopefully the netcdf no-fsync fix along with this fix 2013 // results in negligible impact on runtime with more syncs. 2014 2015 // Need to be able to handle a flushInterval == 1 to force flush 2016 // every time step even in a serial run. 2017 // The default setting for flushInterval is 1, but in the past, 2018 // it was not checked for serial runs. Now, set the default to -1 2019 // and if that is the value and serial, then do the time-based 2020 // check; otherwise, use flushInterval setting... 2021 2022 bool do_flush = true; 2023 if (flushInterval == 1) { 2024 do_flush = true; 2025 } 2026 else if (flushInterval == 0) { 2027 do_flush = false; 2028 } 2029 else if (dbUsage == Ioss::WRITE_HISTORY || !isParallel) { 2030 assert(myProcessor == 0); 2031 time_t cur_time = time(nullptr); 2032 if (cur_time - timeLastFlush >= 10) { 2033 timeLastFlush = cur_time; 2034 do_flush = true; 2035 } 2036 else { 2037 do_flush = false; 2038 } 2039 } 2040 2041 if (!do_flush && flushInterval > 0) { 2042 if (state % flushInterval == 0) { 2043 do_flush = true; 2044 } 2045 } 2046 2047 if (do_flush) { 2048 flush_database__(); 2049 } 2050 } 2051 2052 // common add_attribute_fields(ex_entity_type entity_type,Ioss::GroupingEntity * block,int attribute_count,const std::string & type)2053 void Ioex::BaseDatabaseIO::add_attribute_fields(ex_entity_type entity_type, 2054 Ioss::GroupingEntity *block, int attribute_count, 2055 const std::string &type) 2056 { 2057 /// \todo REFACTOR Some of the attribute knowledge should be at 2058 /// the Ioss::ElementTopology level instead of here. That would 2059 /// make it easier for an application to register a new element 2060 /// type and its attributes. 2061 2062 // Attribute "Conventions" to be used if there are no attribute names on the database: 2063 // from Table 1 in ExodusII manual 2064 // 2065 // Circle 1 Radius [Volume] 2066 // Sphere 1 Radius [Volume] 2067 // Truss 1 Area 2068 // 2D Beam 3 Area, I, J 2069 // 3D Beam 7 Area, I1, I2, J, V1, V2, V3 (V will be a 3D vector named "reference_axis") 2070 // Shell 1 Thickness 2071 // 2072 // Additional conventions not defined in ExodusII manual: 2073 // * If a "beam" has 1 attribute, call it "area" 2074 // * Treat "bar" and "rod" as aliases for "truss" 2075 // * Treat "trishell" as alias for "shell" 2076 // * All "shell" or "trishell" elements -- If #attributes == #node/element, the 2077 // attribute is "nodal_thickness" 2078 // 2079 // If there are attribute names on the database, use those names. 2080 // Always create a variable "attribute" which contains a single 2081 // field for all attributes... 2082 2083 assert(block != nullptr); 2084 if (attribute_count > 0) { 2085 size_t my_element_count = block->entity_count(); 2086 2087 // Get the attribute names. May not exist or may be blank... 2088 char ** names = Ioss::Utils::get_name_array(attribute_count, maximumNameLength); 2089 int64_t id = block->get_property("id").get_int(); 2090 2091 // Some older applications do not want to used named 2092 // attributes; in this case, just create a field for each 2093 // attribute named "attribute_1", "attribute_2", ..., "attribute_#" 2094 // This is controlled by the database property 2095 // "IGNORE_ATTRIBUTE_NAMES" 2096 char field_suffix_separator = get_field_separator(); 2097 bool attributes_named = true; // Possibly reset below; note that even if ignoring 2098 // attribute names, they are still 'named' 2099 2100 if (properties.exists("IGNORE_ATTRIBUTE_NAMES")) { 2101 field_suffix_separator = ' '; // Do not combine into a 2102 // higher-order storage type. 2103 2104 for (int i = 0; i < attribute_count; i++) { 2105 std::string tmp = fmt::format("attribute_{}", i + 1); 2106 Ioss::Utils::copy_string(names[i], tmp, maximumNameLength + 1); 2107 } 2108 } 2109 else { 2110 // Use attribute names if they exist. 2111 { 2112 Ioss::SerializeIO serializeIO__(this); 2113 if (block->entity_count() != 0) { 2114 int ierr = ex_get_attr_names(get_file_pointer(), entity_type, id, &names[0]); 2115 if (ierr < 0) { 2116 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2117 } 2118 } 2119 } 2120 2121 // Sync names across processors... 2122 if (isParallel) { 2123 std::vector<char> cname(attribute_count * (maximumNameLength + 1)); 2124 if (block->entity_count() != 0) { 2125 for (int i = 0; i < attribute_count; i++) { 2126 std::memcpy(&cname[i * (maximumNameLength + 1)], names[i], maximumNameLength + 1); 2127 } 2128 } 2129 util().attribute_reduction(attribute_count * (maximumNameLength + 1), cname.data()); 2130 for (int i = 0; i < attribute_count; i++) { 2131 std::memcpy(names[i], &cname[i * (maximumNameLength + 1)], maximumNameLength + 1); 2132 } 2133 } 2134 2135 // Convert to lowercase. 2136 attributes_named = true; 2137 for (int i = 0; i < attribute_count; i++) { 2138 fix_bad_name(names[i]); 2139 Ioss::Utils::fixup_name(names[i]); 2140 if (names[i][0] == '\0' || (!(std::isalnum(names[i][0]) || names[i][0] == '_'))) { 2141 attributes_named = false; 2142 } 2143 } 2144 } 2145 2146 if (attributes_named) { 2147 std::vector<Ioss::Field> attributes; 2148 Ioss::Utils::get_fields(my_element_count, names, attribute_count, Ioss::Field::ATTRIBUTE, 2149 get_field_recognition(), field_suffix_separator, nullptr, 2150 attributes); 2151 int offset = 1; 2152 for (const auto &field : attributes) { 2153 if (block->field_exists(field.get_name())) { 2154 std::ostringstream errmsg; 2155 fmt::print(errmsg, 2156 "ERROR: In block '{}', attribute '{}' is defined multiple times which is " 2157 "not allowed.\n", 2158 block->name(), field.get_name()); 2159 IOSS_ERROR(errmsg); 2160 } 2161 block->field_add(field); 2162 const Ioss::Field &tmp_field = block->get_fieldref(field.get_name()); 2163 tmp_field.set_index(offset); 2164 offset += field.raw_storage()->component_count(); 2165 } 2166 } 2167 else { 2168 // Attributes are not named.... 2169 // Try to assign some meaningful names based on conventions... 2170 std::string att_name = "attribute"; // Default 2171 int unknown_attributes = 0; 2172 2173 if (type_match(type, "shell") || type_match(type, "trishell")) { 2174 if (attribute_count == block->get_property("topology_node_count").get_int()) { 2175 2176 att_name = "nodal_thickness"; 2177 2178 std::string storage = fmt::format("Real[{}]", attribute_count); 2179 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, storage, 2180 Ioss::Field::ATTRIBUTE, my_element_count, 1)); 2181 } 2182 else { 2183 att_name = "thickness"; 2184 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, IOSS_SCALAR(), 2185 Ioss::Field::ATTRIBUTE, my_element_count, 1)); 2186 unknown_attributes = attribute_count - 1; 2187 } 2188 } 2189 2190 // NOTE: This must appear before the "sphere" check since 2191 // sphere is substring of "sphere-mass" 2192 // Want an exact match here, not substring match... 2193 else if (Ioss::Utils::str_equal(type, "sphere-mass")) { 2194 if (attribute_count != 10) { 2195 if (myProcessor == 0) { 2196 fmt::print(Ioss::WARNING(), 2197 "For element block '{}' of type '{}' there were {} attributes instead of " 2198 "the expected 10 attributes " 2199 "known to the IO Subsystem. " 2200 " The attributes can be accessed as the field named 'attribute'", 2201 block->name(), type, attribute_count); 2202 } 2203 } 2204 else { 2205 // First attribute is concentrated mass... 2206 size_t offset = 1; 2207 block->field_add(Ioss::Field("mass", Ioss::Field::REAL, IOSS_SCALAR(), 2208 Ioss::Field::ATTRIBUTE, my_element_count, offset)); 2209 offset += 1; 2210 2211 // Next six attributes are moment of inertia -- symmetric tensor 2212 block->field_add(Ioss::Field("inertia", Ioss::Field::REAL, IOSS_SYM_TENSOR(), 2213 Ioss::Field::ATTRIBUTE, my_element_count, offset)); 2214 offset += 6; 2215 2216 // Next three attributes are offset from node to CG 2217 block->field_add(Ioss::Field("offset", Ioss::Field::REAL, IOSS_VECTOR_3D(), 2218 Ioss::Field::ATTRIBUTE, my_element_count, offset)); 2219 } 2220 } 2221 2222 else if (type_match(type, "circle") || type_match(type, "sphere")) { 2223 att_name = "radius"; 2224 size_t offset = 1; 2225 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, IOSS_SCALAR(), 2226 Ioss::Field::ATTRIBUTE, my_element_count, offset++)); 2227 if (attribute_count > 1) { 2228 // Default second attribute (from sphgen3d) is "volume" 2229 // which is the volume of the cube which would surround a 2230 // sphere of the given radius. 2231 att_name = "volume"; 2232 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, IOSS_SCALAR(), 2233 Ioss::Field::ATTRIBUTE, my_element_count, offset++)); 2234 } 2235 unknown_attributes = attribute_count - 2; 2236 } 2237 2238 else if (type_match(type, "truss") || type_match(type, "bar") || type_match(type, "beam") || 2239 type_match(type, "rod")) { 2240 // Technically, truss, bar, rod should all only have 1 attribute; however, 2241 // there are some mesh generation codes that treat all of these types the 2242 // same and put "beam-type" attributes on bars... 2243 int index = 1; 2244 att_name = "area"; 2245 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, IOSS_SCALAR(), 2246 Ioss::Field::ATTRIBUTE, my_element_count, index++)); 2247 2248 if (spatialDimension == 2 && attribute_count >= 3) { 2249 block->field_add(Ioss::Field("i", Ioss::Field::REAL, IOSS_SCALAR(), 2250 Ioss::Field::ATTRIBUTE, my_element_count, index++)); 2251 block->field_add(Ioss::Field("j", Ioss::Field::REAL, IOSS_SCALAR(), 2252 Ioss::Field::ATTRIBUTE, my_element_count, index++)); 2253 } 2254 else if (spatialDimension == 3 && attribute_count >= 7) { 2255 block->field_add(Ioss::Field("i1", Ioss::Field::REAL, IOSS_SCALAR(), 2256 Ioss::Field::ATTRIBUTE, my_element_count, index++)); 2257 block->field_add(Ioss::Field("i2", Ioss::Field::REAL, IOSS_SCALAR(), 2258 Ioss::Field::ATTRIBUTE, my_element_count, index++)); 2259 block->field_add(Ioss::Field("j", Ioss::Field::REAL, IOSS_SCALAR(), 2260 Ioss::Field::ATTRIBUTE, my_element_count, index++)); 2261 block->field_add(Ioss::Field("reference_axis", Ioss::Field::REAL, IOSS_VECTOR_3D(), 2262 Ioss::Field::ATTRIBUTE, my_element_count, index)); 2263 index += 3; 2264 if (attribute_count >= 10) { 2265 // Next three attributes would (hopefully) be offset vector... 2266 // This is typically from a NASGEN model. 2267 block->field_add(Ioss::Field("offset", Ioss::Field::REAL, IOSS_VECTOR_3D(), 2268 Ioss::Field::ATTRIBUTE, my_element_count, index)); 2269 index += 3; 2270 } 2271 } 2272 unknown_attributes = attribute_count - (index - 1); 2273 } 2274 2275 else { 2276 unknown_attributes = attribute_count; 2277 } 2278 2279 if (unknown_attributes > 0) { 2280 att_name = "extra_attribute_"; 2281 att_name += std::to_string(unknown_attributes); 2282 std::string storage = fmt::format("Real[{}]", unknown_attributes); 2283 size_t index = attribute_count - unknown_attributes + 1; 2284 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, storage, Ioss::Field::ATTRIBUTE, 2285 my_element_count, index)); 2286 } 2287 } 2288 2289 // Always create a field called "attribute" containing data 2290 // for all attributes on the mesh 2291 std::string att_name = "attribute"; // Default 2292 std::string storage = fmt::format("Real[{}]", attribute_count); 2293 2294 block->field_add(Ioss::Field(att_name, Ioss::Field::REAL, storage, Ioss::Field::ATTRIBUTE, 2295 my_element_count, 1)); 2296 2297 // Release memory... 2298 Ioss::Utils::delete_name_array(names, attribute_count); 2299 } 2300 } 2301 common_write_meta_data(Ioss::IfDatabaseExistsBehavior behavior)2302 void BaseDatabaseIO::common_write_meta_data(Ioss::IfDatabaseExistsBehavior behavior) 2303 { 2304 Ioss::Region *region = get_region(); 2305 2306 // Verify that exodus supports the mesh_type... 2307 if (region->mesh_type() != Ioss::MeshType::UNSTRUCTURED) { 2308 std::ostringstream errmsg; 2309 fmt::print(errmsg, 2310 "ERROR: The mesh type is '{}' which Exodus does not support.\n" 2311 " Only 'Unstructured' is supported at this time.\n", 2312 region->mesh_type_string()); 2313 IOSS_ERROR(errmsg); 2314 } 2315 2316 const Ioss::NodeBlockContainer &node_blocks = region->get_node_blocks(); 2317 assert(node_blocks.size() <= 1); 2318 if (!node_blocks.empty()) { 2319 Ioex::get_id(node_blocks[0], EX_NODE_BLOCK, &ids_); 2320 nodeCount = node_blocks[0]->entity_count(); 2321 spatialDimension = node_blocks[0]->get_property("component_degree").get_int(); 2322 } 2323 else { 2324 spatialDimension = 1; 2325 } 2326 2327 // Assemblies -- 2328 { 2329 const auto &assemblies = region->get_assemblies(); 2330 if (behavior != Ioss::DB_MODIFY) { 2331 // Set ids of all entities that have "id" property... 2332 for (auto &assem : assemblies) { 2333 Ioex::set_id(assem, EX_ASSEMBLY, &ids_); 2334 } 2335 2336 for (auto &assem : assemblies) { 2337 Ioex::get_id(assem, EX_ASSEMBLY, &ids_); 2338 } 2339 } 2340 m_groupCount[EX_ASSEMBLY] = assemblies.size(); 2341 } 2342 2343 // Blobs -- 2344 { 2345 const auto &blobs = region->get_blobs(); 2346 // Set ids of all entities that have "id" property... 2347 if (behavior != Ioss::DB_MODIFY) { 2348 for (auto &blob : blobs) { 2349 Ioex::set_id(blob, EX_BLOB, &ids_); 2350 } 2351 2352 for (auto &blob : blobs) { 2353 Ioex::get_id(blob, EX_BLOB, &ids_); 2354 } 2355 } 2356 m_groupCount[EX_BLOB] = blobs.size(); 2357 } 2358 2359 // Edge Blocks -- 2360 { 2361 const Ioss::EdgeBlockContainer &edge_blocks = region->get_edge_blocks(); 2362 assert(Ioss::Utils::check_block_order(edge_blocks)); 2363 // Set ids of all entities that have "id" property... 2364 if (behavior != Ioss::DB_MODIFY) { 2365 for (auto &edge_block : edge_blocks) { 2366 Ioex::set_id(edge_block, EX_EDGE_BLOCK, &ids_); 2367 } 2368 2369 edgeCount = 0; 2370 for (auto &edge_block : edge_blocks) { 2371 edgeCount += edge_block->entity_count(); 2372 // Set ids of all entities that do not have "id" property... 2373 Ioex::get_id(edge_block, EX_EDGE_BLOCK, &ids_); 2374 } 2375 } 2376 m_groupCount[EX_EDGE_BLOCK] = edge_blocks.size(); 2377 } 2378 2379 // Face Blocks -- 2380 { 2381 const Ioss::FaceBlockContainer &face_blocks = region->get_face_blocks(); 2382 assert(Ioss::Utils::check_block_order(face_blocks)); 2383 // Set ids of all entities that have "id" property... 2384 if (behavior != Ioss::DB_MODIFY) { 2385 for (auto &face_block : face_blocks) { 2386 Ioex::set_id(face_block, EX_FACE_BLOCK, &ids_); 2387 } 2388 2389 faceCount = 0; 2390 for (auto &face_block : face_blocks) { 2391 faceCount += face_block->entity_count(); 2392 // Set ids of all entities that do not have "id" property... 2393 Ioex::get_id(face_block, EX_FACE_BLOCK, &ids_); 2394 } 2395 } 2396 m_groupCount[EX_FACE_BLOCK] = face_blocks.size(); 2397 } 2398 2399 // Element Blocks -- 2400 { 2401 const Ioss::ElementBlockContainer &element_blocks = region->get_element_blocks(); 2402 assert(Ioss::Utils::check_block_order(element_blocks)); 2403 if (behavior != Ioss::DB_MODIFY) { 2404 // Set ids of all entities that have "id" property... 2405 for (auto &element_block : element_blocks) { 2406 Ioex::set_id(element_block, EX_ELEM_BLOCK, &ids_); 2407 } 2408 } 2409 elementCount = 0; 2410 Ioss::Int64Vector element_counts; 2411 element_counts.reserve(element_blocks.size()); 2412 for (auto &element_block : element_blocks) { 2413 elementCount += element_block->entity_count(); 2414 element_counts.push_back(element_block->entity_count()); 2415 // Set ids of all entities that do not have "id" property... 2416 if (behavior != Ioss::DB_MODIFY) { 2417 Ioex::get_id(element_block, EX_ELEM_BLOCK, &ids_); 2418 } 2419 } 2420 m_groupCount[EX_ELEM_BLOCK] = element_blocks.size(); 2421 2422 if (isParallel) { 2423 // Set "global_entity_count" property on all blocks. 2424 // Used to skip output on "globally" empty blocks. 2425 Ioss::Int64Vector global_counts(element_counts.size()); 2426 util().global_count(element_counts, global_counts); 2427 size_t idx = 0; 2428 for (auto &element_block : element_blocks) { 2429 element_block->property_add(Ioss::Property("global_entity_count", global_counts[idx++])); 2430 } 2431 } 2432 } 2433 2434 // NodeSets ... 2435 { 2436 const Ioss::NodeSetContainer &nodesets = region->get_nodesets(); 2437 if (behavior != Ioss::DB_MODIFY) { 2438 for (auto &set : nodesets) { 2439 Ioex::set_id(set, EX_NODE_SET, &ids_); 2440 } 2441 2442 for (auto &set : nodesets) { 2443 Ioex::get_id(set, EX_NODE_SET, &ids_); 2444 } 2445 } 2446 m_groupCount[EX_NODE_SET] = nodesets.size(); 2447 } 2448 2449 // EdgeSets ... 2450 { 2451 const Ioss::EdgeSetContainer &edgesets = region->get_edgesets(); 2452 if (behavior != Ioss::DB_MODIFY) { 2453 for (auto &set : edgesets) { 2454 Ioex::set_id(set, EX_EDGE_SET, &ids_); 2455 } 2456 2457 for (auto &set : edgesets) { 2458 Ioex::get_id(set, EX_EDGE_SET, &ids_); 2459 } 2460 } 2461 m_groupCount[EX_EDGE_SET] = edgesets.size(); 2462 } 2463 2464 // FaceSets ... 2465 { 2466 const Ioss::FaceSetContainer &facesets = region->get_facesets(); 2467 if (behavior != Ioss::DB_MODIFY) { 2468 for (auto &set : facesets) { 2469 Ioex::set_id(set, EX_FACE_SET, &ids_); 2470 } 2471 2472 for (auto &set : facesets) { 2473 Ioex::get_id(set, EX_FACE_SET, &ids_); 2474 } 2475 } 2476 m_groupCount[EX_FACE_SET] = facesets.size(); 2477 } 2478 2479 // ElementSets ... 2480 { 2481 const Ioss::ElementSetContainer &elementsets = region->get_elementsets(); 2482 if (behavior != Ioss::DB_MODIFY) { 2483 for (auto &set : elementsets) { 2484 Ioex::set_id(set, EX_ELEM_SET, &ids_); 2485 } 2486 2487 for (auto &set : elementsets) { 2488 Ioex::get_id(set, EX_ELEM_SET, &ids_); 2489 } 2490 } 2491 m_groupCount[EX_ELEM_SET] = elementsets.size(); 2492 } 2493 2494 // SideSets ... 2495 { 2496 const Ioss::SideSetContainer &ssets = region->get_sidesets(); 2497 if (behavior != Ioss::DB_MODIFY) { 2498 for (auto &set : ssets) { 2499 Ioex::set_id(set, EX_SIDE_SET, &ids_); 2500 } 2501 } 2502 // Get entity counts for all face sets... Create SideSets. 2503 for (auto &set : ssets) { 2504 if (behavior != Ioss::DB_MODIFY) { 2505 Ioex::get_id(set, EX_SIDE_SET, &ids_); 2506 } 2507 int64_t id = set->get_property("id").get_int(); 2508 int64_t entity_count = 0; 2509 int64_t df_count = 0; 2510 2511 const Ioss::SideBlockContainer &side_blocks = set->get_side_blocks(); 2512 for (auto &block : side_blocks) { 2513 // Add "*_offset" properties to specify at what offset 2514 // the data for this block appears in the containing set. 2515 auto *new_block = const_cast<Ioss::SideBlock *>(block); 2516 new_block->property_add(Ioss::Property("set_offset", entity_count)); 2517 new_block->property_add(Ioss::Property("set_df_offset", df_count)); 2518 2519 // If combining sideblocks into sidesets on output, then 2520 // the id of the sideblock must be the same as the sideset 2521 // id. 2522 new_block->property_update("id", id); 2523 new_block->property_update("guid", util().generate_guid(id)); 2524 2525 entity_count += block->entity_count(); 2526 df_count += block->get_property("distribution_factor_count").get_int(); 2527 } 2528 auto *new_entity = const_cast<Ioss::SideSet *>(set); 2529 new_entity->property_add(Ioss::Property("entity_count", entity_count)); 2530 new_entity->property_add(Ioss::Property("distribution_factor_count", df_count)); 2531 } 2532 m_groupCount[EX_SIDE_SET] = ssets.size(); 2533 } 2534 } 2535 output_other_meta_data()2536 void BaseDatabaseIO::output_other_meta_data() 2537 { 2538 // Write attribute names (if any)... 2539 char field_suffix_separator = get_field_separator(); 2540 write_attribute_names(get_file_pointer(), EX_NODE_SET, get_region()->get_nodesets(), 2541 field_suffix_separator); 2542 write_attribute_names(get_file_pointer(), EX_EDGE_SET, get_region()->get_edgesets(), 2543 field_suffix_separator); 2544 write_attribute_names(get_file_pointer(), EX_FACE_SET, get_region()->get_facesets(), 2545 field_suffix_separator); 2546 write_attribute_names(get_file_pointer(), EX_ELEM_SET, get_region()->get_elementsets(), 2547 field_suffix_separator); 2548 write_attribute_names(get_file_pointer(), EX_NODE_BLOCK, get_region()->get_node_blocks(), 2549 field_suffix_separator); 2550 write_attribute_names(get_file_pointer(), EX_EDGE_BLOCK, get_region()->get_edge_blocks(), 2551 field_suffix_separator); 2552 write_attribute_names(get_file_pointer(), EX_FACE_BLOCK, get_region()->get_face_blocks(), 2553 field_suffix_separator); 2554 write_attribute_names(get_file_pointer(), EX_ELEM_BLOCK, get_region()->get_element_blocks(), 2555 field_suffix_separator); 2556 write_attribute_names(get_file_pointer(), EX_ASSEMBLY, get_region()->get_assemblies(), 2557 field_suffix_separator); 2558 write_attribute_names(get_file_pointer(), EX_BLOB, get_region()->get_blobs(), 2559 field_suffix_separator); 2560 2561 // Write "reduction" attributes... 2562 std::vector<Ioss::Region *> regions; 2563 regions.push_back(get_region()); 2564 Ioex::write_reduction_attributes(get_file_pointer(), regions); 2565 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_nodesets()); 2566 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_nodesets()); 2567 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_edgesets()); 2568 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_facesets()); 2569 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_elementsets()); 2570 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_node_blocks()); 2571 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_edge_blocks()); 2572 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_face_blocks()); 2573 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_element_blocks()); 2574 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_assemblies()); 2575 Ioex::write_reduction_attributes(get_file_pointer(), get_region()->get_blobs()); 2576 2577 // Write coordinate names... 2578 if (!get_region()->get_node_blocks().empty()) { 2579 char const *labels[3]; 2580 labels[0] = "x"; 2581 labels[1] = "y"; 2582 labels[2] = "z"; 2583 int ierr = ex_put_coord_names(get_file_pointer(), (char **)labels); 2584 if (ierr < 0) { 2585 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2586 } 2587 } 2588 2589 // Write coordinate frame data... 2590 write_coordinate_frames(get_file_pointer(), get_region()->get_coordinate_frames()); 2591 } 2592 } // namespace Ioex 2593 2594 namespace { 2595 template <typename T> write_attribute_names(int exoid,ex_entity_type type,const std::vector<T * > & entities,const char suffix_separator)2596 void write_attribute_names(int exoid, ex_entity_type type, const std::vector<T *> &entities, 2597 const char suffix_separator) 2598 { 2599 // For the entity, determine the attribute fields and the correct 2600 // order. Write the names of these fields. However, be aware that 2601 // the field "attribute" always exists to contain all attributes 2602 // and its name should not be used even if it is the only 2603 // attribute field. 2604 for (const auto &ge : entities) { 2605 int attribute_count = ge->get_property("attribute_count").get_int(); 2606 if (attribute_count > 0) { 2607 2608 check_attribute_index_order(ge); 2609 2610 std::vector<char *> names(attribute_count); 2611 std::vector<std::string> names_str(attribute_count); 2612 2613 // Get the attribute fields... 2614 Ioss::NameList results_fields; 2615 ge->field_describe(Ioss::Field::ATTRIBUTE, &results_fields); 2616 2617 for (const auto &field_name : results_fields) { 2618 const Ioss::Field &field = ge->get_fieldref(field_name); 2619 assert(field.get_index() != 0); 2620 2621 if (field_name == "attribute") { 2622 field.set_index(1); 2623 continue; 2624 } 2625 2626 const Ioss::VariableType *vtype = field.raw_storage(); 2627 int comp_count = vtype->component_count(); 2628 int field_offset = field.get_index(); 2629 for (int i = 0; i < comp_count; i++) { 2630 names_str[field_offset - 1 + i] = 2631 vtype->label_name(field_name, i + 1, suffix_separator); 2632 names[field_offset - 1 + i] = 2633 const_cast<char *>(names_str[field_offset - 1 + i].c_str()); 2634 } 2635 } 2636 size_t ge_id = ge->get_property("id").get_int(); 2637 int ierr = ex_put_attr_names(exoid, type, ge_id, names.data()); 2638 if (ierr < 0) { 2639 Ioex::exodus_error(exoid, __LINE__, __func__, __FILE__); 2640 } 2641 } 2642 } 2643 } 2644 2645 // common check_attribute_index_order(Ioss::GroupingEntity * block)2646 void check_attribute_index_order(Ioss::GroupingEntity *block) 2647 { 2648 int attribute_count = block->get_property("attribute_count").get_int(); 2649 if (attribute_count == 0) { 2650 return; 2651 } 2652 int component_sum = 0; 2653 2654 std::vector<int> attributes(attribute_count + 1); 2655 2656 // Get the attribute fields... 2657 Ioss::NameList results_fields; 2658 block->field_describe(Ioss::Field::ATTRIBUTE, &results_fields); 2659 2660 bool all_attributes_indexed = true; 2661 bool some_attributes_indexed = false; 2662 2663 for (const auto &field_name : results_fields) { 2664 const Ioss::Field &field = block->get_fieldref(field_name); 2665 2666 if (field_name == "attribute") { 2667 field.set_index(1); 2668 if (results_fields.size() == 1) { 2669 return; 2670 } 2671 continue; 2672 } 2673 2674 int field_offset = field.get_index(); 2675 if (field_offset == 0) { 2676 all_attributes_indexed = false; 2677 } 2678 else { 2679 some_attributes_indexed = true; 2680 } 2681 2682 const Ioss::VariableType *type = field.raw_storage(); 2683 int comp_count = type->component_count(); 2684 component_sum += comp_count; 2685 2686 if (field_offset == 0) { 2687 continue; 2688 } 2689 2690 if (field_offset + comp_count - 1 > attribute_count) { 2691 std::ostringstream errmsg; 2692 fmt::print( 2693 errmsg, 2694 "INTERNAL ERROR: For block '{}', attribute '{}', the indexing is incorrect.\n" 2695 "Something is wrong in the Ioex::BaseDatabaseIO class, function {}. Please report.\n", 2696 block->name(), field_name, __func__); 2697 IOSS_ERROR(errmsg); 2698 } 2699 2700 for (int i = field_offset; i < field_offset + comp_count; i++) { 2701 if (attributes[i] != 0) { 2702 std::ostringstream errmsg; 2703 fmt::print( 2704 errmsg, 2705 "INTERNAL ERROR: For block '{}', attribute '{}', indexes into the same location as a " 2706 "previous attribute.\n" 2707 "Something is wrong in the Ioex::BaseDatabaseIO class, function {}. Please report.\n", 2708 block->name(), field_name, __func__); 2709 IOSS_ERROR(errmsg); 2710 } 2711 attributes[i] = 1; 2712 } 2713 } 2714 2715 if (component_sum > attribute_count) { 2716 std::ostringstream errmsg; 2717 fmt::print( 2718 errmsg, 2719 "INTERNAL ERROR: Block '{}' is supposed to have {} attributes, but {} attributes " 2720 "were counted.\n" 2721 "Something is wrong in the Ioex::BaseDatabaseIO class, function {}. Please report.\n", 2722 block->name(), attribute_count, component_sum, __func__); 2723 IOSS_ERROR(errmsg); 2724 } 2725 2726 // Take care of the easy cases first... 2727 if (all_attributes_indexed) { 2728 // Check that all attributes are defined. This should have 2729 // caught above in the duplicate index check. 2730 for (int i = 1; i <= attribute_count; i++) { 2731 if (attributes[i] == 0) { 2732 std::ostringstream errmsg; 2733 fmt::print( 2734 errmsg, 2735 "INTERNAL ERROR: Block '{}' has an incomplete set of attributes.\n" 2736 "Something is wrong in the Ioex::BaseDatabaseIO class, function {}. Please report.\n", 2737 block->name(), __func__); 2738 IOSS_ERROR(errmsg); 2739 } 2740 } 2741 return; 2742 } 2743 2744 if (!some_attributes_indexed) { 2745 // Index was not set for any of the attributes; set them all... 2746 size_t offset = 1; 2747 for (const auto &field_name : results_fields) { 2748 const Ioss::Field &field = block->get_fieldref(field_name); 2749 2750 if (field_name == "attribute") { 2751 field.set_index(1); 2752 continue; 2753 } 2754 2755 const Ioss::VariableType *type = field.raw_storage(); 2756 int comp_count = type->component_count(); 2757 2758 assert(field.get_index() == 0); 2759 field.set_index(offset); 2760 offset += comp_count; 2761 } 2762 assert((int)offset == attribute_count + 1); 2763 return; 2764 } 2765 2766 // At this point, we have a partially indexed set of attributes. Some have an index and some 2767 // don't 2768 // The easy case is if the missing indices are at the end of the list... 2769 assert(!all_attributes_indexed && some_attributes_indexed); 2770 int last_defined = 0; 2771 for (int i = 1; i < attribute_count + 1; i++) { 2772 if (attributes[i] != 0) { 2773 last_defined = i; 2774 } 2775 } 2776 int first_undefined = attribute_count; 2777 for (int i = attribute_count; i > 0; i--) { 2778 if (attributes[i] == 0) { 2779 first_undefined = i; 2780 } 2781 } 2782 if (last_defined < first_undefined) { 2783 for (const auto &field_name : results_fields) { 2784 const Ioss::Field &field = block->get_fieldref(field_name); 2785 2786 if (field_name == "attribute") { 2787 field.set_index(1); 2788 continue; 2789 } 2790 2791 if (field.get_index() == 0) { 2792 field.set_index(first_undefined); 2793 const Ioss::VariableType *type = field.raw_storage(); 2794 int comp_count = type->component_count(); 2795 first_undefined += comp_count; 2796 } 2797 } 2798 assert(first_undefined == attribute_count + 1); 2799 return; 2800 } 2801 2802 // Take the easy way out... Just reindex all attributes. 2803 size_t offset = 1; 2804 for (const auto &field_name : results_fields) { 2805 const Ioss::Field &field = block->get_fieldref(field_name); 2806 2807 if (field_name == "attribute") { 2808 field.set_index(1); 2809 continue; 2810 } 2811 2812 const Ioss::VariableType *type = field.raw_storage(); 2813 int comp_count = type->component_count(); 2814 2815 assert(field.get_index() == 0); 2816 field.set_index(offset); 2817 offset += comp_count; 2818 } 2819 assert((int)offset == attribute_count + 1); 2820 } 2821 check_variable_consistency(const ex_var_params & exo_params,int my_processor,const std::string & filename,const Ioss::ParallelUtils & util)2822 void check_variable_consistency(const ex_var_params &exo_params, int my_processor, 2823 const std::string &filename, const Ioss::ParallelUtils &util) 2824 { 2825 PAR_UNUSED(exo_params); 2826 PAR_UNUSED(my_processor); 2827 PAR_UNUSED(filename); 2828 PAR_UNUSED(util); 2829 #ifdef SEACAS_HAVE_MPI 2830 const int num_types = 10; 2831 std::vector<int> var_counts(num_types); 2832 var_counts[0] = exo_params.num_glob; 2833 var_counts[1] = exo_params.num_node; 2834 var_counts[2] = exo_params.num_edge; 2835 var_counts[3] = exo_params.num_face; 2836 var_counts[4] = exo_params.num_elem; 2837 var_counts[5] = exo_params.num_nset; 2838 var_counts[6] = exo_params.num_eset; 2839 var_counts[7] = exo_params.num_fset; 2840 var_counts[8] = exo_params.num_sset; 2841 var_counts[9] = exo_params.num_elset; 2842 2843 Ioss::IntVector all_counts; 2844 util.gather(var_counts, all_counts); 2845 2846 bool any_diff = false; 2847 std::ostringstream errmsg; 2848 if (my_processor == 0) { 2849 bool diff[num_types]; 2850 // See if any differ... 2851 for (int iv = 0; iv < 10; iv++) { 2852 diff[iv] = false; 2853 std::string type; 2854 switch (iv) { 2855 case 0: type = "global"; break; 2856 case 1: type = "nodal"; break; 2857 case 2: type = "edge"; break; 2858 case 3: type = "face"; break; 2859 case 4: type = "element"; break; 2860 case 5: type = "nodeset"; break; 2861 case 6: type = "edgeset"; break; 2862 case 7: type = "faceset"; break; 2863 case 8: type = "sideset"; break; 2864 case 9: type = "elementset"; break; 2865 } 2866 2867 for (int ip = 1; ip < util.parallel_size(); ip++) { 2868 if (var_counts[iv] != all_counts[ip * num_types + iv]) { 2869 any_diff = true; 2870 if (!diff[iv]) { 2871 Ioss::FileInfo db(filename); 2872 diff[iv] = true; 2873 fmt::print(errmsg, 2874 "\nERROR: Number of {} variables is not consistent on all processors.\n" 2875 " Database: '{}'\n" 2876 "\tProcessor 0 count = {}\n", 2877 type, db.tailname(), var_counts[iv]); 2878 } 2879 fmt::print(errmsg, "\tProcessor {} count = {}\n", ip, all_counts[ip * num_types + iv]); 2880 } 2881 } 2882 } 2883 } 2884 else { 2885 // Give the other processors something to say... 2886 fmt::print(errmsg, "ERROR: Variable type counts are inconsistent. See processor 0 output for " 2887 "more details.\n"); 2888 } 2889 int idiff = any_diff ? 1 : 0; 2890 MPI_Bcast(&idiff, 1, MPI_INT, 0, util.communicator()); 2891 any_diff = idiff == 1; 2892 2893 if (any_diff) { 2894 std::runtime_error x(errmsg.str()); 2895 throw x; 2896 } 2897 #endif 2898 } 2899 } // namespace 2900