1 // ISSUES: 2 // 1. Does not handle unconnected nodes (not connected to any element) 3 // 4 // 2. SideSet distribution factors are klugy and may not fully work in 5 // strange cases 6 // 7 // 8 // Copyright(C) 1999-2021 National Technology & Engineering Solutions 9 // of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with 10 // NTESS, the U.S. Government retains certain rights in this software. 11 // 12 // See packages/seacas/LICENSE for details 13 14 #include <exodus/Ioex_ParallelDatabaseIO.h> 15 #if defined PARALLEL_AWARE_EXODUS 16 #include <algorithm> 17 #include <cassert> 18 #include <cctype> 19 #include <cfloat> 20 #include <cstddef> 21 #include <cstdlib> 22 #include <cstring> 23 #include <ctime> 24 #include <fmt/ostream.h> 25 #include <functional> 26 #include <iostream> 27 #include <limits> 28 #include <map> 29 #include <numeric> 30 #include <set> 31 #include <string> 32 #include <tokenize.h> 33 #ifndef _MSC_VER 34 #include <unistd.h> 35 #endif 36 #include <utility> 37 #include <vector> 38 39 #include <Ioss_CodeTypes.h> 40 41 #include <exodus/Ioex_DecompositionData.h> 42 #include <exodus/Ioex_Internals.h> 43 #include <exodus/Ioex_Utils.h> 44 #include <vtk_exodusII.h> 45 46 #include <Ioss_Assembly.h> 47 #include <Ioss_Blob.h> 48 #include <Ioss_CommSet.h> 49 #include <Ioss_CoordinateFrame.h> 50 #include <Ioss_DBUsage.h> 51 #include <Ioss_DatabaseIO.h> 52 #include <Ioss_EdgeBlock.h> 53 #include <Ioss_EdgeSet.h> 54 #include <Ioss_ElementBlock.h> 55 #include <Ioss_ElementSet.h> 56 #include <Ioss_ElementTopology.h> 57 #include <Ioss_EntityBlock.h> 58 #include <Ioss_EntitySet.h> 59 #include <Ioss_EntityType.h> 60 #include <Ioss_FaceBlock.h> 61 #include <Ioss_FaceSet.h> 62 #include <Ioss_Field.h> 63 #include <Ioss_FileInfo.h> 64 #include <Ioss_GroupingEntity.h> 65 #include <Ioss_Map.h> 66 #include <Ioss_NodeBlock.h> 67 #include <Ioss_NodeSet.h> 68 #include <Ioss_ParallelUtils.h> 69 #include <Ioss_Property.h> 70 #include <Ioss_Region.h> 71 #include <Ioss_SideBlock.h> 72 #include <Ioss_SideSet.h> 73 #include <Ioss_State.h> 74 #include <Ioss_SurfaceSplit.h> 75 #include <Ioss_Utils.h> 76 #include <Ioss_VariableType.h> 77 78 #include <Ioss_FileInfo.h> 79 #undef MPICPP 80 81 // ======================================================================== 82 // Static internal helper functions 83 // ======================================================================== 84 namespace { 85 const size_t max_line_length = MAX_LINE_LENGTH; 86 SEP()87 const std::string SEP() { return std::string("@"); } // Separator for attribute offset storage 88 const char * complex_suffix[] = {".re", ".im"}; 89 check_node_owning_processor_data(const Ioss::IntVector & nop,size_t file_node_count)90 void check_node_owning_processor_data(const Ioss::IntVector &nop, size_t file_node_count) 91 { 92 // Verify that the nop (NodeOwningProcessor) vector is not empty and is of the correct size. 93 // This vector specifies which rank owns each node on this rank 94 // Throws error if problem, otherwise returns quietly. 95 if (file_node_count == 0) { 96 return; 97 } 98 if (nop.empty()) { 99 std::ostringstream errmsg; 100 fmt::print(errmsg, "ERROR: The use of the 'compose' output option requires the definition of " 101 "the 'owning_processor'" 102 " field prior to the output of nodal data. This field has not yet been " 103 "defined so output is not possible." 104 " For more information, contact gdsjaar@sandia.gov.\n"); 105 IOSS_ERROR(errmsg); 106 } 107 else if (nop.size() < file_node_count) { 108 std::ostringstream errmsg; 109 fmt::print(errmsg, 110 "ERROR: The 'owning_processor' data was defined, but it is not the correct size." 111 " Its size is {}, but it must be at least this size {}." 112 " For more information, contact gdsjaar@sandia.gov.\n", 113 nop.size(), file_node_count); 114 IOSS_ERROR(errmsg); 115 } 116 } get_connectivity_data(int exoid,void * data,ex_entity_type type,ex_entity_id id,int position,int int_size_api)117 void get_connectivity_data(int exoid, void *data, ex_entity_type type, ex_entity_id id, 118 int position, int int_size_api) 119 { 120 int ierr = 0; 121 if (int_size_api == 8) { 122 int64_t *conn[3]; 123 conn[0] = nullptr; 124 conn[1] = nullptr; 125 conn[2] = nullptr; 126 conn[position] = static_cast<int64_t *>(data); 127 assert(1 == 0 && "Unimplemented fixme"); 128 ierr = ex_get_conn(exoid, type, id, conn[0], conn[1], conn[2]); // FIXME 129 } 130 else { 131 int *conn[3]; 132 conn[0] = nullptr; 133 conn[1] = nullptr; 134 conn[2] = nullptr; 135 conn[position] = static_cast<int *>(data); 136 assert(1 == 0 && "Unimplemented fixme"); 137 ierr = ex_get_conn(exoid, type, id, conn[0], conn[1], conn[2]); // FIXME 138 } 139 if (ierr < 0) { 140 Ioex::exodus_error(exoid, __LINE__, __func__, __FILE__); 141 } 142 } 143 144 template <typename T> compute_internal_border_maps(T * entities,T * internal,size_t count,size_t entity_count)145 void compute_internal_border_maps(T *entities, T *internal, size_t count, size_t entity_count) 146 { 147 for (size_t ij = 0; ij < count; ij++) { 148 internal[ij] = 1; 149 } 150 for (size_t J = 0; J < entity_count; J++) { 151 internal[entities[J] - 1] = 0; 152 } 153 154 size_t b = 0; 155 for (size_t ij = 0; ij < count; ij++) { 156 if (internal[ij] == 0) { 157 entities[b++] = ij + 1; 158 } 159 } 160 161 size_t k = 0; 162 for (size_t ij = 0; ij < count; ij++) { 163 if (internal[ij] == 1) { 164 internal[k++] = ij + 1; 165 } 166 } 167 } 168 169 template <typename INT> map_nodeset_id_data(const Ioss::IntVector & owning_processor,Ioss::Int64Vector & owned_nodes,int this_processor,const INT * ids,size_t ids_size,std::vector<INT> & file_data)170 void map_nodeset_id_data(const Ioss::IntVector &owning_processor, Ioss::Int64Vector &owned_nodes, 171 int this_processor, const INT *ids, size_t ids_size, 172 std::vector<INT> &file_data) 173 { 174 // Determine which nodes in this nodeset are owned by this processor. 175 // Save this mapping in the "owned_nodes" vector for use in 176 // mapping nodeset field data (df, transient, attributes, ...) 177 for (size_t i = 0; i < ids_size; i++) { 178 INT node = ids[i]; 179 if (owning_processor[node - 1] == this_processor) { 180 file_data.push_back(ids[i]); 181 owned_nodes.push_back(i); 182 } 183 } 184 } 185 186 template <typename T, typename U> 187 void map_nodeset_data(Ioss::Int64Vector &owned_nodes, const T *data, std::vector<U> &file_data, 188 size_t offset = 0, size_t stride = 1) 189 { 190 // Pull out the locally owned nodeset data 191 for (auto owned_node : owned_nodes) { 192 file_data.push_back(data[stride * owned_node + offset]); 193 } 194 } 195 196 template <typename T> extract_data(std::vector<double> & local_data,T * data,size_t num_entity,size_t offset,size_t comp_count)197 void extract_data(std::vector<double> &local_data, T *data, size_t num_entity, size_t offset, 198 size_t comp_count) 199 { 200 local_data.resize(num_entity); 201 if (comp_count == 1 && offset == 0) { 202 for (size_t j = 0; j < num_entity; j++) { 203 local_data[j] = data[j]; 204 } 205 } 206 else { 207 for (size_t j = 0; j < num_entity; j++) { 208 local_data[j] = data[offset]; 209 offset += comp_count; 210 } 211 } 212 } 213 214 // Ideally, there should only be a single data type for in and out 215 // data, but in the node id map mapping, we have an int64_t coming 216 // in and either an int or int64_t going out... 217 template <typename T, typename U> 218 void filter_owned_nodes(const Ioss::IntVector &owning_processor, int this_processor, 219 const T *data, std::vector<U> &file_data, size_t offset = 0, 220 size_t stride = 1) 221 { 222 size_t index = offset; 223 for (auto op : owning_processor) { 224 if (op == this_processor) { 225 file_data.push_back(data[index]); 226 } 227 index += stride; 228 } 229 } 230 231 // This version can be used *if* the input and output types are the same *and* the 232 // input `data` can be modified / overwritten. 233 template <typename T> filter_owned_nodes(const Ioss::IntVector & owning_processor,int this_processor,T * data)234 void filter_owned_nodes(const Ioss::IntVector &owning_processor, int this_processor, T *data) 235 { 236 size_t index = 0; 237 size_t entry = 0; 238 for (auto op : owning_processor) { 239 if (op == this_processor) { 240 data[entry++] = data[index]; 241 } 242 index++; 243 } 244 } 245 246 template <typename INT> map_local_to_global_implicit(INT * data,size_t count,const std::vector<int64_t> & global_implicit_map)247 void map_local_to_global_implicit(INT *data, size_t count, 248 const std::vector<int64_t> &global_implicit_map) 249 { 250 for (size_t i = 0; i < count; i++) { 251 data[i] = global_implicit_map[data[i] - 1]; 252 } 253 } 254 update_processor_offset_property(Ioss::Region * region,const Ioex::Mesh & mesh)255 void update_processor_offset_property(Ioss::Region *region, const Ioex::Mesh &mesh) 256 { 257 const Ioss::NodeBlockContainer &node_blocks = region->get_node_blocks(); 258 if (!node_blocks.empty()) { 259 node_blocks[0]->property_add( 260 Ioss::Property("_processor_offset", mesh.nodeblocks[0].procOffset)); 261 } 262 const Ioss::EdgeBlockContainer &edge_blocks = region->get_edge_blocks(); 263 for (size_t i = 0; i < edge_blocks.size(); i++) { 264 edge_blocks[i]->property_add( 265 Ioss::Property("_processor_offset", mesh.edgeblocks[i].procOffset)); 266 } 267 const Ioss::FaceBlockContainer &face_blocks = region->get_face_blocks(); 268 for (size_t i = 0; i < face_blocks.size(); i++) { 269 face_blocks[i]->property_add( 270 Ioss::Property("_processor_offset", mesh.faceblocks[i].procOffset)); 271 } 272 273 int64_t offset = 0; // Offset into global element map... 274 const Ioss::ElementBlockContainer &element_blocks = region->get_element_blocks(); 275 for (size_t i = 0; i < element_blocks.size(); i++) { 276 element_blocks[i]->property_add(Ioss::Property("global_map_offset", offset)); 277 offset += mesh.elemblocks[i].entityCount; 278 element_blocks[i]->property_add( 279 Ioss::Property("_processor_offset", mesh.elemblocks[i].procOffset)); 280 } 281 282 const Ioss::NodeSetContainer &nodesets = region->get_nodesets(); 283 for (size_t i = 0; i < nodesets.size(); i++) { 284 nodesets[i]->property_add(Ioss::Property("_processor_offset", mesh.nodesets[i].procOffset)); 285 } 286 const Ioss::EdgeSetContainer &edgesets = region->get_edgesets(); 287 for (size_t i = 0; i < edgesets.size(); i++) { 288 edgesets[i]->property_add(Ioss::Property("_processor_offset", mesh.edgesets[i].procOffset)); 289 } 290 const Ioss::FaceSetContainer &facesets = region->get_facesets(); 291 for (size_t i = 0; i < facesets.size(); i++) { 292 facesets[i]->property_add(Ioss::Property("_processor_offset", mesh.facesets[i].procOffset)); 293 } 294 const Ioss::ElementSetContainer &elementsets = region->get_elementsets(); 295 for (size_t i = 0; i < facesets.size(); i++) { 296 elementsets[i]->property_add( 297 Ioss::Property("_processor_offset", mesh.elemsets[i].procOffset)); 298 } 299 300 const Ioss::SideSetContainer &ssets = region->get_sidesets(); 301 for (size_t i = 0; i < ssets.size(); i++) { 302 ssets[i]->property_add(Ioss::Property("_processor_offset", mesh.sidesets[i].procOffset)); 303 ssets[i]->property_add(Ioss::Property("processor_df_offset", mesh.sidesets[i].dfProcOffset)); 304 305 // Propagate down to owned sideblocks... 306 const Ioss::SideBlockContainer &side_blocks = ssets[i]->get_side_blocks(); 307 for (auto &block : side_blocks) { 308 block->property_add(Ioss::Property("_processor_offset", mesh.sidesets[i].procOffset)); 309 block->property_add(Ioss::Property("processor_df_offset", mesh.sidesets[i].dfProcOffset)); 310 } 311 } 312 const auto &blobs = region->get_blobs(); 313 for (size_t i = 0; i < blobs.size(); i++) { 314 blobs[i]->property_add(Ioss::Property("_processor_offset", mesh.blobs[i].procOffset)); 315 } 316 } 317 } // namespace 318 319 namespace Ioex { ParallelDatabaseIO(Ioss::Region * region,const std::string & filename,Ioss::DatabaseUsage db_usage,MPI_Comm communicator,const Ioss::PropertyManager & props)320 ParallelDatabaseIO::ParallelDatabaseIO(Ioss::Region *region, const std::string &filename, 321 Ioss::DatabaseUsage db_usage, MPI_Comm communicator, 322 const Ioss::PropertyManager &props) 323 : Ioex::BaseDatabaseIO(region, filename, db_usage, communicator, props) 324 { 325 usingParallelIO = true; 326 if (!is_parallel_consistent()) { 327 std::ostringstream errmsg; 328 fmt::print(errmsg, 329 "ERROR: Parallel IO cannot be used in an application that is not guaranteeing " 330 "parallel consistent calls of the get and put field data functions.\n" 331 "The application created this database with a 'false' setting for the " 332 "isParallelConsistent property."); 333 IOSS_ERROR(errmsg); 334 } 335 336 if (!is_input()) { 337 // Check whether appending to or modifying existing file... 338 if (open_create_behavior() == Ioss::DB_APPEND || 339 open_create_behavior() == Ioss::DB_APPEND_GROUP || 340 open_create_behavior() == Ioss::DB_MODIFY) { 341 // Append to file if it already exists -- See if the file exists. 342 Ioss::FileInfo file = Ioss::FileInfo(get_filename()); 343 fileExists = file.exists(); 344 if (fileExists && myProcessor == 0) { 345 fmt::print(Ioss::WARNING(), 346 "Appending to existing database in parallel single-file " 347 "output mode is a new capability; please check results carefully. File '{}'", 348 get_filename()); 349 } 350 } 351 } 352 } 353 release_memory__()354 void ParallelDatabaseIO::release_memory__() 355 { 356 free_file_pointer(); 357 nodeMap.release_memory(); 358 edgeMap.release_memory(); 359 faceMap.release_memory(); 360 elemMap.release_memory(); 361 Ioss::Utils::clear(nodeOwningProcessor); 362 Ioss::Utils::clear(nodeGlobalImplicitMap); 363 Ioss::Utils::clear(elemGlobalImplicitMap); 364 nodeGlobalImplicitMapDefined = false; 365 elemGlobalImplicitMapDefined = false; 366 nodesetOwnedNodes.clear(); 367 try { 368 decomp.reset(); 369 } 370 catch (...) { 371 } 372 } 373 check_valid_file_ptr(bool write_message,std::string * error_msg,int * bad_count,bool abort_if_error)374 bool ParallelDatabaseIO::check_valid_file_ptr(bool write_message, std::string *error_msg, 375 int *bad_count, bool abort_if_error) const 376 { 377 // Check for valid exodus_file_ptr (valid >= 0; invalid < 0) 378 assert(isParallel); 379 int global_file_ptr = util().global_minmax(m_exodusFilePtr, Ioss::ParallelUtils::DO_MIN); 380 381 if (global_file_ptr < 0) { 382 if (write_message || error_msg != nullptr || bad_count != nullptr) { 383 Ioss::IntVector status; 384 util().all_gather(m_exodusFilePtr, status); 385 386 std::string open_create = is_input() ? "open input" : "create output"; 387 if (write_message || error_msg != nullptr) { 388 std::vector<size_t> procs; 389 for (int i = 0; i < util().parallel_size(); i++) { 390 if (status[i] < 0) { 391 procs.push_back(i); 392 } 393 } 394 std::string error_list = Ioss::Utils::format_id_list(procs, "--"); 395 // See which processors could not open/create the file... 396 std::ostringstream errmsg; 397 fmt::print(errmsg, "ERROR: Unable to {} exodus database file '{}' on processors:\n\t{}", 398 open_create, get_filename(), error_list); 399 fmt::print(errmsg, "\n"); 400 if (error_msg != nullptr) { 401 *error_msg = errmsg.str(); 402 } 403 if (write_message && myProcessor == 0) { 404 fmt::print(Ioss::OUTPUT(), "{}", errmsg.str()); 405 } 406 } 407 if (bad_count != nullptr) { 408 *bad_count = std::count_if(status.begin(), status.end(), [](int i) { return i < 0; }); 409 } 410 if (abort_if_error) { 411 std::ostringstream errmsg; 412 fmt::print(errmsg, "ERROR: Cannot {} file '{}'", open_create, get_filename()); 413 IOSS_ERROR(errmsg); 414 } 415 } 416 return false; 417 } 418 return true; 419 } 420 open_input_file(bool write_message,std::string * error_msg,int * bad_count,bool abort_if_error)421 bool ParallelDatabaseIO::open_input_file(bool write_message, std::string *error_msg, 422 int *bad_count, bool abort_if_error) const 423 { 424 int cpu_word_size = sizeof(double); 425 int io_word_size = 0; 426 float version; 427 428 int mode = exodusMode; 429 if (int_byte_size_api() == 8) { 430 mode |= EX_ALL_INT64_API; 431 } 432 433 #if defined EX_DISKLESS 434 // Experimental -- in memory read by netcdf library 435 if (properties.exists("MEMORY_READ")) { 436 mode |= EX_DISKLESS; 437 } 438 #endif 439 440 MPI_Info info = MPI_INFO_NULL; 441 std::string filename = get_filename(); 442 443 // See bug description in thread at 444 // https://www.open-mpi.org/community/lists/users/2015/01/26167.php and 445 // https://prod.sandia.gov/sierra-trac/ticket/14679 446 // Kluge is to set cwd to pathname, open file, then set cwd back to original. 447 // 448 // Since several different mpi implementations are based on the 449 // mpich code which introduced this bug, it has been difficult to 450 // create an ifdef'd version of the fix which is only applied to the 451 // buggy mpiio code. Therefore, we always do chdir call. Maybe in several 452 // years, we can remove this code and everything will work... 453 454 #ifndef _WIN32 455 Ioss::FileInfo file(filename); 456 std::string path = file.pathname(); 457 filename = file.tailname(); 458 char *current_cwd = getcwd(nullptr, 0); 459 chdir(path.c_str()); 460 #endif 461 462 bool do_timer = false; 463 Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer); 464 double t_begin = (do_timer ? Ioss::Utils::timer() : 0); 465 466 int app_opt_val = ex_opts(EX_VERBOSE); 467 m_exodusFilePtr = ex_open_par(filename.c_str(), EX_READ | mode, &cpu_word_size, &io_word_size, 468 &version, util().communicator(), info); 469 470 if (do_timer) { 471 double t_end = Ioss::Utils::timer(); 472 double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX); 473 if (myProcessor == 0) { 474 fmt::print(Ioss::DEBUG(), "File Open Time = {}\n", duration); 475 } 476 } 477 478 #ifndef _WIN32 479 chdir(current_cwd); 480 std::free(current_cwd); 481 #endif 482 483 bool is_ok = check_valid_file_ptr(write_message, error_msg, bad_count, abort_if_error); 484 485 if (is_ok) { 486 finalize_file_open(); 487 } 488 ex_opts(app_opt_val); // Reset back to what it was. 489 return is_ok; 490 } 491 handle_output_file(bool write_message,std::string * error_msg,int * bad_count,bool overwrite,bool abort_if_error)492 bool ParallelDatabaseIO::handle_output_file(bool write_message, std::string *error_msg, 493 int *bad_count, bool overwrite, 494 bool abort_if_error) const 495 { 496 // If 'overwrite' is false, we do not want to overwrite or clobber 497 // the output file if it already exists since the app might be 498 // reading the restart data from this file and then later 499 // clobbering it and then writing restart data to the same 500 // file. So, for output, we first check whether the file exists 501 // and if it it and is writable, assume that we can later create a 502 // new or append to existing file. 503 504 // if 'overwrite' is true, then clobber/append 505 506 if (!overwrite) { 507 // check if file exists and is writeable. If so, return true. 508 // Only need to check on processor 0 509 int int_is_ok = 0; 510 if (myProcessor == 0) { 511 Ioss::FileInfo file(get_filename()); 512 int_is_ok = file.exists() && file.is_writable() ? 1 : 0; 513 } 514 MPI_Bcast(&int_is_ok, 1, MPI_INT, 0, util().communicator()); 515 516 if (int_is_ok == 1) { 517 // Note that at this point, we cannot totally guarantee that 518 // we will be able to create the file when needed, but we have 519 // a pretty good chance. We can't guarantee creation without 520 // creating and the app (or calling function) doesn't want us to overwrite... 521 return true; 522 } 523 // File doesn't exist, so fall through and try to 524 // create file since we won't be overwriting anything... 525 } 526 527 int cpu_word_size = sizeof(double); 528 int io_word_size = 0; 529 float version; 530 531 int mode = exodusMode; 532 if (int_byte_size_api() == 8) { 533 mode |= EX_ALL_INT64_API; 534 } 535 536 #if defined EX_DISKLESS 537 // Experimental -- in memory write by netcdf library 538 if (properties.exists("MEMORY_WRITE")) { 539 mode |= EX_DISKLESS; 540 } 541 #endif 542 543 MPI_Info info = MPI_INFO_NULL; 544 int app_opt_val = ex_opts(EX_VERBOSE); 545 Ioss::DatabaseIO::openDatabase__(); 546 547 std::string filename = get_dwname(); 548 549 #ifndef _WIN32 550 Ioss::FileInfo file(filename); 551 std::string path = file.pathname(); 552 filename = file.tailname(); 553 char *current_cwd = getcwd(nullptr, 0); 554 chdir(path.c_str()); 555 #endif 556 557 bool do_timer = false; 558 Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer); 559 double t_begin = (do_timer ? Ioss::Utils::timer() : 0); 560 561 if (fileExists) { 562 m_exodusFilePtr = ex_open_par(filename.c_str(), EX_WRITE | mode, &cpu_word_size, 563 &io_word_size, &version, util().communicator(), info); 564 } 565 else { 566 // If the first write for this file, create it... 567 if (int_byte_size_api() == 8) { 568 // Check whether client actually wants 4-byte output on db 569 // - If they specified INTEGER_SIZE_DB and the size isn't 8, 570 // then don't change mode and use the default 4-byte output. 571 if (properties.exists("INTEGER_SIZE_DB")) { 572 if (properties.get("INTEGER_SIZE_DB").get_int() == 8) { 573 mode |= EX_ALL_INT64_DB; 574 } 575 } 576 else { 577 mode |= EX_ALL_INT64_DB; 578 } 579 } 580 m_exodusFilePtr = ex_create_par(filename.c_str(), mode, &cpu_word_size, &dbRealWordSize, 581 util().communicator(), info); 582 } 583 584 if (do_timer) { 585 double t_end = Ioss::Utils::timer(); 586 double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX); 587 std::string open_create = fileExists ? "Open" : "Create"; 588 if (myProcessor == 0) { 589 fmt::print(Ioss::DEBUG(), "File {} Time = {}\n", open_create, duration); 590 } 591 } 592 593 #ifndef _WIN32 594 chdir(current_cwd); 595 std::free(current_cwd); 596 #endif 597 598 bool is_ok = check_valid_file_ptr(write_message, error_msg, bad_count, abort_if_error); 599 600 if (is_ok) { 601 ex_set_max_name_length(m_exodusFilePtr, maximumNameLength); 602 603 // Check properties handled post-create/open... 604 if (properties.exists("COMPRESSION_METHOD")) { 605 auto method = properties.get("COMPRESSION_METHOD").get_string(); 606 method = Ioss::Utils::lowercase(method); 607 ex_compression_type exo_method = EX_COMPRESS_ZLIB; 608 if (method == "zlib" || method == "libz" || method == "gzip") { 609 exo_method = EX_COMPRESS_ZLIB; 610 } 611 else if (method == "szip") { 612 #if !defined(NC_HAS_SZIP_WRITE) 613 #define NC_HAS_SZIP_WRITE 0 614 #endif 615 #if NC_HAS_SZIP_WRITE 616 exo_method = EX_COMPRESS_SZIP; 617 #else 618 if (myProcessor == 0) { 619 fmt::print(Ioss::WARNING(), "The NetCDF library does not have SZip compression enabled." 620 " 'zlib' will be used instead.\n\n"); 621 } 622 #endif 623 } 624 else { 625 if (myProcessor == 0) { 626 fmt::print(Ioss::WARNING(), 627 "Unrecognized compression method specified: '{}'." 628 " 'zlib' will be used instead.\n\n", 629 method); 630 } 631 } 632 ex_set_option(m_exodusFilePtr, EX_OPT_COMPRESSION_TYPE, exo_method); 633 } 634 if (properties.exists("COMPRESSION_LEVEL")) { 635 int comp_level = properties.get("COMPRESSION_LEVEL").get_int(); 636 ex_set_option(m_exodusFilePtr, EX_OPT_COMPRESSION_LEVEL, comp_level); 637 } 638 if (properties.exists("COMPRESSION_SHUFFLE")) { 639 int shuffle = properties.get("COMPRESSION_SHUFFLE").get_int(); 640 ex_set_option(m_exodusFilePtr, EX_OPT_COMPRESSION_SHUFFLE, shuffle); 641 } 642 } 643 ex_opts(app_opt_val); // Reset back to what it was. 644 return is_ok; 645 } 646 get_file_pointer()647 int ParallelDatabaseIO::get_file_pointer() const 648 { 649 return Ioex::BaseDatabaseIO::get_file_pointer(); 650 } 651 free_file_pointer()652 int ParallelDatabaseIO::free_file_pointer() const 653 { 654 int flag; 655 MPI_Initialized(&flag); 656 if (flag == 0) { 657 std::ostringstream errmsg; 658 fmt::print(errmsg, "ERROR: MPI is not initialized."); 659 IOSS_ERROR(errmsg); 660 } 661 662 // Make sure all file pointers are valid... 663 int fp_min = util().global_minmax(m_exodusFilePtr, Ioss::ParallelUtils::DO_MIN); 664 int fp_max = util().global_minmax(m_exodusFilePtr, Ioss::ParallelUtils::DO_MAX); 665 if (fp_min != fp_max && fp_min < 0) { 666 std::ostringstream errmsg; 667 fmt::print(errmsg, "ERROR: Inconsistent file pointer values."); 668 IOSS_ERROR(errmsg); 669 } 670 return Ioex::BaseDatabaseIO::free_file_pointer(); 671 } 672 read_meta_data__()673 void ParallelDatabaseIO::read_meta_data__() 674 { 675 int exoid = get_file_pointer(); // get_file_pointer() must be called first. 676 677 // APPENDING: 678 // If parallel (single file, not fpp), we have assumptions 679 // that the writing process (ranks, mesh, decomp, vars) is the 680 // same for the original run that created this database and 681 // for this run which is appending to the database so the 682 // defining of the output database should be the same except 683 // we don't write anything since it is already there. We do 684 // need the number of steps though... 685 if (open_create_behavior() == Ioss::DB_APPEND) { 686 get_step_times__(); 687 return; 688 } 689 690 if (int_byte_size_api() == 8) { 691 decomp = std::unique_ptr<DecompositionDataBase>( 692 new DecompositionData<int64_t>(properties, util().communicator())); 693 } 694 else { 695 decomp = std::unique_ptr<DecompositionDataBase>( 696 new DecompositionData<int>(properties, util().communicator())); 697 } 698 assert(decomp != nullptr); 699 decomp->decompose_model(exoid); 700 701 read_region(); 702 get_elemblocks(); 703 704 get_step_times__(); 705 706 get_nodeblocks(); 707 get_edgeblocks(); 708 get_faceblocks(); 709 710 check_side_topology(); 711 712 get_nodesets(); 713 get_sidesets(); 714 #if 0 715 get_edgesets(); 716 get_facesets(); 717 get_elemsets(); 718 #endif 719 720 get_commsets(); 721 722 // Add assemblies now that all entities should be defined... consistent across processors 723 // (metadata) 724 get_assemblies(); 725 726 get_blobs(); 727 728 handle_groups(); 729 730 add_region_fields(); 731 732 if (!is_input() && open_create_behavior() == Ioss::DB_APPEND) { 733 get_map(EX_NODE_BLOCK); 734 get_map(EX_ELEM_BLOCK); 735 } 736 } 737 read_region()738 void ParallelDatabaseIO::read_region() 739 { 740 // Add properties and fields to the 'owning' region. 741 // Also defines member variables of this class... 742 ex_init_params info{}; 743 int error = ex_get_init_ext(get_file_pointer(), &info); 744 if (error < 0) { 745 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 746 } 747 748 spatialDimension = decomp->spatial_dimension(); 749 nodeCount = decomp->ioss_node_count(); 750 edgeCount = 0; 751 faceCount = 0; 752 elementCount = decomp->ioss_elem_count(); 753 754 m_groupCount[EX_NODE_BLOCK] = 1; 755 m_groupCount[EX_EDGE_BLOCK] = info.num_edge_blk; 756 m_groupCount[EX_FACE_BLOCK] = info.num_face_blk; 757 m_groupCount[EX_ELEM_BLOCK] = info.num_elem_blk; 758 759 m_groupCount[EX_NODE_SET] = info.num_node_sets; 760 m_groupCount[EX_EDGE_SET] = info.num_edge_sets; 761 m_groupCount[EX_FACE_SET] = info.num_face_sets; 762 m_groupCount[EX_ELEM_SET] = info.num_elem_sets; 763 764 m_groupCount[EX_SIDE_SET] = info.num_side_sets; 765 m_groupCount[EX_ASSEMBLY] = info.num_assembly; 766 m_groupCount[EX_BLOB] = info.num_blob; 767 768 // Checks: node, element, blocks > 0; warning if == 0; error if < 0 769 check_valid_values(); 770 771 Ioss::Region *this_region = get_region(); 772 773 // See if any coordinate frames exist on mesh. If so, define them on region. 774 Ioex::add_coordinate_frames(get_file_pointer(), this_region); 775 776 this_region->property_add( 777 Ioss::Property("global_node_count", static_cast<int64_t>(decomp->global_node_count()))); 778 this_region->property_add( 779 Ioss::Property("global_element_count", static_cast<int64_t>(decomp->global_elem_count()))); 780 781 this_region->property_add(Ioss::Property(std::string("title"), info.title)); 782 783 // Get QA records from database and add to qaRecords... 784 int num_qa = ex_inquire_int(get_file_pointer(), EX_INQ_QA); 785 if (num_qa > 0) { 786 struct qa_element 787 { 788 char *qa_record[1][4]; 789 }; 790 791 auto qa = new qa_element[num_qa]; 792 for (int i = 0; i < num_qa; i++) { 793 for (int j = 0; j < 4; j++) { 794 qa[i].qa_record[0][j] = new char[MAX_STR_LENGTH + 1]; 795 } 796 } 797 798 ex_get_qa(get_file_pointer(), qa[0].qa_record); 799 for (int i = 0; i < num_qa; i++) { 800 add_qa_record(qa[i].qa_record[0][0], qa[i].qa_record[0][1], qa[i].qa_record[0][2], 801 qa[i].qa_record[0][3]); 802 } 803 for (int i = 0; i < num_qa; i++) { 804 for (int j = 0; j < 4; j++) { 805 delete[] qa[i].qa_record[0][j]; 806 } 807 } 808 delete[] qa; 809 } 810 811 // Get information records from database and add to informationRecords... 812 int num_info = ex_inquire_int(get_file_pointer(), EX_INQ_INFO); 813 if (num_info > 0) { 814 char **info_rec = Ioss::Utils::get_name_array( 815 num_info, max_line_length); // 'total_lines' pointers to char buffers 816 ex_get_info(get_file_pointer(), info_rec); 817 for (int i = 0; i < num_info; i++) { 818 add_information_record(info_rec[i]); 819 } 820 Ioss::Utils::delete_name_array(info_rec, num_info); 821 } 822 } 823 get_step_times__()824 void ParallelDatabaseIO::get_step_times__() 825 { 826 double last_time = DBL_MAX; 827 int timestep_count = 0; 828 std::vector<double> tsteps(0); 829 830 { 831 timestep_count = ex_inquire_int(get_file_pointer(), EX_INQ_TIME); 832 if (timestep_count <= 0) { 833 return; 834 } 835 836 // For an exodusII file, timesteps are global and are stored in the region. 837 // Read the timesteps and add to the region 838 tsteps.resize(timestep_count); 839 int error = ex_get_all_times(get_file_pointer(), tsteps.data()); 840 if (error < 0) { 841 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 842 } 843 844 // See if the "last_written_time" attribute exists and if it 845 // does, check that it matches the largest time in 'tsteps'. 846 Ioex::read_last_time_attribute(get_file_pointer(), &last_time); 847 } 848 849 // Only add states that are less than or equal to the 850 // 'last_time' value which is either DBL_MAX or the value of 851 // the last time successfully written to the database and 852 // flushed to disk. This is used to avoid corrupt data arising 853 // from a job that crashed during the writing of the last step 854 // on the database. Output a warning message if there is 855 // potentially corrupt data on the database... 856 857 // Check whether user or application wants to limit the times even further... 858 // One use case is that job is restarting at a time prior to what has been 859 // written to the results file, so want to start appending after 860 // restart time instead of at end time on database. 861 int max_step = timestep_count; 862 if (properties.exists("APPEND_OUTPUT_AFTER_STEP")) { 863 max_step = properties.get("APPEND_OUTPUT_AFTER_STEP").get_int(); 864 } 865 if (max_step > timestep_count) { 866 max_step = timestep_count; 867 } 868 869 double max_time = std::numeric_limits<double>::max(); 870 if (properties.exists("APPEND_OUTPUT_AFTER_TIME")) { 871 max_time = properties.get("APPEND_OUTPUT_AFTER_TIME").get_real(); 872 } 873 if (last_time > max_time) { 874 last_time = max_time; 875 } 876 877 Ioss::Region *this_region = get_region(); 878 for (int i = 0; i < max_step; i++) { 879 if (tsteps[i] <= last_time) { 880 this_region->add_state(tsteps[i] * timeScaleFactor); 881 } 882 else { 883 if (myProcessor == 0 && max_time == std::numeric_limits<double>::max()) { 884 // NOTE: Don't want to warn on all processors if there are 885 // corrupt steps on all databases, but this will only print 886 // a warning if there is a corrupt step on processor 887 // 0... Need better warnings which won't overload in the 888 // worst case... 889 fmt::print( 890 Ioss::WARNING(), 891 "Skipping step {:L} at time {} in database file\n\t{}.\nThe data for that step " 892 "is possibly corrupt.\n", 893 i + 1, tsteps[i], get_filename()); 894 } 895 } 896 } 897 } 898 get_map(ex_entity_type type)899 const Ioss::Map &ParallelDatabaseIO::get_map(ex_entity_type type) const 900 { 901 switch (type) { 902 case EX_NODE_BLOCK: 903 case EX_NODE_SET: { 904 assert(decomp != nullptr); 905 size_t offset = decomp->decomp_node_offset(); 906 size_t count = decomp->decomp_node_count(); 907 return get_map(nodeMap, nodeCount, offset, count, EX_NODE_MAP, EX_INQ_NODE_MAP); 908 } 909 case EX_ELEM_BLOCK: 910 case EX_ELEM_SET: { 911 assert(decomp != nullptr); 912 size_t offset = decomp->decomp_elem_offset(); 913 size_t count = decomp->decomp_elem_count(); 914 return get_map(elemMap, elementCount, offset, count, EX_ELEM_MAP, EX_INQ_ELEM_MAP); 915 } 916 917 case EX_FACE_BLOCK: 918 case EX_FACE_SET: return get_map(faceMap, faceCount, 0, 0, EX_FACE_MAP, EX_INQ_FACE_MAP); 919 920 case EX_EDGE_BLOCK: 921 case EX_EDGE_SET: return get_map(edgeMap, edgeCount, 0, 0, EX_EDGE_MAP, EX_INQ_EDGE_MAP); 922 923 default: 924 std::ostringstream errmsg; 925 fmt::print(errmsg, "INTERNAL ERROR: Invalid map type. " 926 "Something is wrong in the Ioex::ParallelDatabaseIO::get_map() function. " 927 "Please report.\n"); 928 IOSS_ERROR(errmsg); 929 } 930 } 931 get_map(Ioss::Map & entity_map,int64_t entity_count,int64_t file_offset,int64_t file_count,ex_entity_type entity_type,ex_inquiry inquiry_type)932 const Ioss::Map &ParallelDatabaseIO::get_map(Ioss::Map &entity_map, int64_t entity_count, 933 int64_t file_offset, int64_t file_count, 934 ex_entity_type entity_type, 935 ex_inquiry inquiry_type) const 936 { 937 // Allocate space for node number map and read it in... 938 // Can be called multiple times, allocate 1 time only 939 if (entity_map.map().empty()) { 940 entity_map.set_size(entity_count); 941 942 if (is_input()) { 943 Ioss::MapContainer file_data(file_count); 944 int error = 0; 945 // Check whether there is a "original_global_id_map" map on 946 // the database. If so, use it instead of the "node_num_map". 947 bool map_read = false; 948 int map_count = ex_inquire_int(get_file_pointer(), inquiry_type); 949 if (map_count > 0) { 950 char **names = Ioss::Utils::get_name_array(map_count, maximumNameLength); 951 int ierr = ex_get_names(get_file_pointer(), entity_type, names); 952 if (ierr < 0) { 953 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 954 } 955 956 if (map_count == 1 && Ioss::Utils::str_equal(names[0], "original_global_id_map")) { 957 if (int_byte_size_api() == 8) { 958 error = ex_get_partial_num_map(get_file_pointer(), entity_type, 1, file_offset + 1, 959 file_count, file_data.data()); 960 } 961 else { 962 // Ioss stores as 64-bit, read as 32-bit and copy over... 963 Ioss::IntVector tmp_map(file_count); 964 error = ex_get_partial_num_map(get_file_pointer(), entity_type, 1, file_offset + 1, 965 file_count, tmp_map.data()); 966 std::copy(tmp_map.begin(), tmp_map.end(), file_data.begin()); 967 } 968 if (error >= 0) { 969 map_read = true; 970 } 971 } 972 Ioss::Utils::delete_name_array(names, map_count); 973 } 974 975 if (!map_read) { 976 if (int_byte_size_api() == 8) { 977 error = ex_get_partial_id_map(get_file_pointer(), entity_type, file_offset + 1, 978 file_count, file_data.data()); 979 } 980 else { 981 // Ioss stores as 64-bit, read as 32-bit and copy over... 982 Ioss::IntVector tmp_map(file_count); 983 error = ex_get_partial_id_map(get_file_pointer(), entity_type, file_offset + 1, 984 file_count, tmp_map.data()); 985 std::copy(tmp_map.begin(), tmp_map.end(), file_data.begin()); 986 } 987 } 988 989 if (error >= 0) { 990 if (entity_type == EX_NODE_MAP) { 991 decomp->communicate_node_data(file_data.data(), &entity_map.map()[1], 1); 992 } 993 else if (entity_type == EX_ELEM_MAP) { 994 decomp->communicate_element_data(file_data.data(), &entity_map.map()[1], 1); 995 } 996 } 997 else { 998 // Clear out the vector... 999 Ioss::MapContainer().swap(entity_map.map()); 1000 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1001 } 1002 1003 // Check for sequential node map. 1004 // If not, build the reverse G2L node map... 1005 entity_map.is_sequential(true); 1006 entity_map.build_reverse_map(); 1007 } 1008 else { 1009 // Output database; entity_map.map() not set yet... Build a default map. 1010 entity_map.set_default(entity_count); 1011 } 1012 } 1013 return entity_map; 1014 } 1015 get_elemblocks()1016 void ParallelDatabaseIO::get_elemblocks() { get_blocks(EX_ELEM_BLOCK, 0, "block"); } 1017 get_faceblocks()1018 void ParallelDatabaseIO::get_faceblocks() 1019 { 1020 // get_blocks(EX_FACE_BLOCK, 1, "faceblock"); 1021 } 1022 get_edgeblocks()1023 void ParallelDatabaseIO::get_edgeblocks() 1024 { 1025 // get_blocks(EX_EDGE_BLOCK, 2, "edgeblock"); 1026 } 1027 get_blocks(ex_entity_type entity_type,int rank_offset,const std::string & basename)1028 void ParallelDatabaseIO::get_blocks(ex_entity_type entity_type, int rank_offset, 1029 const std::string &basename) 1030 { 1031 // Attributes of an X block are: (X = element, face, or edge) 1032 // -- id 1033 // -- name 1034 // -- X type 1035 // -- number of Xs 1036 // -- number of attributes per X 1037 // -- number of nodes per X (derivable from type) 1038 // -- number of faces per X (derivable from type) 1039 // -- number of edges per X (derivable from type) 1040 1041 // In a parallel execution, it is possible that an X block will have 1042 // no Xs on a particular processor... 1043 1044 // NOTE: This routine may be called multiple times on a single database. 1045 // make sure it is not dependent on being called one time only... 1046 1047 // Get exodusII X block metadata 1048 if (m_groupCount[entity_type] == 0) { 1049 return; 1050 } 1051 1052 assert(entity_type == EX_ELEM_BLOCK); 1053 1054 Ioss::Int64Vector X_block_ids(m_groupCount[entity_type]); 1055 int used_blocks = 0; 1056 1057 int error; 1058 if ((ex_int64_status(get_file_pointer()) & EX_IDS_INT64_API) != 0) { 1059 error = ex_get_ids(get_file_pointer(), entity_type, X_block_ids.data()); 1060 } 1061 else { 1062 Ioss::IntVector tmp_set_ids(X_block_ids.size()); 1063 error = ex_get_ids(get_file_pointer(), entity_type, tmp_set_ids.data()); 1064 if (error >= 0) { 1065 std::copy(tmp_set_ids.begin(), tmp_set_ids.end(), X_block_ids.begin()); 1066 } 1067 } 1068 if (error < 0) { 1069 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1070 } 1071 1072 int nvar = std::numeric_limits<int>::max(); // Number of 'block' vars on database. Used to 1073 // skip querying if none. 1074 int nmap = std::numeric_limits<int>::max(); // Number of 'block' maps on database. Used to 1075 // skip querying if none. 1076 for (int iblk = 0; iblk < m_groupCount[entity_type]; iblk++) { 1077 if (decomp->el_blocks[iblk].global_count() == 0) { 1078 continue; 1079 } 1080 1081 int64_t id = decomp->el_blocks[iblk].id(); 1082 1083 bool db_has_name = false; 1084 std::string alias = Ioss::Utils::encode_entity_name(basename, id); 1085 std::string block_name; 1086 if (ignore_database_names()) { 1087 block_name = alias; 1088 } 1089 else { 1090 block_name = Ioex::get_entity_name(get_file_pointer(), entity_type, id, basename, 1091 maximumNameLength, db_has_name); 1092 } 1093 if (get_use_generic_canonical_name()) { 1094 std::swap(block_name, alias); 1095 } 1096 1097 std::string save_type = decomp->el_blocks[iblk].topologyType; 1098 std::string type = Ioss::Utils::fixup_type(decomp->el_blocks[iblk].topologyType, 1099 decomp->el_blocks[iblk].nodesPerEntity, 1100 spatialDimension - rank_offset); 1101 1102 Ioss::EntityBlock *io_block = nullptr; 1103 if (entity_type == EX_ELEM_BLOCK) { 1104 auto *eblock = 1105 new Ioss::ElementBlock(this, block_name, type, decomp->el_blocks[iblk].ioss_count()); 1106 io_block = eblock; 1107 io_block->property_add(Ioss::Property("id", id)); 1108 io_block->property_add(Ioss::Property("guid", util().generate_guid(id))); 1109 io_block->property_add(Ioss::Property("iblk", iblk)); // Sequence in decomp. 1110 1111 if (db_has_name) { 1112 std::string *db_name = &block_name; 1113 if (get_use_generic_canonical_name()) { 1114 db_name = &alias; 1115 } 1116 io_block->property_add(Ioss::Property("db_name", *db_name)); 1117 } 1118 get_region()->add(eblock); 1119 #if 0 1120 } else if (entity_type == EX_FACE_BLOCK) { 1121 Ioss::FaceBlock *fblock = new Ioss::FaceBlock(this, block_name, type, block.num_entry); 1122 io_block = fblock; 1123 io_block->property_add(Ioss::Property("id", id)); 1124 io_block->property_add(Ioss::Property("guid", util().generate_guid(id))); 1125 if (db_has_name) { 1126 std::string *db_name = &block_name; 1127 if (get_use_generic_canonical_name()) { 1128 db_name = &alias; 1129 } 1130 io_block->property_add(Ioss::Property("db_name", *db_name)); 1131 } 1132 get_region()->add(fblock); 1133 } else if (entity_type == EX_EDGE_BLOCK) { 1134 Ioss::EdgeBlock *eblock = new Ioss::EdgeBlock(this, block_name, type, block.num_entry); 1135 io_block = eblock; 1136 io_block->property_add(Ioss::Property("id", id)); 1137 io_block->property_add(Ioss::Property("guid", util().generate_guid(id))); 1138 if (db_has_name) { 1139 std::string *db_name = &block_name; 1140 if (get_use_generic_canonical_name()) { 1141 db_name = &alias; 1142 } 1143 io_block->property_add(Ioss::Property("db_name", *db_name)); 1144 } 1145 get_region()->add(eblock); 1146 #endif 1147 } 1148 else { 1149 std::ostringstream errmsg; 1150 fmt::print(errmsg, "ERROR: Invalid type in get_blocks()"); 1151 IOSS_ERROR(errmsg); 1152 } 1153 1154 #if 0 1155 // See which connectivity options were defined for this block. 1156 // X -> Node is always defined. 1157 // X -> Face? 1158 if (block.num_faces_per_entry > 0 && rank_offset < 1) { 1159 std::string storage = "Real["+std::to_string(block.num_faces_per_entry)+"]"; 1160 io_block->field_add(Ioss::Field("connectivity_face", 1161 io_block->field_int_type(), storage, Ioss::Field::MESH)); 1162 } 1163 // X -> Edge? 1164 if (block.num_edges_per_entry > 0 && rank_offset < 2) { 1165 std::string storage = "Real["+std::to_string(block.num_edges_per_entry)+"]"; 1166 io_block->field_add(Ioss::Field("connectivity_edge", 1167 io_block->field_int_type(), storage, Ioss::Field::MESH)); 1168 } 1169 #endif 1170 1171 // Maintain block order on output database... 1172 io_block->property_add(Ioss::Property("original_block_order", used_blocks++)); 1173 1174 if (save_type != "null" && save_type != "") { 1175 io_block->property_update("original_topology_type", save_type); 1176 } 1177 1178 io_block->property_add(Ioss::Property( 1179 "global_entity_count", static_cast<int64_t>(decomp->el_blocks[iblk].ioss_count()))); 1180 1181 if (block_name != alias) { 1182 get_region()->add_alias(block_name, alias); 1183 } 1184 1185 // Check for additional variables. 1186 add_attribute_fields(entity_type, io_block, decomp->el_blocks[iblk].attributeCount, type); 1187 if (nvar > 0) { 1188 nvar = add_results_fields(entity_type, io_block, iblk); 1189 } 1190 1191 if (entity_type == EX_ELEM_BLOCK) { 1192 if (nmap > 0) { 1193 nmap = 1194 Ioex::add_map_fields(get_file_pointer(), dynamic_cast<Ioss::ElementBlock *>(io_block), 1195 decomp->el_blocks[iblk].ioss_count(), maximumNameLength); 1196 } 1197 assert(blockOmissions.empty() || blockInclusions.empty()); // Only one can be non-empty 1198 1199 // Handle all block omissions or inclusions... 1200 // This only affects the generation of surfaces... 1201 if (!blockOmissions.empty()) { 1202 for (const auto &name : blockOmissions) { 1203 auto block = get_region()->get_element_block(name); 1204 if (block) { 1205 block->property_add(Ioss::Property(std::string("omitted"), 1)); 1206 } 1207 } 1208 } 1209 1210 if (!blockInclusions.empty()) { 1211 const auto &blocks = get_region()->get_element_blocks(); 1212 for (auto &block : blocks) { 1213 block->property_add(Ioss::Property(std::string("omitted"), 1)); 1214 } 1215 1216 // Now, erase the property on any blocks in the inclusion list... 1217 for (const auto &name : blockInclusions) { 1218 auto block = get_region()->get_element_block(name); 1219 if (block != nullptr) { 1220 block->property_erase("omitted"); 1221 } 1222 } 1223 } 1224 } 1225 } 1226 } 1227 compute_node_status()1228 void ParallelDatabaseIO::compute_node_status() const 1229 { 1230 // Create a field for all nodes in the model indicating 1231 // the connectivity 'status' of the node. The status values are: 1232 // 0 -- node not connected to any elements 1233 // 1 -- node only connected to omitted elements 1234 // 2 -- node only connected to active elements 1235 // 3 -- node at border of active and omitted elements. 1236 1237 /// \todo Get working for parallel... 1238 1239 if (nodeConnectivityStatusCalculated) { 1240 return; 1241 } 1242 1243 nodeConnectivityStatus.resize(nodeCount); 1244 1245 const Ioss::ElementBlockContainer &element_blocks = get_region()->get_element_blocks(); 1246 assert(Ioss::Utils::check_block_order(element_blocks)); 1247 1248 for (Ioss::ElementBlock *block : element_blocks) { 1249 unsigned char status = 2; 1250 if (Ioss::Utils::block_is_omitted(block)) { 1251 status = 1; 1252 } 1253 1254 int64_t id = block->get_property("id").get_int(); 1255 int element_nodes = block->topology()->number_nodes(); 1256 int64_t my_element_count = block->entity_count(); 1257 int order = block->get_property("iblk").get_int(); 1258 if (int_byte_size_api() == 8) { 1259 std::vector<int64_t> conn(my_element_count * element_nodes); 1260 decomp->get_block_connectivity(get_file_pointer(), conn.data(), id, order, element_nodes); 1261 for (auto node : conn) { 1262 nodeConnectivityStatus[node - 1] |= status; 1263 } 1264 } 1265 else { 1266 std::vector<int> conn(my_element_count * element_nodes); 1267 decomp->get_block_connectivity(get_file_pointer(), conn.data(), id, order, element_nodes); 1268 for (auto node : conn) { 1269 nodeConnectivityStatus[node - 1] |= status; 1270 } 1271 } 1272 } 1273 nodeConnectivityStatusCalculated = true; 1274 } 1275 get_sidesets()1276 void ParallelDatabaseIO::get_sidesets() 1277 { 1278 // This function creates all sidesets (surfaces) for a 1279 // model. Note that a sideset contains 1 or more sideblocks 1280 // which are homogeneous (same topology). In serial execution, 1281 // this is fairly straightforward since there are no null sets and 1282 // we have all the information we need. (...except see below for 1283 // surface evolution). 1284 // 1285 // However, in a parallel execution, we have the possibility that a 1286 // side set will have no sides or distribution factors on 1287 // a particular processor. We then don't know the block topology of 1288 // the block(s) contained in this set. We could do some 1289 // communication and get a good idea of the topologies that are in 1290 // the set. 1291 1292 if (m_groupCount[EX_SIDE_SET] > 0) { 1293 check_side_topology(); 1294 1295 // Get exodusII sideset metadata 1296 1297 // Get the names (may not exist) of all sidesets and see if they are actually 1298 // side "blocks" (perhaps written by IO system for a restart). In that case, 1299 // they were split by a previous run and we need to reconstruct the side "set" 1300 // that may contain one or more of them. 1301 Ioex::SideSetMap fs_map; 1302 Ioex::SideSetSet fs_set; 1303 1304 { 1305 int error; 1306 1307 for (const auto &ss : decomp->side_sets) { 1308 int64_t id = ss.id(); 1309 std::vector<char> ss_name(maximumNameLength + 1); 1310 error = ex_get_name(get_file_pointer(), EX_SIDE_SET, id, ss_name.data()); 1311 if (error < 0) { 1312 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1313 } 1314 if (ss_name[0] != '\0') { 1315 Ioss::Utils::fixup_name(ss_name.data()); 1316 Ioex::decode_surface_name(fs_map, fs_set, ss_name.data()); 1317 } 1318 } 1319 } 1320 1321 // Create sidesets for each entry in the fs_set... These are the 1322 // sidesets which were probably written by a previous run of the 1323 // IO system and are already split into homogeneous pieces... 1324 { 1325 for (auto &fs_name : fs_set) { 1326 auto *side_set = new Ioss::SideSet(this, fs_name); 1327 get_region()->add(side_set); 1328 int64_t id = Ioex::extract_id(fs_name); 1329 if (id > 0) { 1330 side_set->property_add(Ioss::Property("id", id)); 1331 side_set->property_add(Ioss::Property("guid", util().generate_guid(id))); 1332 } 1333 } 1334 } 1335 1336 for (int iss = 0; iss < m_groupCount[EX_SIDE_SET]; iss++) { 1337 int64_t id = decomp->side_sets[iss].id(); 1338 std::string sid; 1339 Ioex::TopologyMap topo_map; 1340 Ioex::TopologyMap side_map; // Used to determine side consistency 1341 1342 Ioss::SurfaceSplitType split_type = splitType; 1343 std::string side_set_name; 1344 Ioss::SideSet * side_set = nullptr; 1345 1346 { 1347 bool db_has_name = false; 1348 std::string alias = Ioss::Utils::encode_entity_name("surface", id); 1349 if (ignore_database_names()) { 1350 side_set_name = alias; 1351 } 1352 else { 1353 side_set_name = Ioex::get_entity_name(get_file_pointer(), EX_SIDE_SET, id, "surface", 1354 maximumNameLength, db_has_name); 1355 } 1356 1357 if (side_set_name == "universal_sideset") { 1358 split_type = Ioss::SPLIT_BY_DONT_SPLIT; 1359 } 1360 1361 bool in_fs_map = false; 1362 auto FSM = fs_map.find(side_set_name); 1363 if (FSM != fs_map.end()) { 1364 in_fs_map = true; 1365 std::string efs_name = (*FSM).second; 1366 side_set = get_region()->get_sideset(efs_name); 1367 Ioss::Utils::check_non_null(side_set, "sideset", efs_name, __func__); 1368 } 1369 else { 1370 if (get_use_generic_canonical_name()) { 1371 std::swap(side_set_name, alias); 1372 } 1373 side_set = new Ioss::SideSet(this, side_set_name); 1374 side_set->property_add(Ioss::Property("id", id)); 1375 side_set->property_add(Ioss::Property("guid", util().generate_guid(id))); 1376 if (db_has_name) { 1377 std::string *db_name = &side_set_name; 1378 if (get_use_generic_canonical_name()) { 1379 db_name = &alias; 1380 } 1381 side_set->property_add(Ioss::Property("db_name", *db_name)); 1382 } 1383 get_region()->add(side_set); 1384 1385 get_region()->add_alias(side_set_name, alias); 1386 get_region()->add_alias(side_set_name, Ioss::Utils::encode_entity_name("sideset", id)); 1387 } 1388 1389 // split_type = SPLIT_BY_ELEMENT_BLOCK; 1390 // split_type = SPLIT_BY_TOPOLOGIES; 1391 // split_type = SPLIT_BY_DONT_SPLIT; 1392 1393 // Determine how many side blocks compose this side set. 1394 1395 int64_t number_sides = decomp->side_sets[iss].ioss_count(); 1396 // FIXME: Support- number_distribution_factors = decomp->side_sets[iss].df_count(); 1397 1398 Ioss::Int64Vector element(number_sides); 1399 Ioss::Int64Vector sides(number_sides); 1400 1401 // Easier below here if the element and sides are a known 64-bit size... 1402 // Kluge here to do that... 1403 if (int_byte_size_api() == 4) { 1404 Ioss::Field side_field("sides", Ioss::Field::INTEGER, IOSS_SCALAR(), Ioss::Field::MESH, 1405 number_sides); 1406 Ioss::Field elem_field("ids_raw", Ioss::Field::INTEGER, IOSS_SCALAR(), 1407 Ioss::Field::MESH, number_sides); 1408 1409 Ioss::IntVector e32(number_sides); 1410 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, e32.data()); 1411 std::copy(e32.begin(), e32.end(), element.begin()); 1412 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, e32.data()); 1413 std::copy(e32.begin(), e32.end(), sides.begin()); 1414 } 1415 else { 1416 Ioss::Field side_field("sides", Ioss::Field::INT64, IOSS_SCALAR(), Ioss::Field::MESH, 1417 number_sides); 1418 Ioss::Field elem_field("ids_raw", Ioss::Field::INT64, IOSS_SCALAR(), Ioss::Field::MESH, 1419 number_sides); 1420 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 1421 element.data()); 1422 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, sides.data()); 1423 } 1424 1425 if (!blockOmissions.empty() || !blockInclusions.empty()) { 1426 Ioex::filter_element_list(get_region(), element, sides, true); 1427 number_sides = element.size(); 1428 assert(element.size() == sides.size()); 1429 } 1430 1431 if (split_type == Ioss::SPLIT_BY_TOPOLOGIES && sideTopology.size() == 1) { 1432 // There is only one side type for all elements in the model 1433 topo_map[std::make_pair(sideTopology[0].first->name(), sideTopology[0].second)] = 1434 number_sides; 1435 } 1436 else if (split_type == Ioss::SPLIT_BY_DONT_SPLIT) { 1437 const Ioss::ElementTopology *mixed_topo = Ioss::ElementTopology::factory("unknown"); 1438 topo_map[std::make_pair(std::string("unknown"), mixed_topo)] = number_sides; 1439 } 1440 else if (in_fs_map) { 1441 std::vector<std::string> tokens = Ioss::tokenize(side_set_name, "_"); 1442 assert(tokens.size() >= 4); 1443 // The sideset should have only a single topology which is 1444 // given by the sideset name... 1445 const Ioss::ElementTopology *side_topo = 1446 Ioss::ElementTopology::factory(tokens[tokens.size() - 2]); 1447 assert(side_topo != nullptr); 1448 const Ioss::ElementTopology *element_topo = 1449 Ioss::ElementTopology::factory(tokens[tokens.size() - 3], true); 1450 std::string name; 1451 if (element_topo != nullptr) { 1452 name = element_topo->name(); 1453 } 1454 else { 1455 // -4 -3 -2 -1 1456 // Name is of the form name_block_id_sidetopo_id 1457 name = tokens[tokens.size() - 4] + "_" + tokens[tokens.size() - 3]; 1458 } 1459 1460 topo_map[std::make_pair(name, side_topo)] = number_sides; 1461 1462 // We want the id to match the id on the sideset in this 1463 // case so that the generated name will match the current 1464 // name. Instead of converting from string to int back to 1465 // string, we just set a variable to query later. 1466 sid = tokens[tokens.size() - 1]; 1467 } 1468 else if (split_type == Ioss::SPLIT_BY_TOPOLOGIES) { 1469 // There are multiple side types in the model. 1470 // Iterate through the elements in the sideset, determine 1471 // their parent element block using the blocks element 1472 // topology and the side number, determine the side 1473 // type. 1474 1475 for (auto side_topo : sideTopology) { 1476 topo_map[std::make_pair(side_topo.first->name(), side_topo.second)] = 0; 1477 side_map[std::make_pair(side_topo.first->name(), side_topo.second)] = 0; 1478 } 1479 1480 Ioex::separate_surface_element_sides(element, sides, get_region(), topo_map, side_map, 1481 split_type, side_set_name); 1482 } 1483 else if (split_type == Ioss::SPLIT_BY_ELEMENT_BLOCK) { 1484 // There are multiple side types in the model. Iterate 1485 // through the elements in the sideset, determine their parent 1486 // element block using blocks element topology and the side 1487 // number, determine the side type. 1488 1489 // Seed the topo_map map with <block->name, side_topo> 1490 // pairs so we are sure that all processors have the same 1491 // starting topo_map (size and order). 1492 const Ioss::ElementBlockContainer &element_blocks = get_region()->get_element_blocks(); 1493 assert(Ioss::Utils::check_block_order(element_blocks)); 1494 1495 for (Ioss::ElementBlock *block : element_blocks) { 1496 if (!Ioss::Utils::block_is_omitted(block)) { 1497 const std::string & name = block->name(); 1498 const Ioss::ElementTopology *common_ftopo = block->topology()->boundary_type(0); 1499 if (common_ftopo != nullptr) { 1500 // All sides of this element block's topology have the same topology 1501 topo_map[std::make_pair(name, common_ftopo)] = 0; 1502 side_map[std::make_pair(name, common_ftopo)] = 0; 1503 } 1504 else { 1505 // The sides have different topology, iterate over 1506 // them and create an entry for the unique side 1507 // topology types 1508 int par_dim = block->topology()->parametric_dimension(); 1509 if (par_dim == 2 || par_dim == 3) { 1510 int64_t my_side_count = block->topology()->number_boundaries(); 1511 for (int64_t ii = 0; ii < my_side_count; ii++) { 1512 const Ioss::ElementTopology *topo = block->topology()->boundary_type(ii + 1); 1513 topo_map[std::make_pair(name, topo)] = 0; 1514 side_map[std::make_pair(name, topo)] = 0; 1515 } 1516 } 1517 } 1518 } 1519 } 1520 Ioex::separate_surface_element_sides(element, sides, get_region(), topo_map, side_map, 1521 split_type, side_set_name); 1522 } 1523 } 1524 1525 // End of first step in splitting. Check among all processors 1526 // to see which potential splits have sides in them... 1527 Ioss::Int64Vector global_side_counts(topo_map.size()); 1528 { 1529 int64_t i = 0; 1530 for (auto &topo : topo_map) { 1531 global_side_counts[i++] = topo.second; 1532 } 1533 1534 // If splitting by element block, also sync the side_map 1535 // information which specifies whether the sideset has 1536 // consistent sides for all elements. Only really used for 1537 // shells, but easier to just set the value on all surfaces 1538 // in the element block split case. 1539 if (side_map.size() == topo_map.size()) { 1540 global_side_counts.resize(topo_map.size() + side_map.size()); 1541 1542 for (auto &side : side_map) { 1543 global_side_counts[i++] = side.second; 1544 } 1545 } 1546 1547 // See if any processor has non-zero count for the topo_map counts 1548 // For the side_map, need the max value. 1549 util().global_array_minmax(global_side_counts, Ioss::ParallelUtils::DO_MAX); 1550 } 1551 1552 // Create Side Blocks 1553 1554 int64_t i = 0; 1555 for (auto &topo : topo_map) { 1556 if (global_side_counts[i++] > 0) { 1557 const std::string topo_or_block_name = topo.first.first; 1558 const Ioss::ElementTopology *side_topo = topo.first.second; 1559 assert(side_topo != nullptr); 1560 #if 0 1561 if (side_topo->parametric_dimension() == topology_dimension-1 || 1562 split_type == Ioss::SPLIT_BY_DONT_SPLIT ) { 1563 #else 1564 if (true) { 1565 #endif 1566 int64_t my_side_count = topo.second; 1567 1568 std::string side_block_name = "surface_" + topo_or_block_name + "_" + side_topo->name(); 1569 if (side_set_name == "universal_sideset") { 1570 side_block_name = side_set_name; 1571 } 1572 else { 1573 if (sid == "") { 1574 side_block_name = Ioss::Utils::encode_entity_name(side_block_name, id); 1575 } 1576 else { 1577 side_block_name += "_"; 1578 side_block_name += sid; 1579 } 1580 } 1581 1582 Ioss::ElementBlock *block = nullptr; 1583 // Need to get elem_topo.... 1584 const Ioss::ElementTopology *elem_topo = nullptr; 1585 if (split_type == Ioss::SPLIT_BY_TOPOLOGIES) { 1586 elem_topo = Ioss::ElementTopology::factory(topo_or_block_name); 1587 } 1588 else if (split_type == Ioss::SPLIT_BY_ELEMENT_BLOCK) { 1589 block = get_region()->get_element_block(topo_or_block_name); 1590 if (block == nullptr || Ioss::Utils::block_is_omitted(block)) { 1591 std::ostringstream errmsg; 1592 fmt::print(errmsg, 1593 "INTERNAL ERROR: Could not find element block '{}'. Something is wrong " 1594 "in the Ioex::ParallelDatabaseIO class. Please report.\n", 1595 topo_or_block_name); 1596 IOSS_ERROR(errmsg); 1597 } 1598 elem_topo = block->topology(); 1599 } 1600 if (split_type == Ioss::SPLIT_BY_DONT_SPLIT) { 1601 // Most likely this is "unknown", but can be a true 1602 // topology if there is only a single element block in 1603 // the model. 1604 elem_topo = Ioss::ElementTopology::factory(topo_or_block_name); 1605 } 1606 assert(elem_topo != nullptr); 1607 1608 auto *side_block = new Ioss::SideBlock(this, side_block_name, side_topo->name(), 1609 elem_topo->name(), my_side_count); 1610 assert(side_block != nullptr); 1611 side_block->property_add(Ioss::Property("id", id)); 1612 side_block->property_add(Ioss::Property("guid", util().generate_guid(id))); 1613 side_set->add(side_block); 1614 1615 // Note that all sideblocks within a specific 1616 // sideset might have the same id. 1617 1618 // If splitting by element block, need to set the 1619 // element block member on this side block. 1620 if (split_type == Ioss::SPLIT_BY_ELEMENT_BLOCK) { 1621 side_block->set_parent_element_block(block); 1622 } 1623 1624 // If we calculated whether the element side is 1625 // consistent for all sides in this block, then 1626 // tell the block which side it is, or that they are 1627 // inconsistent. If it wasn't calculated above, then it 1628 // will be calculated on the fly when/if requested. 1629 // This is to avoid reading the sideset bulk data in 1630 // cases where we don't need to read it, but if we are 1631 // already reading it (to split the sidesets), then use 1632 // the data when we have it. 1633 if (!side_map.empty()) { 1634 // Set a property indicating which element side 1635 // (1-based) all sides in this block are applied to. 1636 // If they are not all assigned to the same element 1637 // side, indicate this with a side equal to 0. 1638 // 1639 // (note: 'i' has already been incremented earlier in 1640 // the loop. We need previous value here...) 1641 int side = global_side_counts[i - 1 + topo_map.size()]; 1642 if (side == 999) { 1643 side = 0; 1644 } 1645 assert(side <= elem_topo->number_boundaries()); 1646 side_block->set_consistent_side_number(side); 1647 } 1648 1649 // Add an alias... 1650 get_region()->add_alias(side_block); 1651 1652 if (split_type != Ioss::SPLIT_BY_DONT_SPLIT && side_set_name != "universal_sideset") { 1653 std::string storage = "Real["; 1654 storage += std::to_string(side_topo->number_nodes()); 1655 storage += "]"; 1656 side_block->field_add(Ioss::Field("distribution_factors", Ioss::Field::REAL, storage, 1657 Ioss::Field::MESH)); 1658 } 1659 1660 if (side_set_name == "universal_sideset") { 1661 side_block->field_add(Ioss::Field("side_ids", side_block->field_int_type(), "scalar", 1662 Ioss::Field::MESH)); 1663 } 1664 1665 int num_attr = 0; 1666 { 1667 int ierr = ex_get_attr_param(get_file_pointer(), EX_SIDE_SET, 1, &num_attr); 1668 if (ierr < 0) { 1669 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1670 } 1671 } 1672 // Add additional fields 1673 add_attribute_fields(EX_SIDE_SET, side_block, num_attr, ""); 1674 add_results_fields(EX_SIDE_SET, side_block, iss); 1675 } 1676 } 1677 } 1678 } 1679 } 1680 } // namespace Ioex 1681 1682 template <typename T> 1683 void ParallelDatabaseIO::get_sets(ex_entity_type type, int64_t count, const std::string &base, 1684 const T * /* set_type */) 1685 { 1686 // Attributes of a Xset are: 1687 // -- id 1688 // -- name 1689 // -- number of nodes 1690 // -- number of distribution factors (see next comment) 1691 // ----the #distribution factors should equal #Xs or 0, any 1692 // other value does not make sense. If it is 0, then a substitute 1693 // list will be created returning 1.0 for the factor 1694 1695 // In a parallel execution, it is possible that a Xset will have 1696 // no Xs or distribution factors on a particular processor... 1697 1698 // Get exodusII Xset metadata 1699 for (int64_t ins = 0; ins < count; ins++) { 1700 int64_t id = decomp->node_sets[ins].id(); 1701 1702 int num_attr = 0; 1703 int ierr = ex_get_attr_param(get_file_pointer(), type, id, &num_attr); 1704 if (ierr < 0) { 1705 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 1706 } 1707 1708 bool db_has_name = false; 1709 std::string Xset_name; 1710 std::string alias = Ioss::Utils::encode_entity_name(base + "list", id); 1711 if (ignore_database_names()) { 1712 Xset_name = alias; 1713 } 1714 else { 1715 Xset_name = Ioex::get_entity_name(get_file_pointer(), type, id, base + "list", 1716 maximumNameLength, db_has_name); 1717 } 1718 1719 if (get_use_generic_canonical_name()) { 1720 std::swap(Xset_name, alias); 1721 } 1722 1723 T *Xset = new T(this, Xset_name, decomp->node_sets[ins].ioss_count()); 1724 Xset->property_add(Ioss::Property("id", id)); 1725 Xset->property_add(Ioss::Property("guid", util().generate_guid(id))); 1726 if (db_has_name) { 1727 std::string *db_name = &Xset_name; 1728 if (get_use_generic_canonical_name()) { 1729 db_name = &alias; 1730 } 1731 Xset->property_add(Ioss::Property("db_name", *db_name)); 1732 } 1733 get_region()->add(Xset); 1734 get_region()->add_alias(Xset_name, alias); 1735 get_region()->add_alias(Xset_name, Ioss::Utils::encode_entity_name(base + "set", id)); 1736 add_attribute_fields(type, Xset, num_attr, ""); 1737 add_results_fields(type, Xset, ins); 1738 } 1739 } 1740 1741 void ParallelDatabaseIO::get_nodesets() 1742 { 1743 get_sets(EX_NODE_SET, m_groupCount[EX_NODE_SET], "node", (Ioss::NodeSet *)nullptr); 1744 } 1745 1746 void ParallelDatabaseIO::get_edgesets() 1747 { 1748 get_sets(EX_EDGE_SET, m_groupCount[EX_EDGE_SET], "edge", (Ioss::EdgeSet *)nullptr); 1749 } 1750 1751 void ParallelDatabaseIO::get_facesets() 1752 { 1753 get_sets(EX_FACE_SET, m_groupCount[EX_FACE_SET], "face", (Ioss::FaceSet *)nullptr); 1754 } 1755 1756 void ParallelDatabaseIO::get_elemsets() 1757 { 1758 get_sets(EX_ELEM_SET, m_groupCount[EX_ELEM_SET], "element", (Ioss::ElementSet *)nullptr); 1759 } 1760 1761 void ParallelDatabaseIO::get_commsets() 1762 { 1763 // Attributes of a commset are: 1764 // -- id (property) 1765 // -- name (property) 1766 // -- number of node--CPU pairs (field) 1767 1768 // In a parallel execution, it is possible that a commset will have 1769 // no nodes on a particular processor... 1770 1771 // If this is a serial execution, there will be no communication 1772 // nodesets, just return an empty container. 1773 1774 if (isParallel) { 1775 1776 // Create a single node commset 1777 Ioss::CommSet *commset = 1778 new Ioss::CommSet(this, "commset_node", "node", decomp->get_commset_node_size()); 1779 commset->property_add(Ioss::Property("id", 1)); 1780 commset->property_add(Ioss::Property("guid", util().generate_guid(1))); 1781 get_region()->add(commset); 1782 } 1783 } 1784 1785 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::Region *reg, const Ioss::Field &field, 1786 void *data, size_t data_size) const 1787 { 1788 return Ioex::BaseDatabaseIO::get_field_internal(reg, field, data, data_size); 1789 } 1790 1791 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::NodeBlock *nb, const Ioss::Field &field, 1792 void *data, size_t data_size) const 1793 { 1794 size_t num_to_get = field.verify(data_size); 1795 1796 #ifndef NDEBUG 1797 int64_t my_node_count = field.raw_count(); 1798 assert(my_node_count == nodeCount); 1799 #endif 1800 1801 Ioss::Field::RoleType role = field.get_role(); 1802 if (role == Ioss::Field::MESH) { 1803 if (field.get_name() == "mesh_model_coordinates_x" || 1804 field.get_name() == "mesh_model_coordinates_y" || 1805 field.get_name() == "mesh_model_coordinates_z" || 1806 field.get_name() == "mesh_model_coordinates") { 1807 decomp->get_node_coordinates(get_file_pointer(), reinterpret_cast<double *>(data), field); 1808 } 1809 1810 else if (field.get_name() == "ids") { 1811 // Map the local ids in this node block 1812 // (1...node_count) to global node ids. 1813 get_map(EX_NODE_BLOCK).map_implicit_data(data, field, num_to_get, 0); 1814 } 1815 1816 // The 1..global_node_count id. In a parallel-decomposed run, 1817 // it maps the node back to its implicit position in the serial 1818 // undecomposed mesh file. This is ONLY provided for backward- 1819 // compatibility and should not be used unless absolutely required. 1820 else if (field.get_name() == "implicit_ids") { 1821 size_t offset = decomp->decomp_node_offset(); 1822 size_t count = decomp->decomp_node_count(); 1823 if (int_byte_size_api() == 4) { 1824 std::vector<int> file_ids(count); 1825 std::iota(file_ids.begin(), file_ids.end(), offset + 1); 1826 decomp->communicate_node_data(file_ids.data(), reinterpret_cast<int *>(data), 1); 1827 } 1828 else { 1829 std::vector<int64_t> file_ids(count); 1830 std::iota(file_ids.begin(), file_ids.end(), offset + 1); 1831 decomp->communicate_node_data(file_ids.data(), reinterpret_cast<int64_t *>(data), 1); 1832 } 1833 } 1834 1835 else if (field.get_name() == "connectivity") { 1836 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1837 } 1838 else if (field.get_name() == "connectivity_raw") { 1839 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1840 } 1841 else if (field.get_name() == "node_connectivity_status") { 1842 compute_node_status(); 1843 char *status = static_cast<char *>(data); 1844 std::copy(nodeConnectivityStatus.begin(), nodeConnectivityStatus.end(), status); 1845 } 1846 1847 else if (field.get_name() == "owning_processor") { 1848 // If parallel, then set the "locally_owned" property on the nodeblocks. 1849 Ioss::CommSet *css = get_region()->get_commset("commset_node"); 1850 int *idata = static_cast<int *>(data); // Owning processor field is 4-byte int always 1851 for (int64_t i = 0; i < nodeCount; i++) { 1852 idata[i] = myProcessor; 1853 } 1854 1855 if (int_byte_size_api() == 8) { 1856 // Cannot call: 1857 // `css->get_field_data("entity_processor_raw", ent_proc);` 1858 // directly since it will cause a deadlock (in threaded code), 1859 // expand out into corresponding `get_field_internal` call. 1860 Ioss::Field ep_field = css->get_field("entity_processor_raw"); 1861 std::vector<int64_t> ent_proc(ep_field.raw_count() * 1862 ep_field.raw_storage()->component_count()); 1863 size_t ep_data_size = ent_proc.size() * sizeof(int64_t); 1864 get_field_internal(css, ep_field, ent_proc.data(), ep_data_size); 1865 for (size_t i = 0; i < ent_proc.size(); i += 2) { 1866 int64_t node = ent_proc[i + 0]; 1867 int64_t proc = ent_proc[i + 1]; 1868 if (proc < idata[node - 1]) { 1869 idata[node - 1] = proc; 1870 } 1871 } 1872 } 1873 else { 1874 Ioss::Field ep_field = css->get_field("entity_processor_raw"); 1875 std::vector<int> ent_proc(ep_field.raw_count() * ep_field.raw_storage()->component_count()); 1876 size_t ep_data_size = ent_proc.size() * sizeof(int); 1877 get_field_internal(css, ep_field, ent_proc.data(), ep_data_size); 1878 for (size_t i = 0; i < ent_proc.size(); i += 2) { 1879 int node = ent_proc[i + 0]; 1880 int proc = ent_proc[i + 1]; 1881 if (proc < idata[node - 1]) { 1882 idata[node - 1] = proc; 1883 } 1884 } 1885 } 1886 } 1887 1888 else { 1889 num_to_get = Ioss::Utils::field_warning(nb, field, "input"); 1890 } 1891 } 1892 else if (role == Ioss::Field::TRANSIENT) { 1893 // Check if the specified field exists on this node block. 1894 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 1895 // exist on the database as scalars with the appropriate 1896 // extensions. 1897 1898 // Read in each component of the variable and transfer into 1899 // 'data'. Need temporary storage area of size 'number of 1900 // nodes in this block. 1901 num_to_get = read_transient_field(EX_NODE_BLOCK, m_variables[EX_NODE_BLOCK], field, nb, data); 1902 } 1903 else if (role == Ioss::Field::ATTRIBUTE) { 1904 num_to_get = read_attribute_field(EX_NODE_BLOCK, field, nb, data); 1905 } 1906 return num_to_get; 1907 } 1908 1909 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::Blob *blob, const Ioss::Field &field, 1910 void *data, size_t data_size) const 1911 { 1912 { 1913 Ioss::SerializeIO serializeIO__(this); 1914 1915 size_t num_to_get = field.verify(data_size); 1916 if (num_to_get > 0) { 1917 1918 Ioss::Field::RoleType role = field.get_role(); 1919 if (role == Ioss::Field::MESH) { 1920 if (field.get_name() == "ids") { 1921 // Map the local ids in this node block 1922 // (1...node_count) to global node ids. 1923 // get_map(EX_BLOB).map_implicit_data(data, field, num_to_get, 0); 1924 } 1925 1926 else if (field.get_name() == "connectivity") { 1927 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1928 } 1929 else if (field.get_name() == "connectivity_raw") { 1930 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1931 } 1932 else { 1933 num_to_get = Ioss::Utils::field_warning(blob, field, "input"); 1934 } 1935 } 1936 else if (role == Ioss::Field::TRANSIENT) { 1937 // Check if the specified field exists on this blob. 1938 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 1939 // exist on the database as scalars with the appropriate 1940 // extensions. 1941 1942 // Read in each component of the variable and transfer into 1943 // 'data'. Need temporary storage area of size 'number of 1944 // items in this blob. 1945 num_to_get = read_transient_field(EX_BLOB, m_variables[EX_BLOB], field, blob, data); 1946 } 1947 else if (role == Ioss::Field::ATTRIBUTE) { 1948 num_to_get = read_attribute_field(EX_BLOB, field, blob, data); 1949 } 1950 else if (role == Ioss::Field::REDUCTION) { 1951 get_reduction_field(EX_BLOB, field, blob, data); 1952 } 1953 } 1954 return num_to_get; 1955 } 1956 } 1957 1958 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::Assembly *assembly, 1959 const Ioss::Field &field, void *data, 1960 size_t data_size) const 1961 { 1962 { 1963 Ioss::SerializeIO serializeIO__(this); 1964 1965 size_t num_to_get = field.verify(data_size); 1966 if (num_to_get > 0) { 1967 1968 Ioss::Field::RoleType role = field.get_role(); 1969 if (role == Ioss::Field::MESH) { 1970 if (field.get_name() == "ids") { 1971 // Map the local ids in this node block 1972 // (1...node_count) to global node ids. 1973 // get_map(EX_ASSEMBLY).map_implicit_data(data, field, num_to_get, 0); 1974 } 1975 1976 else if (field.get_name() == "connectivity") { 1977 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1978 } 1979 else if (field.get_name() == "connectivity_raw") { 1980 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1981 } 1982 else { 1983 num_to_get = Ioss::Utils::field_warning(assembly, field, "input"); 1984 } 1985 } 1986 else if (role == Ioss::Field::TRANSIENT) { 1987 // Check if the specified field exists on this assembly. 1988 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 1989 // exist on the database as scalars with the appropriate 1990 // extensions. 1991 1992 // Read in each component of the variable and transfer into 1993 // 'data'. Need temporary storage area of size 'number of 1994 // items in this assembly. 1995 num_to_get = 1996 read_transient_field(EX_ASSEMBLY, m_variables[EX_ASSEMBLY], field, assembly, data); 1997 } 1998 else if (role == Ioss::Field::ATTRIBUTE) { 1999 num_to_get = read_attribute_field(EX_ASSEMBLY, field, assembly, data); 2000 } 2001 } 2002 return num_to_get; 2003 } 2004 } 2005 2006 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::ElementBlock *eb, 2007 const Ioss::Field &field, void *data, 2008 size_t data_size) const 2009 { 2010 2011 size_t num_to_get = field.verify(data_size); 2012 2013 int64_t id = Ioex::get_id(eb, EX_ELEM_BLOCK, &ids_); 2014 size_t my_element_count = eb->entity_count(); 2015 Ioss::Field::RoleType role = field.get_role(); 2016 2017 if (role == Ioss::Field::MESH) { 2018 // Handle the MESH fields required for an ExodusII file model. 2019 // (The 'genesis' portion) 2020 2021 if (field.get_name() == "connectivity" || field.get_name() == "connectivity_raw") { 2022 int element_nodes = eb->topology()->number_nodes(); 2023 assert(field.raw_storage()->component_count() == element_nodes); 2024 2025 int order = eb->get_property("iblk").get_int(); 2026 // The connectivity is stored in a 1D array. 2027 // The element_node index varies fastet 2028 2029 decomp->get_block_connectivity(get_file_pointer(), data, id, order, element_nodes); 2030 if (field.get_name() == "connectivity") { 2031 get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get * element_nodes); 2032 } 2033 } 2034 #if 0 2035 else if (field.get_name() == "connectivity_face") { 2036 int face_count = field.raw_storage()->component_count(); 2037 2038 // The connectivity is stored in a 1D array. 2039 // The element_face index varies fastest 2040 if (my_element_count > 0) { 2041 get_connectivity_data(get_file_pointer(), data, EX_ELEM_BLOCK, id, 2, int_byte_size_api()); 2042 get_map(EX_FACE_BLOCK).map_data(data, field, num_to_get*face_count); 2043 } 2044 } 2045 else if (field.get_name() == "connectivity_edge") { 2046 int edge_count = field.raw_storage()->component_count(); 2047 2048 // The connectivity is stored in a 1D array. 2049 // The element_edge index varies fastest 2050 if (my_element_count > 0) { 2051 get_connectivity_data(get_file_pointer(), data, EX_ELEM_BLOCK, id, 1, int_byte_size_api()); 2052 get_map(EX_EDGE_BLOCK).map_data(data, field, num_to_get*edge_count); 2053 } 2054 } 2055 #endif 2056 else if (field.get_name() == "ids" || field.get_name() == "implicit_ids") { 2057 // Map the local ids in this element block 2058 // (eb_offset+1...eb_offset+1+my_element_count) to global element ids. 2059 get_map(EX_ELEM_BLOCK).map_implicit_data(data, field, num_to_get, eb->get_offset()); 2060 } 2061 else if (field.get_name() == "skin") { 2062 // This is (currently) for the skinned body. It maps the 2063 // side element on the skin to the original element/local 2064 // side number. It is a two component field, the first 2065 // component is the global id of the underlying element in 2066 // the initial mesh and its local side number (1-based). 2067 2068 if (field.is_type(Ioss::Field::INTEGER)) { 2069 Ioss::IntVector element(my_element_count); 2070 Ioss::IntVector side(my_element_count); 2071 int * el_side = reinterpret_cast<int *>(data); 2072 2073 // FIX: Hardwired map ids.... 2074 size_t eb_offset = eb->get_offset(); 2075 assert(1 == 0 && "Unimplemented FIXME"); 2076 ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, 1, eb_offset + 1, my_element_count, 2077 element.data()); // FIXME 2078 ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, 2, eb_offset + 1, my_element_count, 2079 side.data()); // FIXME 2080 2081 int index = 0; 2082 for (size_t i = 0; i < my_element_count; i++) { 2083 el_side[index++] = element[i]; 2084 el_side[index++] = side[i]; 2085 } 2086 } 2087 else { 2088 Ioss::Int64Vector element(my_element_count); 2089 Ioss::Int64Vector side(my_element_count); 2090 int64_t * el_side = reinterpret_cast<int64_t *>(data); 2091 2092 // FIX: Hardwired map ids.... 2093 size_t eb_offset = eb->get_offset(); 2094 assert(1 == 0 && "Unimplemented FIXME"); 2095 ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, 1, eb_offset + 1, my_element_count, 2096 element.data()); // FIXME 2097 ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, 2, eb_offset + 1, my_element_count, 2098 side.data()); // FIXME 2099 2100 size_t index = 0; 2101 for (size_t i = 0; i < my_element_count; i++) { 2102 el_side[index++] = element[i]; 2103 el_side[index++] = side[i]; 2104 } 2105 } 2106 } 2107 else { 2108 num_to_get = Ioss::Utils::field_warning(eb, field, "input"); 2109 } 2110 } 2111 else if (role == Ioss::Field::ATTRIBUTE) { 2112 num_to_get = read_attribute_field(EX_ELEM_BLOCK, field, eb, data); 2113 } 2114 else if (role == Ioss::Field::TRANSIENT) { 2115 // Check if the specified field exists on this element block. 2116 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 2117 // exist on the database as scalars with the appropriate 2118 // extensions. 2119 2120 // Read in each component of the variable and transfer into 2121 // 'data'. Need temporary storage area of size 'number of 2122 // elements in this block. 2123 num_to_get = read_transient_field(EX_ELEM_BLOCK, m_variables[EX_ELEM_BLOCK], field, eb, data); 2124 } 2125 else if (role == Ioss::Field::REDUCTION) { 2126 num_to_get = Ioss::Utils::field_warning(eb, field, "input reduction"); 2127 } 2128 return num_to_get; 2129 } 2130 2131 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::FaceBlock *eb, const Ioss::Field &field, 2132 void *data, size_t data_size) const 2133 { 2134 size_t num_to_get = field.verify(data_size); 2135 2136 int64_t id = Ioex::get_id(eb, EX_FACE_BLOCK, &ids_); 2137 size_t my_face_count = eb->entity_count(); 2138 Ioss::Field::RoleType role = field.get_role(); 2139 2140 if (role == Ioss::Field::MESH) { 2141 // Handle the MESH fields required for an ExodusII file model. 2142 // (The 'genesis' portion) 2143 2144 if (field.get_name() == "connectivity") { 2145 int face_nodes = eb->topology()->number_nodes(); 2146 assert(field.raw_storage()->component_count() == face_nodes); 2147 2148 // The connectivity is stored in a 1D array. 2149 // The face_node index varies fastet 2150 if (my_face_count > 0) { 2151 get_connectivity_data(get_file_pointer(), data, EX_FACE_BLOCK, id, 0, int_byte_size_api()); 2152 get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get * face_nodes); 2153 } 2154 } 2155 else if (field.get_name() == "connectivity_edge") { 2156 int edge_count = field.raw_storage()->component_count(); 2157 2158 // The connectivity is stored in a 1D array. 2159 // The face_edge index varies fastest 2160 if (my_face_count > 0) { 2161 get_connectivity_data(get_file_pointer(), data, EX_FACE_BLOCK, id, 1, int_byte_size_api()); 2162 get_map(EX_EDGE_BLOCK).map_data(data, field, num_to_get * edge_count); 2163 } 2164 } 2165 else if (field.get_name() == "connectivity_raw") { 2166 // "connectivity_raw" has nodes in local id space (1-based) 2167 assert(field.raw_storage()->component_count() == eb->topology()->number_nodes()); 2168 2169 // The connectivity is stored in a 1D array. 2170 // The face_node index varies fastet 2171 if (my_face_count > 0) { 2172 get_connectivity_data(get_file_pointer(), data, EX_FACE_BLOCK, id, 0, int_byte_size_api()); 2173 } 2174 } 2175 else if (field.get_name() == "ids") { 2176 // Map the local ids in this face block 2177 // (eb_offset+1...eb_offset+1+my_face_count) to global face ids. 2178 get_map(EX_FACE_BLOCK).map_implicit_data(data, field, num_to_get, eb->get_offset()); 2179 } 2180 else { 2181 num_to_get = Ioss::Utils::field_warning(eb, field, "input"); 2182 } 2183 } 2184 else if (role == Ioss::Field::ATTRIBUTE) { 2185 num_to_get = read_attribute_field(EX_FACE_BLOCK, field, eb, data); 2186 } 2187 else if (role == Ioss::Field::TRANSIENT) { 2188 // Check if the specified field exists on this element block. 2189 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 2190 // exist on the database as scalars with the appropriate 2191 // extensions. 2192 2193 // Read in each component of the variable and transfer into 2194 // 'data'. Need temporary storage area of size 'number of 2195 // elements in this block. 2196 num_to_get = read_transient_field(EX_FACE_BLOCK, m_variables[EX_FACE_BLOCK], field, eb, data); 2197 } 2198 else if (role == Ioss::Field::REDUCTION) { 2199 num_to_get = Ioss::Utils::field_warning(eb, field, "input reduction"); 2200 } 2201 return num_to_get; 2202 } 2203 2204 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::EdgeBlock *eb, const Ioss::Field &field, 2205 void *data, size_t data_size) const 2206 { 2207 size_t num_to_get = field.verify(data_size); 2208 2209 int64_t id = Ioex::get_id(eb, EX_EDGE_BLOCK, &ids_); 2210 int64_t my_edge_count = eb->entity_count(); 2211 Ioss::Field::RoleType role = field.get_role(); 2212 2213 if (role == Ioss::Field::MESH) { 2214 // Handle the MESH fields required for an ExodusII file model. 2215 // (The 'genesis' portion) 2216 2217 if (field.get_name() == "connectivity") { 2218 int edge_nodes = eb->topology()->number_nodes(); 2219 assert(field.raw_storage()->component_count() == edge_nodes); 2220 2221 // The connectivity is stored in a 1D array. 2222 // The edge_node index varies fastet 2223 if (my_edge_count > 0) { 2224 get_connectivity_data(get_file_pointer(), data, EX_EDGE_BLOCK, id, 0, int_byte_size_api()); 2225 get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get * edge_nodes); 2226 } 2227 } 2228 else if (field.get_name() == "connectivity_raw") { 2229 // "connectivity_raw" has nodes in local id space (1-based) 2230 assert(field.raw_storage()->component_count() == eb->topology()->number_nodes()); 2231 2232 // The connectivity is stored in a 1D array. 2233 // The edge_node index varies fastet 2234 if (my_edge_count > 0) { 2235 get_connectivity_data(get_file_pointer(), data, EX_EDGE_BLOCK, id, 0, int_byte_size_api()); 2236 } 2237 } 2238 else if (field.get_name() == "ids") { 2239 // Map the local ids in this edge block 2240 // (eb_offset+1...eb_offset+1+my_edge_count) to global edge ids. 2241 get_map(EX_EDGE_BLOCK).map_implicit_data(data, field, num_to_get, eb->get_offset()); 2242 } 2243 else { 2244 num_to_get = Ioss::Utils::field_warning(eb, field, "input"); 2245 } 2246 } 2247 else if (role == Ioss::Field::ATTRIBUTE) { 2248 num_to_get = read_attribute_field(EX_EDGE_BLOCK, field, eb, data); 2249 } 2250 else if (role == Ioss::Field::TRANSIENT) { 2251 // Check if the specified field exists on this element block. 2252 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 2253 // exist on the database as scalars with the appropriate 2254 // extensions. 2255 2256 // Read in each component of the variable and transfer into 2257 // 'data'. Need temporary storage area of size 'number of 2258 // elements in this block. 2259 num_to_get = read_transient_field(EX_EDGE_BLOCK, m_variables[EX_EDGE_BLOCK], field, eb, data); 2260 } 2261 else if (role == Ioss::Field::REDUCTION) { 2262 num_to_get = Ioss::Utils::field_warning(eb, field, "input reduction"); 2263 } 2264 return num_to_get; 2265 } 2266 2267 int64_t ParallelDatabaseIO::get_Xset_field_internal(ex_entity_type type, const Ioss::EntitySet *ns, 2268 const Ioss::Field &field, void *data, 2269 size_t data_size) const 2270 { 2271 int ierr; 2272 size_t num_to_get = field.verify(data_size); 2273 Ioss::Field::RoleType role = field.get_role(); 2274 2275 // Find corresponding set in file decomp class... 2276 if (role == Ioss::Field::MESH) { 2277 int64_t id = Ioex::get_id(ns, type, &ids_); 2278 2279 if (field.get_name() == "ids" || field.get_name() == "ids_raw") { 2280 if (field.get_type() == Ioss::Field::INTEGER) { 2281 ierr = decomp->get_set_mesh_var(get_file_pointer(), EX_NODE_SET, id, field, 2282 static_cast<int *>(data)); 2283 } 2284 else { 2285 ierr = decomp->get_set_mesh_var(get_file_pointer(), EX_NODE_SET, id, field, 2286 static_cast<int64_t *>(data)); 2287 } 2288 if (ierr < 0) { 2289 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2290 } 2291 2292 if (field.get_name() == "ids") { 2293 // Convert the local node ids to global ids 2294 get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get); 2295 } 2296 } 2297 else if (field.get_name() == "orientation") { 2298 if (field.get_type() == Ioss::Field::INTEGER) { 2299 ierr = decomp->get_set_mesh_var(get_file_pointer(), EX_NODE_SET, id, field, 2300 static_cast<int *>(data)); 2301 } 2302 else { 2303 ierr = decomp->get_set_mesh_var(get_file_pointer(), EX_NODE_SET, id, field, 2304 static_cast<int64_t *>(data)); 2305 } 2306 if (ierr < 0) { 2307 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2308 } 2309 } 2310 else if (field.get_name() == "distribution_factors") { 2311 ierr = decomp->get_set_mesh_double(get_file_pointer(), EX_NODE_SET, id, field, 2312 static_cast<double *>(data)); 2313 if (ierr < 0) { 2314 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2315 } 2316 } 2317 } 2318 else if (role == Ioss::Field::ATTRIBUTE) { 2319 num_to_get = read_attribute_field(type, field, ns, data); 2320 } 2321 else if (role == Ioss::Field::TRANSIENT) { 2322 // Check if the specified field exists on this node block. 2323 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 2324 // exist on the database as scalars with the appropriate 2325 // extensions. 2326 2327 // Read in each component of the variable and transfer into 2328 // 'data'. Need temporary storage area of size 'number of 2329 // nodes in this block. 2330 num_to_get = read_transient_field(type, m_variables[type], field, ns, data); 2331 } 2332 return num_to_get; 2333 } 2334 2335 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::NodeSet *ns, const Ioss::Field &field, 2336 void *data, size_t data_size) const 2337 { 2338 return get_Xset_field_internal(EX_NODE_SET, ns, field, data, data_size); 2339 } 2340 2341 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::EdgeSet *ns, const Ioss::Field &field, 2342 void *data, size_t data_size) const 2343 { 2344 return get_Xset_field_internal(EX_EDGE_SET, ns, field, data, data_size); 2345 } 2346 2347 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::FaceSet *ns, const Ioss::Field &field, 2348 void *data, size_t data_size) const 2349 { 2350 return get_Xset_field_internal(EX_FACE_SET, ns, field, data, data_size); 2351 } 2352 2353 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::ElementSet *ns, const Ioss::Field &field, 2354 void *data, size_t data_size) const 2355 { 2356 return get_Xset_field_internal(EX_ELEM_SET, ns, field, data, data_size); 2357 } 2358 2359 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::SideSet *fs, const Ioss::Field &field, 2360 void *data, size_t data_size) const 2361 { 2362 size_t num_to_get = field.verify(data_size); 2363 if (field.get_name() == "ids") { 2364 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 2365 // However, make sure that the caller gets a consistent answer, i.e., don't leave the buffer 2366 // full of junk 2367 memset(data, 0x00, data_size); 2368 } 2369 else { 2370 num_to_get = Ioss::Utils::field_warning(fs, field, "input"); 2371 } 2372 return num_to_get; 2373 } 2374 2375 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::CommSet *cs, const Ioss::Field &field, 2376 void *data, size_t data_size) const 2377 { 2378 size_t num_to_get = field.verify(data_size); 2379 2380 // Return the <entity (node or side), processor> pair 2381 if (field.get_name() == "entity_processor" || field.get_name() == "entity_processor_raw") { 2382 2383 // Check type -- node or side 2384 std::string type = cs->get_property("entity_type").get_string(); 2385 2386 if (type == "node") { 2387 2388 bool do_map = field.get_name() == "entity_processor"; 2389 // Convert local node id to global node id and store in 'data' 2390 const Ioss::MapContainer &map = get_map(EX_NODE_BLOCK).map(); 2391 if (int_byte_size_api() == 4) { 2392 decomp->get_node_entity_proc_data(static_cast<int *>(data), map, do_map); 2393 } 2394 else { 2395 decomp->get_node_entity_proc_data(static_cast<int64_t *>(data), map, do_map); 2396 } 2397 } 2398 else { 2399 std::ostringstream errmsg; 2400 fmt::print(errmsg, "ERROR: Invalid commset type {}", type); 2401 IOSS_ERROR(errmsg); 2402 } 2403 } 2404 else if (field.get_name() == "ids") { 2405 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 2406 } 2407 else { 2408 num_to_get = Ioss::Utils::field_warning(cs, field, "input"); 2409 } 2410 return num_to_get; 2411 } 2412 2413 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::SideBlock *fb, const Ioss::Field &field, 2414 void *data, size_t data_size) const 2415 { 2416 ssize_t num_to_get = field.verify(data_size); 2417 int ierr = 0; 2418 2419 int64_t id = Ioex::get_id(fb, EX_SIDE_SET, &ids_); 2420 int64_t entity_count = fb->entity_count(); 2421 if (num_to_get != entity_count) { 2422 std::ostringstream errmsg; 2423 fmt::print(errmsg, "ERROR: Partial field input not yet implemented for side blocks"); 2424 IOSS_ERROR(errmsg); 2425 } 2426 2427 auto &set = decomp->get_decomp_set(EX_SIDE_SET, id); 2428 2429 int64_t number_sides = set.ioss_count(); 2430 int64_t number_distribution_factors = set.df_count(); 2431 2432 Ioss::Field::RoleType role = field.get_role(); 2433 if (role == Ioss::Field::MESH) { 2434 2435 // In exodusII, we may have split the sideset into multiple 2436 // side blocks if there are multiple side topologies in the 2437 // sideset. Because of this, the passed in 'data' may not be 2438 // large enough to hold the data residing in the sideset and we 2439 // may need to allocate a temporary array... This can be checked 2440 // by comparing the size of the sideset with the 'side_count' of 2441 // the side block. 2442 2443 // Get size of data stored on the file... 2444 // FIX 64: FIX THIS -- STORING INT IN DOUBLE WON'T WORK 2445 if (field.get_name() == "side_ids" && fb->name() == "universal_sideset") { 2446 // The side ids are being stored as the distribution factor 2447 // field on the universal sideset. There should be no other 2448 // side sets that request this field... (Eventually, 2449 // create an id field to store this info. 2450 2451 if (number_distribution_factors == num_to_get) { 2452 std::vector<double> real_ids(num_to_get); 2453 Ioss::Field df_field("distribution_factor", Ioss::Field::REAL, "scalar", Ioss::Field::MESH, 2454 num_to_get); 2455 decomp->get_set_mesh_double(get_file_pointer(), EX_SIDE_SET, id, df_field, real_ids.data()); 2456 2457 if (field.get_type() == Ioss::Field::INTEGER) { 2458 // Need to convert 'double' to 'int' for Sierra use... 2459 int *ids = static_cast<int *>(data); 2460 for (ssize_t i = 0; i < num_to_get; i++) { 2461 ids[i] = static_cast<int>(real_ids[i]); 2462 } 2463 } 2464 else { 2465 // Need to convert 'double' to 'int' for Sierra use... 2466 int64_t *ids = static_cast<int64_t *>(data); 2467 for (ssize_t i = 0; i < num_to_get; i++) { 2468 ids[i] = static_cast<int64_t>(real_ids[i]); 2469 } 2470 } 2471 } 2472 } 2473 2474 else if (field.get_name() == "side_ids") { 2475 } 2476 2477 else if (field.get_name() == "ids") { 2478 // In exodusII, the 'side set' is stored as a sideset. A 2479 // sideset has a list of elements and a corresponding local 2480 // element side (1-based) The side id is: side_id = 2481 // 10*element_id + local_side_number This assumes that all 2482 // sides in a sideset are boundary sides. Since we 2483 // only have a single array, we need to allocate an extra array 2484 // to store all of the data. Note also that the element_id is 2485 // the global id but only the local id is stored so we need to 2486 // map from local_to_global prior to generating the side id... 2487 2488 Ioss::Field el_side = fb->get_field("element_side"); 2489 std::vector<char> element_side(2 * number_sides * int_byte_size_api()); 2490 get_field_internal(fb, el_side, element_side.data(), element_side.size()); 2491 2492 // At this point, have the 'element_side' data containing 2493 // the global element ids and the sides... Iterate 2494 // through to generate the ids... 2495 if (int_byte_size_api() == 4) { 2496 int64_t int_max = std::numeric_limits<int>::max(); 2497 int * ids = static_cast<int *>(data); 2498 int * els = reinterpret_cast<int *>(element_side.data()); 2499 size_t idx = 0; 2500 for (int64_t iel = 0; iel < 2 * entity_count; iel += 2) { 2501 int64_t new_id = static_cast<int64_t>(10) * els[iel] + els[iel + 1]; 2502 if (new_id > int_max) { 2503 std::ostringstream errmsg; 2504 fmt::print(errmsg, 2505 "ERROR: accessing the sideset field 'ids'\n" 2506 "\t\thas exceeded the integer bounds for entity {}, local side id {}.\n" 2507 "\t\tTry using 64-bit mode to read the file '{}'.\n", 2508 els[iel], els[iel + 1], get_filename()); 2509 IOSS_ERROR(errmsg); 2510 } 2511 2512 ids[idx++] = static_cast<int>(new_id); 2513 } 2514 } 2515 else { 2516 int64_t *ids = static_cast<int64_t *>(data); 2517 int64_t *els = reinterpret_cast<int64_t *>(element_side.data()); 2518 size_t idx = 0; 2519 for (int64_t iel = 0; iel < 2 * entity_count; iel += 2) { 2520 int64_t new_id = 10 * els[iel] + els[iel + 1]; 2521 ids[idx++] = new_id; 2522 } 2523 } 2524 } 2525 else if (field.get_name() == "element_side" || field.get_name() == "element_side_raw") { 2526 // In exodusII, the 'side set' is stored as a sideset. A sideset 2527 // has a list of elements and a corresponding local element side 2528 // (1-based) 2529 2530 // Since we only have a single array, we need to allocate an extra 2531 // array to store all of the data. Note also that the element_id 2532 // is the global id but only the local id is stored so we need to 2533 // map from local_to_global prior to generating the side id... 2534 2535 // Get the element number map (1-based)... 2536 const Ioss::MapContainer &map = get_map(EX_ELEM_BLOCK).map(); 2537 2538 // Allocate space for local side number and element numbers 2539 // numbers. 2540 2541 // See if edges or faces... 2542 int64_t side_offset = Ioss::Utils::get_side_offset(fb); 2543 2544 if (fb->owner()->block_count() == 1 && number_sides == entity_count) { 2545 2546 if (int_byte_size_api() == 4) { 2547 int *element_side = static_cast<int *>(data); 2548 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, field, element_side); 2549 for (int64_t iel = 1; iel < 2 * entity_count; iel += 2) { 2550 element_side[iel] = element_side[iel] - side_offset; 2551 } 2552 } 2553 else { 2554 int64_t *element_side = static_cast<int64_t *>(data); 2555 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, field, element_side); 2556 for (int64_t iel = 1; iel < 2 * entity_count; iel += 2) { 2557 element_side[iel] = element_side[iel] - side_offset; 2558 } 2559 } 2560 } 2561 else { 2562 // Need a larger vector to get the entire sideset and then filter down to the correct 2563 // size... 2564 std::vector<char> element(number_sides * int_byte_size_api()); 2565 std::vector<char> sides(number_sides * int_byte_size_api()); 2566 if (int_byte_size_api() == 4) { 2567 Ioss::Field elem_field("ids", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, 2568 number_sides); 2569 Ioss::Field side_field("sides", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, 2570 number_sides); 2571 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 2572 reinterpret_cast<int *>(element.data())); 2573 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 2574 reinterpret_cast<int *>(sides.data())); 2575 } 2576 else { 2577 Ioss::Field elem_field("ids", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, 2578 number_sides); 2579 Ioss::Field side_field("sides", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, 2580 number_sides); 2581 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 2582 reinterpret_cast<int64_t *>(element.data())); 2583 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 2584 reinterpret_cast<int64_t *>(sides.data())); 2585 } 2586 2587 Ioss::IntVector is_valid_side; 2588 Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(), 2589 element.data(), sides.data(), number_sides, 2590 get_region()); 2591 2592 ssize_t index = 0; 2593 if (int_byte_size_api() == 4) { 2594 int *element_side = static_cast<int *>(data); 2595 int *element32 = reinterpret_cast<int *>(element.data()); 2596 int *sides32 = reinterpret_cast<int *>(sides.data()); 2597 for (int64_t iel = 0; iel < number_sides; iel++) { 2598 if (is_valid_side[iel] == 1) { 2599 // This side belongs in the side block 2600 element_side[index++] = element32[iel]; 2601 element_side[index++] = sides32[iel] - side_offset; 2602 } 2603 } 2604 } 2605 else { 2606 int64_t *element_side = static_cast<int64_t *>(data); 2607 int64_t *element64 = reinterpret_cast<int64_t *>(element.data()); 2608 int64_t *sides64 = reinterpret_cast<int64_t *>(sides.data()); 2609 for (int64_t iel = 0; iel < number_sides; iel++) { 2610 if (is_valid_side[iel] == 1) { 2611 // This side belongs in the side block 2612 element_side[index++] = element64[iel]; 2613 element_side[index++] = sides64[iel] - side_offset; 2614 } 2615 } 2616 } 2617 assert(index / 2 == entity_count); 2618 } 2619 if (field.get_name() == "element_side") { 2620 if (int_byte_size_api() == 4) { 2621 int *element_side = static_cast<int *>(data); 2622 for (int64_t iel = 0; iel < 2 * entity_count; iel += 2) { 2623 element_side[iel] = map[element_side[iel]]; 2624 } 2625 } 2626 else { 2627 int64_t *element_side = static_cast<int64_t *>(data); 2628 for (int64_t iel = 0; iel < 2 * entity_count; iel += 2) { 2629 element_side[iel] = map[element_side[iel]]; 2630 } 2631 } 2632 } 2633 } 2634 else if (field.get_name() == "connectivity") { 2635 // The side connectivity needs to be generated 'on-the-fly' from 2636 // the element number and local side of that element. A sideset 2637 // can span multiple element blocks, and contain multiple side 2638 // types; the side block contains side of similar topology. 2639 ierr = get_side_connectivity(fb, id, entity_count, data, true); 2640 if (ierr < 0) { 2641 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2642 } 2643 } 2644 else if (field.get_name() == "connectivity_raw") { 2645 // The side connectivity needs to be generated 'on-the-fly' from 2646 // the element number and local side of that element. A sideset 2647 // can span multiple element blocks, and contain multiple side 2648 // types; the side block contains side of similar topology. 2649 ierr = get_side_connectivity(fb, id, entity_count, data, false); 2650 if (ierr < 0) { 2651 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2652 } 2653 } 2654 else if (field.get_name() == "distribution_factors") { 2655 ierr = get_side_distributions(fb, id, entity_count, static_cast<double *>(data), 2656 data_size / sizeof(double)); 2657 if (ierr < 0) { 2658 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2659 } 2660 } 2661 else { 2662 num_to_get = Ioss::Utils::field_warning(fb, field, "input"); 2663 } 2664 } 2665 else if (role == Ioss::Field::TRANSIENT) { 2666 // Check if the specified field exists on this block. 2667 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 2668 // exist on the database as scalars with the appropriate 2669 // extensions. 2670 2671 if (fb->owner()->block_count() == 1 && number_sides == entity_count) { 2672 num_to_get = read_transient_field(EX_SIDE_SET, m_variables[EX_SIDE_SET], field, fb, data); 2673 } 2674 else { 2675 // Need to read all values for the specified field and then 2676 // filter down to the elements actually in this side block. 2677 2678 Ioss::IntVector is_valid_side; 2679 // Need a larger vector to get the entire sideset and then filter down to the correct size... 2680 std::vector<char> element(number_sides * int_byte_size_api()); 2681 std::vector<char> sides(number_sides * int_byte_size_api()); 2682 if (int_byte_size_api() == 4) { 2683 Ioss::Field elem_field("ids", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, 2684 number_sides); 2685 Ioss::Field side_field("sides", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, 2686 number_sides); 2687 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 2688 reinterpret_cast<int *>(element.data())); 2689 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 2690 reinterpret_cast<int *>(sides.data())); 2691 } 2692 else { 2693 Ioss::Field elem_field("ids", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, 2694 number_sides); 2695 Ioss::Field side_field("sides", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, 2696 number_sides); 2697 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 2698 reinterpret_cast<int64_t *>(element.data())); 2699 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 2700 reinterpret_cast<int64_t *>(sides.data())); 2701 } 2702 Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(), 2703 element.data(), sides.data(), number_sides, 2704 get_region()); 2705 2706 num_to_get = read_ss_transient_field(field, id, data, is_valid_side); 2707 } 2708 } 2709 return num_to_get; 2710 } 2711 2712 int64_t ParallelDatabaseIO::write_attribute_field(ex_entity_type type, const Ioss::Field &field, 2713 const Ioss::GroupingEntity *ge, void *data) const 2714 { 2715 std::string att_name = ge->name() + SEP() + field.get_name(); 2716 ssize_t num_entity = ge->entity_count(); 2717 ssize_t offset = field.get_index(); 2718 2719 int64_t id = Ioex::get_id(ge, type, &ids_); 2720 assert(offset > 0); 2721 assert(offset - 1 + field.raw_storage()->component_count() <= 2722 ge->get_property("attribute_count").get_int()); 2723 2724 size_t proc_offset = ge->get_optional_property("_processor_offset", 0); 2725 ssize_t file_count = ge->get_optional_property("locally_owned_count", num_entity); 2726 2727 Ioss::Field::BasicType ioss_type = field.get_type(); 2728 assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER || 2729 ioss_type == Ioss::Field::INT64); 2730 2731 int comp_count = field.raw_storage()->component_count(); 2732 2733 if (type == EX_NODAL) { 2734 for (int i = 0; i < comp_count; i++) { 2735 std::vector<double> file_data; 2736 file_data.reserve(file_count); 2737 check_node_owning_processor_data(nodeOwningProcessor, file_count); 2738 if (ioss_type == Ioss::Field::REAL) { 2739 filter_owned_nodes(nodeOwningProcessor, myProcessor, static_cast<double *>(data), file_data, 2740 i, comp_count); 2741 } 2742 else if (ioss_type == Ioss::Field::INTEGER) { 2743 filter_owned_nodes(nodeOwningProcessor, myProcessor, static_cast<int *>(data), file_data, i, 2744 comp_count); 2745 } 2746 else if (ioss_type == Ioss::Field::INT64) { 2747 filter_owned_nodes(nodeOwningProcessor, myProcessor, static_cast<int64_t *>(data), 2748 file_data, i, comp_count); 2749 } 2750 int ierr = ex_put_partial_one_attr(get_file_pointer(), type, id, proc_offset + 1, file_count, 2751 offset + i, file_data.data()); 2752 if (ierr < 0) { 2753 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2754 } 2755 } 2756 } 2757 else if (type == EX_NODE_SET) { 2758 for (int i = 0; i < comp_count; i++) { 2759 std::vector<double> file_data; 2760 file_data.reserve(file_count); 2761 if (ioss_type == Ioss::Field::REAL) { 2762 map_nodeset_data(nodesetOwnedNodes[ge], static_cast<double *>(data), file_data, i, 2763 comp_count); 2764 } 2765 else if (ioss_type == Ioss::Field::INTEGER) { 2766 map_nodeset_data(nodesetOwnedNodes[ge], static_cast<int *>(data), file_data, i, comp_count); 2767 } 2768 else if (ioss_type == Ioss::Field::INT64) { 2769 map_nodeset_data(nodesetOwnedNodes[ge], static_cast<int64_t *>(data), file_data, i, 2770 comp_count); 2771 } 2772 2773 int ierr = ex_put_partial_one_attr(get_file_pointer(), type, id, proc_offset + 1, file_count, 2774 offset + i, file_data.data()); 2775 if (ierr < 0) { 2776 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2777 } 2778 } 2779 } 2780 else { 2781 assert(file_count == num_entity); 2782 std::vector<double> file_data(file_count); 2783 for (int i = 0; i < comp_count; i++) { 2784 if (ioss_type == Ioss::Field::REAL) { 2785 extract_data(file_data, static_cast<double *>(data), num_entity, i, comp_count); 2786 } 2787 else if (ioss_type == Ioss::Field::INTEGER) { 2788 extract_data(file_data, static_cast<int *>(data), num_entity, i, comp_count); 2789 } 2790 else if (ioss_type == Ioss::Field::INT64) { 2791 extract_data(file_data, static_cast<int64_t *>(data), num_entity, i, comp_count); 2792 } 2793 2794 int ierr = ex_put_partial_one_attr(get_file_pointer(), type, id, proc_offset + 1, file_count, 2795 offset + i, file_data.data()); 2796 if (ierr < 0) { 2797 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2798 } 2799 } 2800 } 2801 return num_entity; 2802 } 2803 2804 int64_t ParallelDatabaseIO::read_attribute_field(ex_entity_type type, const Ioss::Field &field, 2805 const Ioss::GroupingEntity *ge, void *data) const 2806 { 2807 int64_t num_entity = ge->entity_count(); 2808 2809 int attribute_count = ge->get_property("attribute_count").get_int(); 2810 int64_t id = Ioex::get_id(ge, type, &ids_); 2811 2812 Ioss::Field::BasicType ioss_type = field.get_type(); 2813 if (ioss_type == Ioss::Field::INTEGER || ioss_type == Ioss::Field::INT64) { 2814 std::ostringstream errmsg; 2815 fmt::print(errmsg, "INTERNAL ERROR: Integer attribute fields are not yet handled for read. " 2816 "Please report.\n"); 2817 IOSS_ERROR(errmsg); 2818 } 2819 2820 std::string att_name = ge->name() + SEP() + field.get_name(); 2821 ssize_t offset = field.get_index(); 2822 assert(offset - 1 + field.raw_storage()->component_count() <= attribute_count); 2823 if (offset == 1 && field.raw_storage()->component_count() == attribute_count) { 2824 // Read all attributes in one big chunk... 2825 int ierr = decomp->get_attr(get_file_pointer(), type, id, attribute_count, 2826 static_cast<double *>(data)); 2827 if (ierr < 0) { 2828 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2829 } 2830 } 2831 else { 2832 // Read a subset of the attributes. If scalar, read one; 2833 // if higher-order (vector3d, ..) read each component and 2834 // put into correct location... 2835 if (field.raw_storage()->component_count() == 1) { 2836 int ierr = 2837 decomp->get_one_attr(get_file_pointer(), type, id, offset, static_cast<double *>(data)); 2838 if (ierr < 0) { 2839 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2840 } 2841 } 2842 else { 2843 // Multi-component... 2844 // Need a local memory space to read data into and 2845 // then push that into the user-supplied data block... 2846 std::vector<double> local_data(num_entity); 2847 int comp_count = field.raw_storage()->component_count(); 2848 double * rdata = static_cast<double *>(data); 2849 for (int i = 0; i < comp_count; i++) { 2850 int ierr = 2851 decomp->get_one_attr(get_file_pointer(), type, id, offset + i, local_data.data()); 2852 if (ierr < 0) { 2853 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2854 } 2855 2856 size_t k = i; 2857 for (int64_t j = 0; j < num_entity; j++) { 2858 rdata[k] = local_data[j]; 2859 k += comp_count; 2860 } 2861 } 2862 } 2863 } 2864 return num_entity; 2865 } 2866 2867 int64_t ParallelDatabaseIO::read_transient_field(ex_entity_type type, 2868 const Ioex::VariableNameMap &variables, 2869 const Ioss::Field & field, 2870 const Ioss::GroupingEntity *ge, void *data) const 2871 { 2872 const Ioss::VariableType *var_type = field.raw_storage(); 2873 2874 // Read into a double variable since that is all ExodusII can store... 2875 size_t num_entity = ge->entity_count(); 2876 std::vector<double> temp(num_entity); 2877 2878 size_t step = get_current_state(); 2879 2880 // get number of components, cycle through each component 2881 // and add suffix to base 'field_name'. Look up index 2882 // of this name in 'nodeVariables' map 2883 size_t comp_count = var_type->component_count(); 2884 2885 char field_suffix_separator = get_field_separator(); 2886 for (size_t i = 0; i < comp_count; i++) { 2887 std::string var_name = var_type->label_name(field.get_name(), i + 1, field_suffix_separator); 2888 2889 // Read the variable... 2890 int64_t id = Ioex::get_id(ge, type, &ids_); 2891 int ierr = 0; 2892 auto var_iter = variables.find(var_name); 2893 if (var_iter == variables.end()) { 2894 std::ostringstream errmsg; 2895 fmt::print(errmsg, "ERROR: Could not find field '{}'\n", var_name); 2896 IOSS_ERROR(errmsg); 2897 } 2898 size_t var_index = var_iter->second; 2899 assert(var_index > 0); 2900 if (type == EX_BLOB) { 2901 size_t offset = ge->get_property("_processor_offset").get_int(); 2902 ierr = ex_get_partial_var(get_file_pointer(), step, type, var_index, id, offset + 1, 2903 num_entity, temp.data()); 2904 } 2905 else { 2906 ierr = decomp->get_var(get_file_pointer(), step, type, var_index, id, num_entity, temp); 2907 } 2908 if (ierr < 0) { 2909 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2910 } 2911 2912 // Transfer to 'data' array. 2913 size_t k = 0; 2914 if (field.get_type() == Ioss::Field::INTEGER) { 2915 int *ivar = static_cast<int *>(data); 2916 for (size_t j = i; j < num_entity * comp_count; j += comp_count) { 2917 ivar[j] = static_cast<int>(temp[k++]); 2918 } 2919 } 2920 else if (field.get_type() == Ioss::Field::INT64) { // FIX 64 UNSAFE 2921 int64_t *ivar = static_cast<int64_t *>(data); 2922 for (size_t j = i; j < num_entity * comp_count; j += comp_count) { 2923 ivar[j] = static_cast<int64_t>(temp[k++]); 2924 } 2925 } 2926 else if (field.get_type() == Ioss::Field::REAL) { 2927 double *rvar = static_cast<double *>(data); 2928 for (size_t j = i; j < num_entity * comp_count; j += comp_count) { 2929 rvar[j] = temp[k++]; 2930 } 2931 } 2932 else { 2933 std::ostringstream errmsg; 2934 fmt::print(errmsg, 2935 "IOSS_ERROR: Field storage type must be either integer or double.\n" 2936 " Field '{}' is invalid.\n", 2937 field.get_name()); 2938 IOSS_ERROR(errmsg); 2939 } 2940 assert(k == num_entity); 2941 } 2942 return num_entity; 2943 } 2944 2945 int64_t ParallelDatabaseIO::read_ss_transient_field(const Ioss::Field &field, int64_t id, 2946 void * variables, 2947 Ioss::IntVector &is_valid_side) const 2948 { 2949 size_t num_valid_sides = 0; 2950 const Ioss::VariableType *var_type = field.raw_storage(); 2951 size_t my_side_count = is_valid_side.size(); 2952 std::vector<double> temp(my_side_count); 2953 2954 size_t step = get_current_state(); 2955 2956 // get number of components, cycle through each component 2957 // and add suffix to base 'field_name'. Look up index 2958 // of this name in 'nodeVariables' map 2959 size_t comp_count = var_type->component_count(); 2960 2961 char field_suffix_separator = get_field_separator(); 2962 for (size_t i = 0; i < comp_count; i++) { 2963 std::string var_name = var_type->label_name(field.get_name(), i + 1, field_suffix_separator); 2964 2965 // Read the variable... 2966 int ierr = 0; 2967 auto var_iter = m_variables[EX_SIDE_SET].find(var_name); 2968 if (var_iter == m_variables[EX_SIDE_SET].end()) { 2969 std::ostringstream errmsg; 2970 fmt::print(errmsg, "ERROR: Could not find Sideset field '{}'\n", var_name); 2971 IOSS_ERROR(errmsg); 2972 } 2973 size_t var_index = var_iter->second; 2974 assert(var_index > 0); 2975 ierr = 2976 decomp->get_var(get_file_pointer(), step, EX_SIDE_SET, var_index, id, my_side_count, temp); 2977 if (ierr < 0) { 2978 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 2979 } 2980 2981 // Transfer to 'variables' array. 2982 size_t j = i; 2983 if (field.get_type() == Ioss::Field::INTEGER) { 2984 int *ivar = static_cast<int *>(variables); 2985 for (size_t k = 0; k < my_side_count; k++) { 2986 if (is_valid_side[k] == 1) { 2987 ivar[j] = static_cast<int>(temp[k]); 2988 j += comp_count; 2989 } 2990 } 2991 } 2992 else if (field.get_type() == Ioss::Field::INT64) { // FIX 64 UNSAFE 2993 int64_t *ivar = static_cast<int64_t *>(variables); 2994 for (size_t k = 0; k < my_side_count; k++) { 2995 if (is_valid_side[k] == 1) { 2996 ivar[j] = static_cast<int64_t>(temp[k]); 2997 j += comp_count; 2998 } 2999 } 3000 } 3001 else if (field.get_type() == Ioss::Field::REAL) { 3002 double *rvar = static_cast<double *>(variables); 3003 for (size_t k = 0; k < my_side_count; k++) { 3004 if (is_valid_side[k] == 1) { 3005 rvar[j] = temp[k]; 3006 j += comp_count; 3007 } 3008 } 3009 } 3010 else { 3011 std::ostringstream errmsg; 3012 fmt::print(errmsg, 3013 "IOSS_ERROR: Field storage type must be either integer or double.\n" 3014 " Field '{}' is invalid.\n", 3015 field.get_name()); 3016 IOSS_ERROR(errmsg); 3017 } 3018 if (i + 1 == comp_count) { 3019 num_valid_sides = j / comp_count; 3020 } 3021 } 3022 return num_valid_sides; 3023 } 3024 3025 int64_t ParallelDatabaseIO::get_side_connectivity(const Ioss::SideBlock *fb, int64_t id, 3026 int64_t /*unused*/, void * fconnect, 3027 bool map_ids) const 3028 { 3029 // Get size of data stored on the file... 3030 ex_set set_param[1]; 3031 set_param[0].id = id; 3032 set_param[0].type = EX_SIDE_SET; 3033 set_param[0].entry_list = nullptr; 3034 set_param[0].extra_list = nullptr; 3035 set_param[0].distribution_factor_list = nullptr; 3036 int ierr = ex_get_sets(get_file_pointer(), 1, set_param); 3037 if (ierr < 0) { 3038 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3039 } 3040 3041 int64_t number_sides = set_param[0].num_entry; 3042 3043 // Allocate space for element and local side number 3044 assert(number_sides > 0); 3045 3046 // Need a larger vector to get the entire sideset and then filter down to the correct size... 3047 std::vector<char> element(number_sides * int_byte_size_api()); 3048 std::vector<char> side(number_sides * int_byte_size_api()); 3049 if (int_byte_size_api() == 4) { 3050 Ioss::Field elem_field("ids", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, number_sides); 3051 Ioss::Field side_field("sides", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, 3052 number_sides); 3053 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 3054 reinterpret_cast<int *>(element.data())); 3055 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 3056 reinterpret_cast<int *>(side.data())); 3057 } 3058 else { 3059 Ioss::Field elem_field("ids", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, number_sides); 3060 Ioss::Field side_field("sides", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, number_sides); 3061 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 3062 reinterpret_cast<int64_t *>(element.data())); 3063 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 3064 reinterpret_cast<int64_t *>(side.data())); 3065 } 3066 3067 Ioss::IntVector is_valid_side; 3068 Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(), 3069 (void *)element.data(), (void *)side.data(), 3070 number_sides, get_region()); 3071 3072 std::vector<char> elconnect; 3073 int64_t elconsize = 0; // Size of currently allocated connectivity block 3074 Ioss::ElementBlock *conn_block = nullptr; // Block that we currently 3075 // have connectivity for 3076 3077 Ioss::ElementBlock *block = nullptr; 3078 3079 int * element32 = nullptr; 3080 int64_t *element64 = nullptr; 3081 int * side32 = nullptr; 3082 int64_t *side64 = nullptr; 3083 3084 int * elconn32 = nullptr; 3085 int64_t *elconn64 = nullptr; 3086 int * fconn32 = nullptr; 3087 int64_t *fconn64 = nullptr; 3088 3089 if (int_byte_size_api() == 4) { 3090 element32 = reinterpret_cast<int *>(element.data()); 3091 side32 = reinterpret_cast<int *>(side.data()); 3092 fconn32 = reinterpret_cast<int *>(fconnect); 3093 } 3094 else { 3095 element64 = reinterpret_cast<int64_t *>(element.data()); 3096 side64 = reinterpret_cast<int64_t *>(side.data()); 3097 fconn64 = reinterpret_cast<int64_t *>(fconnect); 3098 } 3099 3100 Ioss::IntVector side_elem_map; // Maps the side into the elements 3101 // connectivity array 3102 int64_t current_side = -1; 3103 int nelnode = 0; 3104 int nfnodes = 0; 3105 int ieb = 0; 3106 size_t offset = 0; 3107 for (int64_t iel = 0; iel < number_sides; iel++) { 3108 if (is_valid_side[iel] == 1) { 3109 3110 int64_t elem_id = 0; 3111 if (int_byte_size_api() == 4) { 3112 elem_id = element32[iel]; 3113 } 3114 else { 3115 elem_id = element64[iel]; 3116 } 3117 3118 // ensure we have correct connectivity 3119 block = get_region()->get_element_block(elem_id); 3120 if (conn_block != block) { 3121 ssize_t nelem = block->entity_count(); 3122 nelnode = block->topology()->number_nodes(); 3123 // Used to map element number into position in connectivity array. 3124 // E.g., element 97 is the (97-offset)th element in this block and 3125 // is stored in array index (97-offset-1). 3126 offset = block->get_offset() + 1; 3127 if (elconsize < nelem * nelnode) { 3128 elconsize = nelem * nelnode; 3129 elconnect.resize(elconsize * int_byte_size_api()); 3130 if (int_byte_size_api() == 4) { 3131 elconn32 = reinterpret_cast<int *>(elconnect.data()); 3132 } 3133 else { 3134 elconn64 = reinterpret_cast<int64_t *>(elconnect.data()); 3135 } 3136 } 3137 if (map_ids) { 3138 get_field_internal(block, block->get_field("connectivity"), elconnect.data(), 3139 nelem * nelnode * int_byte_size_api()); 3140 } 3141 else { 3142 get_field_internal(block, block->get_field("connectivity_raw"), elconnect.data(), 3143 nelem * nelnode * int_byte_size_api()); 3144 } 3145 conn_block = block; 3146 current_side = -1; 3147 } 3148 3149 // NOTE: Element connectivity is returned with nodes in global id space if "map_ids" false, 3150 // otherwise it is in local space. 3151 int64_t side_id = 0; 3152 if (int_byte_size_api() == 4) { 3153 side_id = side32[iel]; 3154 } 3155 else { 3156 side_id = side64[iel]; 3157 } 3158 3159 if (current_side != side_id) { 3160 side_elem_map = block->topology()->boundary_connectivity(side_id); 3161 current_side = side_id; 3162 nfnodes = block->topology()->boundary_type(side_id)->number_nodes(); 3163 } 3164 for (int inode = 0; inode < nfnodes; inode++) { 3165 size_t index = (elem_id - offset) * nelnode + side_elem_map[inode]; 3166 if (int_byte_size_api() == 4) { 3167 fconn32[ieb++] = elconn32[index]; 3168 } 3169 else { 3170 fconn64[ieb++] = elconn64[index]; 3171 } 3172 } 3173 } 3174 } 3175 return ierr; 3176 } 3177 3178 // Get distribution factors for the specified side block 3179 int64_t ParallelDatabaseIO::get_side_distributions(const Ioss::SideBlock *fb, int64_t id, 3180 int64_t my_side_count, double *dist_fact, 3181 size_t /* data_size */) const 3182 { 3183 // Allocate space for elements and local side numbers 3184 // Get size of data stored on the file... 3185 3186 auto & set = decomp->get_decomp_set(EX_SIDE_SET, id); 3187 int64_t number_sides = set.ioss_count(); 3188 int64_t number_distribution_factors = set.df_count(); 3189 3190 const Ioss::ElementTopology *ftopo = fb->topology(); 3191 int nfnodes = ftopo->number_nodes(); 3192 3193 if (set.distributionFactorConstant) { 3194 // Fill in the array with the constant value... 3195 for (int64_t i = 0; i < nfnodes * my_side_count; i++) { 3196 dist_fact[i] = set.distributionFactorValue; 3197 } 3198 return 0; 3199 } 3200 3201 // Take care of the easy situation -- If 'side_count' == 3202 // 'number_sides' then the sideset is stored in a single sideblock 3203 // and all distribution factors on the database are transferred 3204 // 1-to-1 into 'dist_fact' array. 3205 int64_t entity_count = fb->entity_count(); 3206 if (fb->owner()->block_count() == 1 && number_sides == entity_count) { 3207 assert(number_sides == 0 || number_distribution_factors % number_sides == 0); 3208 assert(number_sides == 0 || number_distribution_factors / number_sides == nfnodes); 3209 if (number_sides * nfnodes != number_distribution_factors && 3210 number_sides != number_distribution_factors) { 3211 std::ostringstream errmsg; 3212 fmt::print(errmsg, 3213 "ERROR: SideBlock '{}' has incorrect distribution factor count.\n" 3214 "\tThere are {} '{}' sides with {} nodes per side, but there are {}" 3215 " distribution factors which is not correct.\n" 3216 "\tThere should be either {} or {} distribution factors.\n", 3217 fb->name(), number_sides, ftopo->name(), nfnodes, number_distribution_factors, 3218 number_sides, number_sides * nfnodes); 3219 IOSS_ERROR(errmsg); 3220 } 3221 std::string storage = "Real[" + std::to_string(nfnodes) + "]"; 3222 Ioss::Field dist("distribution_factors", Ioss::Field::REAL, storage, Ioss::Field::MESH, 3223 number_sides); 3224 decomp->get_set_mesh_double(get_file_pointer(), EX_SIDE_SET, id, dist, dist_fact); 3225 return 0; 3226 } 3227 3228 std::string storage = "Real[" + std::to_string(nfnodes) + "]"; 3229 Ioss::Field field("distribution_factors", Ioss::Field::REAL, storage, Ioss::Field::MESH, 3230 number_distribution_factors / nfnodes); 3231 std::vector<double> dist(number_distribution_factors); 3232 decomp->get_set_mesh_double(get_file_pointer(), EX_SIDE_SET, id, field, dist.data()); 3233 3234 // Another easy situation (and common for exodusII) is if the input 3235 // distribution factors are all the same value (typically 1). In 3236 // that case, we only have to fill in the output array with that 3237 // value. 3238 { 3239 double value = number_distribution_factors > 0 ? dist[0] : 0.0; 3240 bool constant = true; 3241 for (auto df : dist) { 3242 if (df != value) { 3243 constant = false; 3244 break; 3245 } 3246 } 3247 3248 constant = util().global_minmax(constant ? 1 : 0, Ioss::ParallelUtils::DO_MIN); 3249 3250 if (constant) { 3251 if (value == 0.0) { 3252 value = 1.0; // Take care of some buggy mesh generators 3253 } 3254 for (ssize_t j = 0; j < my_side_count * nfnodes; j++) { 3255 dist_fact[j] = value; 3256 } 3257 return 0; 3258 } 3259 } 3260 3261 // If we get to here, the underlying sideset contains multiple side 3262 // topologies and the distribution factors are non-constant. Need to 3263 // allocate space to store all distribution factors and then pull 3264 // out those that are applied to sides with the correct topology. 3265 3266 // Allocate space for element and local side number (this is bulk 3267 // data...) 3268 //---- 3269 std::vector<char> element(number_sides * int_byte_size_api()); 3270 std::vector<char> sides(number_sides * int_byte_size_api()); 3271 if (int_byte_size_api() == 4) { 3272 Ioss::Field elem_field("ids", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, number_sides); 3273 Ioss::Field side_field("sides", Ioss::Field::INTEGER, "scalar", Ioss::Field::MESH, 3274 number_sides); 3275 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 3276 reinterpret_cast<int *>(element.data())); 3277 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 3278 reinterpret_cast<int *>(sides.data())); 3279 } 3280 else { 3281 Ioss::Field elem_field("ids", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, number_sides); 3282 Ioss::Field side_field("sides", Ioss::Field::INT64, "scalar", Ioss::Field::MESH, number_sides); 3283 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, elem_field, 3284 reinterpret_cast<int64_t *>(element.data())); 3285 decomp->get_set_mesh_var(get_file_pointer(), EX_SIDE_SET, id, side_field, 3286 reinterpret_cast<int64_t *>(sides.data())); 3287 } 3288 //---- 3289 3290 Ioss::IntVector is_valid_side; 3291 Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(), 3292 element.data(), sides.data(), number_sides, 3293 get_region()); 3294 3295 int64_t ieb = 0; // counter for distribution factors in this sideblock 3296 int64_t idb = 0; // counter for distribution factors read from database 3297 Ioss::ElementBlock *block = nullptr; 3298 3299 int * element32 = nullptr; 3300 int64_t *element64 = nullptr; 3301 int * side32 = nullptr; 3302 int64_t *side64 = nullptr; 3303 3304 if (int_byte_size_api() == 4) { 3305 element32 = reinterpret_cast<int *>(element.data()); 3306 side32 = reinterpret_cast<int *>(sides.data()); 3307 } 3308 else { 3309 element64 = reinterpret_cast<int64_t *>(element.data()); 3310 side64 = reinterpret_cast<int64_t *>(sides.data()); 3311 } 3312 3313 for (int64_t iel = 0; iel < number_sides; iel++) { 3314 int64_t elem_id = 0; 3315 int64_t side_id = 0; 3316 if (int_byte_size_api() == 4) { 3317 elem_id = element32[iel]; 3318 side_id = side32[iel]; 3319 } 3320 else { 3321 elem_id = element64[iel]; 3322 side_id = side64[iel]; 3323 } 3324 3325 if (block == nullptr || !block->contains(elem_id)) { 3326 block = get_region()->get_element_block(elem_id); 3327 } 3328 3329 if (block == nullptr) { 3330 std::ostringstream errmsg; 3331 fmt::print(errmsg, 3332 "INTERNAL ERROR: Could not find element block containing element with id {}." 3333 " Something is wrong in the Ioex::ParallelDatabaseIO class. Please report.\n", 3334 elem_id); 3335 IOSS_ERROR(errmsg); 3336 } 3337 3338 const Ioss::ElementTopology *topo = block->topology()->boundary_type(side_id); 3339 3340 if (topo == nullptr) { 3341 std::ostringstream errmsg; 3342 fmt::print(errmsg, 3343 "INTERNAL ERROR: Could not find topology of element block boundary. " 3344 "Something is wrong in the Ioex::ParallelDatabaseIO class. Please report.\n"); 3345 IOSS_ERROR(errmsg); 3346 } 3347 3348 int nside_nodes = topo->number_nodes(); 3349 3350 if (is_valid_side[iel] == 1) { 3351 // This side belongs in the sideblock 3352 for (int64_t i = 0; i < nside_nodes; i++) { 3353 dist_fact[ieb++] = dist[idb++]; 3354 } 3355 } 3356 else { 3357 // Skip over unused 'dist' factors 3358 idb += topo->number_nodes(); 3359 } 3360 } 3361 3362 assert(ieb == my_side_count * nfnodes); 3363 // If the following assert fails, it may be due to bug in Patran 3364 // which writes too many distribution factors to the database in a 3365 // mixed element case. Note that this is checked earlier also with a 3366 // better error message. 3367 assert(idb == number_distribution_factors); 3368 return 0; 3369 } 3370 3371 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::Region *reg, const Ioss::Field &field, 3372 void *data, size_t data_size) const 3373 { 3374 return Ioex::BaseDatabaseIO::put_field_internal(reg, field, data, data_size); 3375 } 3376 3377 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::NodeBlock *nb, const Ioss::Field &field, 3378 void *data, size_t data_size) const 3379 { 3380 size_t num_to_get = field.verify(data_size); 3381 3382 size_t proc_offset = nb->get_optional_property("_processor_offset", 0); 3383 size_t file_count = nb->get_optional_property("locally_owned_count", num_to_get); 3384 3385 Ioss::Field::RoleType role = field.get_role(); 3386 3387 if (role == Ioss::Field::MESH) { 3388 if (field.get_name() == "owning_processor") { 3389 // Set the nodeOwningProcessor vector for all nodes on this processor. 3390 // Value is the processor that owns the node. 3391 3392 // NOTE: The owning_processor field is always int32 3393 nodeOwningProcessor.reserve(num_to_get); 3394 int *owned = (int *)data; 3395 for (size_t i = 0; i < num_to_get; i++) { 3396 nodeOwningProcessor.push_back(owned[i]); 3397 } 3398 3399 // Now create the "implicit local" to "implicit global" 3400 // map which maps data from its local implicit position 3401 // to its implicit (1..num_global_node) position in the 3402 // global file. This is needed for the global-to-local 3403 // mapping of element connectivity and nodeset nodelists. 3404 create_implicit_global_map(); 3405 } 3406 3407 else if (field.get_name() == "mesh_model_coordinates_x") { 3408 double * rdata = static_cast<double *>(data); 3409 std::vector<double> file_data; 3410 file_data.reserve(file_count); 3411 check_node_owning_processor_data(nodeOwningProcessor, file_count); 3412 filter_owned_nodes(nodeOwningProcessor, myProcessor, rdata, file_data); 3413 3414 int ierr = ex_put_partial_coord_component(get_file_pointer(), proc_offset + 1, file_count, 1, 3415 file_data.data()); 3416 if (ierr < 0) { 3417 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3418 } 3419 } 3420 3421 else if (field.get_name() == "mesh_model_coordinates_y") { 3422 double * rdata = static_cast<double *>(data); 3423 std::vector<double> file_data; 3424 file_data.reserve(file_count); 3425 check_node_owning_processor_data(nodeOwningProcessor, file_count); 3426 filter_owned_nodes(nodeOwningProcessor, myProcessor, rdata, file_data); 3427 int ierr = ex_put_partial_coord_component(get_file_pointer(), proc_offset + 1, file_count, 2, 3428 file_data.data()); 3429 if (ierr < 0) { 3430 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3431 } 3432 } 3433 3434 else if (field.get_name() == "mesh_model_coordinates_z") { 3435 double * rdata = static_cast<double *>(data); 3436 std::vector<double> file_data; 3437 file_data.reserve(file_count); 3438 check_node_owning_processor_data(nodeOwningProcessor, file_count); 3439 filter_owned_nodes(nodeOwningProcessor, myProcessor, rdata, file_data); 3440 int ierr = ex_put_partial_coord_component(get_file_pointer(), proc_offset + 1, file_count, 3, 3441 file_data.data()); 3442 if (ierr < 0) { 3443 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3444 } 3445 } 3446 3447 else if (field.get_name() == "mesh_model_coordinates") { 3448 // Data required by upper classes store x0, y0, z0, ... xn, yn, zn 3449 // Data stored in exodusII file is x0, ..., xn, y0, ..., yn, z0, ..., zn 3450 // so we have to allocate some scratch memory to read in the data 3451 // and then map into supplied 'data' 3452 std::vector<double> x; 3453 std::vector<double> y; 3454 std::vector<double> z; 3455 3456 x.reserve(file_count > 0 ? file_count : 1); 3457 if (spatialDimension > 1) { 3458 y.reserve(file_count > 0 ? file_count : 1); 3459 } 3460 if (spatialDimension == 3) { 3461 z.reserve(file_count > 0 ? file_count : 1); 3462 } 3463 3464 // Cast 'data' to correct size -- double 3465 double *rdata = static_cast<double *>(data); 3466 check_node_owning_processor_data(nodeOwningProcessor, file_count); 3467 filter_owned_nodes(nodeOwningProcessor, myProcessor, rdata, x, 0, spatialDimension); 3468 if (spatialDimension > 1) { 3469 filter_owned_nodes(nodeOwningProcessor, myProcessor, rdata, y, 1, spatialDimension); 3470 } 3471 if (spatialDimension == 3) { 3472 filter_owned_nodes(nodeOwningProcessor, myProcessor, rdata, z, 2, spatialDimension); 3473 } 3474 3475 int ierr = ex_put_partial_coord(get_file_pointer(), proc_offset + 1, file_count, x.data(), 3476 y.data(), z.data()); 3477 if (ierr < 0) { 3478 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3479 } 3480 } 3481 else if (field.get_name() == "ids") { 3482 // The ids coming in are the global ids; their position is the 3483 // local id -1 (That is, data[0] contains the global id of local 3484 // node 1) 3485 3486 // Another 'const-cast' since we are modifying the database just 3487 // for efficiency; which the client does not see... 3488 handle_node_ids(data, num_to_get, proc_offset, file_count); 3489 } 3490 else if (field.get_name() == "connectivity") { 3491 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 3492 } 3493 else if (field.get_name() == "connectivity_raw") { 3494 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 3495 } 3496 else if (field.get_name() == "node_connectivity_status") { 3497 // Do nothing, input only field. 3498 } 3499 else if (field.get_name() == "implicit_ids") { 3500 // Do nothing, input only field. 3501 } 3502 else { 3503 return Ioss::Utils::field_warning(nb, field, "mesh output"); 3504 } 3505 } 3506 else if (role == Ioss::Field::TRANSIENT) { 3507 // Check if the specified field exists on this node block. 3508 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 3509 // exist on the database as scalars with the appropriate 3510 // extensions. 3511 3512 // Transfer each component of the variable into 'data' and then 3513 // output. Need temporary storage area of size 'number of 3514 // nodes in this block. 3515 write_nodal_transient_field(EX_NODE_BLOCK, field, nb, num_to_get, data); 3516 } 3517 else if (role == Ioss::Field::REDUCTION) { 3518 store_reduction_field(EX_NODE_BLOCK, field, nb, data); 3519 } 3520 return num_to_get; 3521 } 3522 3523 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::Blob *blob, const Ioss::Field &field, 3524 void *data, size_t data_size) const 3525 { 3526 { 3527 Ioss::SerializeIO serializeIO__(this); 3528 3529 size_t num_to_get = field.verify(data_size); 3530 if (num_to_get > 0) { 3531 3532 Ioss::Field::RoleType role = field.get_role(); 3533 3534 if (role == Ioss::Field::MESH) { 3535 if (field.get_name() == "ids") { 3536 // The ids coming in are the global ids; their position is the 3537 // local id -1 (That is, data[0] contains the global id of local 3538 // node 1) 3539 // handle_node_ids(data, num_to_get); 3540 } 3541 else if (field.get_name() == "connectivity") { 3542 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 3543 } 3544 else if (field.get_name() == "connectivity_raw") { 3545 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 3546 } 3547 else if (field.get_name() == "node_connectivity_status") { 3548 // Do nothing, input only field. 3549 } 3550 else if (field.get_name() == "implicit_ids") { 3551 // Do nothing, input only field. 3552 } 3553 else { 3554 return Ioss::Utils::field_warning(blob, field, "mesh output"); 3555 } 3556 } 3557 else if (role == Ioss::Field::TRANSIENT) { 3558 // Check if the specified field exists on this node block. 3559 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 3560 // exist on the database as scalars with the appropriate 3561 // extensions. 3562 3563 // Transfer each component of the variable into 'data' and then 3564 // output. Need temporary storage area of size 'number of 3565 // nodes in this block. 3566 write_entity_transient_field(EX_BLOB, field, blob, num_to_get, data); 3567 } 3568 else if (role == Ioss::Field::REDUCTION) { 3569 store_reduction_field(EX_BLOB, field, blob, data); 3570 } 3571 else if (role == Ioss::Field::ATTRIBUTE) { 3572 num_to_get = write_attribute_field(EX_BLOB, field, blob, data); 3573 } 3574 } 3575 return num_to_get; 3576 } 3577 } 3578 3579 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::Assembly *assembly, 3580 const Ioss::Field &field, void *data, 3581 size_t data_size) const 3582 { 3583 { 3584 Ioss::SerializeIO serializeIO__(this); 3585 3586 size_t num_to_get = field.verify(data_size); 3587 if (num_to_get > 0) { 3588 3589 Ioss::Field::RoleType role = field.get_role(); 3590 3591 if (role == Ioss::Field::MESH) { 3592 if (field.get_name() == "ids") { 3593 // The ids coming in are the global ids; their position is the 3594 // local id -1 (That is, data[0] contains the global id of local 3595 // node 1) 3596 // handle_node_ids(data, num_to_get); 3597 } 3598 else if (field.get_name() == "connectivity") { 3599 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 3600 } 3601 else if (field.get_name() == "connectivity_raw") { 3602 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 3603 } 3604 else if (field.get_name() == "node_connectivity_status") { 3605 // Do nothing, input only field. 3606 } 3607 else if (field.get_name() == "implicit_ids") { 3608 // Do nothing, input only field. 3609 } 3610 else { 3611 return Ioss::Utils::field_warning(assembly, field, "mesh output"); 3612 } 3613 } 3614 else if (role == Ioss::Field::TRANSIENT) { 3615 // Check if the specified field exists on this node block. 3616 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 3617 // exist on the database as scalars with the appropriate 3618 // extensions. 3619 3620 // Transfer each component of the variable into 'data' and then 3621 // output. Need temporary storage area of size 'number of 3622 // nodes in this block. 3623 write_entity_transient_field(EX_ASSEMBLY, field, assembly, num_to_get, data); 3624 } 3625 else if (role == Ioss::Field::REDUCTION) { 3626 store_reduction_field(EX_ASSEMBLY, field, assembly, data); 3627 } 3628 else if (role == Ioss::Field::ATTRIBUTE) { 3629 num_to_get = write_attribute_field(EX_ASSEMBLY, field, assembly, data); 3630 } 3631 } 3632 return num_to_get; 3633 } 3634 } 3635 3636 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::ElementBlock *eb, 3637 const Ioss::Field &field, void *data, 3638 size_t data_size) const 3639 { 3640 size_t num_to_get = field.verify(data_size); 3641 3642 int ierr = 0; 3643 3644 // Get the element block id and element count 3645 int64_t id = Ioex::get_id(eb, EX_ELEM_BLOCK, &ids_); 3646 int64_t my_element_count = eb->entity_count(); 3647 Ioss::Field::RoleType role = field.get_role(); 3648 3649 size_t proc_offset = eb->get_optional_property("_processor_offset", 0); 3650 size_t file_count = eb->get_optional_property("locally_owned_count", num_to_get); 3651 3652 if (role == Ioss::Field::MESH) { 3653 // Handle the MESH fields required for an ExodusII file model. 3654 // (The 'genesis' portion) 3655 if (field.get_name() == "connectivity") { 3656 // Map element connectivity from global node id to local node id. 3657 int element_nodes = eb->topology()->number_nodes(); 3658 3659 // Maps global to local 3660 nodeMap.reverse_map_data(data, field, num_to_get * element_nodes); 3661 3662 // Maps local to "global_implicit" 3663 if (int_byte_size_api() == 4) { 3664 map_local_to_global_implicit(reinterpret_cast<int *>(data), num_to_get * element_nodes, 3665 nodeGlobalImplicitMap); 3666 } 3667 else { 3668 map_local_to_global_implicit(reinterpret_cast<int64_t *>(data), num_to_get * element_nodes, 3669 nodeGlobalImplicitMap); 3670 } 3671 3672 ierr = ex_put_partial_conn(get_file_pointer(), EX_ELEM_BLOCK, id, proc_offset + 1, file_count, 3673 data, nullptr, nullptr); 3674 if (ierr < 0) { 3675 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3676 } 3677 } 3678 else if (field.get_name() == "connectivity_edge") { 3679 // Map element connectivity from global edge id to local edge id. 3680 int element_edges = field.transformed_storage()->component_count(); 3681 edgeMap.reverse_map_data(data, field, num_to_get * element_edges); 3682 ierr = ex_put_conn(get_file_pointer(), EX_ELEM_BLOCK, id, nullptr, data, nullptr); 3683 if (ierr < 0) { 3684 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3685 } 3686 } 3687 else if (field.get_name() == "connectivity_face") { 3688 // Map element connectivity from global face id to local face id. 3689 int element_faces = field.transformed_storage()->component_count(); 3690 faceMap.reverse_map_data(data, field, num_to_get * element_faces); 3691 ierr = ex_put_conn(get_file_pointer(), EX_ELEM_BLOCK, id, nullptr, nullptr, data); 3692 if (ierr < 0) { 3693 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3694 } 3695 } 3696 else if (field.get_name() == "connectivity_raw") { 3697 // Element connectivity is already in local node id, map local to "global_implicit" 3698 int element_nodes = eb->topology()->number_nodes(); 3699 if (int_byte_size_api() == 4) { 3700 map_local_to_global_implicit(reinterpret_cast<int *>(data), num_to_get * element_nodes, 3701 nodeGlobalImplicitMap); 3702 } 3703 else { 3704 map_local_to_global_implicit(reinterpret_cast<int64_t *>(data), num_to_get * element_nodes, 3705 nodeGlobalImplicitMap); 3706 } 3707 3708 ierr = ex_put_partial_conn(get_file_pointer(), EX_ELEM_BLOCK, id, proc_offset + 1, file_count, 3709 data, nullptr, nullptr); 3710 if (ierr < 0) { 3711 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3712 } 3713 } 3714 else if (field.get_name() == "ids") { 3715 size_t glob_map_offset = eb->get_property("global_map_offset").get_int(); 3716 handle_element_ids(eb, data, num_to_get, glob_map_offset + proc_offset, file_count); 3717 } 3718 else if (field.get_name() == "implicit_ids") { 3719 // Do nothing, input only field. 3720 } 3721 else if (field.get_name() == "skin") { 3722 // This is (currently) for the skinned body. It maps the 3723 // side element on the skin to the original element/local 3724 // side number. It is a two component field, the first 3725 // component is the global id of the underlying element in 3726 // the initial mesh and its local side number (1-based). 3727 3728 // FIX: Hardwired map ids.... 3729 int map_count = ex_inquire_int(get_file_pointer(), EX_INQ_ELEM_MAP); 3730 if (map_count == 0) { 3731 // This needs to be fixed... Currently hardwired.... 3732 ierr = ex_put_map_param(get_file_pointer(), 0, 2); 3733 if (ierr < 0) { 3734 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3735 } 3736 } 3737 3738 std::vector<char> element(my_element_count * int_byte_size_api()); 3739 std::vector<char> side(my_element_count * int_byte_size_api()); 3740 3741 if (int_byte_size_api() == 4) { 3742 int *el_side = reinterpret_cast<int *>(data); 3743 int *element32 = reinterpret_cast<int *>(element.data()); 3744 int *side32 = reinterpret_cast<int *>(side.data()); 3745 3746 int index = 0; 3747 for (int i = 0; i < my_element_count; i++) { 3748 element32[i] = el_side[index++]; 3749 side32[i] = el_side[index++]; 3750 } 3751 } 3752 else { 3753 int64_t *el_side = reinterpret_cast<int64_t *>(data); 3754 int64_t *element64 = reinterpret_cast<int64_t *>(element.data()); 3755 int64_t *side64 = reinterpret_cast<int64_t *>(side.data()); 3756 3757 int64_t index = 0; 3758 for (int64_t i = 0; i < my_element_count; i++) { 3759 element64[i] = el_side[index++]; 3760 side64[i] = el_side[index++]; 3761 } 3762 } 3763 3764 size_t eb_offset = eb->get_offset() + proc_offset; 3765 ierr = ex_put_partial_num_map(get_file_pointer(), EX_ELEM_MAP, 1, eb_offset + 1, file_count, 3766 element.data()); 3767 if (ierr < 0) { 3768 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3769 } 3770 3771 ierr = ex_put_partial_num_map(get_file_pointer(), EX_ELEM_MAP, 2, eb_offset + 1, file_count, 3772 side.data()); 3773 if (ierr < 0) { 3774 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3775 } 3776 3777 if (map_count == 0) { 3778 // NOTE: ex_put_*num_map must be called prior to defining the name... 3779 ierr = ex_put_name(get_file_pointer(), EX_ELEM_MAP, 1, "skin:parent_element_id"); 3780 if (ierr < 0) { 3781 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3782 } 3783 3784 ierr = ex_put_name(get_file_pointer(), EX_ELEM_MAP, 2, "skin:parent_element_side_number"); 3785 if (ierr < 0) { 3786 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3787 } 3788 } 3789 } 3790 else { 3791 num_to_get = Ioss::Utils::field_warning(eb, field, "mesh output"); 3792 } 3793 } 3794 else if (role == Ioss::Field::ATTRIBUTE) { 3795 num_to_get = write_attribute_field(EX_ELEM_BLOCK, field, eb, data); 3796 } 3797 else if (role == Ioss::Field::TRANSIENT) { 3798 // Check if the specified field exists on this element block. 3799 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 3800 // exist on the database as scalars with the appropriate 3801 // extensions. 3802 3803 // Transfer each component of the variable into 'data' and then 3804 // output. Need temporary storage area of size 'number of 3805 // elements in this block. 3806 auto global_entity_count = eb->get_property("global_entity_count").get_int(); 3807 if (global_entity_count > 0) { 3808 write_entity_transient_field(EX_ELEM_BLOCK, field, eb, my_element_count, data); 3809 } 3810 } 3811 else if (role == Ioss::Field::REDUCTION) { 3812 store_reduction_field(EX_ELEM_BLOCK, field, eb, data); 3813 } 3814 return num_to_get; 3815 } 3816 3817 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::FaceBlock *eb, const Ioss::Field &field, 3818 void *data, size_t data_size) const 3819 { 3820 size_t num_to_get = field.verify(data_size); 3821 3822 int ierr = 0; 3823 3824 // Get the face block id and face count 3825 int64_t id = Ioex::get_id(eb, EX_FACE_BLOCK, &ids_); 3826 int64_t my_face_count = eb->entity_count(); 3827 Ioss::Field::RoleType role = field.get_role(); 3828 3829 if (role == Ioss::Field::MESH) { 3830 // Handle the MESH fields required for an ExodusII file model. 3831 // (The 'genesis' portion) 3832 if (field.get_name() == "connectivity") { 3833 if (my_face_count > 0) { 3834 // Map face connectivity from global node id to local node id. 3835 int face_nodes = eb->topology()->number_nodes(); 3836 nodeMap.reverse_map_data(data, field, num_to_get * face_nodes); 3837 ierr = ex_put_conn(get_file_pointer(), EX_FACE_BLOCK, id, data, nullptr, nullptr); 3838 if (ierr < 0) { 3839 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3840 } 3841 } 3842 } 3843 else if (field.get_name() == "connectivity_edge") { 3844 if (my_face_count > 0) { 3845 // Map face connectivity from global edge id to local edge id. 3846 // Do it in 'data' ... 3847 int face_edges = field.transformed_storage()->component_count(); 3848 edgeMap.reverse_map_data(data, field, num_to_get * face_edges); 3849 ierr = ex_put_conn(get_file_pointer(), EX_FACE_BLOCK, id, nullptr, data, nullptr); 3850 if (ierr < 0) { 3851 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3852 } 3853 } 3854 } 3855 else if (field.get_name() == "connectivity_raw") { 3856 // Do nothing, input only field. 3857 } 3858 else if (field.get_name() == "ids") { 3859 handle_face_ids(eb, data, num_to_get); 3860 } 3861 else { 3862 num_to_get = Ioss::Utils::field_warning(eb, field, "mesh output"); 3863 } 3864 } 3865 else if (role == Ioss::Field::ATTRIBUTE) { 3866 num_to_get = write_attribute_field(EX_FACE_BLOCK, field, eb, data); 3867 } 3868 else if (role == Ioss::Field::TRANSIENT) { 3869 // Check if the specified field exists on this face block. 3870 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 3871 // exist on the database as scalars with the appropriate 3872 // extensions. 3873 3874 // Transfer each component of the variable into 'data' and then 3875 // output. Need temporary storage area of size 'number of 3876 // faces in this block. 3877 write_entity_transient_field(EX_FACE_BLOCK, field, eb, my_face_count, data); 3878 } 3879 else if (role == Ioss::Field::REDUCTION) { 3880 store_reduction_field(EX_FACE_BLOCK, field, eb, data); 3881 } 3882 return num_to_get; 3883 } 3884 3885 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::EdgeBlock *eb, const Ioss::Field &field, 3886 void *data, size_t data_size) const 3887 { 3888 size_t num_to_get = field.verify(data_size); 3889 3890 int ierr = 0; 3891 3892 // Get the edge block id and edge count 3893 int64_t id = Ioex::get_id(eb, EX_EDGE_BLOCK, &ids_); 3894 int64_t my_edge_count = eb->entity_count(); 3895 Ioss::Field::RoleType role = field.get_role(); 3896 3897 if (role == Ioss::Field::MESH) { 3898 // Handle the MESH fields required for an ExodusII file model. (The 'genesis' portion) 3899 if (field.get_name() == "connectivity") { 3900 if (my_edge_count > 0) { 3901 // Map edge connectivity from global node id to local node id. 3902 int edge_nodes = eb->topology()->number_nodes(); 3903 nodeMap.reverse_map_data(data, field, num_to_get * edge_nodes); 3904 ierr = ex_put_conn(get_file_pointer(), EX_EDGE_BLOCK, id, data, nullptr, nullptr); 3905 if (ierr < 0) { 3906 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 3907 } 3908 } 3909 } 3910 else if (field.get_name() == "connectivity_raw") { 3911 // Do nothing, input only field. 3912 } 3913 else if (field.get_name() == "ids") { 3914 handle_edge_ids(eb, data, num_to_get); 3915 } 3916 else { 3917 num_to_get = Ioss::Utils::field_warning(eb, field, "mesh output"); 3918 } 3919 } 3920 else if (role == Ioss::Field::ATTRIBUTE) { 3921 num_to_get = write_attribute_field(EX_EDGE_BLOCK, field, eb, data); 3922 } 3923 else if (role == Ioss::Field::TRANSIENT) { 3924 // Check if the specified field exists on this edge block. 3925 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 3926 // exist on the database as scalars with the appropriate 3927 // extensions. 3928 3929 // Transfer each component of the variable into 'data' and then 3930 // output. Need temporary storage area of size 'number of 3931 // edges in this block. 3932 write_entity_transient_field(EX_EDGE_BLOCK, field, eb, my_edge_count, data); 3933 } 3934 else if (role == Ioss::Field::REDUCTION) { 3935 store_reduction_field(EX_EDGE_BLOCK, field, eb, data); 3936 } 3937 return num_to_get; 3938 } 3939 3940 int64_t ParallelDatabaseIO::handle_node_ids(void *ids, int64_t num_to_get, size_t /* offset */, 3941 size_t /*count*/) const 3942 { 3943 /*! 3944 * There are two modes we need to support in this routine: 3945 * 1. Initial definition of node map (local->global) and 3946 * nodeMap.reverse (global->local). 3947 * 2. Redefinition of node map via 'reordering' of the original 3948 * map when the nodes on this processor are the same, but their 3949 * order is changed (or count because of ghosting) 3950 * 3951 * So, there will be two maps the 'nodeMap.map' map is a 'direct lookup' 3952 * map which maps current local position to global id and the 3953 * 'nodeMap.reverse' is an associative lookup which maps the 3954 * global id to 'original local'. There is also a 3955 * 'nodeMap.reorder' which is direct lookup and maps current local 3956 * position to original local. 3957 3958 * The ids coming in are the global ids; their position is the 3959 * "local id-1" (That is, data[0] contains the global id of local 3960 * node 1 in this node block). 3961 * 3962 * int local_position = nodeMap.reverse[NodeMap[i+1]] 3963 * (the nodeMap.map and nodeMap.reverse are 1-based) 3964 * 3965 * To determine which map to update on a call to this function, we 3966 * use the following hueristics: 3967 * -- If the database state is 'STATE_MODEL:', then update the 3968 * 'nodeMap.reverse' and 'nodeMap.map' 3969 * 3970 * -- If the database state is not STATE_MODEL, then leave the 3971 * 'nodeMap.reverse' and 'nodeMap.map' alone since they correspond to the 3972 * information already written to the database. [May want to add a 3973 * STATE_REDEFINE_MODEL] 3974 * 3975 * -- In both cases, update the nodeMap.reorder 3976 * 3977 * NOTE: The mapping is done on TRANSIENT fields only; MODEL fields 3978 * should be in the original order... 3979 */ 3980 nodeMap.set_size(num_to_get); 3981 3982 bool in_define = (dbState == Ioss::STATE_MODEL) || (dbState == Ioss::STATE_DEFINE_MODEL); 3983 if (int_byte_size_api() == 4) { 3984 nodeMap.set_map(static_cast<int *>(ids), num_to_get, 0, in_define); 3985 } 3986 else { 3987 nodeMap.set_map(static_cast<int64_t *>(ids), num_to_get, 0, in_define); 3988 } 3989 3990 nodeMap.set_defined(true); 3991 return num_to_get; 3992 } 3993 3994 int64_t ParallelDatabaseIO::handle_element_ids(const Ioss::ElementBlock *eb, void *ids, 3995 size_t num_to_get, size_t offset, size_t count) const 3996 { 3997 if (dbState == Ioss::STATE_MODEL) { 3998 if (elemGlobalImplicitMap.empty()) { 3999 elemGlobalImplicitMap.resize(elementCount); 4000 } 4001 // Build the implicit_global map used to map an elements 4002 // local-implicit position to the global-implicit 4003 // position. Primarily used for sideset elements. 'count' 4004 // Elements starting at 'eb_offset' map to the global implicit 4005 // position of 'offset' 4006 int64_t eb_offset = eb->get_offset(); 4007 for (size_t i = 0; i < count; i++) { 4008 elemGlobalImplicitMap[eb_offset + i] = offset + i + 1; 4009 } 4010 elemGlobalImplicitMapDefined = true; 4011 } 4012 4013 elemMap.set_size(elementCount); 4014 return handle_block_ids(eb, EX_ELEM_MAP, elemMap, ids, num_to_get, offset); 4015 } 4016 4017 int64_t ParallelDatabaseIO::handle_face_ids(const Ioss::FaceBlock *eb, void *ids, 4018 size_t num_to_get) const 4019 { 4020 faceMap.set_size(faceCount); 4021 return handle_block_ids(eb, EX_FACE_MAP, faceMap, ids, num_to_get, 0); 4022 } 4023 4024 int64_t ParallelDatabaseIO::handle_edge_ids(const Ioss::EdgeBlock *eb, void *ids, 4025 size_t num_to_get) const 4026 { 4027 edgeMap.set_size(edgeCount); 4028 return handle_block_ids(eb, EX_EDGE_MAP, edgeMap, ids, num_to_get, 0); 4029 } 4030 4031 void ParallelDatabaseIO::write_nodal_transient_field(ex_entity_type /* type */, 4032 const Ioss::Field & field, 4033 const Ioss::NodeBlock *nb, int64_t count, 4034 void *variables) const 4035 { 4036 Ioss::Field::BasicType ioss_type = field.get_type(); 4037 assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER || 4038 ioss_type == Ioss::Field::INT64 || ioss_type == Ioss::Field::COMPLEX); 4039 4040 // Note that if the field's basic type is COMPLEX, then each component of 4041 // the VariableType is a complex variable consisting of a real and 4042 // imaginary part. Since exodus cannot handle complex variables, 4043 // we have to output a (real and imaginary) X (number of 4044 // components) fields. For example, if V is a 3d vector of complex 4045 // data, the data in the 'variables' array are v_x, v.im_x, v_y, 4046 // v.im_y, v_z, v.im_z which need to be output in six separate 4047 // exodus fields. These fields were already defined in 4048 // "write_results_metadata". 4049 4050 const Ioss::VariableType *var_type = field.transformed_storage(); 4051 std::vector<double> temp(count); 4052 4053 int step = get_current_state(); 4054 step = get_database_step(step); 4055 4056 // get number of components, cycle through each component 4057 // and add suffix to base 'field_name'. Look up index 4058 // of this name in 'm_variables[EX_NODE_BLOCK]' map 4059 int comp_count = var_type->component_count(); 4060 4061 int re_im = 1; 4062 if (ioss_type == Ioss::Field::COMPLEX) { 4063 re_im = 2; 4064 } 4065 for (int complex_comp = 0; complex_comp < re_im; complex_comp++) { 4066 std::string field_name = field.get_name(); 4067 if (re_im == 2) { 4068 field_name += complex_suffix[complex_comp]; 4069 } 4070 4071 char field_suffix_separator = get_field_separator(); 4072 for (int i = 0; i < comp_count; i++) { 4073 std::string var_name = var_type->label_name(field_name, i + 1, field_suffix_separator); 4074 4075 auto var_iter = m_variables[EX_NODE_BLOCK].find(var_name); 4076 if (var_iter == m_variables[EX_NODE_BLOCK].end()) { 4077 std::ostringstream errmsg; 4078 fmt::print(errmsg, "ERROR: Could not find nodal variable '{}'\n", var_name); 4079 IOSS_ERROR(errmsg); 4080 } 4081 4082 int var_index = var_iter->second; 4083 4084 size_t begin_offset = (re_im * i) + complex_comp; 4085 size_t stride = re_im * comp_count; 4086 size_t num_out = 0; 4087 4088 if (ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::COMPLEX) { 4089 num_out = nodeMap.map_field_to_db_scalar_order(static_cast<double *>(variables), temp, 4090 begin_offset, count, stride, 0); 4091 } 4092 else if (ioss_type == Ioss::Field::INTEGER) { 4093 num_out = nodeMap.map_field_to_db_scalar_order(static_cast<int *>(variables), temp, 4094 begin_offset, count, stride, 0); 4095 } 4096 else if (ioss_type == Ioss::Field::INT64) { 4097 num_out = nodeMap.map_field_to_db_scalar_order(static_cast<int64_t *>(variables), temp, 4098 begin_offset, count, stride, 0); 4099 } 4100 4101 if (num_out != static_cast<size_t>(nodeCount)) { 4102 std::ostringstream errmsg; 4103 fmt::print(errmsg, 4104 "ERROR: Problem outputting nodal variable '{}' with index = {} to file '{}' on " 4105 "processor {}\n" 4106 "\tShould have output {:L} values, but instead only output {:L} values.\n", 4107 var_name, var_index, get_filename(), myProcessor, nodeCount, num_out); 4108 IOSS_ERROR(errmsg); 4109 } 4110 4111 // Write the variable... 4112 size_t proc_offset = nb->get_optional_property("_processor_offset", 0); 4113 size_t file_count = nb->get_optional_property("locally_owned_count", num_out); 4114 4115 check_node_owning_processor_data(nodeOwningProcessor, file_count); 4116 filter_owned_nodes(nodeOwningProcessor, myProcessor, temp.data()); 4117 int ierr = ex_put_partial_var(get_file_pointer(), step, EX_NODE_BLOCK, var_index, 0, 4118 proc_offset + 1, file_count, temp.data()); 4119 if (ierr < 0) { 4120 std::ostringstream errmsg; 4121 fmt::print(errmsg, 4122 "Problem outputting nodal variable '{}' with index = {} on processor {}\n", 4123 var_name, var_index, myProcessor); 4124 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__, errmsg.str()); 4125 } 4126 } 4127 } 4128 } 4129 4130 void ParallelDatabaseIO::write_entity_transient_field(ex_entity_type type, const Ioss::Field &field, 4131 const Ioss::GroupingEntity *ge, int64_t count, 4132 void *variables) const 4133 { 4134 static Ioss::Map non_element_map; // Used as an empty map for ge->type() != element block. 4135 const Ioss::VariableType *var_type = field.transformed_storage(); 4136 std::vector<double> temp(count); 4137 4138 int step = get_current_state(); 4139 step = get_database_step(step); 4140 4141 Ioss::Map *map = nullptr; 4142 ssize_t eb_offset = 0; 4143 if (ge->type() == Ioss::ELEMENTBLOCK) { 4144 const Ioss::ElementBlock *elb = dynamic_cast<const Ioss::ElementBlock *>(ge); 4145 Ioss::Utils::check_dynamic_cast(elb); 4146 eb_offset = elb->get_offset(); 4147 map = &elemMap; 4148 } 4149 else { 4150 map = &non_element_map; 4151 } 4152 4153 Ioss::Field::BasicType ioss_type = field.get_type(); 4154 assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER || 4155 ioss_type == Ioss::Field::INT64 || ioss_type == Ioss::Field::COMPLEX); 4156 4157 // Note that if the field's basic type is COMPLEX, then each component of 4158 // the VariableType is a complex variable consisting of a real and 4159 // imaginary part. Since exodus cannot handle complex variables, 4160 // we have to output a (real and imaginary) X (number of 4161 // components) fields. For example, if V is a 3d vector of complex 4162 // data, the data in the 'variables' array are v_x, v.im_x, v_y, 4163 // v.im_y, v_z, v.im_z which need to be output in six separate 4164 // exodus fields. These fields were already defined in 4165 // "write_results_metadata". 4166 4167 // get number of components, cycle through each component 4168 // and add suffix to base 'field_name'. Look up index 4169 // of this name in 'm_variables[type]' map 4170 int comp_count = var_type->component_count(); 4171 4172 int re_im = 1; 4173 if (ioss_type == Ioss::Field::COMPLEX) { 4174 re_im = 2; 4175 } 4176 for (int complex_comp = 0; complex_comp < re_im; complex_comp++) { 4177 std::string field_name = field.get_name(); 4178 if (re_im == 2) { 4179 field_name += complex_suffix[complex_comp]; 4180 } 4181 4182 char field_suffix_separator = get_field_separator(); 4183 for (int i = 0; i < comp_count; i++) { 4184 std::string var_name = var_type->label_name(field_name, i + 1, field_suffix_separator); 4185 4186 auto var_iter = m_variables[type].find(var_name); 4187 if (var_iter == m_variables[type].end()) { 4188 std::ostringstream errmsg; 4189 fmt::print(errmsg, "ERROR: Could not find field '{}'\n", var_name); 4190 IOSS_ERROR(errmsg); 4191 } 4192 int var_index = var_iter->second; 4193 assert(var_index > 0); 4194 4195 // var is a [count,comp,re_im] array; re_im = 1(real) or 2(complex) 4196 // beg_offset = (re_im*i)+complex_comp 4197 // number_values = count 4198 // stride = re_im*comp_count 4199 ssize_t begin_offset = (re_im * i) + complex_comp; 4200 ssize_t stride = re_im * comp_count; 4201 4202 if (ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::COMPLEX) { 4203 map->map_field_to_db_scalar_order(static_cast<double *>(variables), temp, begin_offset, 4204 count, stride, eb_offset); 4205 } 4206 else if (ioss_type == Ioss::Field::INTEGER) { 4207 map->map_field_to_db_scalar_order(static_cast<int *>(variables), temp, begin_offset, count, 4208 stride, eb_offset); 4209 } 4210 else if (ioss_type == Ioss::Field::INT64) { 4211 map->map_field_to_db_scalar_order(static_cast<int64_t *>(variables), temp, begin_offset, 4212 count, stride, eb_offset); 4213 } 4214 4215 // Write the variable... 4216 size_t proc_offset = ge->get_optional_property("_processor_offset", 0); 4217 size_t file_count = ge->get_optional_property("locally_owned_count", count); 4218 4219 int64_t id = Ioex::get_id(ge, type, &ids_); 4220 int ierr; 4221 if (type == EX_SIDE_SET) { 4222 size_t offset = ge->get_property("set_offset").get_int(); 4223 ierr = ex_put_partial_var(get_file_pointer(), step, type, var_index, id, 4224 proc_offset + offset + 1, count, temp.data()); 4225 } 4226 else if (type == EX_NODE_SET) { 4227 std::vector<double> file_data; 4228 file_data.reserve(file_count); 4229 map_nodeset_data(nodesetOwnedNodes[ge], temp.data(), file_data); 4230 ierr = ex_put_partial_var(get_file_pointer(), step, type, var_index, id, proc_offset + 1, 4231 file_count, file_data.data()); 4232 } 4233 else { 4234 ierr = ex_put_partial_var(get_file_pointer(), step, type, var_index, id, proc_offset + 1, 4235 file_count, temp.data()); 4236 } 4237 4238 if (ierr < 0) { 4239 std::ostringstream extra_info; 4240 fmt::print(extra_info, "Outputting component {} of field '{}' at step {:L} on {} '{}'.", i, 4241 field_name, step, ge->type_string(), ge->name()); 4242 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__, extra_info.str()); 4243 } 4244 } 4245 } 4246 } 4247 4248 int64_t ParallelDatabaseIO::put_Xset_field_internal(ex_entity_type type, const Ioss::EntitySet *ns, 4249 const Ioss::Field &field, void *data, 4250 size_t data_size) const 4251 { 4252 size_t entity_count = ns->entity_count(); 4253 size_t num_to_get = field.verify(data_size); 4254 4255 int64_t id = Ioex::get_id(ns, type, &ids_); 4256 Ioss::Field::RoleType role = field.get_role(); 4257 4258 if (role == Ioss::Field::MESH) { 4259 4260 void * out_data = data; 4261 std::vector<int> i32data; 4262 std::vector<int64_t> i64data; 4263 std::vector<double> dbldata; 4264 4265 size_t proc_offset = ns->get_optional_property("_processor_offset", 0); 4266 size_t file_count = ns->get_optional_property("locally_owned_count", num_to_get); 4267 4268 if (field.get_name() == "ids" || field.get_name() == "ids_raw") { 4269 // Map node id from global node id to local node id. 4270 // Do it in 'data' ... 4271 4272 if (field.get_name() == "ids") { 4273 nodeMap.reverse_map_data(data, field, num_to_get); 4274 } 4275 4276 if (type == EX_NODE_SET) { 4277 nodesetOwnedNodes[ns].reserve(file_count); 4278 if (int_byte_size_api() == 4) { 4279 i32data.reserve(file_count); 4280 check_node_owning_processor_data(nodeOwningProcessor, file_count); 4281 map_nodeset_id_data(nodeOwningProcessor, nodesetOwnedNodes[ns], myProcessor, 4282 reinterpret_cast<int *>(data), num_to_get, i32data); 4283 assert(i32data.size() == file_count); 4284 // Maps local to "global_implicit" 4285 map_local_to_global_implicit(i32data.data(), file_count, nodeGlobalImplicitMap); 4286 out_data = i32data.data(); 4287 } 4288 else { 4289 i64data.reserve(file_count); 4290 check_node_owning_processor_data(nodeOwningProcessor, file_count); 4291 map_nodeset_id_data(nodeOwningProcessor, nodesetOwnedNodes[ns], myProcessor, 4292 reinterpret_cast<int64_t *>(data), num_to_get, i64data); 4293 assert(i64data.size() == file_count); 4294 map_local_to_global_implicit(i64data.data(), file_count, nodeGlobalImplicitMap); 4295 out_data = i64data.data(); 4296 } 4297 } 4298 int ierr = ex_put_partial_set(get_file_pointer(), type, id, proc_offset + 1, file_count, 4299 out_data, nullptr); 4300 if (ierr < 0) { 4301 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4302 } 4303 } 4304 else if (field.get_name() == "orientation") { 4305 int ierr = ex_put_partial_set(get_file_pointer(), type, id, proc_offset + 1, file_count, 4306 nullptr, out_data); 4307 if (ierr < 0) { 4308 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4309 } 4310 } 4311 else if (field.get_name() == "distribution_factors") { 4312 int ierr = 0; 4313 if (type == EX_NODE_SET) { 4314 map_nodeset_data(nodesetOwnedNodes[ns], reinterpret_cast<double *>(data), dbldata); 4315 ierr = ex_put_partial_set_dist_fact(get_file_pointer(), type, id, proc_offset + 1, 4316 file_count, dbldata.data()); 4317 } 4318 else { 4319 ierr = ex_put_partial_set_dist_fact(get_file_pointer(), type, id, proc_offset + 1, 4320 num_to_get, static_cast<double *>(out_data)); 4321 } 4322 if (ierr < 0) { 4323 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4324 } 4325 } 4326 else { 4327 num_to_get = Ioss::Utils::field_warning(ns, field, "output"); 4328 } 4329 } 4330 else if (role == Ioss::Field::TRANSIENT) { 4331 // Check if the specified field exists on this element block. 4332 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 4333 // exist on the database as scalars with the appropriate 4334 // extensions. 4335 4336 // Transfer each component of the variable into 'data' and then 4337 // output. Need temporary storage area of size 'number of 4338 // elements in this block. 4339 write_entity_transient_field(type, field, ns, entity_count, data); 4340 } 4341 else if (role == Ioss::Field::ATTRIBUTE) { 4342 num_to_get = write_attribute_field(type, field, ns, data); 4343 } 4344 else if (role == Ioss::Field::REDUCTION) { 4345 store_reduction_field(type, field, ns, data); 4346 } 4347 return num_to_get; 4348 } 4349 4350 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::NodeSet *ns, const Ioss::Field &field, 4351 void *data, size_t data_size) const 4352 { 4353 return put_Xset_field_internal(EX_NODE_SET, ns, field, data, data_size); 4354 } 4355 4356 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::EdgeSet *ns, const Ioss::Field &field, 4357 void *data, size_t data_size) const 4358 { 4359 return put_Xset_field_internal(EX_EDGE_SET, ns, field, data, data_size); 4360 } 4361 4362 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::FaceSet *ns, const Ioss::Field &field, 4363 void *data, size_t data_size) const 4364 { 4365 return put_Xset_field_internal(EX_FACE_SET, ns, field, data, data_size); 4366 } 4367 4368 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::ElementSet *ns, const Ioss::Field &field, 4369 void *data, size_t data_size) const 4370 { 4371 return put_Xset_field_internal(EX_ELEM_SET, ns, field, data, data_size); 4372 } 4373 4374 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::SideSet *fs, const Ioss::Field &field, 4375 void * /* data */, size_t data_size) const 4376 { 4377 size_t num_to_get = field.verify(data_size); 4378 if (field.get_name() == "ids") { 4379 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 4380 } 4381 else { 4382 num_to_get = Ioss::Utils::field_warning(fs, field, "output"); 4383 } 4384 return num_to_get; 4385 } 4386 4387 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::CommSet * /*cs*/, 4388 const Ioss::Field &field, void * /*data*/, 4389 size_t data_size) const 4390 { 4391 size_t num_to_get = field.verify(data_size); 4392 return num_to_get; 4393 } 4394 4395 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::SideBlock *fb, const Ioss::Field &field, 4396 void *data, size_t data_size) const 4397 { 4398 size_t num_to_get = field.verify(data_size); 4399 int64_t id = Ioex::get_id(fb, EX_SIDE_SET, &ids_); 4400 4401 size_t entity_count = fb->entity_count(); 4402 size_t offset = fb->get_property("set_offset").get_int(); 4403 4404 Ioss::Field::RoleType role = field.get_role(); 4405 4406 if (role == Ioss::Field::MESH) { 4407 if (field.get_name() == "side_ids" && fb->name() == "universal_sideset") { 4408 // The side ids are being stored as the distribution factor 4409 // field on the universal sideset. There should be no other 4410 // side sets that request this field... (Eventually, 4411 // create an id field to store this info. 4412 4413 // Need to convert 'ints' to 'double' for storage on mesh... 4414 // FIX 64 4415 if (field.get_type() == Ioss::Field::INTEGER) { 4416 int * ids = static_cast<int *>(data); 4417 std::vector<double> real_ids(num_to_get); 4418 for (size_t i = 0; i < num_to_get; i++) { 4419 real_ids[i] = static_cast<double>(ids[i]); 4420 } 4421 int ierr = ex_put_partial_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, offset + 1, 4422 entity_count, real_ids.data()); 4423 if (ierr < 0) { 4424 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4425 } 4426 } 4427 else { 4428 int64_t * ids = static_cast<int64_t *>(data); 4429 std::vector<double> real_ids(num_to_get); 4430 for (size_t i = 0; i < num_to_get; i++) { 4431 real_ids[i] = static_cast<double>(ids[i]); 4432 } 4433 int ierr = ex_put_partial_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, offset + 1, 4434 entity_count, real_ids.data()); 4435 if (ierr < 0) { 4436 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4437 } 4438 } 4439 } 4440 4441 else if (field.get_name() == "side_ids") { 4442 } 4443 4444 else if (field.get_name() == "ids") { 4445 // ============================================================= 4446 // NOTE: Code is currently commented out since we have 4447 // redundant ways of getting the data (element/side) out to 4448 // the database. The 'ids' field method relies on a numbering 4449 // kluge, so for now trying the 'element_side' field... 4450 // ============================================================= 4451 } 4452 4453 else if (field.get_name() == "distribution_factors") { 4454 int ierr; 4455 size_t df_offset = fb->get_property("set_df_offset").get_int(); 4456 size_t proc_df_offset = fb->get_property("processor_df_offset").get_int(); 4457 size_t df_count = fb->get_property("distribution_factor_count").get_int(); 4458 ierr = ex_put_partial_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, 4459 proc_df_offset + df_offset + 1, df_count, 4460 static_cast<double *>(data)); 4461 if (ierr < 0) { 4462 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4463 } 4464 } 4465 else if (field.get_name() == "element_side") { 4466 // In exodusII, the 'side block' is stored as a sideset. A 4467 // sideset has a list of elements and a corresponding local 4468 // element side (1-based) 4469 4470 // The 'data' passed into the function is stored as a 4471 // 2D vector e0,f0,e1,f1,... (e=element, f=side) 4472 4473 // To avoid overwriting the passed in data, we allocate 4474 // two arrays to store the data for this sideset. 4475 4476 // The element_id passed in is the global id; we need to 4477 // output the local id. 4478 4479 // Allocate space for local side number and element numbers 4480 // numbers. 4481 // See if edges or faces... 4482 size_t side_offset = Ioss::Utils::get_side_offset(fb); 4483 4484 size_t index = 0; 4485 4486 size_t proc_offset = fb->get_optional_property("_processor_offset", 0); 4487 4488 if (field.get_type() == Ioss::Field::INTEGER) { 4489 Ioss::IntVector element(num_to_get); 4490 Ioss::IntVector side(num_to_get); 4491 int * el_side = reinterpret_cast<int *>(data); 4492 4493 for (size_t i = 0; i < num_to_get; i++) { 4494 element[i] = elemMap.global_to_local(el_side[index++]); 4495 side[i] = el_side[index++] + side_offset; 4496 } 4497 4498 map_local_to_global_implicit(element.data(), num_to_get, elemGlobalImplicitMap); 4499 int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, proc_offset + offset + 1, 4500 num_to_get, element.data(), side.data()); 4501 if (ierr < 0) { 4502 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4503 } 4504 } 4505 else { 4506 Ioss::Int64Vector element(num_to_get); 4507 Ioss::Int64Vector side(num_to_get); 4508 int64_t * el_side = reinterpret_cast<int64_t *>(data); 4509 4510 for (size_t i = 0; i < num_to_get; i++) { 4511 element[i] = elemMap.global_to_local(el_side[index++]); 4512 side[i] = el_side[index++] + side_offset; 4513 } 4514 4515 map_local_to_global_implicit(element.data(), num_to_get, elemGlobalImplicitMap); 4516 int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, proc_offset + offset + 1, 4517 num_to_get, element.data(), side.data()); 4518 if (ierr < 0) { 4519 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4520 } 4521 } 4522 } 4523 else if (field.get_name() == "element_side_raw") { 4524 // In exodusII, the 'side block' is stored as a sideset. A 4525 // sideset has a list of elements and a corresponding local 4526 // element side (1-based) 4527 4528 // The 'data' passed into the function is stored as a 4529 // 2D vector e0,f0,e1,f1,... (e=element, f=side) 4530 4531 // To avoid overwriting the passed in data, we allocate 4532 // two arrays to store the data for this sideset. 4533 4534 // The element_id passed in is the local id. 4535 4536 // See if edges or faces... 4537 size_t side_offset = Ioss::Utils::get_side_offset(fb); 4538 4539 size_t index = 0; 4540 if (field.get_type() == Ioss::Field::INTEGER) { 4541 Ioss::IntVector element(num_to_get); 4542 Ioss::IntVector side(num_to_get); 4543 int * el_side = reinterpret_cast<int *>(data); 4544 4545 for (size_t i = 0; i < num_to_get; i++) { 4546 element[i] = el_side[index++]; 4547 side[i] = el_side[index++] + side_offset; 4548 } 4549 4550 int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, offset + 1, entity_count, 4551 element.data(), side.data()); 4552 if (ierr < 0) { 4553 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4554 } 4555 } 4556 else { 4557 Ioss::Int64Vector element(num_to_get); 4558 Ioss::Int64Vector side(num_to_get); 4559 int64_t * el_side = reinterpret_cast<int64_t *>(data); 4560 4561 for (size_t i = 0; i < num_to_get; i++) { 4562 element[i] = el_side[index++]; 4563 side[i] = el_side[index++] + side_offset; 4564 } 4565 4566 int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, offset + 1, entity_count, 4567 element.data(), side.data()); 4568 if (ierr < 0) { 4569 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4570 } 4571 } 4572 } 4573 else if (field.get_name() == "connectivity") { 4574 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 4575 } 4576 else if (field.get_name() == "connectivity_raw") { 4577 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 4578 } 4579 else { 4580 num_to_get = Ioss::Utils::field_warning(fb, field, "output"); 4581 } 4582 } 4583 else if (role == Ioss::Field::TRANSIENT) { 4584 // Check if the specified field exists on this block. 4585 // Note that 'higher-order' storage types (e.g. SYM_TENSOR) 4586 // exist on the database as scalars with the appropriate 4587 // extensions. 4588 4589 // Transfer each component of the variable into 'data' and then 4590 // output. Need temporary storage area of size 'number of 4591 // entities in this block. 4592 write_entity_transient_field(EX_SIDE_SET, field, fb, entity_count, data); 4593 } 4594 else if (role == Ioss::Field::ATTRIBUTE) { 4595 num_to_get = write_attribute_field(EX_SIDE_SET, field, fb, data); 4596 } 4597 else if (role == Ioss::Field::REDUCTION) { 4598 store_reduction_field(EX_SIDE_SET, field, fb, data); 4599 } 4600 return num_to_get; 4601 } 4602 4603 void ParallelDatabaseIO::write_meta_data(Ioss::IfDatabaseExistsBehavior behavior) 4604 { 4605 Ioss::Region *region = get_region(); 4606 common_write_meta_data(behavior); 4607 4608 char the_title[max_line_length + 1]; 4609 4610 // Title... 4611 if (region->property_exists("title")) { 4612 std::string title_str = region->get_property("title").get_string(); 4613 Ioss::Utils::copy_string(the_title, title_str); 4614 } 4615 else { 4616 Ioss::Utils::copy_string(the_title, "IOSS Default Output Title"); 4617 } 4618 4619 bool file_per_processor = false; 4620 Ioex::Mesh mesh(spatialDimension, the_title, util(), file_per_processor); 4621 mesh.populate(region); 4622 4623 if (behavior != Ioss::DB_APPEND && behavior != Ioss::DB_MODIFY) { 4624 if (!properties.exists("OMIT_QA_RECORDS")) { 4625 put_qa(); 4626 } 4627 if (!properties.exists("OMIT_INFO_RECORDS")) { 4628 put_info(); 4629 } 4630 4631 // Write the metadata to the exodusII file... 4632 Ioex::Internals data(get_file_pointer(), maximumNameLength, util()); 4633 mesh.comm.outputNemesis = false; 4634 4635 int ierr = data.write_meta_data(mesh); 4636 4637 if (ierr < 0) { 4638 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4639 } 4640 } 4641 4642 metaDataWritten = true; 4643 4644 // Set the processor offset property. Specifies where in the global list, the data from this 4645 // processor begins... 4646 update_processor_offset_property(region, mesh); 4647 4648 if (behavior != Ioss::DB_APPEND && behavior != Ioss::DB_MODIFY) { 4649 output_node_map(); 4650 output_other_meta_data(); 4651 } 4652 } 4653 4654 void ParallelDatabaseIO::create_implicit_global_map() const 4655 { 4656 // If the node is locally owned, then its position is basically 4657 // determined by removing all shared nodes from the list and 4658 // then compressing the list. This location plus the proc_offset 4659 // gives its location in the global-implicit file. 4660 // 4661 // Do this over in the DecompositionData class since it has 4662 // several utilities in place for MPI communication. 4663 4664 DecompositionData<int64_t> compose(Ioss::PropertyManager(), util().communicator()); 4665 int64_t locally_owned_count = 0; 4666 int64_t processor_offset = 0; 4667 compose.create_implicit_global_map(nodeOwningProcessor, nodeGlobalImplicitMap, nodeMap, 4668 &locally_owned_count, &processor_offset); 4669 4670 nodeGlobalImplicitMapDefined = true; 4671 const Ioss::NodeBlockContainer &node_blocks = get_region()->get_node_blocks(); 4672 if (!node_blocks[0]->property_exists("locally_owned_count")) { 4673 node_blocks[0]->property_add(Ioss::Property("locally_owned_count", locally_owned_count)); 4674 } 4675 if (!node_blocks[0]->property_exists("_processor_offset")) { 4676 node_blocks[0]->property_add(Ioss::Property("_processor_offset", processor_offset)); 4677 } 4678 4679 output_node_map(); 4680 } 4681 4682 void ParallelDatabaseIO::output_node_map() const 4683 { 4684 // Write the partial nodemap to the database... This is called 4685 // two times -- once from create_implicit_global_map() and once 4686 // from write_meta_data(). It will only output the map if 4687 // the metadata has been written to the output database AND if 4688 // the nodeMap.map and nodeGlobalImplicitMap are defined. 4689 4690 if (metaDataWritten) { 4691 const Ioss::NodeBlockContainer &node_blocks = get_region()->get_node_blocks(); 4692 if (node_blocks.empty()) { 4693 return; 4694 } 4695 assert(node_blocks[0]->property_exists("_processor_offset")); 4696 assert(node_blocks[0]->property_exists("locally_owned_count")); 4697 size_t processor_offset = node_blocks[0]->get_property("_processor_offset").get_int(); 4698 size_t locally_owned_count = node_blocks[0]->get_property("locally_owned_count").get_int(); 4699 4700 int ierr = 0; 4701 if (nodeMap.defined() && nodeGlobalImplicitMapDefined) { 4702 4703 if (int_byte_size_api() == 4) { 4704 std::vector<int> file_ids; 4705 file_ids.reserve(locally_owned_count); 4706 check_node_owning_processor_data(nodeOwningProcessor, locally_owned_count); 4707 filter_owned_nodes(nodeOwningProcessor, myProcessor, &nodeMap.map()[1], file_ids); 4708 ierr = ex_put_partial_id_map(get_file_pointer(), EX_NODE_MAP, processor_offset + 1, 4709 locally_owned_count, file_ids.data()); 4710 } 4711 else { 4712 std::vector<int64_t> file_ids; 4713 file_ids.reserve(locally_owned_count); 4714 check_node_owning_processor_data(nodeOwningProcessor, locally_owned_count); 4715 filter_owned_nodes(nodeOwningProcessor, myProcessor, &nodeMap.map()[1], file_ids); 4716 ierr = ex_put_partial_id_map(get_file_pointer(), EX_NODE_MAP, processor_offset + 1, 4717 locally_owned_count, file_ids.data()); 4718 } 4719 } 4720 if (ierr < 0) { 4721 Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__); 4722 } 4723 } 4724 } 4725 4726 void ParallelDatabaseIO::check_valid_values() const 4727 { 4728 std::vector<int64_t> counts{nodeCount, elementCount, m_groupCount[EX_ELEM_BLOCK]}; 4729 std::vector<int64_t> all_counts; 4730 util().all_gather(counts, all_counts); 4731 // Get minimum value in `all_counts`. If >0, then don't need to check further... 4732 auto min_val = *std::min_element(all_counts.begin(), all_counts.end()); 4733 4734 if (myProcessor == 0) { 4735 size_t proc_count = all_counts.size() / 3; 4736 4737 if (min_val < 0) { 4738 static std::array<std::string, 3> label{"node", "element", "element block"}; 4739 // Error on one or more of the counts... 4740 for (size_t j = 0; j < 3; j++) { 4741 std::vector<size_t> bad_proc; 4742 for (size_t i = 0; i < proc_count; i++) { 4743 if (all_counts[3 * i + j] < 0) { 4744 bad_proc.push_back(i); 4745 } 4746 } 4747 4748 if (!bad_proc.empty()) { 4749 std::ostringstream errmsg; 4750 fmt::print(errmsg, "ERROR: Negative {} count on {} processor{}:\n\t{}\n\n", label[j], 4751 bad_proc.size(), bad_proc.size() > 1 ? "s" : "", 4752 Ioss::Utils::format_id_list(bad_proc, ":")); 4753 IOSS_ERROR(errmsg); 4754 } 4755 } 4756 } 4757 4758 // Now check for warning (count == 0) 4759 if (min_val <= 0) { 4760 static std::array<std::string, 3> label{"nodes or elements", "elements", "element blocks"}; 4761 // Possible warning on one or more of the counts... 4762 // Note that it is possible to have nodes on a processor with no elements, 4763 // but not possible to have elements if no nodes... 4764 for (size_t j = 0; j < 3; j++) { 4765 std::vector<size_t> bad_proc; 4766 for (size_t i = 0; i < proc_count; i++) { 4767 if (all_counts[3 * i + j] == 0) { 4768 bad_proc.push_back(i); 4769 } 4770 } 4771 4772 if (!bad_proc.empty()) { 4773 fmt::print(Ioss::WARNING(), "No {} on processor{}:\n\t{}\n\n", label[j], 4774 bad_proc.size() > 1 ? "s" : "", Ioss::Utils::format_id_list(bad_proc, ":")); 4775 if (j == 0) { 4776 break; 4777 } 4778 } 4779 } 4780 } 4781 } 4782 else { // All other processors; need to abort if negative count 4783 if (min_val < 0) { 4784 std::ostringstream errmsg; 4785 IOSS_ERROR(errmsg); 4786 } 4787 } 4788 } 4789 } // namespace Ioex 4790 #endif 4791