// Copyright(C) 1999-2024 National Technology & Engineering Solutions
// of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
// NTESS, the U.S. Government retains certain rights in this software.
//
// See packages/seacas/LICENSE for details

#include "Ioss_CodeTypes.h"
#include "Ioss_FileInfo.h"
#include "Ioss_ParallelUtils.h"
#include "Ioss_SerializeIO.h"
#include "Ioss_SmartAssert.h"
#include "Ioss_SurfaceSplit.h"
#include "Ioss_Utils.h"
#include "exodus/Ioex_DatabaseIO.h"
#include "exodus/Ioex_Internals.h"
#include "exodus/Ioex_Utils.h"
#include <array>
#include <cassert>
#include <cfloat>
#include <cstring>
#include <ctime>
#include <vtk_exodusII.h>
#include "vtk_fmt.h"
#include VTK_FMT(fmt/format.h)
#include VTK_FMT(fmt/ostream.h)
#include <limits>
#include <map>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <tokenize.h>
#include <vector>

#include "Ioex_BaseDatabaseIO.h"
#include "Ioss_Assembly.h"
#include "Ioss_Blob.h"
#include "Ioss_CommSet.h"
#include "Ioss_DBUsage.h"
#include "Ioss_EdgeBlock.h"
#include "Ioss_EdgeSet.h"
#include "Ioss_ElementBlock.h"
#include "Ioss_ElementSet.h"
#include "Ioss_ElementTopology.h"
#include "Ioss_EntityBlock.h"
#include "Ioss_EntitySet.h"
#include "Ioss_EntityType.h"
#include "Ioss_FaceBlock.h"
#include "Ioss_FaceSet.h"
#include "Ioss_Field.h"
#include "Ioss_GroupingEntity.h"
#include "Ioss_Map.h"
#include "Ioss_NodeBlock.h"
#include "Ioss_NodeSet.h"
#include "Ioss_Property.h"
#include "Ioss_PropertyManager.h"
#include "Ioss_Region.h"
#include "Ioss_SideBlock.h"
#include "Ioss_SideSet.h"
#include "Ioss_State.h"
#include "Ioss_VariableType.h"
#include "hopscotch_hash.h"

// ========================================================================
// Static internal helper functions
// ========================================================================
namespace {
  const size_t max_line_length = MAX_LINE_LENGTH;

  std::string SEP() { return std::string("@"); } // Separator for attribute offset storage
  std::array<std::string, 2> complex_suffix{".re", ".im"};

  void get_connectivity_data(int exoid, void *data, ex_entity_type type, ex_entity_id id,
                             int position)
  {
    int ierr = 0;
    if ((ex_int64_status(exoid) & EX_BULK_INT64_API) != 0) {
      std::array<int64_t *, 3> conn{nullptr, nullptr, nullptr};
      conn[position] = static_cast<int64_t *>(data);
      ierr           = ex_get_conn(exoid, type, id, conn[0], conn[1], conn[2]);
    }
    else {
      std::array<int *, 3> conn{nullptr, nullptr, nullptr};
      conn[position] = static_cast<int *>(data);
      ierr           = ex_get_conn(exoid, type, id, conn[0], conn[1], conn[2]);
    }
    if (ierr < 0) {
      Ioex::exodus_error(exoid, __LINE__, __func__, __FILE__);
    }
  }

  template <typename T>
  void compute_internal_border_maps(T *entities, T *internal, size_t count, size_t entity_count)
  {
    // Construct the node/element map (internal vs. border).
    // Border nodes/elements are those in the communication map (use entities array)
    // Internal nodes/elements are the rest.  Allocate array to hold all nodes/elements,
    // initialize all to '1', then zero out the nodes/elements in 'entities'.
    // Iterate through array again and consolidate all '1's
    for (size_t ij = 0; ij < count; ij++) {
      internal[ij] = 1;
    }
    for (size_t J = 0; J < entity_count; J++) {
      internal[entities[J] - 1] = 0;
    }

    size_t b = 0;
    for (size_t ij = 0; ij < count; ij++) {
      if (internal[ij] == 0) {
        entities[b++] = ij + 1;
      }
    }

    size_t k = 0;
    for (size_t ij = 0; ij < count; ij++) {
      if (internal[ij] == 1) {
        internal[k++] = ij + 1;
      }
    }
  }

  template <typename T>
  void extract_data(std::vector<double> &local_data, T *data, size_t num_entity, size_t comp_count,
                    size_t offset)
  {
    local_data.resize(num_entity);
    if (comp_count == 1 && offset == 0) {
      for (size_t j = 0; j < num_entity; j++) {
        local_data[j] = data[j];
      }
    }
    else {
      for (size_t j = 0; j < num_entity; j++) {
        local_data[j] = data[offset];
        offset += comp_count;
      }
    }
  }

} // namespace

namespace Ioex {
  DatabaseIO::DatabaseIO(Ioss::Region *region, const std::string &filename,
                         Ioss::DatabaseUsage db_usage, Ioss_MPI_Comm communicator,
                         const Ioss::PropertyManager &props)
      : Ioex::BaseDatabaseIO(region, filename, db_usage, communicator, props)
  {
    if (!is_input()) {
      // Check whether appending to or modify existing file...
      if (open_create_behavior() == Ioss::DB_APPEND ||
          open_create_behavior() == Ioss::DB_APPEND_GROUP ||
          open_create_behavior() == Ioss::DB_MODIFY) {
        // Append to file if it already exists -- See if the file exists.
        auto file  = Ioss::FileInfo(decoded_filename());
        fileExists = file.exists();
      }
    }

    if (properties.exists("processor_count") && properties.exists("my_processor")) {
      if (!isParallel) {
        isSerialParallel = true;
      }
      else {
        std::ostringstream errmsg;
        fmt::print(
            errmsg,
            "ERROR: Processor id and processor count are specified via the "
            "'processor_count' and 'processor_id' properties which indicates that this "
            "database is "
            "being run in 'serial-parallel' mode, but the database constructor was passed an "
            "mpi communicator which has more than 1 processor. This is not allowed.\n");
        IOSS_ERROR(errmsg);
      }
    }
  }

  bool DatabaseIO::check_valid_file_ptr(bool write_message, std::string *error_msg, int *bad_count,
                                        bool abort_if_error) const
  {
    bool no_collective_calls = Ioss::SerializeIO::isEnabled();
    if (isParallel && no_collective_calls) {
      // Can't output a nice error message on processor 0 and throw a consistent error.
      // Have to just write message on processors that have issue and throw exception.
      if (m_exodusFilePtr < 0) {
        std::ostringstream errmsg;
        std::string        open_create = is_input() ? "open input" : "create output";
        fmt::print(errmsg, "ERROR: Unable to {} exodus decomposed database file '{}'\n",
                   open_create, decoded_filename());

        if (abort_if_error) {
          IOSS_ERROR(errmsg);
        }
        Ioss::WarnOut() << errmsg.str();
        return false;
      }
      return true; // At least on this processor...
    }

    // Check for valid exodus_file_ptr (valid >= 0; invalid < 0)
    int global_file_ptr = m_exodusFilePtr;
    if (isParallel) {
      global_file_ptr = util().global_minmax(m_exodusFilePtr, Ioss::ParallelUtils::DO_MIN);
    }

    if (global_file_ptr < 0) {
      if (write_message || error_msg != nullptr || bad_count != nullptr) {
        Ioss::IntVector status;
        if (isParallel) {
          util().all_gather(m_exodusFilePtr, status);
        }
        else {
          status.push_back(m_exodusFilePtr);
        }

        std::string open_create = is_input() ? "open input" : "create output";
        if (write_message || error_msg != nullptr) {
          // See which processors could not open/create the file...
          std::ostringstream errmsg;
          if (isParallel) {
            fmt::print(errmsg, "ERROR: Unable to {} exodus decomposed database files:\n",
                       open_create);
            for (int i = 0; i < util().parallel_size(); i++) {
              if (status[i] < 0) {
                fmt::print(errmsg, "\t{}\n",
                           Ioss::Utils::decode_filename(get_filename(), i, util().parallel_size()));
              }
            }
          }
          else {
            fmt::print(errmsg, "ERROR: Unable to {} database '{}' of type 'exodusII'", open_create,
                       get_filename());
          }
          fmt::print(errmsg, "\n");
          if (error_msg != nullptr) {
            *error_msg = errmsg.str();
          }
          if (write_message && myProcessor == 0) {
            Ioss::WarnOut() << errmsg.str();
          }
        }
        if (bad_count != nullptr) {
          *bad_count = std::count_if(status.begin(), status.end(), [](int i) { return i < 0; });
        }
        if (abort_if_error) {
          std::ostringstream errmsg;
          fmt::print(errmsg, "ERROR: Cannot {} file '{}'\n", open_create, get_filename());
          IOSS_ERROR(errmsg);
        }
      }
      return false;
    }
    return true;
  }

  bool DatabaseIO::open_input_file(bool write_message, std::string *error_msg, int *bad_count,
                                   bool abort_if_error) const
  {
    int   cpu_word_size = sizeof(double);
    int   io_word_size  = 0;
    float version;

    int mode = exodusMode;
    if (int_byte_size_api() == 8) {
      mode |= EX_ALL_INT64_API;
    }

#if defined EX_DISKLESS
    // Experimental -- in memory read by netcdf library
    if (properties.exists("MEMORY_READ")) {
      mode |= EX_DISKLESS;
    }
#endif

    bool do_timer = false;
    Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer);
    double t_begin = ((do_timer && isParallel) ? Ioss::Utils::timer() : 0);

    int app_opt_val = ex_opts(EX_VERBOSE);
    m_exodusFilePtr = ex_open(decoded_filename().c_str(), EX_READ | mode, &cpu_word_size,
                              &io_word_size, &version);

    if (do_timer && isParallel) {
      double t_end    = Ioss::Utils::timer();
      double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX);
      if (myProcessor == 0) {
        fmt::print(Ioss::DebugOut(), "Input File Open Time = {} ({})\n", duration,
                   decoded_filename());
      }
    }

    bool is_ok = check_valid_file_ptr(write_message, error_msg, bad_count, abort_if_error);

    if (is_ok) {
      finalize_file_open();
    }
    ex_opts(app_opt_val); // Reset back to what it was.
    return is_ok;
  }

  bool DatabaseIO::handle_output_file(bool write_message, std::string *error_msg, int *bad_count,
                                      bool overwrite, bool abort_if_error) const
  {
    // If 'overwrite' is false, we do not want to overwrite or clobber
    // the output file if it already exists since the app might be
    // reading the restart data from this file and then later
    // clobbering it and then writing restart data to the same
    // file. So, for output, we first check whether the file exists
    // and if it it and is writable, assume that we can later create a
    // new or append to existing file.

    // if 'overwrite' is true, then clobber/append
    bool is_ok = false;

    if (!overwrite) {
      // check if file exists and is writeable. If so, return true.
      Ioss::FileInfo file(decoded_filename());
      int            int_is_ok = file.exists() && file.is_writable() ? 1 : 0;

      // Check for consistency among all processors.
      // OK if *all* 0 or *all* 1
      int sum = util().global_minmax(int_is_ok, Ioss::ParallelUtils::DO_SUM);
      if (sum == util().parallel_size()) {
        // Note that at this point, we cannot totally guarantee that
        // we will be able to create the file when needed, but we have
        // a pretty good chance.  We can't guarantee creation without
        // creating and the app (or calling function) doesn't want us to overwrite...
        return true;
      }
      // File doesn't exist on any or all processors, so fall through and try to
      // create file since we won't be overwriting anything...
    }

    int   cpu_word_size = sizeof(double);
    int   io_word_size  = 0;
    float version;

    int mode = exodusMode;
    if (int_byte_size_api() == 8) {
      mode |= EX_ALL_INT64_API;
    }

#if defined EX_DISKLESS
    // Experimental -- in memory write by netcdf library
    if (properties.exists("MEMORY_WRITE")) {
      mode |= EX_DISKLESS;
    }
#endif
    int app_opt_val = ex_opts(EX_VERBOSE);
    if (fileExists) {
      m_exodusFilePtr = ex_open(decoded_filename().c_str(), EX_WRITE | mode, &cpu_word_size,
                                &io_word_size, &version);
    }
    else {
      // If the first write for this file, create it...
      if (int_byte_size_api() == 8) {
        // Check whether client actually wants 4-byte output on db
        // - If they specified INTEGER_SIZE_DB and the size isn't 8,
        //   then don't change mode and use the default 4-byte output.
        if (properties.exists("INTEGER_SIZE_DB")) {
          if (properties.get("INTEGER_SIZE_DB").get_int() == 8) {
            mode |= EX_ALL_INT64_DB;
          }
        }
        else {
          mode |= EX_ALL_INT64_DB;
        }
      }
      m_exodusFilePtr =
          ex_create(decoded_filename().c_str(), mode, &cpu_word_size, &dbRealWordSize);
    }

    is_ok = check_valid_file_ptr(write_message, error_msg, bad_count, abort_if_error);

    if (is_ok) {
      ex_set_max_name_length(m_exodusFilePtr, maximumNameLength);

      // Check properties handled post-create/open...
      if (properties.exists("COMPRESSION_METHOD")) {
        auto method                    = properties.get("COMPRESSION_METHOD").get_string();
        method                         = Ioss::Utils::lowercase(method);
        ex_compression_type exo_method = EX_COMPRESS_ZLIB;
        if (method == "zlib" || method == "libz" || method == "gzip") {
          exo_method = EX_COMPRESS_ZLIB;
        }
        else if (method == "szip") {
#if !defined(NC_HAS_SZIP_WRITE)
#define NC_HAS_SZIP_WRITE 0
#endif
#if NC_HAS_SZIP_WRITE
          exo_method = EX_COMPRESS_SZIP;
#else
          fmt::print(Ioss::WarnOut(), "The NetCDF library does not have SZip compression enabled."
                                      " 'zlib' will be used instead.\n\n");
#endif
        }
        else {
          fmt::print(Ioss::WarnOut(),
                     "Unrecognized compression method specified: '{}'."
                     " 'zlib' will be used instead.\n\n",
                     method);
        }
        ex_set_option(m_exodusFilePtr, EX_OPT_COMPRESSION_TYPE, exo_method);
      }

      if (properties.exists("COMPRESSION_LEVEL")) {
        int comp_level = properties.get("COMPRESSION_LEVEL").get_int();
        ex_set_option(m_exodusFilePtr, EX_OPT_COMPRESSION_LEVEL, comp_level);
      }

      if (properties.exists("COMPRESSION_SHUFFLE")) {
        int shuffle = properties.get("COMPRESSION_SHUFFLE").get_int();
        ex_set_option(m_exodusFilePtr, EX_OPT_COMPRESSION_SHUFFLE, shuffle);
      }
    }
    ex_opts(app_opt_val); // Reset back to what it was.
    return is_ok;
  }

  int DatabaseIO::get_file_pointer() const
  {
    // Returns the file_pointer used to access the file on disk.
    // Checks that the file is open and if not, opens it first.
    if (Ioss::SerializeIO::isEnabled()) {
      if (!Ioss::SerializeIO::inBarrier()) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Process {} is attempting to do I/O without serialized I/O",
                   Ioss::SerializeIO::getRank());
        IOSS_ERROR(errmsg);
      }

      if (!Ioss::SerializeIO::inMyGroup()) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Process {} is attempting to do I/O while {} owns the token",
                   Ioss::SerializeIO::getRank(), Ioss::SerializeIO::getOwner());
        IOSS_ERROR(errmsg);
      }
    }

    return Ioex::BaseDatabaseIO::get_file_pointer();
  }

  void DatabaseIO::read_meta_data_nl()
  {
    // If this is a HISTORY file, there isn't really any metadata
    // Other than a single node and single element.  Just hardwire
    // it here (needed when appending to existing history file)
    if (dbUsage == Ioss::WRITE_HISTORY) {
      if (myProcessor == 0) {
        nodeCount    = 1;
        elementCount = 1;
        auto *nb     = new Ioss::NodeBlock(this, "nodeblock_1", 1, 3);
        get_region()->add(nb);

        // Element Block
        auto *eb = new Ioss::ElementBlock(this, "e1", "sphere", 1);
        eb->property_add(Ioss::Property("id", 1));
        eb->property_add(Ioss::Property("guid", util().generate_guid(1)));
        get_region()->add(eb);
        get_step_times_nl();
        add_region_fields();
      }
      return;
    }

    // APPENDING:
    // There is an assumption that the writing process (mesh, vars) is
    // the same for the original run that created this database and
    // for this run which is appending to the database so the defining
    // of the output database should be the same except we don't write
    // anything since it is already there.  We do need the number of
    // steps though...
    if (open_create_behavior() == Ioss::DB_APPEND || dbUsage == Ioss::QUERY_TIMESTEPS_ONLY) {
      get_step_times_nl();
      return;
    }

    {
      Ioss::SerializeIO serializeIO_(this);

      if (isParallel) {
        Ioex::check_processor_info(decoded_filename(), get_file_pointer(), util().parallel_size(),
                                   myProcessor);
      }

      read_region();
      read_communication_metadata();

      Ioex::read_exodus_basis(get_file_pointer());
      Ioex::read_exodus_quadrature(get_file_pointer());
    }

    get_step_times_nl();

    get_nodeblocks();
    get_edgeblocks();
    get_faceblocks();
    get_elemblocks();

    check_side_topology();

    get_sidesets();
    get_nodesets();
    get_edgesets();
    get_facesets();
    get_elemsets();

    get_commsets();

    // Add assemblies now that all entities should be defined... consistent across processors
    // (metadata)
    get_assemblies();
    get_blobs();

    handle_groups();

    add_region_fields();

    if (!is_input() && open_create_behavior() == Ioss::DB_APPEND) {
      (void)get_map(EX_NODE_BLOCK);
      (void)get_map(EX_EDGE_BLOCK);
      (void)get_map(EX_FACE_BLOCK);
      (void)get_map(EX_ELEM_BLOCK);
    }
  }

  void DatabaseIO::read_region()
  {
    // Add properties and fields to the 'owning' region.
    // Also defines member variables of this class...
    ex_init_params info{};
    int            error = ex_get_init_ext(get_file_pointer(), &info);
    if (error < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }

    spatialDimension = info.num_dim;
    nodeCount        = info.num_nodes;
    edgeCount        = info.num_edge;
    faceCount        = info.num_face;
    elementCount     = info.num_elem;

    m_groupCount[EX_NODE_BLOCK] = 1;
    m_groupCount[EX_EDGE_BLOCK] = info.num_edge_blk;
    m_groupCount[EX_FACE_BLOCK] = info.num_face_blk;
    m_groupCount[EX_ELEM_BLOCK] = info.num_elem_blk;

    m_groupCount[EX_NODE_SET] = info.num_node_sets;
    m_groupCount[EX_EDGE_SET] = info.num_edge_sets;
    m_groupCount[EX_FACE_SET] = info.num_face_sets;
    m_groupCount[EX_ELEM_SET] = info.num_elem_sets;

    m_groupCount[EX_SIDE_SET] = info.num_side_sets;
    m_groupCount[EX_ASSEMBLY] = info.num_assembly;
    m_groupCount[EX_BLOB]     = info.num_blob;

    if (nodeCount == 0 && info.num_blob == 0) {
      fmt::print(Ioss::WarnOut(), "No nodes were found in the model, file '{}'\n",
                 decoded_filename());
    }
    else if (nodeCount < 0) {
      // NOTE: Code will not continue past this call...
      std::ostringstream errmsg;
      fmt::print(errmsg,
                 "ERROR: Negative node count was found in the model\n"
                 "       File: '{}'.\n",
                 decoded_filename());
      IOSS_ERROR(errmsg);
    }

    if (elementCount == 0 && info.num_blob == 0) {
      fmt::print(Ioss::WarnOut(), "No elements were found in the model, file '{}'\n",
                 decoded_filename());
    }

    if (elementCount < 0) {
      // NOTE: Code will not continue past this call...
      std::ostringstream errmsg;
      fmt::print(errmsg,
                 "ERROR: Negative element count was found in the model\n"
                 "       File: '{}'.\n",
                 decoded_filename());
      IOSS_ERROR(errmsg);
    }

    if (elementCount > 0 && m_groupCount[EX_ELEM_BLOCK] <= 0) {
      // NOTE: Code will not continue past this call...
      std::ostringstream errmsg;
      fmt::print(errmsg,
                 "ERROR: No element blocks were found in the model\n"
                 "       File: '{}'.\n",
                 decoded_filename());
      IOSS_ERROR(errmsg);
    }

    Ioss::Region *this_region = get_region();

    // See if any coordinate frames exist on mesh.  If so, define them on region.
    Ioex::add_coordinate_frames(get_file_pointer(), this_region);

    this_region->property_add(Ioss::Property(std::string("title"), info.title));

    // Get QA records from database and add to qaRecords...
    int num_qa = ex_inquire_int(get_file_pointer(), EX_INQ_QA);
    if (num_qa > 0) {
      struct qa_element
      {
        char *qa_record[1][4];
      };

      std::vector<qa_element> qa(num_qa);
      for (int i = 0; i < num_qa; i++) {
        for (int j = 0; j < 4; j++) {
          qa[i].qa_record[0][j] = new char[MAX_STR_LENGTH + 1];
        }
      }

      ex_get_qa(get_file_pointer(), qa[0].qa_record);
      for (int i = 0; i < num_qa; i++) {
        add_qa_record(qa[i].qa_record[0][0], qa[i].qa_record[0][1], qa[i].qa_record[0][2],
                      qa[i].qa_record[0][3]);
      }
      for (int i = 0; i < num_qa; i++) {
        for (int j = 0; j < 4; j++) {
          delete[] qa[i].qa_record[0][j];
        }
      }
    }

    // Get information records from database and add to informationRecords...
    int num_info = ex_inquire_int(get_file_pointer(), EX_INQ_INFO);
    if (num_info > 0) {
      char **info_rec =
          Ioex::get_name_array(num_info, max_line_length); // 'total_lines' pointers to char buffers
      ex_get_info(get_file_pointer(), info_rec);
      for (int i = 0; i < num_info; i++) {
        add_information_record(info_rec[i]);
      }
      Ioex::delete_name_array(info_rec, num_info);
    }
  }

  void DatabaseIO::get_step_times_nl()
  {
    bool                exists         = false;
    double              last_time      = DBL_MAX;
    int                 timestep_count = 0;
    std::vector<double> tsteps(0);

    if (dbUsage == Ioss::WRITE_HISTORY) {
      if (myProcessor == 0) {
        timestep_count = ex_inquire_int(get_file_pointer(), EX_INQ_TIME);
        if (timestep_count <= 0) {
          return;
        }

        // For an exodus file, timesteps are global and are stored in the region.
        // A history file only stores that last time / step
        // Read the timesteps and add them to the region.
        // Since we can't access the Region's stateCount directly, we just add
        // all of the steps and assume the Region is dealing with them directly...
        tsteps.resize(timestep_count);

        int error = ex_get_all_times(get_file_pointer(), Data(tsteps));
        if (error < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        int max_step = properties.get_optional("APPEND_OUTPUT_AFTER_STEP", timestep_count);
        max_step     = std::min(max_step, timestep_count);

        double max_time =
            properties.get_optional("APPEND_OUTPUT_AFTER_TIME", std::numeric_limits<double>::max());

        Ioss::Region *this_region = get_region();
        for (int i = 0; i < max_step; i++) {
          if (tsteps[i] <= max_time) {
            this_region->add_state_nl(tsteps[i] * timeScaleFactor);
          }
        }
      }
    }
    else {
      {
        Ioss::SerializeIO serializeIO_(this);
        timestep_count = ex_inquire_int(get_file_pointer(), EX_INQ_TIME);
        if (timestep_count <= 0) {
          return;
        }

        // For an exodus file, timesteps are global and are stored in the region.
        // Read the timesteps and add to the region
        tsteps.resize(timestep_count, -std::numeric_limits<double>::max());

        // The `EXODUS_CALL_GET_ALL_TIMES=NO` is typically only used in
        // isSerialParallel mode and the client is responsible for
        // making sure that the step times are handled correctly.  All
        // databases will know about the number of timesteps, but if
        // this is skipped, then the times will all be zero.  Use case
        // is that in isSerialParallel, each call to
        // `ex_get_all_times` for all files is performed sequentially,
        // so if you have hundreds to thousands of files, the time for
        // the call is additive and since timesteps are record
        // variables in netCDF, accessing the data for all timesteps
        // involves lseeks throughout the file.
        bool call_ex_get_all_times = true;
        Ioss::Utils::check_set_bool_property(properties, "EXODUS_CALL_GET_ALL_TIMES",
                                             call_ex_get_all_times);
        if (call_ex_get_all_times) {
          int error = ex_get_all_times(get_file_pointer(), Data(tsteps));
          if (error < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }

        // See if the "last_written_time" attribute exists and if it
        // does, check that it matches the largest time in 'tsteps'.
        exists = Ioex::read_last_time_attribute(get_file_pointer(), &last_time);
      }
      if (exists && isParallel) {
        // Assume that if it exists on 1 processor, it exists on
        // all... Sync value among processors since could have a
        // corrupt step on only a single database.
        last_time = util().global_minmax(last_time, Ioss::ParallelUtils::DO_MIN);
      }

      // Only add states that are less than or equal to the
      // 'last_time' value which is either DBL_MAX or the value of
      // the last time successfully written to the database and
      // flushed to disk.  This is used to avoid corrupt data arising
      // from a job that crashed during the writing of the last step
      // on the database.  Output a warning message if there is
      // potentially corrupt data on the database...

      // Check whether user or application wants to limit the times even further...
      // One use case is that job is restarting at a time prior to what has been
      // written to the results file, so want to start appending after
      // restart time instead of at end time on database.
      int max_step = properties.get_optional("APPEND_OUTPUT_AFTER_STEP", timestep_count);
      max_step     = std::min(max_step, timestep_count);

      double max_time =
          properties.get_optional("APPEND_OUTPUT_AFTER_TIME", std::numeric_limits<double>::max());
      last_time = std::min(last_time, max_time);

      Ioss::Region *this_region = get_region();
      for (int i = 0; i < max_step; i++) {
        if (tsteps[i] <= last_time) {
          this_region->add_state_nl(tsteps[i] * timeScaleFactor);
        }
        else {
          if (myProcessor == 0 && max_time == std::numeric_limits<double>::max()) {
            // NOTE: Don't want to warn on all processors if there are
            // corrupt steps on all databases, but this will only print
            // a warning if there is a corrupt step on processor
            // 0... Need better warnings which won't overload in the
            // worst case...
            fmt::print(Ioss::WarnOut(),
                       "Skipping step {} at time {} in database file\n\t{}.\n"
                       "\tThe data for that step is possibly corrupt since the last time written "
                       "successfully was {}.\n",
                       fmt::group_digits(i + 1), tsteps[i], get_filename(), last_time);
          }
        }
      }
    }
  }

  void DatabaseIO::read_communication_metadata()
  {
    // Check that file is nemesis.
    int  num_proc;         // Number of processors file was decomposed for
    int  num_proc_in_file; // Number of processors this file has info for
    char file_type[2];     // "s" for scalar, "p" for parallel

    // Get global data (over all processors)
    int64_t global_nodes    = nodeCount;
    int64_t global_elements = elementCount;
    int64_t global_eblocks  = 0; // unused
    int64_t global_nsets    = 0; // unused
    int64_t global_ssets    = 0; // unused

    int64_t num_external_nodes; // unused
    int64_t num_elem_cmaps     = 0;
    int64_t num_node_cmaps     = 0;
    int64_t num_internal_nodes = nodeCount;
    int64_t num_border_nodes   = 0;
    int64_t num_internal_elems = elementCount;
    int64_t num_border_elems   = 0;

    bool nemesis_file = true;
    int  error = ex_get_init_info(get_file_pointer(), &num_proc, &num_proc_in_file, &file_type[0]);
    if (error < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }

    if (num_proc == 1 && num_proc_in_file == 1) {
      // Not a nemesis file
      nemesis_file = false;
      if (isParallel && util().parallel_size() > 1) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Exodus file '{}' does not contain nemesis information.\n",
                   get_filename());
        IOSS_ERROR(errmsg);
      }
      file_type[0] = 'p';
    }
    else {
      if (!isParallel) {
        // The file contains nemesis parallel information.
        // Even though we are running in serial, make the information
        // available to the application.
        isSerialParallel = true;
        get_region()->property_add(Ioss::Property("processor_count", num_proc));
      }
    }

    if (isParallel && num_proc != util().parallel_size() && util().parallel_size() > 1) {
      std::ostringstream errmsg;
      fmt::print(
          errmsg,
          "ERROR: Exodus file '{}' was decomposed for {} processors; application is currently "
          "being run on {} processors",
          get_filename(), num_proc, util().parallel_size());
      IOSS_ERROR(errmsg);
    }
    if (num_proc_in_file != 1) {
      std::ostringstream errmsg;
      fmt::print(errmsg,
                 "ERROR: Exodus file '{}' contains data for {} processors; application requires 1 "
                 "processor per file.",
                 get_filename(), num_proc_in_file);
      IOSS_ERROR(errmsg);
    }
    if (file_type[0] != 'p') {
      std::ostringstream errmsg;
      fmt::print(
          errmsg,
          "ERROR: Exodus file '{}' contains scalar nemesis data; application requires parallel "
          "nemesis data.",
          get_filename());
      IOSS_ERROR(errmsg);
    }

    bool minimal_nemesis = false;
    if (nemesis_file) {
      // A nemesis file typically separates nodes into multiple
      // communication sets by processor.  (each set specifies
      // nodes/elements that communicate with only a single
      // processor).  For Sierra, we want a single node communication
      // map and a single element communication map specifying all
      // communications so we combine all sets into a single set.

      if (int_byte_size_api() == 4) {
        int gn, ge, geb, gns, gss;
        error           = ex_get_init_global(get_file_pointer(), &gn, &ge, &geb, &gns, &gss);
        global_nodes    = gn;
        global_elements = ge;
        global_eblocks  = geb;
        global_nsets    = gns;
        global_ssets    = gss;
      }
      else {
        error = ex_get_init_global(get_file_pointer(), &global_nodes, &global_elements,
                                   &global_eblocks, &global_nsets, &global_ssets);
      }
      if (error < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }

      // If the file was written with minimal nemesis data, then these will not exist...
      if (int_byte_size_api() == 4) {
        int nin, nbn, nen, nie, nbe, nnc, nec;
        error = ex_get_loadbal_param(get_file_pointer(), &nin, &nbn, &nen, &nie, &nbe, &nnc, &nec,
                                     myProcessor);
        if (error == 0) {
          num_external_nodes = nen;
          num_elem_cmaps     = nec;
          num_node_cmaps     = nnc;
          num_internal_nodes = nin;
          num_border_nodes   = nbn;
          num_internal_elems = nie;
          num_border_elems   = nbe;
        }
        else {
          minimal_nemesis = true;
        }
      }
      else {
        error = ex_get_loadbal_param(get_file_pointer(), &num_internal_nodes, &num_border_nodes,
                                     &num_external_nodes, &num_internal_elems, &num_border_elems,
                                     &num_node_cmaps, &num_elem_cmaps, myProcessor);
        if (error < 0) {
          minimal_nemesis    = true;
          num_external_nodes = 0;
          num_elem_cmaps     = 0;
          num_node_cmaps     = 0;
          num_internal_nodes = 0;
          num_border_nodes   = 0;
          num_internal_elems = 0;
          num_border_elems   = 0;
        }
      }
    }

    commsetNodeCount = num_node_cmaps;
    commsetElemCount = num_elem_cmaps;

    Ioss::Region *region = get_region();
    region->property_add(Ioss::Property("global_node_count", global_nodes));
    region->property_add(Ioss::Property("global_element_count", global_elements));
    region->property_add(Ioss::Property("global_element_block_count", global_eblocks));
    region->property_add(Ioss::Property("global_node_set_count", global_nsets));
    region->property_add(Ioss::Property("global_side_set_count", global_ssets));

    if (!minimal_nemesis) {
      region->property_add(Ioss::Property("internal_node_count", num_internal_nodes));
      region->property_add(Ioss::Property("border_node_count", num_border_nodes));
      region->property_add(Ioss::Property("internal_element_count", num_internal_elems));
      region->property_add(Ioss::Property("border_element_count", num_border_elems));

      // Possibly, the following 4 fields should be nodesets and element
      // sets instead of fields on the region...
      region->field_add(Ioss::Field("internal_nodes", region->field_int_type(), IOSS_SCALAR(),
                                    Ioss::Field::COMMUNICATION, num_internal_nodes));
      region->field_add(Ioss::Field("border_nodes", region->field_int_type(), IOSS_SCALAR(),
                                    Ioss::Field::COMMUNICATION, num_border_nodes));
      region->field_add(Ioss::Field("internal_elements", region->field_int_type(), IOSS_SCALAR(),
                                    Ioss::Field::COMMUNICATION, num_internal_elems));
      region->field_add(Ioss::Field("border_elements", region->field_int_type(), IOSS_SCALAR(),
                                    Ioss::Field::COMMUNICATION, num_border_elems));

      SMART_ASSERT(nodeCount == num_internal_nodes + num_border_nodes)
      (nodeCount)(num_internal_nodes)(num_border_nodes);
      SMART_ASSERT(elementCount == num_internal_elems + num_border_elems)
      (elementCount)(num_internal_elems)(num_border_elems);
    }
  }

  const Ioss::Map &DatabaseIO::get_map(ex_entity_type type) const
  {
    switch (type) {
    case EX_NODE_BLOCK:
    case EX_NODE_SET: return get_map(nodeMap, nodeCount, EX_NODE_MAP, EX_INQ_NODE_MAP);

    case EX_ELEM_BLOCK:
    case EX_ELEM_SET: return get_map(elemMap, elementCount, EX_ELEM_MAP, EX_INQ_ELEM_MAP);

    case EX_FACE_BLOCK:
    case EX_FACE_SET: return get_map(faceMap, faceCount, EX_FACE_MAP, EX_INQ_FACE_MAP);

    case EX_EDGE_BLOCK:
    case EX_EDGE_SET: return get_map(edgeMap, edgeCount, EX_EDGE_MAP, EX_INQ_EDGE_MAP);

    default:
      std::ostringstream errmsg;
      fmt::print(errmsg, "INTERNAL ERROR: Invalid map type. "
                         "Something is wrong in the Ioex::DatabaseIO::get_map() function. "
                         "Please report.\n");
      IOSS_ERROR(errmsg);
    }
  }

  const Ioss::Map &DatabaseIO::get_map(Ioss::Map &entity_map, int64_t entity_count,
                                       ex_entity_type entity_type, ex_inquiry inquiry_type) const
  {
    // Allocate space for node number map and read it in...
    // Can be called multiple times, allocate 1 time only
    bool read_exodus_map = true;
    if (entity_map.map().empty()) {
      switch (entity_type) {
      case EX_NODE_MAP: read_exodus_map = !properties.exists("IGNORE_NODE_MAP"); break;
      case EX_ELEM_MAP: read_exodus_map = !properties.exists("IGNORE_ELEMENT_MAP"); break;
      case EX_FACE_MAP: read_exodus_map = !properties.exists("IGNORE_FACE_MAP"); break;
      case EX_EDGE_MAP: read_exodus_map = !properties.exists("IGNORE_EDGE_MAP"); break;
      default:
        std::ostringstream errmsg;
        fmt::print(errmsg, "INTERNAL ERROR: Invalid map type. "
                           "Something is wrong in the Ioex::DatabaseIO::get_map() function. "
                           "Please report.\n");
        IOSS_ERROR(errmsg);
      }

      if (is_input() || open_create_behavior() == Ioss::DB_APPEND) {
        entity_map.set_size(entity_count);

        Ioss::SerializeIO serializeIO_(this);
        // Check whether there is a "original_global_id_map" map on
        // the database. If so, use it instead of the "node_num_map".
        bool map_read  = false;
        int  map_count = ex_inquire_int(get_file_pointer(), inquiry_type);
        if (read_exodus_map && map_count > 0) {
          char **names = Ioex::get_name_array(map_count, maximumNameLength);
          int    ierr  = ex_get_names(get_file_pointer(), entity_type, names);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }

          for (int i = 0; i < map_count; i++) {
            if (Ioss::Utils::str_equal(names[i], "original_global_id_map")) {
              int error = 0;
              if ((ex_int64_status(get_file_pointer()) & EX_BULK_INT64_API) != 0) {
                Ioss::Int64Vector tmp_map(entity_map.size());
                error = ex_get_num_map(get_file_pointer(), entity_type, i + 1, Data(tmp_map));
                if (error >= 0) {
                  entity_map.set_map(Data(tmp_map), tmp_map.size(), 0, true);
                  map_read = true;
                  break;
                }
              }
              else {
                // Ioss stores as 64-bit, read as 32-bit and copy over...
                Ioss::IntVector tmp_map(entity_map.size());
                error = ex_get_num_map(get_file_pointer(), entity_type, i + 1, Data(tmp_map));
                if (error >= 0) {
                  entity_map.set_map(Data(tmp_map), tmp_map.size(), 0, true);
                  map_read = true;
                  break;
                }
              }
            }
          }
          Ioex::delete_name_array(names, map_count);
        }

        if (!map_read) {
          if (isParallel || read_exodus_map) {
            int error = 0;
            if ((ex_int64_status(get_file_pointer()) & EX_BULK_INT64_API) != 0) {
              Ioss::Int64Vector tmp_map(entity_map.size());
              error = ex_get_id_map(get_file_pointer(), entity_type, Data(tmp_map));
              if (error >= 0) {
                entity_map.set_map(Data(tmp_map), tmp_map.size(), 0, true);
              }
            }
            else {
              // Ioss stores as 64-bit, read as 32-bit and copy over...
              Ioss::IntVector tmp_map(entity_map.size());
              error = ex_get_id_map(get_file_pointer(), entity_type, Data(tmp_map));
              if (error >= 0) {
                entity_map.set_map(Data(tmp_map), tmp_map.size(), 0, true);
              }
            }
            if (error < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else {
            // Ignoring exodus database map (if any)
            entity_map.set_default(entity_map.size());
          }
        }
      }
      else {
        // Output database; entity_map.map not set yet... Build a default map.
        entity_map.set_default(entity_count);
      }
    }
    return entity_map;
  }

  void DatabaseIO::get_elemblocks() { get_blocks(EX_ELEM_BLOCK, 0, "block"); }

  void DatabaseIO::get_faceblocks() { get_blocks(EX_FACE_BLOCK, 1, "faceblock"); }

  void DatabaseIO::get_edgeblocks() { get_blocks(EX_EDGE_BLOCK, 2, "edgeblock"); }

  void DatabaseIO::get_blocks(ex_entity_type entity_type, int rank_offset,
                              const std::string &basename)
  {
    // Attributes of an X block are:  (X = element, face, or edge)
    // -- id
    // -- name
    // -- X type
    // -- number of Xs
    // -- number of attributes per X
    // -- number of nodes per X (derivable from type)
    // -- number of faces per X (derivable from type)
    // -- number of edges per X (derivable from type)

    // In a parallel execution, it is possible that an X block will have
    // no Xs on a particular processor...

    // NOTE: This routine may be called multiple times on a single database.
    //       make sure it is not dependent on being called one time only...

    // Get exodus X block metadata
    if (m_groupCount[entity_type] == 0) {
      return;
    }

    // If the model contains assemblies, we want to retain the empty blocks since the blocks
    // might be in an assembly.  This is typically the case when an application is running
    // in parallel, but is telling IOSS that it is "serial" (MPI_COMM_SELF) and taking care
    // of synchronization at the app level instead of down here...
    bool retain_empty_blocks = m_groupCount[EX_ASSEMBLY] > 0;

    // The application can override this setting using the `RETAIN_EMPTY_BLOCKS` property.
    // This can either set to TRUE or FALSE...  Note that `retain_empty_blocks` will not be
    // changed unless the property exists.
    Ioss::Utils::check_set_bool_property(properties, "RETAIN_EMPTY_BLOCKS", retain_empty_blocks);

    Ioss::Int64Vector X_block_ids(m_groupCount[entity_type]);

    int error;
    {
      Ioss::SerializeIO serializeIO_(this);

      if ((ex_int64_status(get_file_pointer()) & EX_IDS_INT64_API) != 0) {
        error = ex_get_ids(get_file_pointer(), entity_type, Data(X_block_ids));
      }
      else {
        Ioss::IntVector tmp_set_ids(X_block_ids.size());
        error = ex_get_ids(get_file_pointer(), entity_type, Data(tmp_set_ids));
        if (error >= 0) {
          std::copy(tmp_set_ids.begin(), tmp_set_ids.end(), X_block_ids.begin());
        }
      }
      if (error < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }
    }

    size_t            all_X_type_length = m_groupCount[entity_type] * (MAX_STR_LENGTH + 1);
    std::vector<char> all_X_type(all_X_type_length);

    Ioss::Int64Vector counts(m_groupCount[entity_type] * 4);
    Ioss::Int64Vector local_X_count(m_groupCount[entity_type]);
    Ioss::Int64Vector global_X_count(m_groupCount[entity_type]);
    int               iblk;

    {
      Ioss::SerializeIO serializeIO_(this);

      for (iblk = 0; iblk < m_groupCount[entity_type]; iblk++) {
        int     index = 4 * iblk;
        int64_t id    = X_block_ids[iblk];

        char *const X_type = Data(all_X_type) + iblk * (MAX_STR_LENGTH + 1);

        ex_block block{};
        block.id   = id;
        block.type = entity_type;
        error      = ex_get_block_param(get_file_pointer(), &block);
        if (error < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        local_X_count[iblk] = block.num_entry;

        counts[index + 0] = block.num_nodes_per_entry;
        counts[index + 1] = block.num_edges_per_entry;
        counts[index + 2] = block.num_faces_per_entry;
        counts[index + 3] = block.num_attribute;

        if (block.num_entry == 0) {
          std::memset(X_type, 0, MAX_STR_LENGTH + 1);
        }
        else {
          Ioss::Utils::copy_string(X_type, block.topology, MAX_STR_LENGTH + 1);
        }
      }
    }

    // This is a collective call...
    util().attribute_reduction(all_X_type_length, Data(all_X_type));

    // This is a collective call...
    util().global_array_minmax(counts, Ioss::ParallelUtils::DO_MAX);

    // Determine global X count for each X block....
    // Can also get this from an nemesis call, but the data may not always be there in all cases.
    util().global_count(local_X_count, global_X_count);

    // The 'offset' is used to map an X location within an X
    // block to the X 'file descriptor'.  For example, the file
    // descriptor of the 37th X in the 4th block is calculated by:
    // file_descriptor = offset of block 4 + 37 This can also be used to
    // determine which X block an X with a file_descriptor
    // maps into. An particular X block contains all Xs in
    // the range:
    //     offset < file_descriptor <= offset+number_Xs_per_block
    int used_blocks = 0;

    int nvar = std::numeric_limits<int>::max(); // Number of 'block' vars on database. Used to skip
                                                // querying if none.
    int nmap = std::numeric_limits<int>::max(); // Number of 'block' vars on database. Used to skip
                                                // querying if none.
    for (iblk = 0; iblk < m_groupCount[entity_type]; iblk++) {
      int     index       = 4 * iblk;
      int64_t nodes_per_X = counts[index + 0];
      int64_t edges_per_X = counts[index + 1];
      int64_t faces_per_X = counts[index + 2];
      int64_t attributes  = counts[index + 3];

      int64_t     id     = X_block_ids[iblk];
      std::string alias  = Ioss::Utils::encode_entity_name(basename, id);
      char *const X_type = Data(all_X_type) + iblk * (MAX_STR_LENGTH + 1);

      bool        db_has_name = false;
      std::string block_name;
      if (ignore_database_names()) {
        block_name = alias;
      }
      else {
        Ioss::SerializeIO serializeio_(this);
        block_name = Ioex::get_entity_name(get_file_pointer(), entity_type, id, basename,
                                           maximumNameLength, db_has_name);
      }
      if (get_use_generic_canonical_name()) {
        std::swap(block_name, alias);
      }

      std::string save_type = X_type;
      std::string type =
          Ioss::Utils::fixup_type(X_type, nodes_per_X, spatialDimension - rank_offset);
      if (local_X_count[iblk] == 0 && type.empty()) {
        // For an empty block, exodus does not store the X
        // type information and returns "nullptr" If there are no
        // Xs on any processors for this block, it will have
        // an empty type which is invalid and will throw an
        // exception in the XBlock constructor. Try to discern
        // the correct X type based on the block_name.
        auto tokens = Ioss::tokenize(block_name, "_");
        if (tokens.size() >= 2) {
          // Check whether last token names an X topology type...
          const Ioss::ElementTopology *topology =
              Ioss::ElementTopology::factory(tokens.back(), true);
          if (topology != nullptr) {
            type = topology->name();
          }
        }
      }

      if (type == "null" || type.empty()) {
        // If we have no idea what the topology type for an empty
        // X block is, call it "unknown"
        type = "unknown";

        // If there are no Xs on any processor for this block and
        // we have no idea what the topology type is, skip it...
        if (!retain_empty_blocks && global_X_count[iblk] == 0) {
          continue;
        }
      }

      Ioss::EntityBlock *block = nullptr;
      if (entity_type == EX_ELEM_BLOCK) {
        auto eblock = new Ioss::ElementBlock(this, block_name, type, local_X_count[iblk]);
        block       = eblock;
        block->property_add(
            Ioss::Property("id", id)); // Do before adding for better error messages.
        get_region()->add(eblock);
      }
      else if (entity_type == EX_FACE_BLOCK) {
        auto fblock = new Ioss::FaceBlock(this, block_name, type, local_X_count[iblk]);
        block       = fblock;
        block->property_add(
            Ioss::Property("id", id)); // Do before adding for better error messages.
        get_region()->add(fblock);
      }
      else if (entity_type == EX_EDGE_BLOCK) {
        auto eblock = new Ioss::EdgeBlock(this, block_name, type, local_X_count[iblk]);
        block       = eblock;
        block->property_add(
            Ioss::Property("id", id)); // Do before adding for better error messages.
        get_region()->add(eblock);
      }
      else {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Invalid type in get_blocks()");
        IOSS_ERROR(errmsg);
      }

      // See which connectivity options were defined for this block.
      // X -> Node is always defined.
      // X -> Face?
      if (faces_per_X > 0 && rank_offset < 1) {
        std::string storage = "Real[" + std::to_string(faces_per_X) + "]";
        block->field_add(
            Ioss::Field("connectivity_face", block->field_int_type(), storage, Ioss::Field::MESH));
      }
      // X -> Edge?
      if (edges_per_X > 0 && rank_offset < 2) {
        std::string storage = "Real[" + std::to_string(edges_per_X) + "]";
        block->field_add(
            Ioss::Field("connectivity_edge", block->field_int_type(), storage, Ioss::Field::MESH));
      }

      block->property_add(Ioss::Property("guid", util().generate_guid(id)));
      if (db_has_name) {
        std::string *db_name = &block_name;
        if (get_use_generic_canonical_name()) {
          db_name = &alias;
        }
        if (alias != block_name) {
          block->property_add(Ioss::Property("db_name", *db_name));
        }
      }

      // Maintain block order on output database...
      block->property_add(Ioss::Property("original_block_order", used_blocks++));

      if (save_type != "null" && !save_type.empty()) {
        block->property_update("original_topology_type", save_type);
      }

      block->property_add(Ioss::Property("global_entity_count", global_X_count[iblk]));

      get_region()->add_alias(block_name, alias, block->type());

      // Check for additional variables.
      add_attribute_fields(block, attributes, type);
      if (nvar > 0) {
        nvar = add_results_fields(block, iblk);
      }
      add_reduction_results_fields(block);
      add_mesh_reduction_fields(id, block);

      if (entity_type == EX_ELEM_BLOCK) {
        Ioss::SerializeIO serializeio_(this);
        if (nmap > 0) {
          auto *elb = dynamic_cast<Ioss::ElementBlock *>(block);
          Ioss::Utils::check_dynamic_cast(elb);
          nmap =
              Ioex::add_map_fields(get_file_pointer(), elb, local_X_count[iblk], maximumNameLength);
        }
      }
    }
    m_groupCount[entity_type] = used_blocks;

    if (entity_type == EX_ELEM_BLOCK) {
      if (!assemblyOmissions.empty() || !assemblyInclusions.empty()) {
        update_block_omissions_from_assemblies();
      }

      assert(blockOmissions.empty() || blockInclusions.empty()); // Only one can be non-empty

      // Handle all block omissions or inclusions...
      // This only affects the generation of surfaces...
      if (!blockOmissions.empty()) {
        for (const auto &name : blockOmissions) {
          auto block = get_region()->get_element_block(name);
          if (block != nullptr) {
            block->property_add(Ioss::Property(std::string("omitted"), 1));
          }
        }
      }

      if (!blockInclusions.empty()) {
        const auto &blocks = get_region()->get_element_blocks();
        for (auto &block : blocks) {
          block->property_add(Ioss::Property(std::string("omitted"), 1));
        }

        // Now, erase the property on any blocks in the inclusion list...
        for (const auto &name : blockInclusions) {
          auto block = get_region()->get_element_block(name);
          if (block != nullptr) {
            block->property_erase("omitted");
          }
        }
      }
    }
  }

  void DatabaseIO::compute_node_status() const
  {
    // Create a field for all nodes in the model indicating
    // the connectivity 'status' of the node.  The status values are:
    // 0 -- node not connected to any elements
    // 1 -- node only connected to omitted elements
    // 2 -- node only connected to active elements
    // 3 -- node at border of active and omitted elements.

    /// \todo Get working for parallel...

    if (nodeConnectivityStatusCalculated) {
      return;
    }

    nodeConnectivityStatus.resize(nodeCount);

    const Ioss::ElementBlockContainer &element_blocks = get_region()->get_element_blocks();
    assert(Ioss::Utils::check_block_order(element_blocks));

    for (Ioss::ElementBlock *block : element_blocks) {
      unsigned char status = 2;
      if (Ioss::Utils::block_is_omitted(block)) {
        status = 1;
      }

      int64_t id               = block->get_property("id").get_int();
      int     element_nodes    = block->topology()->number_nodes();
      int64_t my_element_count = block->entity_count();
      if (my_element_count > 0) {
        if ((ex_int64_status(get_file_pointer()) & EX_BULK_INT64_API) != 0) {
          std::vector<int64_t> conn(my_element_count * element_nodes);
          ex_get_conn(get_file_pointer(), EX_ELEM_BLOCK, id, Data(conn), nullptr, nullptr);
          for (int64_t j = 0; j < my_element_count * element_nodes; j++) {
            nodeConnectivityStatus[conn[j] - 1] |= status;
          }
        }
        else {
          std::vector<int> conn(my_element_count * element_nodes);
          ex_get_conn(get_file_pointer(), EX_ELEM_BLOCK, id, Data(conn), nullptr, nullptr);
          for (int64_t j = 0; j < my_element_count * element_nodes; j++) {
            nodeConnectivityStatus[conn[j] - 1] |= status;
          }
        }
      }
    }
    nodeConnectivityStatusCalculated = true;
  }

  namespace {
    void get_element_sides_lists(int exoid, int64_t id, int int_byte_size, int64_t number_sides,
                                 Ioss::Int64Vector &element, Ioss::Int64Vector &sides)
    {
      // Check whether we already populated the element/sides vectors.
      if (element.empty() && sides.empty() && number_sides > 0) {
        element.resize(number_sides);
        sides.resize(number_sides);
        // Easier below here if the element and sides are a known 64-bit size...
        // Kluge here to do that...
        if (int_byte_size == 4) {
          Ioss::IntVector e32(number_sides);
          Ioss::IntVector s32(number_sides);
          int             ierr = ex_get_set(exoid, EX_SIDE_SET, id, Data(e32), Data(s32));
          if (ierr < 0) {
            Ioex::exodus_error(exoid, __LINE__, __func__, __FILE__);
          }
          std::copy(e32.begin(), e32.end(), element.begin());
          std::copy(s32.begin(), s32.end(), sides.begin());
        }
        else {
          int ierr = ex_get_set(exoid, EX_SIDE_SET, id, Data(element), Data(sides));
          if (ierr < 0) {
            Ioex::exodus_error(exoid, __LINE__, __func__, __FILE__);
          }
        }
      }
    }
  } // namespace

  void DatabaseIO::get_sidesets()
  {
    // This function creates all sidesets (surfaces) for a
    // model.  Note that a sideset contains 1 or more sideblocks
    // which are homogeneous (same topology). In serial execution,
    // this is fairly straightforward since there are no null sets and
    // we have all the information we need. (...except see below for
    // surface evolution).
    //
    // However, in a parallel execution, we have the possibility that a
    // side set will have no sides or distribution factors on
    // a particular processor.  We then don't know the block topology of
    // the block(s) contained in this set. We could do some
    // communication and get a good idea of the topologies that are in
    // the set.

    if (m_groupCount[EX_SIDE_SET] > 0) {
      check_side_topology();

      // Get exodus sideset metadata

      // Get the names (may not exist) of all sidesets and see if they are actually
      // side "blocks" (perhaps written by IO system for a restart).  In that case,
      // they were split by a previous run and we need to reconstruct the side "set"
      // that may contain one or more of them.
      Ioex::SideSetMap fs_map;
      Ioex::SideSetSet fs_set;

      Ioss::Int64Vector side_set_ids(m_groupCount[EX_SIDE_SET]);
      {
        Ioss::SerializeIO serializeio_(this);
        int               error;
        if ((ex_int64_status(get_file_pointer()) & EX_IDS_INT64_API) != 0) {
          error = ex_get_ids(get_file_pointer(), EX_SIDE_SET, Data(side_set_ids));
        }
        else {
          Ioss::IntVector tmp_set_ids(side_set_ids.size());
          error = ex_get_ids(get_file_pointer(), EX_SIDE_SET, Data(tmp_set_ids));
          if (error >= 0) {
            std::copy(tmp_set_ids.begin(), tmp_set_ids.end(), side_set_ids.begin());
          }
        }
        if (error < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        for (const auto &id : side_set_ids) {
          std::vector<char> ss_name(maximumNameLength + 1);
          error = ex_get_name(get_file_pointer(), EX_SIDE_SET, id, Data(ss_name));
          if (error < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
          if (ss_name[0] != '\0') {
            Ioss::Utils::fixup_name(Data(ss_name));
            Ioex::decode_surface_name(fs_map, fs_set, Data(ss_name));
          }
        }
      }

      // Create sidesets for each entry in the fs_set... These are the
      // sidesets which were probably written by a previous run of the
      // IO system and are already split into homogeneous pieces...
      {
        for (const auto &fs_name : fs_set) {
          auto side_set = new Ioss::SideSet(this, fs_name);
          get_region()->add(side_set);
          int64_t id = Ioex::extract_id(fs_name);
          if (id > 0) {
            side_set->property_add(Ioss::Property("id", id));
            side_set->property_add(Ioss::Property("guid", util().generate_guid(id)));
          }
        }
      }

      for (int iss = 0; iss < m_groupCount[EX_SIDE_SET]; iss++) {
        int64_t           id = side_set_ids[iss];
        std::string       sid;
        Ioex::TopologyMap topo_map;
        Ioex::TopologyMap side_map; // Used to determine side consistency

        Ioss::SurfaceSplitType split_type = splitType;
        std::string            side_set_name;
        Ioss::SideSet         *side_set = nullptr;

        bool db_has_name = false;
        {
          Ioss::SerializeIO serializeio_(this);

          std::string alias = Ioss::Utils::encode_entity_name("surface", id);
          if (ignore_database_names()) {
            side_set_name = alias;
          }
          else {
            side_set_name = Ioex::get_entity_name(get_file_pointer(), EX_SIDE_SET, id, "surface",
                                                  maximumNameLength, db_has_name);
          }

          if (side_set_name == "universal_sideset") {
            split_type = Ioss::SPLIT_BY_DONT_SPLIT;
          }

          bool in_fs_map = false;
          auto FSM       = fs_map.find(side_set_name);
          if (FSM != fs_map.end()) {
            in_fs_map            = true;
            std::string efs_name = (*FSM).second;
            side_set             = get_region()->get_sideset(efs_name);
            Ioss::Utils::check_non_null(side_set, "sideset", efs_name, __func__);
          }
          else {
            if (get_use_generic_canonical_name()) {
              std::swap(side_set_name, alias);
            }
            side_set = new Ioss::SideSet(this, side_set_name);
            side_set->property_add(Ioss::Property("id", id));
            side_set->property_add(Ioss::Property("guid", util().generate_guid(id)));
            if (db_has_name) {
              std::string *db_name = &side_set_name;
              if (get_use_generic_canonical_name()) {
                db_name = &alias;
              }
              if (alias != side_set_name) {
                side_set->property_add(Ioss::Property("db_name", *db_name));
              }
            }

            get_region()->add(side_set);

            get_region()->add_alias(side_set_name, alias, Ioss::SIDESET);
            get_region()->add_alias(side_set_name, Ioss::Utils::encode_entity_name("sideset", id),
                                    Ioss::SIDESET);
          }

          //      split_type = SPLIT_BY_ELEMENT_BLOCK;
          //      split_type = SPLIT_BY_TOPOLOGIES;
          //      split_type = SPLIT_BY_DONT_SPLIT;

          // Determine how many side blocks compose this side set.
          ex_set set_param[1];
          set_param[0].id                       = id;
          set_param[0].type                     = EX_SIDE_SET;
          set_param[0].entry_list               = nullptr;
          set_param[0].extra_list               = nullptr;
          set_param[0].distribution_factor_list = nullptr;

          int error = ex_get_sets(get_file_pointer(), 1, set_param);
          if (error < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }

          int64_t number_sides = set_param[0].num_entry;

          Ioss::Int64Vector element;
          Ioss::Int64Vector sides;

          if (!blockOmissions.empty() || !blockInclusions.empty()) {
            get_element_sides_lists(get_file_pointer(), id, int_byte_size_api(), number_sides,
                                    element, sides);
            Ioex::filter_element_list(get_region(), element, sides, true);
            number_sides = element.size();
            SMART_ASSERT(element.size() == sides.size())(element.size())(sides.size());
          }

          if (split_type == Ioss::SPLIT_BY_TOPOLOGIES && sideTopology.size() == 1) {
            // There is only one side type for all elements in the model
            topo_map[std::make_pair(sideTopology[0].first->name(), sideTopology[0].second)] =
                number_sides;
          }
          else if (split_type == Ioss::SPLIT_BY_DONT_SPLIT) {
            const Ioss::ElementTopology *mixed_topo = Ioss::ElementTopology::factory("unknown");
            topo_map[std::make_pair(std::string("unknown"), mixed_topo)] = number_sides;
          }
          else if (in_fs_map) {
            auto tokens = Ioss::tokenize(side_set_name, "_");
            SMART_ASSERT(tokens.size() >= 4)(tokens.size());
            // The sideset should have only a single topology which is
            // given by the sideset name...
            const Ioss::ElementTopology *side_topo =
                Ioss::ElementTopology::factory(tokens[tokens.size() - 2]);
            assert(side_topo != nullptr);
            const Ioss::ElementTopology *element_topo =
                Ioss::ElementTopology::factory(tokens[tokens.size() - 3], true);
            std::string name;
            if (element_topo != nullptr) {
              name = element_topo->name();
            }
            else {
              //                           -4   -3   -2     -1
              // Name is of the form name_block_id_sidetopo_id
              name = tokens[tokens.size() - 4] + "_" + tokens[tokens.size() - 3];
            }

            topo_map[std::make_pair(name, side_topo)] = number_sides;

            // We want the id to match the id on the sideset in this
            // case so that the generated name will match the current
            // name.  Instead of converting from string to int back to
            // string, we just set a variable to query later.
            sid = tokens.back();
          }
          else if (split_type == Ioss::SPLIT_BY_TOPOLOGIES) {
            // There are multiple side types in the model.
            // Iterate through the elements in the sideset, determine
            // their parent element block using the blocks element
            // topology and the side number, determine the side
            // type.

            for (const auto &elem : sideTopology) {
              topo_map[std::make_pair(elem.first->name(), elem.second)] = 0;
              side_map[std::make_pair(elem.first->name(), elem.second)] = 0;
            }

            get_element_sides_lists(get_file_pointer(), id, int_byte_size_api(), number_sides,
                                    element, sides);
            Ioex::separate_surface_element_sides(element, sides, get_region(), topo_map, side_map,
                                                 split_type, side_set_name);
          }
          else if (split_type == Ioss::SPLIT_BY_ELEMENT_BLOCK) {
            // There are multiple side types in the model.  Iterate
            // through the elements in the sideset, determine their parent
            // element block using blocks element topology and the side
            // number, determine the side type.

            // Seed the topo_map map with <block->name, side_topo>
            // pairs so we are sure that all processors have the same
            // starting topo_map (size and order).
            const Ioss::ElementBlockContainer &element_blocks = get_region()->get_element_blocks();
            assert(Ioss::Utils::check_block_order(element_blocks));

            for (Ioss::ElementBlock *block : element_blocks) {
              if (!Ioss::Utils::block_is_omitted(block)) {
                std::string                  name         = block->name();
                const Ioss::ElementTopology *common_ftopo = block->topology()->boundary_type(0);
                if (common_ftopo != nullptr) {
                  // All sides of this element block's topology have the same topology
                  topo_map[std::make_pair(name, common_ftopo)] = 0;
                  side_map[std::make_pair(name, common_ftopo)] = 0;
                }
                else {
                  // The sides have different topology, iterate over
                  // them and create an entry for the unique side
                  // topology types
                  int par_dim = block->topology()->parametric_dimension();
                  if (par_dim == 2 || par_dim == 3) {
                    int64_t my_side_count = block->topology()->number_boundaries();
                    for (int64_t ii = 0; ii < my_side_count; ii++) {
                      const Ioss::ElementTopology *topo = block->topology()->boundary_type(ii + 1);
                      topo_map[std::make_pair(name, topo)] = 0;
                      side_map[std::make_pair(name, topo)] = 0;
                    }
                  }
                }
              }
            }
            get_element_sides_lists(get_file_pointer(), id, int_byte_size_api(), number_sides,
                                    element, sides);
            Ioex::separate_surface_element_sides(element, sides, get_region(), topo_map, side_map,
                                                 split_type, side_set_name);
          }
        }

        // End of first step in splitting.  Check among all processors
        // to see which potential splits have sides in them...
        Ioss::Int64Vector global_side_counts(topo_map.size());
        {
          int64_t i = 0;
          for (const auto &topo : topo_map) {
            global_side_counts[i++] = topo.second;
          }

          // If splitting by element block, also sync the side_map
          // information which specifies whether the sideset has
          // consistent sides for all elements. Only really used for
          // shells, but easier to just set the value on all surfaces
          // in the element block split case.
          if (side_map.size() == topo_map.size()) {
            global_side_counts.resize(topo_map.size() + side_map.size());
            for (const auto &side : side_map) {
              global_side_counts[i++] = side.second;
            }
          }

          // See if any processor has non-zero count for the topo_map counts
          // For the side_map, need the max value.
          util().global_array_minmax(global_side_counts, Ioss::ParallelUtils::DO_MAX);
        }

        // Create Side Blocks

        int64_t i = 0;
        for (const auto &topo : topo_map) {
          if (global_side_counts[i++] > 0) {
            const std::string            topo_or_block_name = topo.first.first;
            const Ioss::ElementTopology *side_topo          = topo.first.second;
            assert(side_topo != nullptr);
#if 0
            if (side_topo->parametric_dimension() == topology_dimension-1 ||
                split_type == Ioss::SPLIT_BY_DONT_SPLIT) {
#else
            if (true) {
#endif
              int64_t my_side_count = topo.second;

              std::string side_block_name =
                  "surface_" + topo_or_block_name + "_" + side_topo->name();
              if (split_type == Ioss::SPLIT_BY_DONT_SPLIT) {
                side_block_name = side_set_name;
              }
              else {
                if (db_has_name) {
                  side_block_name =
                      side_set->name() + "_" + topo_or_block_name + "_" + side_topo->name();
                }
                else {
                  if (sid.empty()) {
                    side_block_name = Ioss::Utils::encode_entity_name(side_block_name, id);
                  }
                  else {
                    side_block_name += "_";
                    side_block_name += sid;
                  }
                }
              }

              Ioss::ElementBlock *block = nullptr;
              // Need to get elem_topo....
              const Ioss::ElementTopology *elem_topo = nullptr;
              if (split_type == Ioss::SPLIT_BY_TOPOLOGIES) {
                elem_topo = Ioss::ElementTopology::factory(topo_or_block_name);
              }
              else if (split_type == Ioss::SPLIT_BY_ELEMENT_BLOCK) {
                block = get_region()->get_element_block(topo_or_block_name);
                if (block == nullptr || Ioss::Utils::block_is_omitted(block)) {
                  std::ostringstream errmsg;
                  fmt::print(errmsg,
                             "INTERNAL ERROR: Could not find element block '{}' Something is wrong "
                             "in the Ioex::DatabaseIO class. Please report.\n",
                             topo_or_block_name);
                  IOSS_ERROR(errmsg);
                }
                elem_topo = block->topology();
              }
              else if (split_type == Ioss::SPLIT_BY_DONT_SPLIT) {
                // Most likely this is "unknown", but can be a true
                // topology if there is only a single element block in
                // the model.
                elem_topo = Ioss::ElementTopology::factory(topo_or_block_name);
              }
              else {
                std::ostringstream errmsg;
                fmt::print(
                    errmsg,
                    "INTERNAL ERROR: Invalid setting for `split_type` {}. Something is wrong "
                    "in the Ioex::DatabaseIO class. Please report.\n",
                    static_cast<int>(split_type));
                IOSS_ERROR(errmsg);
              }
              assert(elem_topo != nullptr);

              auto side_block = new Ioss::SideBlock(this, side_block_name, side_topo->name(),
                                                    elem_topo->name(), my_side_count);
              side_set->add(side_block);

              // Note that all sideblocks within a specific
              // sideset might have the same id.
              assert(side_block != nullptr);
              side_block->property_add(Ioss::Property("id", id));
              side_block->property_add(Ioss::Property("guid", util().generate_guid(id)));

              // If splitting by element block, need to set the
              // element block member on this side block.
              if (split_type == Ioss::SPLIT_BY_ELEMENT_BLOCK) {
                side_block->set_parent_element_block(block);
              }

              // If we calculated whether the element side is
              // consistent for all sides in this block, then
              // tell the block which side it is, or that they are
              // inconsistent. If it wasn't calculated above, then it
              // will be calculated on the fly when/if requested.
              // This is to avoid reading the sideset bulk data in
              // cases where we don't need to read it, but if we are
              // already reading it (to split the sidesets), then use
              // the data when we have it.
              if (!side_map.empty()) {
                // Set a property indicating which element side
                // (1-based) all sides in this block are applied to.
                // If they are not all assigned to the same element
                // side, indicate this with a side equal to 0.
                //
                // (note: 'i' has already been incremented earlier in
                // the loop.  We need previous value here...)
                int side = global_side_counts[i - 1 + topo_map.size()];
                if (side == 999) {
                  side = 0;
                }
                SMART_ASSERT(side <= elem_topo->number_boundaries())
                (side)(elem_topo->number_boundaries());
                side_block->set_consistent_side_number(side);
              }

              // Add an alias...
              get_region()->add_alias(side_block);

              if (split_type != Ioss::SPLIT_BY_DONT_SPLIT && side_set_name != "universal_sideset") {
                std::string storage = "Real[";
                storage += std::to_string(side_topo->number_nodes());
                storage += "]";
                side_block->field_add(Ioss::Field("distribution_factors", Ioss::Field::REAL,
                                                  storage, Ioss::Field::MESH));
              }

              if (side_set_name == "universal_sideset") {
                side_block->field_add(Ioss::Field("side_ids", side_block->field_int_type(),
                                                  "scalar", Ioss::Field::MESH));
              }

              int num_attr = 0;
              {
                Ioss::SerializeIO serializeio_(this);
                int ierr = ex_get_attr_param(get_file_pointer(), EX_SIDE_SET, 1, &num_attr);
                if (ierr < 0) {
                  Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
                }
              }
              // Add additional fields
              add_attribute_fields(side_block, num_attr, "");
              add_results_fields(side_block, iss);
            }
          }
        }
      }
    }
  } // namespace Ioex

  template <typename T>
  void DatabaseIO::get_sets(ex_entity_type type, int64_t count, const std::string &base,
                            const T * /* set_type */)
  {
    // Attributes of a Xset are:
    // -- id
    // -- name
    // -- number of nodes
    // -- number of distribution factors (see next comment)
    // ----the #distribution factors should equal #Xs or 0, any
    //     other value does not make sense. If it is 0, then a substitute
    //     list will be created returning 1.0 for the factor

    // In a parallel execution, it is possible that a Xset will have
    // no Xs or distribution factors on a particular processor...

    // Get exodus Xset metadata
    if (count > 0) {
      Ioss::Int64Vector Xset_ids(count);
      Ioss::IntVector   attributes(count);
      std::vector<T *>  Xsets(count);
      {
        Ioss::SerializeIO serializeio_(this);
        if (ex_int64_status(get_file_pointer()) & EX_IDS_INT64_API) {
          int error = ex_get_ids(get_file_pointer(), type, Data(Xset_ids));
          if (error < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
        else {
          Ioss::IntVector tmp_set_ids(count);
          int             error = ex_get_ids(get_file_pointer(), type, Data(tmp_set_ids));
          if (error < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
          std::copy(tmp_set_ids.begin(), tmp_set_ids.end(), Xset_ids.begin());
        }

        std::vector<ex_set> set_params(count);
        for (int ins = 0; ins < count; ins++) {
          set_params[ins].type                     = type;
          set_params[ins].id                       = Xset_ids[ins];
          set_params[ins].entry_list               = nullptr;
          set_params[ins].extra_list               = nullptr;
          set_params[ins].distribution_factor_list = nullptr;
        }

        int error = ex_get_sets(get_file_pointer(), count, Data(set_params));
        if (error < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        for (int ins = 0; ins < count; ins++) {
          int64_t id       = set_params[ins].id;
          int     num_attr = 0;
          int     ierr     = ex_get_attr_param(get_file_pointer(), type, id, &num_attr);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
          attributes[ins] = num_attr;

          bool        db_has_name = false;
          std::string alias       = Ioss::Utils::encode_entity_name(base + "list", id);
          std::string Xset_name;
          if (ignore_database_names()) {
            Xset_name = alias;
          }
          else {
            Xset_name = Ioex::get_entity_name(get_file_pointer(), type, id, base + "list",
                                              maximumNameLength, db_has_name);
          }

          if (get_use_generic_canonical_name()) {
            std::swap(Xset_name, alias);
          }

          bool              filtered          = false;
          int64_t           original_set_size = set_params[ins].num_entry;
          Ioss::Int64Vector active_node_index;
          if ((!blockOmissions.empty() || !blockInclusions.empty()) && type == EX_NODE_SET) {
            active_node_index.resize(set_params[ins].num_entry);
            set_params[ins].entry_list = Data(active_node_index);

            int old_status = ex_int64_status(get_file_pointer());
            ex_set_int64_status(get_file_pointer(), EX_BULK_INT64_API);
            error = ex_get_sets(get_file_pointer(), 1, &set_params[ins]);
            if (error < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
            ex_set_int64_status(get_file_pointer(), old_status);

            compute_node_status();
            filtered = Ioex::filter_node_list(active_node_index, nodeConnectivityStatus);
            set_params[ins].num_entry = active_node_index.size();
          }
          auto *Xset = new T(this, Xset_name, set_params[ins].num_entry);
          Xsets[ins] = Xset;
          Xset->property_add(Ioss::Property("id", id));
          Xset->property_add(Ioss::Property("guid", util().generate_guid(id)));
          if (db_has_name) {
            std::string *db_name = &Xset_name;
            if (get_use_generic_canonical_name()) {
              db_name = &alias;
            }
            if (alias != Xset_name) {
              Xset->property_add(Ioss::Property("db_name", *db_name));
            }
          }
          if (filtered && type == EX_NODE_SET) {
            Xset->property_add(Ioss::Property("filtered_db_set_size", original_set_size));
            activeNodeSetNodesIndex[Xset_name].swap(active_node_index);
          }
          get_region()->add(Xset);
          get_region()->add_alias(Xset_name, alias, Xset->type());
          get_region()->add_alias(Xset_name, Ioss::Utils::encode_entity_name(base + "set", id),
                                  Xset->type());
        }
      }

      // The attribute count will either be 0 if there are no
      // entities in the grouping entity on this processor, or it will be
      // the number of attributes (> 0). Therefore, if we take the 'max'
      // over all processors, each processor will then have the correct
      // attribute count...
      // This is a collective call...
      util().global_array_minmax(attributes, Ioss::ParallelUtils::DO_MAX);

      for (int ins = 0; ins < count; ins++) {
        add_attribute_fields(Xsets[ins], attributes[ins], "");
        add_results_fields(Xsets[ins], ins);
      }
    }
  }

  void DatabaseIO::get_nodesets()
  {
    get_sets(EX_NODE_SET, m_groupCount[EX_NODE_SET], "node", (Ioss::NodeSet *)nullptr);
  }

  void DatabaseIO::get_edgesets()
  {
    get_sets(EX_EDGE_SET, m_groupCount[EX_EDGE_SET], "edge", (Ioss::EdgeSet *)nullptr);
  }

  void DatabaseIO::get_facesets()
  {
    get_sets(EX_FACE_SET, m_groupCount[EX_FACE_SET], "face", (Ioss::FaceSet *)nullptr);
  }

  void DatabaseIO::get_elemsets()
  {
    get_sets(EX_ELEM_SET, m_groupCount[EX_ELEM_SET], "element", (Ioss::ElementSet *)nullptr);
  }

  void DatabaseIO::get_commsets()
  {
    // Attributes of a commset are:
    // -- id (property)
    // -- name (property)
    // -- number of node--CPU pairs (field)

    // In a parallel execution, it is possible that a commset will have
    // no nodes on a particular processor...

    // If this is a serial execution, there will be no communication
    // nodesets, just return an empty container.

    if (isParallel || isSerialParallel) {
      Ioss::SerializeIO serializeio_(this);
      // This is a parallel run. There should be communications data
      // Get nemesis commset metadata
      int64_t my_node_count = 0;
      int64_t elem_count    = 0;

      // NOTE: It is possible for a parallel run to have no
      // communications maps if the decomposition occurs along contact
      // surfaces.  In this case, we create empty node and element
      // communication maps.
      if (commsetNodeCount > 0 || commsetElemCount > 0) {
        if (commsetNodeCount > 0) {
          nodeCmapIds.resize(commsetNodeCount);
          nodeCmapNodeCnts.resize(commsetNodeCount);
        }
        if (commsetElemCount > 0) {
          elemCmapIds.resize(commsetElemCount);
          elemCmapElemCnts.resize(commsetElemCount);
        }

        int error;
        if (int_byte_size_api() == 4) {
          Ioss::IntVector nci(nodeCmapIds.size());
          Ioss::IntVector ncnc(nodeCmapNodeCnts.size());
          Ioss::IntVector eci(elemCmapIds.size());
          Ioss::IntVector ecec(elemCmapElemCnts.size());
          error = ex_get_cmap_params(get_file_pointer(), Data(nci), Data(ncnc), Data(eci),
                                     Data(ecec), myProcessor);
          if (error >= 0) {
            std::copy(nci.begin(), nci.end(), nodeCmapIds.begin());
            std::copy(ncnc.begin(), ncnc.end(), nodeCmapNodeCnts.begin());
            std::copy(eci.begin(), eci.end(), elemCmapIds.begin());
            std::copy(ecec.begin(), ecec.end(), elemCmapElemCnts.begin());
          }
        }
        else {
          error = ex_get_cmap_params(get_file_pointer(), Data(nodeCmapIds), Data(nodeCmapNodeCnts),
                                     Data(elemCmapIds), Data(elemCmapElemCnts), myProcessor);
        }
        if (error < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        // Count nodes, elements, and convert counts to offsets.
        my_node_count +=
            std::accumulate(nodeCmapNodeCnts.begin(), nodeCmapNodeCnts.end(), int64_t(0));

        elem_count += std::accumulate(elemCmapElemCnts.begin(), elemCmapElemCnts.end(), int64_t(0));
      }
      // Create a single node commset and a single element commset
      auto *commset = new Ioss::CommSet(this, "commset_node", "node", my_node_count);
      commset->property_add(Ioss::Property("id", 1));
      commset->property_add(Ioss::Property("guid", util().generate_guid(1)));
      get_region()->add(commset);

      commset = new Ioss::CommSet(this, "commset_side", "side", elem_count);
      commset->property_add(Ioss::Property("id", 1));
      commset->property_add(Ioss::Property("guid", util().generate_guid(1)));
      get_region()->add(commset);
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::Region *reg, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return Ioex::BaseDatabaseIO::get_field_internal(reg, field, data, data_size);
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::NodeBlock *nb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

#ifndef NDEBUG
        int64_t my_node_count = field.raw_count();
        SMART_ASSERT(my_node_count == nodeCount)(my_node_count)(nodeCount);
#endif

        Ioss::Field::RoleType role = field.get_role();
        if (role == Ioss::Field::MESH) {
          if (field.get_name() == "mesh_model_coordinates_x") {
            auto *rdata = static_cast<double *>(data);
            int   ierr  = ex_get_coord(get_file_pointer(), rdata, nullptr, nullptr);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }

          else if (field.get_name() == "mesh_model_coordinates_y") {
            auto *rdata = static_cast<double *>(data);
            int   ierr  = ex_get_coord(get_file_pointer(), nullptr, rdata, nullptr);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }

          else if (field.get_name() == "mesh_model_coordinates_z") {
            auto *rdata = static_cast<double *>(data);
            int   ierr  = ex_get_coord(get_file_pointer(), nullptr, nullptr, rdata);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }

          else if (field.get_name() == "mesh_model_coordinates") {
            // Data required by upper classes store x0, y0, z0, ... xn,
            // yn, zn. Data stored in exodus file is x0, ..., xn, y0,
            // ..., yn, z0, ..., zn so we have to allocate some scratch
            // memory to read in the data and then map into supplied
            // 'data'
            std::vector<double> x(num_to_get);
            std::vector<double> y;
            if (spatialDimension > 1) {
              y.resize(num_to_get);
            }
            std::vector<double> z;
            if (spatialDimension == 3) {
              z.resize(num_to_get);
            }

            // Cast 'data' to correct size -- double
            auto *rdata = static_cast<double *>(data);

            int ierr = ex_get_coord(get_file_pointer(), Data(x), Data(y), Data(z));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }

            size_t index = 0;
            for (size_t i = 0; i < num_to_get; i++) {
              rdata[index++] = x[i];
              if (spatialDimension > 1) {
                rdata[index++] = y[i];
              }
              if (spatialDimension == 3) {
                rdata[index++] = z[i];
              }
            }
          }

          else if (field.get_name() == "ids") {
            // Map the local ids in this node block
            // (1...node_count) to global node ids.
            get_map(EX_NODE_BLOCK).map_implicit_data(data, field, num_to_get, 0);
          }

          // The 1..global_node_count id.  In a parallel-decomposed run,
          // it maps the node back to its implicit position in the serial
          // undecomposed mesh file.  This is ONLY provided for backward-
          // compatibility and should not be used unless absolutely required.
          else if (field.get_name() == "implicit_ids") {
            // If not parallel, then this is just 1..node_count
            // If parallel, then it is the data in the ex_get_id_map created by nem_spread.
            if (isParallel) {
              int error = ex_get_id_map(get_file_pointer(), EX_NODE_MAP, data);
              if (error < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
            else {
              if ((ex_int64_status(get_file_pointer()) & EX_BULK_INT64_API) != 0) {
                auto *idata = static_cast<int64_t *>(data);
                for (int64_t i = 0; i < nodeCount; i++) {
                  idata[i] = i + 1;
                }
              }
              else {
                auto *idata = static_cast<int *>(data);
                for (int64_t i = 0; i < nodeCount; i++) {
                  idata[i] = i + 1;
                }
              }
            }
          }

          else if (field.get_name() == "connectivity") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "node_connectivity_status") {
            compute_node_status();
            char *status = static_cast<char *>(data);
            std::copy(nodeConnectivityStatus.begin(), nodeConnectivityStatus.end(), status);
          }

          else if (field.get_name() == "owning_processor") {
            // owning_processor field is always 32-bit.
            if (isParallel) {
              Ioss::CommSet *css   = get_region()->get_commset("commset_node");
              int           *idata = static_cast<int *>(data);
              for (int64_t i = 0; i < nodeCount; i++) {
                idata[i] = myProcessor;
              }

              if ((ex_int64_status(get_file_pointer()) & EX_BULK_INT64_API) != 0) {
                Ioss::Field          ep_field = css->get_field("entity_processor_raw");
                std::vector<int64_t> ent_proc(
                    ep_field.raw_count() * ep_field.get_component_count(Ioss::Field::InOut::INPUT));
                size_t ep_data_size = ent_proc.size() * sizeof(int64_t);
                get_field_internal(css, ep_field, Data(ent_proc), ep_data_size);
                for (size_t i = 0; i < ent_proc.size(); i += 2) {
                  int64_t node = ent_proc[i + 0];
                  int64_t proc = ent_proc[i + 1];
                  if (proc < myProcessor) {
                    idata[node - 1] = proc;
                  }
                }
              }
              else {
                Ioss::Field      ep_field = css->get_field("entity_processor_raw");
                std::vector<int> ent_proc(ep_field.raw_count() *
                                          ep_field.get_component_count(Ioss::Field::InOut::INPUT));
                size_t           ep_data_size = ent_proc.size() * sizeof(int);
                get_field_internal(css, ep_field, Data(ent_proc), ep_data_size);
                for (size_t i = 0; i < ent_proc.size(); i += 2) {
                  int node = ent_proc[i + 0];
                  int proc = ent_proc[i + 1];
                  if (proc < myProcessor) {
                    idata[node - 1] = proc;
                  }
                }
              }
            }
            else {
              // Serial case...
              int *idata = static_cast<int *>(data);
              for (int64_t i = 0; i < nodeCount; i++) {
                idata[i] = 0;
              }
            }
          }
          else {
            num_to_get = Ioss::Utils::field_warning(nb, field, "input");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this node block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // nodes in this block.
          num_to_get = read_transient_field(m_variables[EX_NODE_BLOCK], field, nb, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, nb, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, nb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::Blob *blob, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        Ioss::Field::RoleType role = field.get_role();
        if (role == Ioss::Field::MESH) {
          if (field.get_name() == "ids") {
            // Map the local ids in this node block
            // (1...node_count) to global node ids.
            //          get_map(EX_BLOB).map_implicit_data(data, field, num_to_get, 0);
          }

          else if (field.get_name() == "connectivity") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else {
            num_to_get = Ioss::Utils::field_warning(blob, field, "input");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this blob.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // items in this blob.
          num_to_get = read_transient_field(m_variables[EX_BLOB], field, blob, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, blob, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, blob, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::Assembly *assembly, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        Ioss::Field::RoleType role = field.get_role();
        if (role == Ioss::Field::MESH) {
          if (field.get_name() == "ids") {
            // Map the local ids in this node block
            // (1...node_count) to global node ids.
            //          get_map(EX_ASSEMBLY).map_implicit_data(data, field, num_to_get, 0);
          }

          else if (field.get_name() == "connectivity") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else {
            num_to_get = Ioss::Utils::field_warning(assembly, field, "input");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this assembly.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // items in this assembly.
          num_to_get = read_transient_field(m_variables[EX_ASSEMBLY], field, assembly, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, assembly, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, assembly, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::ElementBlock *eb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        int64_t               id               = Ioex::get_id(eb, &ids_);
        size_t                my_element_count = eb->entity_count();
        Ioss::Field::RoleType role             = field.get_role();

        if (role == Ioss::Field::MESH) {
          // Handle the MESH fields required for an Exodus file model.
          // (The 'genesis' portion)

          if (field.get_name() == "connectivity") {
            int element_nodes = eb->topology()->number_nodes();
            assert(field.get_component_count(Ioss::Field::InOut::INPUT) == element_nodes);

            // The connectivity is stored in a 1D array.
            // The element_node index varies fastet
            if (my_element_count > 0) {
              get_connectivity_data(get_file_pointer(), data, EX_ELEM_BLOCK, id, 0);
              get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get * element_nodes);
            }
          }
          else if (field.get_name() == "connectivity_face") {
            // The connectivity is stored in a 1D array.
            // The element_face index varies fastest
            if (my_element_count > 0) {
              int face_count = field.get_component_count(Ioss::Field::InOut::INPUT);
              get_connectivity_data(get_file_pointer(), data, EX_ELEM_BLOCK, id, 2);
              get_map(EX_FACE_BLOCK).map_data(data, field, num_to_get * face_count);
            }
          }
          else if (field.get_name() == "connectivity_edge") {
            // The connectivity is stored in a 1D array.
            // The element_edge index varies fastest
            if (my_element_count > 0) {
              int edge_count = field.get_component_count(Ioss::Field::InOut::INPUT);
              get_connectivity_data(get_file_pointer(), data, EX_ELEM_BLOCK, id, 1);
              get_map(EX_EDGE_BLOCK).map_data(data, field, num_to_get * edge_count);
            }
          }
          else if (field.get_name() == "connectivity_raw") {
            // "connectivity_raw" has nodes in local id space (1-based)
            assert(field.get_component_count(Ioss::Field::InOut::INPUT) ==
                   eb->topology()->number_nodes());

            // The connectivity is stored in a 1D array.
            // The element_node index varies fastest
            if (my_element_count > 0) {
              get_connectivity_data(get_file_pointer(), data, EX_ELEM_BLOCK, id, 0);
            }
          }
          else if (field.get_name() == "ids") {
            // Map the local ids in this element block
            // (eb_offset+1...eb_offset+1+my_element_count) to global element ids.
            get_map(EX_ELEM_BLOCK).map_implicit_data(data, field, num_to_get, eb->get_offset());
          }
          else if (field.get_name() == "implicit_ids") {
            // If not parallel, then this is just one..element_count
            // If parallel, then it is the data in the ex_get_id_map created by nem_spread.
            size_t eb_offset_plus_one = eb->get_offset() + 1;
            if (isParallel) {
              int error = ex_get_partial_id_map(get_file_pointer(), EX_ELEM_MAP, eb_offset_plus_one,
                                                my_element_count, data);
              if (error < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
            else {
              if ((ex_int64_status(get_file_pointer()) & EX_BULK_INT64_API) != 0) {
                auto *idata = static_cast<int64_t *>(data);
                for (size_t i = 0; i < my_element_count; i++) {
                  idata[i] = eb_offset_plus_one + i;
                }
              }
              else {
                auto *idata = static_cast<int *>(data);
                for (size_t i = 0; i < my_element_count; i++) {
                  idata[i] = eb_offset_plus_one + i;
                }
              }
            }
          }
        }
        else if (role == Ioss::Field::MAP) {
          int component_count = field.get_component_count(Ioss::Field::InOut::INPUT);
          if (component_count == 1) {
            if (field.is_type(Ioss::Field::INTEGER)) {
              Ioss::IntVector component(my_element_count);
              int            *data32 = reinterpret_cast<int *>(data);

              size_t eb_offset = eb->get_offset();
              ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, field.get_index(),
                                     eb_offset + 1, my_element_count, Data(component));

              for (size_t i = 0; i < my_element_count; i++) {
                data32[i] = component[i];
              }
            }
            else {
              Ioss::Int64Vector component(my_element_count);
              auto             *data64 = reinterpret_cast<int64_t *>(data);

              size_t eb_offset = eb->get_offset();
              ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, field.get_index(),
                                     eb_offset + 1, my_element_count, Data(component));

              for (size_t i = 0; i < my_element_count; i++) {
                data64[i] = component[i];
              }
            }
          }
          else {
            // Multi-component field...
            if (field.is_type(Ioss::Field::INTEGER)) {
              Ioss::IntVector component(my_element_count);
              int            *data32    = reinterpret_cast<int *>(data);
              size_t          eb_offset = eb->get_offset();

              for (int comp = 0; comp < component_count; comp++) {
                ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, field.get_index() + comp,
                                       eb_offset + 1, my_element_count, Data(component));

                int index = comp;
                for (size_t i = 0; i < my_element_count; i++) {
                  data32[index] = component[i];
                  index += component_count;
                }
              }
            }
            else {
              Ioss::Int64Vector component(my_element_count);
              auto             *data64    = reinterpret_cast<int64_t *>(data);
              size_t            eb_offset = eb->get_offset();

              for (int comp = 0; comp < component_count; comp++) {
                ex_get_partial_num_map(get_file_pointer(), EX_ELEM_MAP, field.get_index() + comp,
                                       eb_offset + 1, my_element_count, Data(component));

                int index = comp;
                for (size_t i = 0; i < my_element_count; i++) {
                  data64[index] = component[i];
                  index += component_count;
                }
              }
            }
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, eb, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this element block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // elements in this block.
          num_to_get = read_transient_field(m_variables[EX_ELEM_BLOCK], field, eb, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, eb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::FaceBlock *eb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        int64_t               id            = Ioex::get_id(eb, &ids_);
        size_t                my_face_count = eb->entity_count();
        Ioss::Field::RoleType role          = field.get_role();

        if (role == Ioss::Field::MESH) {
          // Handle the MESH fields required for an Exodus file model.
          // (The 'genesis' portion)

          if (field.get_name() == "connectivity") {
            // The connectivity is stored in a 1D array.
            // The face_node index varies fastet
            if (my_face_count > 0) {
              int face_nodes = eb->topology()->number_nodes();
              assert(field.get_component_count(Ioss::Field::InOut::INPUT) == face_nodes);

              get_connectivity_data(get_file_pointer(), data, EX_FACE_BLOCK, id, 0);
              get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get * face_nodes);
            }
          }
          else if (field.get_name() == "connectivity_edge") {
            // The connectivity is stored in a 1D array.
            // The face_edge index varies fastest
            if (my_face_count > 0) {
              int edge_count = field.get_component_count(Ioss::Field::InOut::INPUT);
              get_connectivity_data(get_file_pointer(), data, EX_FACE_BLOCK, id, 1);
              get_map(EX_EDGE_BLOCK).map_data(data, field, num_to_get * edge_count);
            }
          }
          else if (field.get_name() == "connectivity_raw") {
            // The connectivity is stored in a 1D array.
            // The face_node index varies fastet
            if (my_face_count > 0) {
              // "connectivity_raw" has nodes in local id space (1-based)
              assert(field.get_component_count(Ioss::Field::InOut::INPUT) ==
                     eb->topology()->number_nodes());

              get_connectivity_data(get_file_pointer(), data, EX_FACE_BLOCK, id, 0);
            }
          }
          else if (field.get_name() == "ids") {
            // Map the local ids in this face block
            // (eb_offset+1...eb_offset+1+my_face_count) to global face ids.
            get_map(EX_FACE_BLOCK).map_implicit_data(data, field, num_to_get, eb->get_offset());
          }
          else {
            num_to_get = Ioss::Utils::field_warning(eb, field, "input");
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, eb, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this element block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // elements in this block.
          num_to_get = read_transient_field(m_variables[EX_FACE_BLOCK], field, eb, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, eb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::EdgeBlock *eb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        int64_t               id            = Ioex::get_id(eb, &ids_);
        int64_t               my_edge_count = eb->entity_count();
        Ioss::Field::RoleType role          = field.get_role();

        if (role == Ioss::Field::MESH) {
          // Handle the MESH fields required for an Exodus file model.
          // (The 'genesis' portion)

          if (field.get_name() == "connectivity") {
            int edge_nodes = eb->topology()->number_nodes();
            assert(field.get_component_count(Ioss::Field::InOut::INPUT) == edge_nodes);

            // The connectivity is stored in a 1D array.
            // The edge_node index varies fastet
            if (my_edge_count > 0) {
              get_connectivity_data(get_file_pointer(), data, EX_EDGE_BLOCK, id, 0);
              get_map(EX_NODE_BLOCK).map_data(data, field, num_to_get * edge_nodes);
            }
          }
          else if (field.get_name() == "connectivity_raw") {
            // "connectivity_raw" has nodes in local id space (1-based)
            assert(field.get_component_count(Ioss::Field::InOut::INPUT) ==
                   eb->topology()->number_nodes());

            // The connectivity is stored in a 1D array.
            // The edge_node index varies fastet
            if (my_edge_count > 0) {
              get_connectivity_data(get_file_pointer(), data, EX_EDGE_BLOCK, id, 0);
            }
          }
          else if (field.get_name() == "ids") {
            // Map the local ids in this edge block
            // (eb_offset+1...eb_offset+1+my_edge_count) to global edge ids.
            get_map(EX_EDGE_BLOCK).map_implicit_data(data, field, num_to_get, eb->get_offset());
          }
          else {
            num_to_get = Ioss::Utils::field_warning(eb, field, "input");
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, eb, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this element block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // elements in this block.
          num_to_get = read_transient_field(m_variables[EX_EDGE_BLOCK], field, eb, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, eb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_Xset_field_internal(const Ioss::EntitySet *ns, const Ioss::Field &field,
                                              void *data, size_t data_size) const
  {
    {
      int               ierr;
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        ex_entity_type        type = Ioex::map_exodus_type(ns->type());
        int64_t               id   = Ioex::get_id(ns, &ids_);
        Ioss::Field::RoleType role = field.get_role();
        if (role == Ioss::Field::MESH) {

          if (field.get_name() == "ids" || field.get_name() == "ids_raw") {
            if (field.get_type() == Ioss::Field::INTEGER) {
              ierr = ex_get_set(get_file_pointer(), type, id, static_cast<int *>(data), nullptr);
            }
            else {
              ierr =
                  ex_get_set(get_file_pointer(), type, id, static_cast<int64_t *>(data), nullptr);
            }
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }

            if (field.get_name() == "ids") {
              // Convert the local node ids to global ids
              get_map(type).map_data(data, field, num_to_get);
            }
          }
          else if (field.get_name() == "orientation") {
            if (field.get_type() == Ioss::Field::INTEGER) {
              ierr = ex_get_set(get_file_pointer(), type, id, nullptr, static_cast<int *>(data));
            }
            else {
              ierr =
                  ex_get_set(get_file_pointer(), type, id, nullptr, static_cast<int64_t *>(data));
            }
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else if (field.get_name() == "distribution_factors") {
            ex_set set_param[1];
            set_param[0].id                       = id;
            set_param[0].type                     = type;
            set_param[0].entry_list               = nullptr;
            set_param[0].extra_list               = nullptr;
            set_param[0].distribution_factor_list = nullptr;
            ierr                                  = ex_get_sets(get_file_pointer(), 1, set_param);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }

            if (set_param[0].num_distribution_factor == 0) {
              auto *rdata = static_cast<double *>(data);
              for (size_t i = 0; i < num_to_get; i++) {
                rdata[i] = 1.0;
              }
            }
            else {
              set_param[0].distribution_factor_list = static_cast<double *>(data);
              ierr                                  = ex_get_sets(get_file_pointer(), 1, set_param);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else {
            num_to_get = Ioss::Utils::field_warning(ns, field, "input");
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = read_attribute_field(field, ns, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          get_reduction_field(field, ns, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this node block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Read in each component of the variable and transfer into
          // 'data'.  Need temporary storage area of size 'number of
          // nodes in this block.
          num_to_get = read_transient_field(m_variables[type], field, ns, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::NodeSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    if (!ns->property_exists("filtered_db_set_size")) {
      return get_Xset_field_internal(ns, field, data, data_size);
    }

    size_t db_size = ns->get_property("filtered_db_set_size").get_int();

    int               ierr;
    Ioss::SerializeIO serializeio_(this);

    size_t num_to_get = field.verify(data_size);
    if (num_to_get > 0) {

      int64_t               id   = Ioex::get_id(ns, &ids_);
      Ioss::Field::RoleType role = field.get_role();
      if (role == Ioss::Field::MESH) {

        if (field.get_name() == "ids" || field.get_name() == "ids_raw") {
          if (field.get_type() == Ioss::Field::INTEGER) {
            Ioss::IntVector dbvals(db_size);
            ierr = ex_get_set(get_file_pointer(), EX_NODE_SET, id, Data(dbvals), nullptr);
            if (ierr >= 0) {
              Ioex::filter_node_list(static_cast<int *>(data), dbvals,
                                     activeNodeSetNodesIndex[ns->name()]);
            }
          }
          else {
            Ioss::Int64Vector dbvals(db_size);
            ierr = ex_get_set(get_file_pointer(), EX_NODE_SET, id, Data(dbvals), nullptr);
            if (ierr >= 0) {
              Ioex::filter_node_list(static_cast<int64_t *>(data), dbvals,
                                     activeNodeSetNodesIndex[ns->name()]);
            }
          }
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }

          if (field.get_name() == "ids") {
            // Convert the local node ids to global ids
            get_map(EX_NODE_SET).map_data(data, field, num_to_get);
          }
        }
        else if (field.get_name() == "distribution_factors") {
          ex_set set_param[1];
          set_param[0].id                       = id;
          set_param[0].type                     = EX_NODE_SET;
          set_param[0].entry_list               = nullptr;
          set_param[0].extra_list               = nullptr;
          set_param[0].distribution_factor_list = nullptr;
          ierr                                  = ex_get_sets(get_file_pointer(), 1, set_param);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }

          if (set_param[0].num_distribution_factor == 0) {
            auto *rdata = static_cast<double *>(data);
            for (size_t i = 0; i < num_to_get; i++) {
              rdata[i] = 1.0;
            }
          }
          else {
            std::vector<double> dbvals(db_size);
            set_param[0].distribution_factor_list = Data(dbvals);
            ierr                                  = ex_get_sets(get_file_pointer(), 1, set_param);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
            Ioex::filter_node_list(static_cast<double *>(data), dbvals,
                                   activeNodeSetNodesIndex[ns->name()]);
            set_param[0].distribution_factor_list = nullptr;
          }
        }
        else {
          num_to_get = Ioss::Utils::field_warning(ns, field, "input");
        }
      }
      else if (role == Ioss::Field::ATTRIBUTE) {
        num_to_get = Ioss::Utils::field_warning(ns, field, "input");
      }
      else if (role == Ioss::Field::REDUCTION) {
        num_to_get = Ioss::Utils::field_warning(ns, field, "input");
      }
      else if (role == Ioss::Field::TRANSIENT) {
        // Filtered not currently implemented for transient or attributes....
      }
    }
    return num_to_get;
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::EdgeSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return get_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::FaceSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return get_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::ElementSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return get_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::SideSet *fs, const Ioss::Field &field,
                                         void * /* data */, size_t data_size) const
  {
    size_t num_to_get = field.verify(data_size);
    if (field.get_name() == "ids") {
      // Do nothing, just handles an idiosyncrasy of the GroupingEntity
    }
    else {
      num_to_get = Ioss::Utils::field_warning(fs, field, "input");
    }
    return num_to_get;
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::CommSet *cs, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);

      if (num_to_get > 0) {
        int64_t entity_count = cs->entity_count();

        // Return the <entity (node or side), processor> pair
        if (field.get_name() == "entity_processor" || field.get_name() == "entity_processor_raw") {

          // Check type -- node or side
          std::string type = cs->get_property("entity_type").get_string();

          // Allocate temporary storage space
          std::vector<char> entities(num_to_get * int_byte_size_api());
          std::vector<char> procs(num_to_get * int_byte_size_api());

          if (type == "node") {
            int64_t cm_offset = 0;

            for (int64_t i = 0; i < commsetNodeCount; i++) {
              int ierr = ex_get_node_cmap(get_file_pointer(), nodeCmapIds[i], &entities[cm_offset],
                                          &procs[cm_offset], myProcessor);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
              cm_offset += (nodeCmapNodeCnts[i] * int_byte_size_api());
            }
            assert(cm_offset == entity_count * int_byte_size_api());

            // Convert local node id to global node id and store in 'data'
            if (int_byte_size_api() == 4) {
              int *entity_proc = static_cast<int *>(data);
              int *ents        = reinterpret_cast<int *>(Data(entities));
              int *pros        = reinterpret_cast<int *>(Data(procs));

              size_t j = 0;
              if (field.get_name() == "entity_processor") {
                const Ioss::MapContainer &map = get_map(EX_NODE_BLOCK).map();

                for (int64_t i = 0; i < entity_count; i++) {
                  int local_id     = ents[i];
                  entity_proc[j++] = map[local_id];
                  entity_proc[j++] = pros[i];
                }
              }
              else { // "entity_processor_raw"
                for (int64_t i = 0; i < entity_count; i++) {
                  entity_proc[j++] = ents[i];
                  entity_proc[j++] = pros[i];
                }
              }
            }
            else {
              auto *entity_proc = static_cast<int64_t *>(data);
              auto *ents        = reinterpret_cast<int64_t *>(Data(entities));
              auto *pros        = reinterpret_cast<int64_t *>(Data(procs));

              size_t j = 0;
              if (field.get_name() == "entity_processor") {
                const Ioss::MapContainer &map = get_map(EX_NODE_BLOCK).map();

                for (int64_t i = 0; i < entity_count; i++) {
                  int64_t local_id = ents[i];
                  entity_proc[j++] = map[local_id];
                  entity_proc[j++] = pros[i];
                }
              }
              else { // "entity_processor_raw"
                for (int64_t i = 0; i < entity_count; i++) {
                  entity_proc[j++] = ents[i];
                  entity_proc[j++] = pros[i];
                }
              }
            }
          }
          else if (type == "side") {
            std::vector<char> sides(entity_count * int_byte_size_api());
            int64_t           cm_offset = 0;
            for (int64_t i = 0; i < commsetElemCount; i++) {
              int ierr = ex_get_elem_cmap(get_file_pointer(), elemCmapIds[i], &entities[cm_offset],
                                          &sides[cm_offset], &procs[cm_offset], myProcessor);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
              cm_offset += (elemCmapElemCnts[i] * int_byte_size_api());
            }
            assert(cm_offset == entity_count * int_byte_size_api());

            if (int_byte_size_api() == 4) {
              int *entity_proc = static_cast<int *>(data);
              int *ents        = reinterpret_cast<int *>(Data(entities));
              int *pros        = reinterpret_cast<int *>(Data(procs));
              int *sids        = reinterpret_cast<int *>(Data(sides));

              size_t j = 0;
              if (field.get_name() == "entity_processor") {
                const Ioss::MapContainer &map = get_map(EX_ELEM_BLOCK).map();

                for (int64_t i = 0; i < entity_count; i++) {
                  entity_proc[j++] = map[ents[i]];
                  entity_proc[j++] = sids[i];
                  entity_proc[j++] = pros[i];
                }
              }
              else { // "entity_processor_raw"
                for (int64_t i = 0; i < entity_count; i++) {
                  entity_proc[j++] = ents[i];
                  entity_proc[j++] = sids[i];
                  entity_proc[j++] = pros[i];
                }
              }
            }
            else {
              auto *entity_proc = static_cast<int64_t *>(data);
              auto *ents        = reinterpret_cast<int64_t *>(Data(entities));
              auto *pros        = reinterpret_cast<int64_t *>(Data(procs));
              auto *sids        = reinterpret_cast<int64_t *>(Data(sides));

              size_t j = 0;
              if (field.get_name() == "entity_processor") {
                const Ioss::MapContainer &map = get_map(EX_ELEM_BLOCK).map();

                for (int64_t i = 0; i < entity_count; i++) {
                  entity_proc[j++] = map[ents[i]];
                  entity_proc[j++] = sids[i];
                  entity_proc[j++] = pros[i];
                }
              }
              else { // "entity_processor_raw"
                for (int64_t i = 0; i < entity_count; i++) {
                  entity_proc[j++] = ents[i];
                  entity_proc[j++] = sids[i];
                  entity_proc[j++] = pros[i];
                }
              }
            }
          }
          else {
            std::ostringstream errmsg;
            fmt::print(errmsg, "ERROR: Invalid commset type {}", type);
            IOSS_ERROR(errmsg);
          }
        }
        else if (field.get_name() == "ids") {
          // Do nothing, just handles an idiosyncrasy of the GroupingEntity
        }
        else {
          num_to_get = Ioss::Utils::field_warning(cs, field, "input");
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::get_field_internal(const Ioss::SideBlock *fb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    Ioss::SerializeIO serializeio_(this);
    int64_t           num_to_get = field.verify(data_size);
    if (num_to_get > 0) {

      int64_t id           = Ioex::get_id(fb, &ids_);
      int64_t entity_count = fb->entity_count();
      if (num_to_get != entity_count) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Partial field input not yet implemented for side blocks");
        IOSS_ERROR(errmsg);
      }

      ex_set set_param[1];
      set_param[0].id                       = id;
      set_param[0].type                     = EX_SIDE_SET;
      set_param[0].entry_list               = nullptr;
      set_param[0].extra_list               = nullptr;
      set_param[0].distribution_factor_list = nullptr;
      int ierr                              = ex_get_sets(get_file_pointer(), 1, set_param);
      if (ierr < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }

      int64_t number_sides                = set_param[0].num_entry;
      int64_t number_distribution_factors = set_param[0].num_distribution_factor;

      Ioss::Field::RoleType role = field.get_role();
      if (role == Ioss::Field::MESH) {

        // In exodus, we may have split the sideset into multiple
        // side blocks if there are multiple side topologies in the
        // sideset.  Because of this, the passed in 'data' may not be
        // large enough to hold the data residing in the sideset and we
        // may need to allocate a temporary array...  This can be checked
        // by comparing the size of the sideset with the 'side_count' of
        // the side block.

        // Get size of data stored on the file...
        // FIX 64: FIX THIS -- STORING INT IN DOUBLE WON'T WORK
        if (field.get_name() == "side_ids" && fb->name() == "universal_sideset") {
          // The side ids are being stored as the distribution factor
          // field on the universal sideset.  There should be no other
          // side sets that request this field...  (Eventually,
          // create an id field to store this info.

          if (number_distribution_factors == num_to_get) {
            std::vector<double> real_ids(num_to_get);
            set_param[0].distribution_factor_list = Data(real_ids);
            ierr                                  = ex_get_sets(get_file_pointer(), 1, set_param);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
            set_param[0].distribution_factor_list = nullptr;

            if (field.get_type() == Ioss::Field::INTEGER) {
              // Need to convert 'double' to 'int' for Sierra use...
              int *ids = static_cast<int *>(data);
              for (int64_t i = 0; i < num_to_get; i++) {
                ids[i] = static_cast<int>(real_ids[i]);
              }
            }
            else {
              // Need to convert 'double' to 'int' for Sierra use...
              auto *ids = static_cast<int64_t *>(data);
              for (int64_t i = 0; i < num_to_get; i++) {
                ids[i] = static_cast<int64_t>(real_ids[i]);
              }
            }
          }
        }

        else if (field.get_name() == "side_ids") {
        }

        else if (field.get_name() == "ids") {
          // In exodus, the 'side set' is stored as a sideset.  A
          // sideset has a list of elements and a corresponding local
          // element side (1-based) The side id is: side_id =
          // 10*element_id + local_side_number This assumes that all
          // sides in a sideset are boundary sides.  Since we
          // only have a single array, we need to allocate an extra array
          // to store all of the data.  Note also that the element_id is
          // the global id but only the local id is stored so we need to
          // map from local_to_global prior to generating the side id...

          Ioss::Field       el_side = fb->get_field("element_side");
          std::vector<char> element_side(2 * number_sides * int_byte_size_api());
          get_field_internal(fb, el_side, Data(element_side), element_side.size());

          // At this point, have the 'element_side' data containing
          // the global element ids and the sides...  Iterate
          // through to generate the ids...
          if (int_byte_size_api() == 4) {
            int64_t int_max = std::numeric_limits<int>::max();
            int    *ids     = static_cast<int *>(data);
            int    *els     = reinterpret_cast<int *>(Data(element_side));
            size_t  idx     = 0;
            for (int64_t iel = 0; iel < 2 * entity_count; iel += 2) {
              int64_t new_id = static_cast<int64_t>(10) * els[iel] + els[iel + 1];
              if (new_id > int_max) {
                std::ostringstream errmsg;
                fmt::print(errmsg,
                           "ERROR: Process {} accessing the sideset field 'ids'\n"
                           "\t\thas exceeded the integer bounds for entity {}, local side id {}"
                           ".\n\t\tTry using 64-bit mode to read the file '{}'.\n",
                           Ioss::SerializeIO::getRank(), els[iel], els[iel + 1],
                           decoded_filename());
                IOSS_ERROR(errmsg);
              }

              ids[idx++] = static_cast<int>(new_id);
            }
          }
          else {
            auto  *ids = static_cast<int64_t *>(data);
            auto  *els = reinterpret_cast<int64_t *>(Data(element_side));
            size_t idx = 0;
            for (int64_t iel = 0; iel < 2 * entity_count; iel += 2) {
              int64_t new_id = 10 * els[iel] + els[iel + 1];
              ids[idx++]     = new_id;
            }
          }
        }
        else if (field.get_name() == "element_side" || field.get_name() == "element_side_raw") {
          // In exodus, the 'side set' is stored as a sideset.  A sideset
          // has a list of elements and a corresponding local element side
          // (1-based)

          // Since we only have a single array, we need to allocate an extra
          // array to store all of the data.  Note also that the element_id
          // is the global id but only the local id is stored so we need to
          // map from local_to_global prior to generating the side  id...

          // Get the element number map (1-based)...
          // Not used for "element_side_raw", but reduces code duplication...
          // Although combining code for raw and regular introduces an if in the
          // middle of a loop, we rely on the branch predictor to learn that
          // `do_map` does not change during the iteration...
          const Ioss::MapContainer &map = get_map(EX_ELEM_BLOCK).map();

          bool do_map = field.get_name() == "element_side";

          // Allocate space for local side number and element numbers
          // numbers.

          // See if edges or faces...
          int64_t side_offset = Ioss::Utils::get_side_offset(fb);

          std::vector<char> element(number_sides * int_byte_size_api());
          std::vector<char> sides(number_sides * int_byte_size_api());

          ierr = ex_get_set(get_file_pointer(), EX_SIDE_SET, id, Data(element), Data(sides));
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }

          if (number_sides == entity_count) {
            int64_t index = 0;
            if (int_byte_size_api() == 4) {
              auto *element_side = static_cast<int *>(data);
              auto *element32    = reinterpret_cast<int *>(Data(element));
              auto *sides32      = reinterpret_cast<int *>(Data(sides));
              for (int64_t iel = 0; iel < entity_count; iel++) {
                if (do_map) {
                  element_side[index++] = map[element32[iel]];
                }
                else {
                  element_side[index++] = element32[iel];
                }
                element_side[index++] = sides32[iel] - side_offset;
              }
            }
            else {
              auto *element_side = static_cast<int64_t *>(data);
              auto *element64    = reinterpret_cast<int64_t *>(Data(element));
              auto *sides64      = reinterpret_cast<int64_t *>(Data(sides));
              for (int64_t iel = 0; iel < entity_count; iel++) {
                if (do_map) {
                  element_side[index++] = map[element64[iel]];
                }
                else {
                  element_side[index++] = element64[iel];
                }
                element_side[index++] = sides64[iel] - side_offset;
              }
            }
            assert(index / 2 == entity_count);
          }
          else {
            Ioss::IntVector is_valid_side;
            Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(),
                                                        Data(element), Data(sides), number_sides,
                                                        get_region());

            int64_t index = 0;
            if (int_byte_size_api() == 4) {
              auto *element_side = static_cast<int *>(data);
              auto *element32    = reinterpret_cast<int *>(Data(element));
              auto *sides32      = reinterpret_cast<int *>(Data(sides));
              for (int64_t iel = 0; iel < number_sides; iel++) {
                if (is_valid_side[iel] == 1) {
                  // This side  belongs in the side block
                  if (do_map) {
                    element_side[index++] = map[element32[iel]];
                  }
                  else {
                    element_side[index++] = element32[iel];
                  }
                  element_side[index++] = sides32[iel] - side_offset;
                }
              }
            }
            else {
              auto *element_side = static_cast<int64_t *>(data);
              auto *element64    = reinterpret_cast<int64_t *>(Data(element));
              auto *sides64      = reinterpret_cast<int64_t *>(Data(sides));
              for (int64_t iel = 0; iel < number_sides; iel++) {
                if (is_valid_side[iel] == 1) {
                  // This side  belongs in the side block
                  if (do_map) {
                    element_side[index++] = map[element64[iel]];
                  }
                  else {
                    element_side[index++] = element64[iel];
                  }
                  element_side[index++] = sides64[iel] - side_offset;
                }
              }
            }
            assert(index / 2 == entity_count);
          }
        }
        else if (field.get_name() == "connectivity") {
          // The side  connectivity needs to be generated 'on-the-fly' from
          // the element number and local side  of that element. A sideset
          // can span multiple element blocks, and contain multiple side
          // types; the side block contains side of similar topology.
          ierr = get_side_connectivity(fb, id, entity_count, data, true);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
        else if (field.get_name() == "connectivity_raw") {
          // The side  connectivity needs to be generated 'on-the-fly' from
          // the element number and local side  of that element. A sideset
          // can span multiple element blocks, and contain multiple side
          // types; the side block contains side of similar topology.
          ierr = get_side_connectivity(fb, id, entity_count, data, false);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
        else if (field.get_name() == "distribution_factors") {
          ierr = get_side_distributions(fb, id, entity_count, static_cast<double *>(data),
                                        data_size / sizeof(double));
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
        else {
          num_to_get = Ioss::Utils::field_warning(fb, field, "input");
        }
      }
      else if (role == Ioss::Field::TRANSIENT) {
        // Check if the specified field exists on this block.
        // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
        // exist on the database as scalars with the appropriate
        // extensions.

        if (number_sides == entity_count) {
          num_to_get = read_transient_field(m_variables[EX_SIDE_SET], field, fb, data);
        }
        else {
          // Need to read all values for the specified field and then
          // filter down to the elements actually in this side block.

          // Determine which sides are member of this block
          Ioss::IntVector is_valid_side;
          {
            //----
            std::vector<char> element(number_sides * int_byte_size_api());
            std::vector<char> sides(number_sides * int_byte_size_api());
            ierr = ex_get_set(get_file_pointer(), EX_SIDE_SET, id, Data(element), Data(sides));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
            //----
            Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(),
                                                        Data(element), Data(sides), number_sides,
                                                        get_region());
          }

          num_to_get = read_ss_transient_field(field, id, data, is_valid_side);
        }
      }
    }
    return num_to_get;
  }

  int64_t DatabaseIO::write_attribute_field(const Ioss::Field          &field,
                                            const Ioss::GroupingEntity *ge, void *data) const
  {
    std::string att_name   = ge->name() + SEP() + field.get_name();
    int64_t     num_entity = ge->entity_count();
    int64_t     fld_offset = field.get_index();

    ex_entity_type type            = Ioex::map_exodus_type(ge->type());
    int64_t        id              = Ioex::get_id(ge, &ids_);
    int            attribute_count = ge->get_property("attribute_count").get_int();
    assert(fld_offset > 0);
    assert(fld_offset - 1 + field.get_component_count(Ioss::Field::InOut::OUTPUT) <=
           attribute_count);

    Ioss::Field::BasicType ioss_type = field.get_type();
    assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER ||
           ioss_type == Ioss::Field::INT64);

    if (ioss_type == Ioss::Field::INT64) {
      Ioss::Utils::check_int_to_real_overflow(field, (int64_t *)data, num_entity);
    }

    if (fld_offset == 1 &&
        field.get_component_count(Ioss::Field::InOut::OUTPUT) == attribute_count) {
      // Write all attributes in one big chunk...
      std::vector<double> temp;
      double             *rdata = nullptr;
      if (ioss_type == Ioss::Field::INTEGER) {
        int *idata = static_cast<int *>(data);
        extract_data(temp, idata, attribute_count * num_entity, 1, 0);
        rdata = Data(temp);
      }
      else if (ioss_type == Ioss::Field::INT64) {
        auto *idata = static_cast<int64_t *>(data);
        extract_data(temp, idata, attribute_count * num_entity, 1, 0);
        rdata = Data(temp);
      }
      else {
        rdata = static_cast<double *>(data);
      }

      int ierr = ex_put_attr(get_file_pointer(), type, id, rdata);
      if (ierr < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }
    }
    else {
      // Write a subset of the attributes.  If scalar, write one;
      // if higher-order (vector3d, ..) write each component.
      if (field.get_component_count(Ioss::Field::InOut::OUTPUT) == 1) {
        std::vector<double> temp;
        double             *rdata = nullptr;
        if (ioss_type == Ioss::Field::INTEGER) {
          int *idata = static_cast<int *>(data);
          extract_data(temp, idata, num_entity, 1, 0);
          rdata = Data(temp);
        }
        else if (ioss_type == Ioss::Field::INT64) {
          auto *idata = static_cast<int64_t *>(data);
          extract_data(temp, idata, num_entity, 1, 0);
          rdata = Data(temp);
        }
        else {
          rdata = static_cast<double *>(data);
        }

        int ierr = ex_put_one_attr(get_file_pointer(), type, id, fld_offset, rdata);
        if (ierr < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }
      }
      else {
        // Multi-component...  Need a local memory space to push
        // data into and then write that out to the file...
        std::vector<double> local_data(num_entity);
        int                 comp_count = field.get_component_count(Ioss::Field::InOut::OUTPUT);
        for (int i = 0; i < comp_count; i++) {
          size_t offset = i;
          if (ioss_type == Ioss::Field::REAL) {
            auto *rdata = static_cast<double *>(data);
            extract_data(local_data, rdata, num_entity, comp_count, offset);
          }
          else if (ioss_type == Ioss::Field::INTEGER) {
            int *idata = static_cast<int *>(data);
            extract_data(local_data, idata, num_entity, comp_count, offset);
          }
          else if (ioss_type == Ioss::Field::INT64) {
            auto *idata = static_cast<int64_t *>(data);
            extract_data(local_data, idata, num_entity, comp_count, offset);
          }

          int ierr =
              ex_put_one_attr(get_file_pointer(), type, id, fld_offset + i, Data(local_data));
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
      }
    }
    return num_entity;
  }

  int64_t DatabaseIO::read_attribute_field(const Ioss::Field &field, const Ioss::GroupingEntity *ge,
                                           void *data) const
  {
    // TODO: Handle INTEGER fields...

    int64_t num_entity = ge->entity_count();
    if (num_entity == 0) {
      return 0;
    }

    Ioss::Field::BasicType ioss_type = field.get_type();
    if (ioss_type == Ioss::Field::INTEGER || ioss_type == Ioss::Field::INT64) {
      std::ostringstream errmsg;
      fmt::print(errmsg, "INTERNAL ERROR: Integer attribute fields are not yet handled for read. "
                         "Please report.\n");
      IOSS_ERROR(errmsg);
    }

    int            attribute_count = ge->get_property("attribute_count").get_int();
    ex_entity_type type            = Ioex::map_exodus_type(ge->type());
    int64_t        id              = Ioex::get_id(ge, &ids_);

    std::string att_name = ge->name() + SEP() + field.get_name();
    int64_t     offset   = field.get_index();
    assert(offset - 1 + field.get_component_count(Ioss::Field::InOut::INPUT) <= attribute_count);
    if (offset == 1 && field.get_component_count(Ioss::Field::InOut::INPUT) == attribute_count) {
      // Read all attributes in one big chunk...
      int ierr = ex_get_attr(get_file_pointer(), type, id, static_cast<double *>(data));
      if (ierr < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }
    }
    else {
      // Read a subset of the attributes.  If scalar, read one;
      // if higher-order (vector3d, ..) read each component and
      // put into correct location...
      if (field.get_component_count(Ioss::Field::InOut::INPUT) == 1) {
        int ierr =
            ex_get_one_attr(get_file_pointer(), type, id, offset, static_cast<double *>(data));
        if (ierr < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }
      }
      else {
        // Multi-component...
        // Need a local memory space to read data into and
        // then push that into the user-supplied data block...
        std::vector<double> local_data(num_entity);
        int                 comp_count = field.get_component_count(Ioss::Field::InOut::INPUT);
        auto               *rdata      = static_cast<double *>(data);
        for (int i = 0; i < comp_count; i++) {
          int ierr = ex_get_one_attr(get_file_pointer(), type, id, offset + i, Data(local_data));
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }

          size_t k = i;
          for (int64_t j = 0; j < num_entity; j++) {
            rdata[k] = local_data[j];
            k += comp_count;
          }
        }
      }
    }
    return num_entity;
  }

  int64_t DatabaseIO::read_transient_field(const Ioex::VariableNameMap &variables,
                                           const Ioss::Field &field, const Ioss::GroupingEntity *ge,
                                           void *data) const
  {
    // Read into a double variable since that is all Exodus can store...
    size_t              num_entity = ge->entity_count();
    std::vector<double> temp(num_entity);

    size_t step = get_current_state();

    // get number of components, cycle through each component
    // and add suffix to base 'field_name'.  Look up index
    // of this name in 'nodeVariables' map
    size_t comp_count = field.get_component_count(Ioss::Field::InOut::INPUT);

    ex_entity_type type = Ioex::map_exodus_type(ge->type());
    if (comp_count == 1 && field.get_type() == Ioss::Field::REAL) {
      std::string var_name = get_component_name(field, Ioss::Field::InOut::INPUT, 1);
      if (lowerCaseVariableNames) {
        Ioss::Utils::fixup_name(var_name);
      }

      // Read the variable...
      int64_t id   = Ioex::get_id(ge, &ids_);
      int     ierr = 0;

      auto var_iter = variables.find(var_name);
      if (var_iter == variables.end()) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Could not find field '{}'\n", var_name);
        IOSS_ERROR(errmsg);
      }
      size_t var_index = var_iter->second;
      assert(var_index > 0);
      ierr = ex_get_var(get_file_pointer(), step, type, var_index, id, num_entity, data);
      if (ierr < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }
    }
    else {
      for (size_t i = 0; i < comp_count; i++) {
        std::string var_name = get_component_name(field, Ioss::Field::InOut::INPUT, i + 1);
        if (lowerCaseVariableNames) {
          Ioss::Utils::fixup_name(var_name);
        }

        // Read the variable...
        int64_t id       = Ioex::get_id(ge, &ids_);
        int     ierr     = 0;
        auto    var_iter = variables.find(var_name);
        if (var_iter == variables.end()) {
          std::ostringstream errmsg;
          fmt::print(errmsg, "ERROR: Could not find field '{}'\n", var_name);
          IOSS_ERROR(errmsg);
        }
        size_t var_index = var_iter->second;
        assert(var_index > 0);
        ierr = ex_get_var(get_file_pointer(), step, type, var_index, id, num_entity, Data(temp));
        if (ierr < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        // Transfer to 'data' array.
        size_t k = 0;
        if (field.get_type() == Ioss::Field::INTEGER) {
          int *ivar = static_cast<int *>(data);
          for (size_t j = i; j < num_entity * comp_count; j += comp_count) {
            ivar[j] = static_cast<int>(temp[k++]);
          }
        }
        else if (field.get_type() == Ioss::Field::INT64) { // FIX 64 UNSAFE
          auto *ivar = static_cast<int64_t *>(data);
          for (size_t j = i; j < num_entity * comp_count; j += comp_count) {
            ivar[j] = static_cast<int64_t>(temp[k++]);
          }
        }
        else if (field.get_type() == Ioss::Field::REAL) {
          auto *rvar = static_cast<double *>(data);
          for (size_t j = i; j < num_entity * comp_count; j += comp_count) {
            rvar[j] = temp[k++];
          }
        }
        else {
          std::ostringstream errmsg;
          fmt::print(errmsg,
                     "IOSS_ERROR: Field storage type must be either integer or double.\n"
                     "       Field '{}' is invalid.\n",
                     field.get_name());
          IOSS_ERROR(errmsg);
        }
        assert(k == num_entity);
      }
    }
    return num_entity;
  }

  int64_t DatabaseIO::read_ss_transient_field(const Ioss::Field &field, int64_t id, void *variables,
                                              Ioss::IntVector &is_valid_side) const
  {
    size_t              num_valid_sides = 0;
    size_t              my_side_count   = is_valid_side.size();
    std::vector<double> temp(my_side_count);

    size_t step = get_current_state();

    // get number of components, cycle through each component
    // and add suffix to base 'field_name'.  Look up index
    // of this name in 'nodeVariables' map
    size_t comp_count = field.get_component_count(Ioss::Field::InOut::INPUT);

    for (size_t i = 0; i < comp_count; i++) {
      std::string var_name = get_component_name(field, Ioss::Field::InOut::INPUT, i + 1);
      if (lowerCaseVariableNames) {
        Ioss::Utils::fixup_name(var_name);
      }

      // Read the variable...
      int  ierr     = 0;
      auto var_iter = m_variables[EX_SIDE_SET].find(var_name);
      if (var_iter == m_variables[EX_SIDE_SET].end()) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Could not find Sideset field '{}'\n", var_name);
        IOSS_ERROR(errmsg);
      }
      size_t var_index = var_iter->second;
      assert(var_index > 0);
      ierr = ex_get_var(get_file_pointer(), step, EX_SIDE_SET, var_index, id, my_side_count,
                        Data(temp));
      if (ierr < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }

      // Transfer to 'variables' array.
      size_t j = i;
      if (field.get_type() == Ioss::Field::INTEGER) {
        int *ivar = static_cast<int *>(variables);
        for (size_t k = 0; k < my_side_count; k++) {
          if (is_valid_side[k] == 1) {
            ivar[j] = static_cast<int>(temp[k]);
            j += comp_count;
          }
        }
      }
      else if (field.get_type() == Ioss::Field::INT64) { // FIX 64 UNSAFE
        auto *ivar = static_cast<int64_t *>(variables);
        for (size_t k = 0; k < my_side_count; k++) {
          if (is_valid_side[k] == 1) {
            ivar[j] = static_cast<int64_t>(temp[k]);
            j += comp_count;
          }
        }
      }
      else if (field.get_type() == Ioss::Field::REAL) {
        auto *rvar = static_cast<double *>(variables);
        for (size_t k = 0; k < my_side_count; k++) {
          if (is_valid_side[k] == 1) {
            rvar[j] = temp[k];
            j += comp_count;
          }
        }
      }
      else {
        std::ostringstream errmsg;
        fmt::print(errmsg,
                   "IOSS_ERROR: Field storage type must be either integer or double.\n"
                   "       Field '{}' is invalid.\n",
                   field.get_name());
        IOSS_ERROR(errmsg);
      }
      if (i + 1 == comp_count) {
        num_valid_sides = j / comp_count;
      }
    }
    return num_valid_sides;
  }

  template <typename INT>
  int64_t DatabaseIO::get_side_connectivity_internal(const Ioss::SideBlock *fb, int64_t id,
                                                     int64_t /*unused*/, INT           *fconnect,
                                                     bool map_ids) const
  {
    // Get size of data stored on the file...
    ex_set set_param[1];
    set_param[0].id                       = id;
    set_param[0].type                     = EX_SIDE_SET;
    set_param[0].entry_list               = nullptr;
    set_param[0].extra_list               = nullptr;
    set_param[0].distribution_factor_list = nullptr;
    int ierr                              = ex_get_sets(get_file_pointer(), 1, set_param);
    if (ierr < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }

    int64_t number_sides = set_param[0].num_entry;

    // Allocate space for element and local side number
    assert(number_sides > 0);
    //----
    std::vector<INT> element(number_sides);
    std::vector<INT> side(number_sides);

    set_param[0].entry_list = Data(element);
    set_param[0].extra_list = Data(side);
    ierr                    = ex_get_sets(get_file_pointer(), 1, set_param);
    if (ierr < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }
    //----

    Ioss::IntVector is_valid_side;
    Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(),
                                                (void *)Data(element), (void *)Data(side),
                                                number_sides, get_region());

    std::vector<INT>    elconnect;
    int64_t             elconsize  = 0;       // Size of currently allocated connectivity block
    Ioss::ElementBlock *conn_block = nullptr; // Block that we currently
    // have connectivity for

    Ioss::ElementBlock *block = nullptr;
    Ioss::IntVector     side_elem_map; // Maps the side into the elements

    // connectivity array
    int64_t current_side = -1;
    int     nelnode      = 0;
    int     nfnodes      = 0;
    int     ieb          = 0;
    size_t  offset       = 0;
    for (int64_t iel = 0; iel < number_sides; iel++) {
      if (is_valid_side[iel] == 1) {

        int64_t elem_id = element[iel];

        // ensure we have correct connectivity
        block = get_region()->get_element_block(elem_id);
        assert(block != nullptr);
        if (conn_block != block) {
          int64_t nelem = block->entity_count();
          nelnode       = block->topology()->number_nodes();
          // Used to map element number into position in connectivity array.
          // E.g., element 97 is the (97-offset)th element in this block and
          // is stored in array index (97-offset-1).
          offset = block->get_offset() + 1;
          if (elconsize < nelem * nelnode) {
            elconsize = nelem * nelnode;
            elconnect.resize(elconsize);
          }
          if (map_ids) {
            get_field_internal(block, block->get_field("connectivity"), Data(elconnect),
                               nelem * nelnode * int_byte_size_api());
          }
          else {
            get_field_internal(block, block->get_field("connectivity_raw"), Data(elconnect),
                               nelem * nelnode * int_byte_size_api());
          }
          conn_block   = block;
          current_side = -1;
        }

        // NOTE: Element connectivity is returned with nodes in global id space if "map_ids" false,
        //       otherwise it is in local space.
        int64_t side_id = side[iel];

        if (current_side != side_id) {
          side_elem_map = block->topology()->boundary_connectivity(side_id);
          current_side  = side_id;
          nfnodes       = block->topology()->boundary_type(side_id)->number_nodes();
        }
        for (int inode = 0; inode < nfnodes; inode++) {
          size_t index    = (elem_id - offset) * nelnode + side_elem_map[inode];
          fconnect[ieb++] = elconnect[index];
        }
      }
    }
    return ierr;
  }

  int64_t DatabaseIO::get_side_connectivity(const Ioss::SideBlock *fb, int64_t id,
                                            int64_t my_side_count, void *fconnect,
                                            bool map_ids) const
  {
    if (int_byte_size_api() == 4) {
      return get_side_connectivity_internal(fb, id, my_side_count,
                                            reinterpret_cast<int *>(fconnect), map_ids);
    }
    return get_side_connectivity_internal(fb, id, my_side_count,
                                          reinterpret_cast<int64_t *>(fconnect), map_ids);
  }

  // Get distribution factors for the specified side block
  int64_t DatabaseIO::get_side_distributions(const Ioss::SideBlock *fb, int64_t id,
                                             int64_t my_side_count, double *dist_fact,
                                             size_t /* data_size */) const
  {
    // Allocate space for elements and local side numbers
    // Get size of data stored on the file...
    ex_set set_param[1];
    set_param[0].id                       = id;
    set_param[0].type                     = EX_SIDE_SET;
    set_param[0].entry_list               = nullptr;
    set_param[0].extra_list               = nullptr;
    set_param[0].distribution_factor_list = nullptr;

    int error = ex_get_sets(get_file_pointer(), 1, set_param);
    if (error < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }
    int64_t number_sides                = set_param[0].num_entry;
    int64_t number_distribution_factors = set_param[0].num_distribution_factor;

    const Ioss::ElementTopology *ftopo   = fb->topology();
    int                          nfnodes = ftopo->number_nodes();

    if (number_distribution_factors == 0) {
      // Fill in the array with '1.0'...
      for (int64_t i = 0; i < nfnodes * my_side_count; i++) {
        dist_fact[i] = 1.0;
      }
      return 0;
    }

    // Take care of the easy situation -- If 'side_count' ==
    // 'number_sides' then the sideset is stored in a single sideblock
    // and all distribution factors on the database are transferred
    // 1-to-1 into 'dist_fact' array.
    if (my_side_count == number_sides) {
      // Verify that number_distribution_factors is sane...
      if (number_sides * nfnodes != number_distribution_factors &&
          number_sides != number_distribution_factors) {
        std::ostringstream errmsg;
        fmt::print(
            errmsg,
            "ERROR: SideBlock '{}' in file '{}'\n"
            "\thas incorrect distribution factor count.\n"
            "\tThere are {} '{}' sides with "
            "{} nodes per side, but there are {} distribution factors which is not correct.\n"
            "\tThere should be either {} or {} distribution factors.\n",
            fb->name(), get_filename(), number_sides, ftopo->name(), nfnodes,
            number_distribution_factors, number_sides, number_sides * nfnodes);
        IOSS_ERROR(errmsg);
      }
      return ex_get_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, dist_fact);
    }

    // Allocate space for distribution factors.
    std::vector<double> dist(number_distribution_factors);
    int ierr = ex_get_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, Data(dist));
    if (ierr < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }

    // Another easy situation (and common for exodus) is if the input
    // distribution factors are all the same value (typically 1).  In
    // that case, we only have to fill in the output array with that
    // value.
    {
      double value    = dist[0];
      bool   constant = true;
      for (int64_t i = 1; i < number_distribution_factors; i++) {
        if (dist[i] != value) {
          constant = false;
          break;
        }
      }
      if (constant) {
        if (value == 0.0) {
          value = 1.0; // Take care of some buggy mesh generators
        }
        for (int64_t j = 0; j < my_side_count * nfnodes; j++) {
          dist_fact[j] = value;
        }
        return 0;
      }
    }

    // If we get to here, the underlying sideset contains multiple side
    // topologies and the distribution factors are non-constant. Need to
    // allocate space to store all distribution factors and then pull
    // out those that are applied to sides with the correct topology.

    // Allocate space for element and local side number (this is bulk
    // data...)
    //----
    std::vector<char> element(number_sides * int_byte_size_api());
    std::vector<char> side(number_sides * int_byte_size_api());

    ierr = ex_get_set(get_file_pointer(), EX_SIDE_SET, id, Data(element), Data(side));
    if (ierr < 0) {
      Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
    }
    //----

    Ioss::IntVector is_valid_side;
    Ioss::Utils::calculate_sideblock_membership(is_valid_side, fb, int_byte_size_api(),
                                                Data(element), Data(side), number_sides,
                                                get_region());

    int64_t             ieb   = 0; // counter for distribution factors in this sideblock
    int64_t             idb   = 0; // counter for distribution factors read from database
    Ioss::ElementBlock *block = nullptr;

    int     *element32 = nullptr;
    int64_t *element64 = nullptr;
    int     *side32    = nullptr;
    int64_t *side64    = nullptr;

    if (int_byte_size_api() == 4) {
      element32 = reinterpret_cast<int *>(Data(element));
      side32    = reinterpret_cast<int *>(Data(side));
    }
    else {
      element64 = reinterpret_cast<int64_t *>(Data(element));
      side64    = reinterpret_cast<int64_t *>(Data(side));
    }

    for (int64_t iel = 0; iel < number_sides; iel++) {
      int64_t elem_id = 0;
      int64_t side_id = 0;
      if (int_byte_size_api() == 4) {
        elem_id = element32[iel];
        side_id = side32[iel];
      }
      else {
        elem_id = element64[iel];
        side_id = side64[iel];
      }

      if (block == nullptr || !block->contains(elem_id)) {
        block = get_region()->get_element_block(elem_id);
      }

      if (block == nullptr) {
        std::ostringstream errmsg;
        fmt::print(errmsg,
                   "INTERNAL ERROR: Could not find element block containing element with id {}. "
                   "Something is wrong in the Ioex::DatabaseIO class. Please report.\n",
                   elem_id);
        IOSS_ERROR(errmsg);
      }

      const Ioss::ElementTopology *topo = block->topology()->boundary_type(side_id);

      if (topo == nullptr) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "INTERNAL ERROR: Could not find topology of element block boundary. "
                           "Something is wrong in the Ioex::DatabaseIO class. Please report.\n");
        IOSS_ERROR(errmsg);
      }

      int nside_nodes = topo->number_nodes();

      if (is_valid_side[iel] == 1) {
        // This side belongs in the sideblock
        for (int64_t i = 0; i < nside_nodes; i++) {
          dist_fact[ieb++] = dist[idb++];
        }
      }
      else {
        // Skip over unused 'dist' factors
        idb += topo->number_nodes();
      }
    }

    assert(ieb == my_side_count * nfnodes);
    // If the following assert fails, it may be due to bug in Patran
    // which writes too many distribution factors to the database in a
    // mixed element case. Note that this is checked earlier also with a
    // better error message.
    assert(idb == number_distribution_factors);
    return 0;
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::Region *reg, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return Ioex::BaseDatabaseIO::put_field_internal(reg, field, data, data_size);
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::NodeBlock *nb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        Ioss::Field::RoleType role = field.get_role();

        if (role == Ioss::Field::MESH) {
          if (field.get_name() == "mesh_model_coordinates_x") {
            auto *rdata = static_cast<double *>(data);
            int   ierr  = ex_put_coord(get_file_pointer(), rdata, nullptr, nullptr);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }

          else if (field.get_name() == "mesh_model_coordinates_y") {
            auto *rdata = static_cast<double *>(data);
            int   ierr  = ex_put_coord(get_file_pointer(), nullptr, rdata, nullptr);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }

          else if (field.get_name() == "mesh_model_coordinates_z") {
            auto *rdata = static_cast<double *>(data);
            int   ierr  = ex_put_coord(get_file_pointer(), nullptr, nullptr, rdata);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }

          else if (field.get_name() == "mesh_model_coordinates") {
            // Data required by upper classes store x0, y0, z0, ... xn, yn, zn
            // Data stored in exodus file is x0, ..., xn, y0, ..., yn, z0, ..., zn
            // so we have to allocate some scratch memory to read in the data
            // and then map into supplied 'data'
            std::vector<double> x;
            x.reserve(num_to_get);
            std::vector<double> y;
            if (spatialDimension > 1) {
              y.reserve(num_to_get);
            }
            std::vector<double> z;
            if (spatialDimension == 3) {
              z.reserve(num_to_get);
            }

            // Cast 'data' to correct size -- double
            auto *rdata = static_cast<double *>(data);

            size_t index = 0;
            for (size_t i = 0; i < num_to_get; i++) {
              x.push_back(rdata[index++]);
              if (spatialDimension > 1) {
                y.push_back(rdata[index++]);
              }
              if (spatialDimension == 3) {
                z.push_back(rdata[index++]);
              }
            }
            int ierr = ex_put_coord(get_file_pointer(), Data(x), Data(y), Data(z));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else if (field.get_name() == "ids") {
            // The ids coming in are the global ids; their position is the
            // local id -1 (That is, data[0] contains the global id of local
            // node 1)
            handle_node_ids(data, num_to_get);
          }
          else if (field.get_name() == "connectivity") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "node_connectivity_status") {
            // Do nothing, input only field.
          }
          else if (field.get_name() == "implicit_ids") {
            // Do nothing, input only field.
          }
          else {
            return Ioss::Utils::field_warning(nb, field, "mesh output");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this node block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // nodes in this block.
          write_nodal_transient_field(field, nb, num_to_get, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, nb, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, nb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::Blob *blob, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        Ioss::Field::RoleType role = field.get_role();

        if (role == Ioss::Field::MESH) {
          if (field.get_name() == "ids") {
            // The ids coming in are the global ids; their position is the
            // local id -1 (That is, data[0] contains the global id of local
            // node 1)
            //          handle_node_ids(data, num_to_get);
          }
          else if (field.get_name() == "connectivity") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "node_connectivity_status") {
            // Do nothing, input only field.
          }
          else if (field.get_name() == "implicit_ids") {
            // Do nothing, input only field.
          }
          else {
            return Ioss::Utils::field_warning(blob, field, "mesh output");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this node block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // nodes in this block.
          write_entity_transient_field(field, blob, num_to_get, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, blob, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, blob, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::Assembly *assembly, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);
      if (num_to_get > 0) {

        Ioss::Field::RoleType role = field.get_role();

        if (role == Ioss::Field::MESH) {
          if (field.get_name() == "ids") {
            // The ids coming in are the global ids; their position is the
            // local id -1 (That is, data[0] contains the global id of local
            // node 1)
            //          handle_node_ids(data, num_to_get);
          }
          else if (field.get_name() == "connectivity") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, just handles an idiosyncrasy of the GroupingEntity
          }
          else if (field.get_name() == "node_connectivity_status") {
            // Do nothing, input only field.
          }
          else if (field.get_name() == "implicit_ids") {
            // Do nothing, input only field.
          }
          else {
            return Ioss::Utils::field_warning(assembly, field, "mesh output");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this node block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // nodes in this block.
          write_entity_transient_field(field, assembly, num_to_get, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, assembly, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, assembly, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::ElementBlock *eb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);

      if (num_to_get > 0) {
        int ierr = 0;

        // Get the element block id and element count
        int64_t               id               = Ioex::get_id(eb, &ids_);
        size_t                my_element_count = eb->entity_count();
        Ioss::Field::RoleType role             = field.get_role();

        if (role == Ioss::Field::MESH) {
          // Handle the MESH fields required for an Exodus file model.
          // (The 'genesis' portion)
          if (field.get_name() == "connectivity") {
            if (my_element_count > 0) {
              // Map element connectivity from global node id to local node id.
              int element_nodes = eb->topology()->number_nodes();
              nodeMap.reverse_map_data(data, field, num_to_get * element_nodes);
              ierr = ex_put_conn(get_file_pointer(), EX_ELEM_BLOCK, id, data, nullptr, nullptr);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "connectivity_edge") {
            if (my_element_count > 0) {
              // Map element connectivity from global edge id to local edge id.
              int element_edges = field.transformed_storage()->component_count();
              edgeMap.reverse_map_data(data, field, num_to_get * element_edges);
              ierr = ex_put_conn(get_file_pointer(), EX_ELEM_BLOCK, id, nullptr, data, nullptr);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "connectivity_face") {
            if (my_element_count > 0) {
              // Map element connectivity from global face id to local face id.
              int element_faces = field.transformed_storage()->component_count();
              faceMap.reverse_map_data(data, field, num_to_get * element_faces);
              ierr = ex_put_conn(get_file_pointer(), EX_ELEM_BLOCK, id, nullptr, nullptr, data);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "connectivity_raw") {
            if (my_element_count > 0) {
              // Element connectivity is already in local node id.
              ierr = ex_put_conn(get_file_pointer(), EX_ELEM_BLOCK, id, data, nullptr, nullptr);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "ids") {
            handle_element_ids(eb, data, num_to_get);
          }
          else if (field.get_name() == "implicit_ids") {
            // Do nothing, input only field.
          }
        }
        else if (role == Ioss::Field::MAP) {
          int comp_count = field.get_component_count(Ioss::Field::InOut::OUTPUT);
          for (int comp = 0; comp < comp_count; comp++) {
            std::vector<char> component(my_element_count * int_byte_size_api());

            if (int_byte_size_api() == 4) {
              auto *data32 = reinterpret_cast<int *>(data);
              auto *comp32 = reinterpret_cast<int *>(Data(component));

              int index = comp;
              for (size_t i = 0; i < my_element_count; i++) {
                comp32[i] = data32[index];
                index += comp_count;
              }
            }
            else {
              auto *data64 = reinterpret_cast<int64_t *>(data);
              auto *comp64 = reinterpret_cast<int64_t *>(Data(component));

              int index = comp;
              for (size_t i = 0; i < my_element_count; i++) {
                comp64[i] = data64[index];
                index += comp_count;
              }
            }
            size_t eb_offset = eb->get_offset();
            int    index     = -1 * (static_cast<int>(field.get_index()) + comp);
            ierr = ex_put_partial_num_map(get_file_pointer(), EX_ELEM_MAP, index, eb_offset + 1,
                                          my_element_count, Data(component));
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, eb, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this element block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // elements in this block.
          write_entity_transient_field(field, eb, my_element_count, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, eb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::FaceBlock *eb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);

      if (num_to_get > 0) {
        int ierr = 0;

        // Get the face block id and face count
        int64_t               id            = Ioex::get_id(eb, &ids_);
        int64_t               my_face_count = eb->entity_count();
        Ioss::Field::RoleType role          = field.get_role();

        if (role == Ioss::Field::MESH) {
          // Handle the MESH fields required for an Exodus file model.
          // (The 'genesis' portion)
          if (field.get_name() == "connectivity") {
            if (my_face_count > 0) {
              // Map face connectivity from global node id to local node id.
              // Do it in 'data' ...
              int face_nodes = eb->topology()->number_nodes();
              nodeMap.reverse_map_data(data, field, num_to_get * face_nodes);
              ierr = ex_put_conn(get_file_pointer(), EX_FACE_BLOCK, id, data, nullptr, nullptr);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "connectivity_edge") {
            if (my_face_count > 0) {
              // Map face connectivity from global edge id to local edge id.
              int face_edges = field.transformed_storage()->component_count();
              edgeMap.reverse_map_data(data, field, num_to_get * face_edges);
              ierr = ex_put_conn(get_file_pointer(), EX_FACE_BLOCK, id, nullptr, data, nullptr);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, input only field.
          }
          else if (field.get_name() == "ids") {
            handle_face_ids(eb, data, num_to_get);
          }
          else {
            num_to_get = Ioss::Utils::field_warning(eb, field, "mesh output");
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, eb, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this face block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // faces in this block.
          write_entity_transient_field(field, eb, my_face_count, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, eb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::EdgeBlock *eb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);

      size_t num_to_get = field.verify(data_size);

      if (num_to_get > 0) {
        int ierr = 0;

        // Get the edge block id and edge count
        int64_t               id            = Ioex::get_id(eb, &ids_);
        int64_t               my_edge_count = eb->entity_count();
        Ioss::Field::RoleType role          = field.get_role();

        if (role == Ioss::Field::MESH) {
          // Handle the MESH fields required for an Exodus file model. (The 'genesis' portion)
          if (field.get_name() == "connectivity") {
            if (my_edge_count > 0) {
              // Map edge connectivity from global node id to local node id.
              // Do it in 'data' ...
              int edge_nodes = eb->topology()->number_nodes();
              nodeMap.reverse_map_data(data, field, num_to_get * edge_nodes);
              ierr = ex_put_conn(get_file_pointer(), EX_EDGE_BLOCK, id, data, nullptr, nullptr);
              if (ierr < 0) {
                Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
              }
            }
          }
          else if (field.get_name() == "connectivity_raw") {
            // Do nothing, input only field.
          }
          else if (field.get_name() == "ids") {
            handle_edge_ids(eb, data, num_to_get);
          }
          else {
            num_to_get = Ioss::Utils::field_warning(eb, field, "mesh output");
          }
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, eb, data);
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this edge block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // edges in this block.
          write_entity_transient_field(field, eb, my_edge_count, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, eb, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::handle_node_ids(void *ids, int64_t num_to_get) const
  {
    /*!
     * There are two modes we need to support in this routine:
     * 1. Initial definition of node map (local->global) and
     * nodeMap.reverse (global->local).
     * 2. Redefinition of node map via 'reordering' of the original
     * map when the nodes on this processor are the same, but their
     * order is changed (or count because of ghosting)
     *
     * So, there will be two maps the 'nodeMap.map' map is a 'direct lookup'
     * map which maps current local position to global id and the
     * 'nodeMap.reverse' is an associative lookup which maps the
     * global id to 'original local'.  There is also a
     * 'nodeMap.reorder' which is direct lookup and maps current local
     * position to original local.

     * The ids coming in are the global ids; their position is the
     * "local id-1" (That is, data[0] contains the global id of local
     * node 1 in this node block).
     *
     * int local_position = nodeMap.reverse[NodeMap[i+1]]
     * (the nodeMap.map and nodeMap.reverse are 1-based)
     *
     * To determine which map to update on a call to this function, we
     * use the following heuristics:
     * -- If the database state is 'STATE_MODEL:', then update the
     *    'nodeMap.reverse' and 'nodeMap.map'
     *
     * -- If the database state is not STATE_MODEL, then leave the
     *    'nodeMap.reverse' and 'nodeMap.map' alone since they correspond to the
     *    information already written to the database. [May want to add a
     *    STATE_REDEFINE_MODEL]
     *
     * -- In both cases, update the nodeMap.reorder
     *
     * NOTE: The mapping is done on TRANSIENT fields only; MODEL fields
     *       should be in the original order...
     */
    SMART_ASSERT(num_to_get == nodeCount)(num_to_get)(nodeCount);

    nodeMap.set_size(nodeCount);

    bool in_define = (dbState == Ioss::STATE_MODEL) || (dbState == Ioss::STATE_DEFINE_MODEL);
    if (int_byte_size_api() == 4) {
      nodeMap.set_map(static_cast<int *>(ids), num_to_get, 0, in_define);
    }
    else {
      nodeMap.set_map(static_cast<int64_t *>(ids), num_to_get, 0, in_define);
    }

    if (in_define) {
      // Only a single nodeblock and all set
      assert(get_region()->get_property("node_block_count").get_int() == 1);

      // Write to the database...
      if (ex_put_id_map(get_file_pointer(), EX_NODE_MAP, ids) < 0) {
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
      }
    }
    return num_to_get;
  }

  int64_t DatabaseIO::handle_element_ids(const Ioss::ElementBlock *eb, void *ids,
                                         size_t num_to_get) const
  {
    elemMap.set_size(elementCount);
    size_t offset = eb->get_offset();
    return handle_block_ids(eb, EX_ELEM_MAP, elemMap, ids, num_to_get, offset);
  }

  int64_t DatabaseIO::handle_face_ids(const Ioss::FaceBlock *eb, void *ids, size_t num_to_get) const
  {
    faceMap.set_size(faceCount);
    size_t offset = eb->get_offset();
    return handle_block_ids(eb, EX_FACE_MAP, faceMap, ids, num_to_get, offset);
  }

  int64_t DatabaseIO::handle_edge_ids(const Ioss::EdgeBlock *eb, void *ids, size_t num_to_get) const
  {
    edgeMap.set_size(edgeCount);
    size_t offset = eb->get_offset();
    return handle_block_ids(eb, EX_EDGE_MAP, edgeMap, ids, num_to_get, offset);
  }

  void DatabaseIO::write_nodal_transient_field(const Ioss::Field &field,
                                               const Ioss::NodeBlock * /* ge */, int64_t count,
                                               void *variables) const
  {
    Ioss::Field::BasicType ioss_type = field.get_type();
    assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER ||
           ioss_type == Ioss::Field::INT64 || ioss_type == Ioss::Field::COMPLEX);

    if (ioss_type == Ioss::Field::INT64) {
      Ioss::Utils::check_int_to_real_overflow(field, (int64_t *)variables, count);
    }

    // Note that if the field's basic type is COMPLEX, then each component of
    // the VariableType is a complex variable consisting of a real and
    // imaginary part.  Since exodus cannot handle complex variables,
    // we have to output a (real and imaginary) X (number of
    // components) fields. For example, if V is a 3d vector of complex
    // data, the data in the 'variables' array are v_x, v.im_x, v_y,
    // v.im_y, v_z, v.im_z which need to be output in six separate
    // exodus fields.  These fields were already defined in
    // "write_results_metadata".

    std::vector<double> temp(count);

    int step = get_current_state();
    step     = get_database_step(step);

    // get number of components, cycle through each component
    // and add suffix to base 'field_name'.  Look up index
    // of this name in 'm_variables[EX_NODE_BLOCK]' map
    int comp_count = field.get_component_count(Ioss::Field::InOut::OUTPUT);
    int var_index  = 0;

    int re_im = 1;
    if (ioss_type == Ioss::Field::COMPLEX) {
      re_im = 2;
    }
    for (int complex_comp = 0; complex_comp < re_im; complex_comp++) {
      std::string field_name = field.get_name();
      if (re_im == 2) {
        field_name += complex_suffix[complex_comp];
      }

      for (int i = 0; i < comp_count; i++) {
        std::string var_name = get_component_name(field, Ioss::Field::InOut::OUTPUT, i + 1);

        auto var_iter = m_variables[EX_NODE_BLOCK].find(var_name);
        if (var_iter == m_variables[EX_NODE_BLOCK].end()) {
          std::ostringstream errmsg;
          fmt::print(errmsg, "ERROR: Could not find nodal variable '{}'\n", var_name);
          IOSS_ERROR(errmsg);
        }
        var_index = var_iter->second;

        size_t  begin_offset = (re_im * i) + complex_comp;
        size_t  stride       = re_im * comp_count;
        int64_t num_out      = 0;

        if (ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::COMPLEX) {
          num_out = nodeMap.map_field_to_db_scalar_order(static_cast<double *>(variables), temp,
                                                         begin_offset, count, stride, 0);
        }
        else if (ioss_type == Ioss::Field::INTEGER) {
          num_out = nodeMap.map_field_to_db_scalar_order(static_cast<int *>(variables), temp,
                                                         begin_offset, count, stride, 0);
        }
        else if (ioss_type == Ioss::Field::INT64) {
          num_out = nodeMap.map_field_to_db_scalar_order(static_cast<int64_t *>(variables), temp,
                                                         begin_offset, count, stride, 0);
        }

        if (num_out != nodeCount) {
          std::ostringstream errmsg;
          fmt::print(errmsg,
                     "ERROR: Problem outputting nodal variable '{}' with index = {} to file '{}'\n"
                     "Should have output {} values, but instead only output {} values.\n",
                     var_name, var_index, decoded_filename(), nodeCount, num_out);
          IOSS_ERROR(errmsg);
        }

        // Write the variable...
        int ierr =
            ex_put_var(get_file_pointer(), step, EX_NODE_BLOCK, var_index, 0, num_out, Data(temp));
        if (ierr < 0) {
          std::ostringstream errmsg;
          fmt::print(errmsg,
                     "Problem outputting nodal variable '{}' with index = {} to file '{}'\n",
                     var_name, var_index, decoded_filename());
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__, errmsg.str());
        }
      }
    }
  }

  void DatabaseIO::write_entity_transient_field(const Ioss::Field          &field,
                                                const Ioss::GroupingEntity *ge, int64_t count,
                                                void *variables) const
  {
    static Ioss::Map    non_element_map; // Used as an empty map for ge->type() != element block.
    std::vector<double> temp(count);

    int step = get_current_state();
    step     = get_database_step(step);

    Ioss::Map *map       = nullptr;
    int64_t    eb_offset = 0;
    if (ge->type() == Ioss::ELEMENTBLOCK) {
      const auto *elb = dynamic_cast<const Ioss::ElementBlock *>(ge);
      Ioss::Utils::check_dynamic_cast(elb);
      eb_offset = elb->get_offset();
      map       = &elemMap;
    }
    else {
      map = &non_element_map;
    }

    Ioss::Field::BasicType ioss_type = field.get_type();
    assert(ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::INTEGER ||
           ioss_type == Ioss::Field::INT64 || ioss_type == Ioss::Field::COMPLEX);

    if (ioss_type == Ioss::Field::INT64) {
      Ioss::Utils::check_int_to_real_overflow(field, (int64_t *)variables, count);
    }

    // Note that if the field's basic type is COMPLEX, then each component of
    // the VariableType is a complex variable consisting of a real and
    // imaginary part.  Since exodus cannot handle complex variables,
    // we have to output a (real and imaginary) X (number of
    // components) fields. For example, if V is a 3d vector of complex
    // data, the data in the 'variables' array are v_x, v.im_x, v_y,
    // v.im_y, v_z, v.im_z which need to be output in six separate
    // exodus fields.  These fields were already defined in
    // "write_results_metadata".

    // get number of components, cycle through each component
    // and add suffix to base 'field_name'.  Look up index
    // of this name in 'm_variables[type]' map
    int comp_count = field.get_component_count(Ioss::Field::InOut::OUTPUT);

    // Handle quick easy, hopefully common case first...
    ex_entity_type type = Ioex::map_exodus_type(ge->type());
    if (comp_count == 1 && ioss_type == Ioss::Field::REAL && type != EX_SIDE_SET &&
        !map->reorders()) {
      // Simply output the variable...
      int64_t     id       = Ioex::get_id(ge, &ids_);
      std::string var_name = get_component_name(field, Ioss::Field::InOut::OUTPUT, 1);
      auto        var_iter = m_variables[type].find(var_name);
      if (var_iter == m_variables[type].end()) {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Could not find field '{}'\n", var_name);
        IOSS_ERROR(errmsg);
      }
      int var_index = var_iter->second;
      assert(var_index > 0);
      int ierr = ex_put_var(get_file_pointer(), step, type, var_index, id, count, variables);

      if (ierr < 0) {
        std::ostringstream extra_info;
        fmt::print(extra_info, "Outputting field {} at step {} on {} {}.", field.get_name(),
                   fmt::group_digits(step), ge->type_string(), ge->name());
        Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__, extra_info.str());
      }
      return;
    }
    int re_im = 1;
    if (ioss_type == Ioss::Field::COMPLEX) {
      re_im = 2;
    }
    for (int complex_comp = 0; complex_comp < re_im; complex_comp++) {
      std::string field_name = field.get_name();
      if (re_im == 2) {
        field_name += complex_suffix[complex_comp];
      }

      for (int i = 0; i < comp_count; i++) {
        std::string var_name = get_component_name(field, Ioss::Field::InOut::OUTPUT, i + 1);

        auto var_iter = m_variables[type].find(var_name);
        if (var_iter == m_variables[type].end()) {
          std::ostringstream errmsg;
          fmt::print(errmsg, "ERROR: Could not find field '{}'\n", var_name);
          IOSS_ERROR(errmsg);
        }
        int var_index = var_iter->second;
        assert(var_index > 0);

        // var is a [count,comp,re_im] array;  re_im = 1(real) or 2(complex)
        // beg_offset = (re_im*i)+complex_comp
        // number_values = count
        // stride = re_im*comp_count
        int64_t begin_offset = (re_im * i) + complex_comp;
        int64_t stride       = re_im * comp_count;

        if (ioss_type == Ioss::Field::REAL || ioss_type == Ioss::Field::COMPLEX) {
          map->map_field_to_db_scalar_order(static_cast<double *>(variables), temp, begin_offset,
                                            count, stride, eb_offset);
        }
        else if (ioss_type == Ioss::Field::INTEGER) {
          map->map_field_to_db_scalar_order(static_cast<int *>(variables), temp, begin_offset,
                                            count, stride, eb_offset);
        }
        else if (ioss_type == Ioss::Field::INT64) {
          map->map_field_to_db_scalar_order(static_cast<int64_t *>(variables), temp, begin_offset,
                                            count, stride, eb_offset);
        }

        // Write the variable...
        int64_t id = Ioex::get_id(ge, &ids_);
        int     ierr;
        if (type == EX_SIDE_SET) {
          size_t offset = ge->get_property("set_offset").get_int();
          ierr = ex_put_partial_var(get_file_pointer(), step, type, var_index, id, offset + 1,
                                    count, Data(temp));
        }
        else {
          ierr = ex_put_var(get_file_pointer(), step, type, var_index, id, count, Data(temp));
        }

        if (ierr < 0) {
          std::ostringstream extra_info;
          fmt::print(extra_info, "Outputting component {} of field {} at step {} on {} {}.", i,
                     field_name, fmt::group_digits(step), ge->type_string(), ge->name());
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__, extra_info.str());
        }
      }
    }
  }

  int64_t DatabaseIO::put_Xset_field_internal(const Ioss::EntitySet *ns, const Ioss::Field &field,
                                              void *data, size_t data_size) const
  {
    {
      Ioss::SerializeIO serializeio_(this);
      //    ex_update(get_file_pointer());

      size_t entity_count = ns->entity_count();
      size_t num_to_get   = field.verify(data_size);
      if (num_to_get > 0) {

        ex_entity_type        type = Ioex::map_exodus_type(ns->type());
        int64_t               id   = Ioex::get_id(ns, &ids_);
        Ioss::Field::RoleType role = field.get_role();

        if (role == Ioss::Field::MESH) {

          if (field.get_name() == "ids" || field.get_name() == "ids_raw") {
            // Map node id from global node id to local node id.
            // Do it in 'data' ...

            if (field.get_name() == "ids") {
              nodeMap.reverse_map_data(data, field, num_to_get);
            }
            int ierr = ex_put_set(get_file_pointer(), type, id, data, nullptr);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else if (field.get_name() == "orientation") {
            int ierr = ex_put_set(get_file_pointer(), type, id, nullptr, data);
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else if (field.get_name() == "distribution_factors") {
            int ierr =
                ex_put_set_dist_fact(get_file_pointer(), type, id, static_cast<double *>(data));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else {
            num_to_get = Ioss::Utils::field_warning(ns, field, "output");
          }
        }
        else if (role == Ioss::Field::TRANSIENT) {
          // Check if the specified field exists on this element block.
          // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
          // exist on the database as scalars with the appropriate
          // extensions.

          // Transfer each component of the variable into 'data' and then
          // output.  Need temporary storage area of size 'number of
          // elements in this block.
          write_entity_transient_field(field, ns, entity_count, data);
        }
        else if (role == Ioss::Field::ATTRIBUTE) {
          num_to_get = write_attribute_field(field, ns, data);
        }
        else if (role == Ioss::Field::REDUCTION) {
          store_reduction_field(field, ns, data);
        }
      }
      return num_to_get;
    }
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::NodeSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return put_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::EdgeSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return put_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::FaceSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return put_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::ElementSet *ns, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    return put_Xset_field_internal(ns, field, data, data_size);
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::CommSet *cs, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    size_t num_to_get   = field.verify(data_size);
    size_t entity_count = cs->entity_count();

    assert(num_to_get == entity_count);
    if (num_to_get == 0) {
      return 0;
    }

    // Return the <entity (node or side), processor> pair
    if (field.get_name() == "entity_processor") {

      // Check type -- node or side
      std::string type = cs->get_property("entity_type").get_string();

      // Allocate temporary storage space
      std::vector<char> entities(entity_count * int_byte_size_api());
      std::vector<char> procs(entity_count * int_byte_size_api());

      if (type == "node") {
        Ioss::SerializeIO serializeio_(this);
        // Convert global node id to local node id and store in 'entities'
        if (int_byte_size_api() == 4) {
          int *entity_proc = static_cast<int *>(data);
          int *ent         = reinterpret_cast<int *>(Data(entities));
          int *pro         = reinterpret_cast<int *>(Data(procs));
          int  j           = 0;
          for (size_t i = 0; i < entity_count; i++) {
            int global_id = entity_proc[j++];
            ent[i]        = nodeMap.global_to_local(global_id, true);
            pro[i]        = entity_proc[j++];
          }
        }
        else {
          auto   *entity_proc = static_cast<int64_t *>(data);
          auto   *ent         = reinterpret_cast<int64_t *>(Data(entities));
          auto   *pro         = reinterpret_cast<int64_t *>(Data(procs));
          int64_t j           = 0;
          for (size_t i = 0; i < entity_count; i++) {
            int64_t global_id = entity_proc[j++];
            ent[i]            = nodeMap.global_to_local(global_id, true);
            pro[i]            = entity_proc[j++];
          }
        }

        if (commsetNodeCount > 0) {
          int ierr = ex_put_node_cmap(get_file_pointer(), Ioex::get_id(cs, &ids_), Data(entities),
                                      Data(procs), myProcessor);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }

        if (commsetNodeCount == 1) {
          // NOTE: The internal and border node maps must be output in one call.
          //       In this routine, we only have one commset at a time and can't
          //       construct the entire map at one time.  This is not really needed,
          //       so for now we just skip if there is more than one commset.  If
          //       this information is really needed, need to cache the information
          //       until all commsets have been processed.  Also need to change
          //       write_communication_metada() [Maybe, unless client sets correct
          //       properties.]

          // Construct the node map (internal vs. border).
          // Border nodes are those in the communication map (use entities array)
          // Internal nodes are the rest.  Allocate array to hold all nodes,
          // initialize all to '1', then zero out the nodes in 'entities'.
          // Iterate through array again and consolidate all '1's

          std::vector<char> internal(nodeCount * int_byte_size_api());
          if (int_byte_size_api() == 4) {
            compute_internal_border_maps(reinterpret_cast<int *>(Data(entities)),
                                         reinterpret_cast<int *>(Data(internal)), nodeCount,
                                         entity_count);
          }
          else {
            compute_internal_border_maps(reinterpret_cast<int64_t *>(Data(entities)),
                                         reinterpret_cast<int64_t *>(Data(internal)), nodeCount,
                                         entity_count);
          }

          int ierr = ex_put_processor_node_maps(get_file_pointer(), Data(internal), Data(entities),
                                                nullptr, myProcessor);
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
      }
      else if (type == "side") {
        Ioss::SerializeIO serializeio_(this);
        std::vector<char> sides(entity_count * int_byte_size_api());
        if (int_byte_size_api() == 4) {
          int *entity_proc = static_cast<int *>(data);
          int *ent         = reinterpret_cast<int *>(Data(entities));
          int *sid         = reinterpret_cast<int *>(Data(sides));
          int *pro         = reinterpret_cast<int *>(Data(procs));
          int  j           = 0;
          for (size_t i = 0; i < entity_count; i++) {
            ent[i] = elemMap.global_to_local(entity_proc[j++]);
            sid[i] = entity_proc[j++];
            pro[i] = entity_proc[j++];
          }
        }
        else {
          auto   *entity_proc = static_cast<int64_t *>(data);
          auto   *ent         = reinterpret_cast<int64_t *>(Data(entities));
          auto   *sid         = reinterpret_cast<int64_t *>(Data(sides));
          auto   *pro         = reinterpret_cast<int64_t *>(Data(procs));
          int64_t j           = 0;
          for (size_t i = 0; i < entity_count; i++) {
            ent[i] = elemMap.global_to_local(entity_proc[j++]);
            sid[i] = entity_proc[j++];
            pro[i] = entity_proc[j++];
          }
        }

        int ierr = ex_put_elem_cmap(get_file_pointer(), Ioex::get_id(cs, &ids_), Data(entities),
                                    Data(sides), Data(procs), myProcessor);
        if (ierr < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }

        // Construct the element map (internal vs. border).
        // Border elements are those in the communication map (use entities array)
        // Internal elements are the rest.  Allocate array to hold all elements,
        // initialize all to '1', then zero out the elements in 'entities'.
        // Iterate through array again and consolidate all '1's
        std::vector<char> internal(elementCount * int_byte_size_api());
        if (int_byte_size_api() == 4) {
          compute_internal_border_maps(reinterpret_cast<int *>(Data(entities)),
                                       reinterpret_cast<int *>(Data(internal)), elementCount,
                                       entity_count);
        }
        else {
          compute_internal_border_maps(reinterpret_cast<int64_t *>(Data(entities)),
                                       reinterpret_cast<int64_t *>(Data(internal)), elementCount,
                                       entity_count);
        }

        ierr = ex_put_processor_elem_maps(get_file_pointer(), Data(internal), Data(entities),
                                          myProcessor);
        if (ierr < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }
      }
      else {
        std::ostringstream errmsg;
        fmt::print(errmsg, "ERROR: Invalid commset type {}", type);
        IOSS_ERROR(errmsg);
      }
    }
    else if (field.get_name() == "ids") {
      // Do nothing, just handles an idiosyncrasy of the GroupingEntity
    }
    else {
      num_to_get = Ioss::Utils::field_warning(cs, field, "output");
    }
    return num_to_get;
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::SideSet *fs, const Ioss::Field &field,
                                         void * /* data */, size_t data_size) const
  {
    size_t num_to_get = field.verify(data_size);
    if (field.get_name() == "ids") {
      // Do nothing, just handles an idiosyncrasy of the GroupingEntity
    }
    else {
      num_to_get = Ioss::Utils::field_warning(fs, field, "output");
    }
    return num_to_get;
  }

  int64_t DatabaseIO::put_field_internal(const Ioss::SideBlock *fb, const Ioss::Field &field,
                                         void *data, size_t data_size) const
  {
    Ioss::SerializeIO serializeio_(this);
    size_t            num_to_get = field.verify(data_size);
    if (num_to_get > 0) {

      int64_t id = Ioex::get_id(fb, &ids_);

      size_t entity_count = fb->entity_count();
      size_t offset       = fb->get_property("set_offset").get_int();

      Ioss::Field::RoleType role = field.get_role();

      if (role == Ioss::Field::MESH) {
        if (field.get_name() == "side_ids" && fb->name() == "universal_sideset") {
          // The side ids are being stored as the distribution factor
          // field on the universal sideset.  There should be no other
          // side sets that request this field...  (Eventually,
          // create an id field to store this info.

          // Need to convert 'ints' to 'double' for storage on mesh...
          // FIX 64
          if (field.get_type() == Ioss::Field::INTEGER) {
            int                *ids = static_cast<int *>(data);
            std::vector<double> real_ids(num_to_get);
            for (size_t i = 0; i < num_to_get; i++) {
              real_ids[i] = static_cast<double>(ids[i]);
            }
            int ierr = ex_put_partial_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, offset + 1,
                                                    entity_count, Data(real_ids));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else {
            auto               *ids = static_cast<int64_t *>(data);
            std::vector<double> real_ids(num_to_get);
            for (size_t i = 0; i < num_to_get; i++) {
              real_ids[i] = static_cast<double>(ids[i]);
            }
            int ierr = ex_put_partial_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, offset + 1,
                                                    entity_count, Data(real_ids));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
        }

        else if (field.get_name() == "side_ids") {
        }

        else if (field.get_name() == "ids") {
          // =============================================================
          // NOTE: Code is currently commented out since we have
          // redundant ways of getting the data (element/side) out to
          // the database.  The 'ids' field method relies on a numbering
          // kluge, so for now trying the 'element_side' field...
          // =============================================================
        }

        else if (field.get_name() == "distribution_factors") {
          int    ierr;
          size_t df_offset = fb->get_property("set_df_offset").get_int();
          size_t df_count  = fb->get_property("distribution_factor_count").get_int();
          ierr = ex_put_partial_set_dist_fact(get_file_pointer(), EX_SIDE_SET, id, df_offset + 1,
                                              df_count, static_cast<double *>(data));
          if (ierr < 0) {
            Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
          }
        }
        else if (field.get_name() == "element_side") {
          // In exodus, the 'side block' is stored as a sideset.  A
          // sideset has a list of elements and a corresponding local
          // element side (1-based)

          // The 'data' passed into the function is stored as a
          // 2D vector e0,f0,e1,f1,... (e=element, f=side)

          // To avoid overwriting the passed in data, we allocate
          // two arrays to store the data for this sideset.

          // The element_id passed in is the global id; we need to
          // output the local id.

          // Allocate space for local side number and element numbers
          // numbers.
          // See if edges or faces...
          size_t side_offset = Ioss::Utils::get_side_offset(fb);

          size_t index = 0;

          if (field.get_type() == Ioss::Field::INTEGER) {
            Ioss::IntVector element(num_to_get);
            Ioss::IntVector side(num_to_get);
            int            *el_side = reinterpret_cast<int *>(data);

            try {
              for (size_t i = 0; i < num_to_get; i++) {
                element[i] = elemMap.global_to_local(el_side[index++]);
                side[i]    = el_side[index++] + side_offset;
              }
            }
            catch (const std::runtime_error &x) {
              std::ostringstream errmsg;
              fmt::print(errmsg, "{}On SideBlock `{}` while outputting field `elem_side`\n",
                         x.what(), fb->name());
              IOSS_ERROR(errmsg);
            }

            int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, offset + 1,
                                          entity_count, Data(element), Data(side));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else {
            Ioss::Int64Vector element(num_to_get);
            Ioss::Int64Vector side(num_to_get);
            auto             *el_side = reinterpret_cast<int64_t *>(data);

            try {
              for (size_t i = 0; i < num_to_get; i++) {
                element[i] = elemMap.global_to_local(el_side[index++]);
                side[i]    = el_side[index++] + side_offset;
              }
            }
            catch (const std::runtime_error &x) {
              std::ostringstream errmsg;
              fmt::print(errmsg, "{}On SideBlock `{}` while outputting field `elem_side`\n",
                         x.what(), fb->name());
              IOSS_ERROR(errmsg);
            }

            int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, offset + 1,
                                          entity_count, Data(element), Data(side));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
        }
        else if (field.get_name() == "element_side_raw") {
          // In exodus, the 'side block' is stored as a sideset.  A
          // sideset has a list of elements and a corresponding local
          // element side (1-based)

          // The 'data' passed into the function is stored as a
          // 2D vector e0,f0,e1,f1,... (e=element, f=side)

          // To avoid overwriting the passed in data, we allocate
          // two arrays to store the data for this sideset.

          // The element_id passed in is the local id.

          // See if edges or faces...
          size_t side_offset = Ioss::Utils::get_side_offset(fb);

          size_t index = 0;
          if (field.get_type() == Ioss::Field::INTEGER) {
            Ioss::IntVector element(num_to_get);
            Ioss::IntVector side(num_to_get);
            int            *el_side = reinterpret_cast<int *>(data);

            for (size_t i = 0; i < num_to_get; i++) {
              element[i] = el_side[index++];
              side[i]    = el_side[index++] + side_offset;
            }

            int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, offset + 1,
                                          entity_count, Data(element), Data(side));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
          else {
            Ioss::Int64Vector element(num_to_get);
            Ioss::Int64Vector side(num_to_get);
            auto             *el_side = reinterpret_cast<int64_t *>(data);

            for (size_t i = 0; i < num_to_get; i++) {
              element[i] = el_side[index++];
              side[i]    = el_side[index++] + side_offset;
            }

            int ierr = ex_put_partial_set(get_file_pointer(), EX_SIDE_SET, id, offset + 1,
                                          entity_count, Data(element), Data(side));
            if (ierr < 0) {
              Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
            }
          }
        }
        else if (field.get_name() == "connectivity") {
          // Do nothing, just handles an idiosyncrasy of the GroupingEntity
        }
        else if (field.get_name() == "connectivity_raw") {
          // Do nothing, just handles an idiosyncrasy of the GroupingEntity
        }
        else {
          num_to_get = Ioss::Utils::field_warning(fb, field, "output");
        }
      }
      else if (role == Ioss::Field::TRANSIENT) {
        // Check if the specified field exists on this block.
        // Note that 'higher-order' storage types (e.g. SYM_TENSOR)
        // exist on the database as scalars with the appropriate
        // extensions.

        // Transfer each component of the variable into 'data' and then
        // output.  Need temporary storage area of size 'number of
        // entities in this block.
        write_entity_transient_field(field, fb, entity_count, data);
      }
      else if (role == Ioss::Field::ATTRIBUTE) {
        num_to_get = write_attribute_field(field, fb, data);
      }
      else if (role == Ioss::Field::REDUCTION) {
        store_reduction_field(field, fb, data);
      }
    }
    return num_to_get;
  }

  void DatabaseIO::write_meta_data(Ioss::IfDatabaseExistsBehavior behavior)
  {
    Ioss::Region *region = get_region();
    common_write_metadata(behavior);

    char the_title[max_line_length + 1];

    // Title...
    if (region->property_exists("title")) {
      std::string title_str = region->get_property("title").get_string();
      Ioss::Utils::copy_string(the_title, title_str);
    }
    else {
      Ioss::Utils::copy_string(the_title, "IOSS Default Output Title");
    }

    bool       file_per_processor = true;
    Ioex::Mesh mesh(spatialDimension, the_title, util(), file_per_processor);
    {
      bool omit_maps = false;
      Ioss::Utils::check_set_bool_property(properties, "OMIT_EXODUS_NUM_MAPS", omit_maps);
      if (omit_maps) {
        // Used for special cases only -- typically very large meshes with *known* 1..count maps
        // and workarounds that avoid calling the "ids" put_field calls.
        mesh.use_node_map = false;
        mesh.use_elem_map = false;
        mesh.use_face_map = false;
        mesh.use_edge_map = false;
      }

      bool minimal_nemesis = false;
      Ioss::Utils::check_set_bool_property(properties, "MINIMAL_NEMESIS_DATA", minimal_nemesis);
      if (minimal_nemesis) {
        // Only output the node communication map data... This is all that stk/sierra needs
        mesh.full_nemesis_data = false;
      }

      Ioss::SerializeIO serializeio_(this);
      mesh.populate(region);
      gather_communication_metadata(&mesh.comm);

      if (behavior != Ioss::DB_APPEND && behavior != Ioss::DB_MODIFY) {
        bool omit_qa = false;
        Ioss::Utils::check_set_bool_property(properties, "OMIT_QA_RECORDS", omit_qa);
        if (!omit_qa) {
          put_qa();
        }

        bool omit_info = false;
        Ioss::Utils::check_set_bool_property(properties, "OMIT_INFO_RECORDS", omit_info);
        if (!omit_info) {
          put_info();
        }

        // Write the metadata to the exodus file...
        Ioex::Internals data(get_file_pointer(), maximumNameLength, util());
        int             ierr = data.write_meta_data(mesh);

        if (ierr < 0) {
          Ioex::exodus_error(get_file_pointer(), __LINE__, __func__, __FILE__);
        }
        output_other_metadata();
      }
    }
  }

  void DatabaseIO::gather_communication_metadata(Ioex::CommunicationMetaData *meta)
  {
    // It's possible that we are a serial program outputting information
    // for later use by a parallel program.

    meta->processorCount = 0;
    meta->processorId    = 0;
    meta->outputNemesis  = false;

    if (isParallel) {
      meta->processorCount = util().parallel_size();
      meta->processorId    = myProcessor;
      meta->outputNemesis  = true;
    }
    else {
      if (properties.exists("processor_count")) {
        meta->processorCount = properties.get("processor_count").get_int();
      }
      else if (get_region()->property_exists("processor_count")) {
        meta->processorCount = get_region()->get_property("processor_count").get_int();
      }

      if (properties.exists("my_processor")) {
        meta->processorId = properties.get("my_processor").get_int();
      }
      else if (get_region()->property_exists("my_processor")) {
        meta->processorId = get_region()->get_property("my_processor").get_int();
      }

      if (!get_region()->get_commsets().empty()) {
        isSerialParallel    = true;
        meta->outputNemesis = true;
      }
    }

    if (isSerialParallel || meta->processorCount > 0) {
      meta->globalNodes         = get_region()->get_optional_property("global_node_count", 1);
      meta->globalElements      = get_region()->get_optional_property("global_element_count", 1);
      meta->globalElementBlocks = get_region()->get_optional_property(
          "global_element_block_count", get_region()->get_element_blocks().size());
      meta->globalNodeSets = get_region()->get_optional_property(
          "global_node_set_count", get_region()->get_nodesets().size());
      meta->globalSideSets = get_region()->get_optional_property(
          "global_side_set_count", get_region()->get_sidesets().size());

      // ========================================================================
      // Load balance parameters (NEMESIS, p15)
      meta->nodesInternal = get_region()->get_optional_property("internal_node_count", nodeCount);
      meta->nodesBorder   = get_region()->get_optional_property("border_node_count", 0);
      meta->nodesExternal = 0; // Shadow nodes == 0 for now
      meta->elementsInternal =
          get_region()->get_optional_property("internal_element_count", elementCount);
      meta->elementsBorder = get_region()->get_optional_property("border_element_count", 0);

      const Ioss::CommSetContainer &comm_sets = get_region()->get_commsets();
      for (auto &cs : comm_sets) {
        std::string type  = cs->get_property("entity_type").get_string();
        size_t      count = cs->entity_count();
        int64_t     id    = Ioex::get_id(cs, &ids_);

        if (type == "node") {
          meta->nodeMap.emplace_back(id, count, 'n');
        }
        else if (type == "side") {
          meta->elementMap.emplace_back(id, count, 'e');
        }
        else {
          std::ostringstream errmsg;
          fmt::print(errmsg, "Internal Program Error...");
          IOSS_ERROR(errmsg);
        }
      }
    }
    commsetNodeCount = meta->nodeMap.size();
    commsetElemCount = meta->elementMap.size();
  }
} // namespace Ioex
