// Copyright(C) 2009-2010-2017 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.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//
//     * Redistributions in binary form must reproduce the above
//       copyright notice, this list of conditions and the following
//       disclaimer in the documentation and/or other materials provided
//       with the distribution.
//     * Neither the name of NTESS nor the names of its
//       contributors may be used to endorse or promote products derived
//       from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <algorithm>
#include <exception>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <numeric>
#include <set>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

#include <cfloat>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <limits>
#include <unistd.h>

#include <cctype>
#include <cstring>
#include <ctime>
#include <sys/times.h>
#include <sys/utsname.h>

#include "add_to_log.h"
#include "adler.h"
#include "smart_assert.h"
#include <exodusII.h>

#if EX_API_VERS_NODOT <= 467
#error "Requires exodusII version 4.68 or later"
#endif

namespace {
  template <typename T> void clear(std::vector<T> &vec)
  {
    vec.clear();
    vec.shrink_to_fit();
    SMART_ASSERT(vec.capacity() == 0);
  }

  template <typename T> bool approx_equal(T v1, T v2)
  {
#if 0
  static const T tolerance = 100.0 * std::numeric_limits<T>::epsilon();
  return std::fabs(v1 - v2) <= std::fabs(v1+v2)*tolerance;
#else
    return (float)v1 == (float)v2;
#endif
  }
} // namespace

struct NodeInfo
{
  NodeInfo() : id(0), x(0.0), y(0.0), z(0.0) {}
  NodeInfo(size_t id_, double x_, double y_, double z_) : id(id_), x(x_), y(y_), z(z_) {}
  size_t id;
  double x, y, z;

  bool operator==(const NodeInfo &other) const
  {
    return id == other.id && approx_equal(x, other.x) && approx_equal(y, other.y) &&
           approx_equal(z, other.z);
  }

  bool operator!=(const NodeInfo &other) const { return !(*this == other); }

  bool operator<(const NodeInfo &other) const
  {
    if (id < other.id) {
      return true;
    }
    if (id > other.id) {
      return false;
    }
    SMART_ASSERT(id == other.id);
    if (!approx_equal(x, other.x) && x < other.x) {
      return true;
    }
    if (!approx_equal(x, other.x) && x > other.x) {
      return false;
    }
    SMART_ASSERT(approx_equal(x, other.x));
    if (!approx_equal(y, other.y) && y < other.y) {
      return true;
    }
    if (!approx_equal(y, other.y) && y > other.y) {
      return false;
    }
    SMART_ASSERT(approx_equal(y, other.y));
    if (!approx_equal(z, other.z) && z < other.z) {
      return true;
    }
    return false;
  }
};

using GlobalMap = std::vector<NodeInfo>;
using GMapIter  = GlobalMap::iterator;

typedef std::vector<std::pair<int, size_t>> GlobalElemMap;
using GElemMapIter = GlobalElemMap::iterator;

#include "CJ_ExodusEntity.h"
#include "CJ_ExodusFile.h"
#include "CJ_Internals.h"
#include "CJ_ObjectType.h"
#include "CJ_SystemInterface.h"
#include "CJ_Variables.h"
#include "CJ_Version.h"

extern double seacas_timer();

namespace {
  template <class T> struct TimeStepMap
  {
    TimeStepMap(size_t part, int step, T time)
        : partNumber(part), localStepNumber(step), timeValue(time)
    {
    }
    size_t partNumber;
    int    localStepNumber;
    T      timeValue;
  };

  std::string time_stamp(const std::string &format);
  std::string format_time(double seconds);
  int         get_width(size_t max_value);

  ex_entity_type exodus_object_type(Excn::ObjectType &conjoin_type)
  {
    switch (conjoin_type) {
    case Excn::EBLK: return EX_ELEM_BLOCK;
    case Excn::SSET: return EX_SIDE_SET;
    case Excn::NSET: return EX_NODE_SET;
    default:
      throw std::runtime_error("Invalid Object Type in exodus_object_type: " +
                               std::to_string(conjoin_type));
    }
  }

  char **get_name_array(int size, size_t length)
  {
    char **names = nullptr;
    if (size > 0) {
      names = new char *[size];
      for (int i = 0; i < size; i++) {
        names[i] = new char[length + 1];
        std::memset(names[i], '\0', length + 1);
      }
    }
    return names;
  }

  void free_name_array(char **names, int size)
  {
    for (int i = 0; i < size; i++) {
      delete[] names[i];
    }
    delete[] names;
    names = nullptr;
  }
} // namespace

std::string tsFormat = "[%H:%M:%S] ";

// prototypes

// The main program templated to permit float/double transfer.
template <typename T, typename INT>
int conjoin(Excn::SystemInterface &interface, T /* dummy */, INT /* dummy int */);

namespace {
  void                         compress_white_space(char *str);
  void                         add_info_record(char *info_record, int size);
  template <typename INT> void put_mesh_summary(const Excn::Mesh<INT> &mesh);

  template <typename T>
  void get_put_qa(int id, int id_out, const std::vector<TimeStepMap<T>> &global_times,
                  Excn::SystemInterface &interface);
  int  get_put_coordinate_names(int in, int out, int dimensionality);

  template <typename T, typename INT>
  int get_put_coordinates(Excn::Mesh<INT> &global, size_t part_count,
                          std::vector<Excn::Mesh<INT>> &local_mesh, T dummy);

  template <typename T, typename INT>
  int get_coordinates(int id, int dimensionality, size_t num_nodes,
                      std::vector<INT> local_node_to_global, size_t part, std::vector<T> &x,
                      std::vector<T> &y, std::vector<T> &z);

  StringVector get_exodus_variable_names(int id, ex_entity_type elType, size_t var_count);

  template <typename T, typename INT>
  void filter_truth_table(int id, Excn::Mesh<INT> &global, std::vector<T> &glob_blocks,
                          Excn::Variables &vars, const StringIdVector &variable_names);

  template <typename T, typename INT>
  void get_truth_table(Excn::Mesh<INT> &global, std::vector<std::vector<T>> &blocks,
                       std::vector<T> &glob_blocks, Excn::Variables &vars, int debug);

  template <typename U>
  void create_output_truth_table(std::vector<U> &global_sets, Excn::Variables &vars,
                                 std::vector<int> &truth_table);

  template <typename T, typename U, typename INT>
  int read_write_master_values(Excn::Variables &vars, const Excn::Mesh<INT> &global,
                               std::vector<U> &              global_sets,
                               std::vector<Excn::Mesh<INT>> &local_mesh,
                               std::vector<std::vector<U>> &local_sets, std::vector<T> &values,
                               size_t p, Excn::ExodusFile &id, int time_step, int time_step_out);

  void get_variable_params(int id, Excn::Variables &vars, const StringIdVector &variable_list);

  void get_put_variable_names(int id, int out, Excn::Variables &vars, Excn::SystemInterface &si,
                              int *combined_status_variable_index = nullptr);

  template <typename INT>
  void build_reverse_node_map(std::vector<Excn::Mesh<INT>> &local_mesh, Excn::Mesh<INT> *global,
                              size_t part_count, GlobalMap &global_node_map);

  template <typename INT>
  void build_reverse_node_map(std::vector<Excn::Mesh<INT>> &local_mesh, Excn::Mesh<INT> *global,
                              size_t part_count, std::vector<INT> &global_node_map);

  template <typename INT>
  void build_reverse_element_map(std::vector<Excn::Mesh<INT>> &         local_mesh,
                                 std::vector<std::vector<Excn::Block>> &blocks,
                                 std::vector<Excn::Block> &glob_blocks, Excn::Mesh<INT> *global,
                                 size_t part_count, GlobalElemMap &global_element_map);

  template <typename INT>
  void get_nodesets(size_t total_node_count, std::vector<Excn::Mesh<INT>> &local_mesh,
                    std::vector<std::vector<Excn::NodeSet<INT>>> &nodesets,
                    std::vector<Excn::NodeSet<INT>> &             glob_sets);

  template <typename INT>
  void get_element_blocks(const std::vector<Excn::Mesh<INT>> &   local_mesh,
                          const Excn::Mesh<INT> &                global,
                          std::vector<std::vector<Excn::Block>> &blocks,
                          std::vector<Excn::Block> &             glob_blocks);
  template <typename T, typename INT>
  void put_element_blocks(std::vector<Excn::Mesh<INT>> &         local_mesh,
                          std::vector<std::vector<Excn::Block>> &blocks,
                          std::vector<Excn::Block> &glob_blocks, T single_or_double);

  template <typename INT> void put_nodesets(std::vector<Excn::NodeSet<INT>> &glob_sets);

  template <typename INT>
  void                         get_sideset_metadata(std::vector<Excn::Mesh<INT>> &                local_mesh,
                                                    std::vector<std::vector<Excn::SideSet<INT>>> &sets,
                                                    std::vector<Excn::SideSet<INT>> &             glob_ssets);
  template <typename INT> void get_put_sidesets(std::vector<Excn::SideSet<INT>> &glob_ssets);

  template <typename T, typename INT>
  void add_status_variable(int id_out, const Excn::Mesh<INT> &global,
                           const std::vector<Excn::Block> &blocks,
                           const std::vector<Excn::Block> &glob_blocks,
                           const std::vector<INT> &local_element_to_global, int step, int variable,
                           T alive, int combined_variable_index);

  template <typename INT>
  size_t find_max_entity_count(size_t part_count, std::vector<Excn::Mesh<INT>> &local_mesh,
                               const Excn::Mesh<INT> &                       global,
                               std::vector<std::vector<Excn::Block>> &       blocks,
                               std::vector<std::vector<Excn::NodeSet<INT>>> &nodesets,
                               std::vector<std::vector<Excn::SideSet<INT>>> &sidesets);

  int case_compare(const std::string &s1, const std::string &s2);

  template <typename T>
  void verify_set_position_mapping(const std::string &type, size_t part_count,
                                   const std::vector<T> &             global_sets,
                                   const std::vector<std::vector<T>> &sets)
  {
    bool problem = false;
    for (size_t i = 0; i < global_sets.size(); i++) {
      size_t loc_pos = global_sets[i].position_;
      for (size_t p = 0; p < part_count; p++) {
        if (global_sets[i].id != sets[p][loc_pos].id ||
            case_compare(global_sets[i].name_, sets[p][loc_pos].name_) != 0 ||
            sets[p][loc_pos].position_ != i) {
          problem = true;
          std::cerr << "\nMismatch for global " << type << " at position " << i << " and local "
                    << type << " in part " << p + 1 << " position " << loc_pos << "\n";
          global_sets[i].dump();
          sets[p][loc_pos].dump();
        }
      }
    }

    if (problem) {
      throw std::runtime_error(type + " mismatch");
    }
  }

  // SEE: http://lemire.me/blog/2017/04/10/removing-duplicates-from-lists-quickly
  template <typename T> size_t unique(std::vector<T> &out)
  {
    if (out.empty())
      return 0;
    size_t pos  = 1;
    T      oldv = out[0];
    for (size_t i = 1; i < out.size(); ++i) {
      T newv   = out[i];
      out[pos] = newv;
      pos += (newv != oldv);
      oldv = newv;
    }
    return pos;
  }

  template <typename T> void uniquify(std::vector<T> &vec)
  {
    std::sort(vec.begin(), vec.end());
    vec.resize(unique(vec));
    vec.shrink_to_fit();
  }
} // namespace

unsigned int debug_level = 0;
const float  FILL_VALUE  = FLT_MAX;

int main(int argc, char *argv[])
{
#if defined(__LIBCATAMOUNT__)
  setlinebuf(stderr);
#endif
  try {
    time_t begin_time = time(nullptr);
    Excn::SystemInterface::show_version();

    Excn::SystemInterface interface;
    bool                  ok = interface.parse_options(argc, argv);

    if (!ok) {
      std::cerr << "\nERROR: Problems parsing command line arguments.\n\n";
      exit(EXIT_FAILURE);
    }

    debug_level = interface.debug();

    if ((debug_level & 64) != 0u) {
      ex_opts(EX_VERBOSE | EX_DEBUG);
    }
    else {
      ex_opts(0);
    }

    int error = 0;

    if (!Excn::ExodusFile::initialize(interface)) {
      std::cerr << "ERROR: Problem initializing input and/or output files.\n";
      exit(EXIT_FAILURE);
    }

    int int_byte_size = 4;
    if (interface.ints_64_bit()) {
      int_byte_size = 8;
    }

    if (Excn::ExodusFile::io_word_size() == 4) {
      if (int_byte_size == 4) {
        error = conjoin(interface, static_cast<float>(0.0), 0);
      }
      else {
        error = conjoin(interface, static_cast<float>(0.0), static_cast<int64_t>(0));
      }
    }
    else {
      if (int_byte_size == 4) {
        error = conjoin(interface, 0.0, 0);
      }
      else {
        error = conjoin(interface, 0.0, static_cast<int64_t>(0));
      }
    }

    Excn::ExodusFile::close_all();

    time_t end_time = time(nullptr);
    add_to_log(argv[0], end_time - begin_time);
    return (error);
  }
  catch (std::exception &e) {
    std::cerr << "ERROR: Standard exception: " << e.what() << '\n';
  }
}

template <typename T, typename INT>
int conjoin(Excn::SystemInterface &interface, T /* dummy */, INT /* dummy int */)
{
  SMART_ASSERT(sizeof(T) == Excn::ExodusFile::io_word_size());

  const T alive      = interface.alive_value();
  size_t  part_count = interface.inputFiles_.size();

  auto mytitle = new char[MAX_LINE_LENGTH + 1];
  memset(mytitle, '\0', MAX_LINE_LENGTH + 1);

  Excn::Mesh<INT> global;

  std::vector<Excn::Mesh<INT>> local_mesh(part_count);

  // ******************************************************************
  // 1. Read global info

  int error = 0;

  if (debug_level & 1) {
    std::cerr << time_stamp(tsFormat);
  }

  std::string title0;

  for (size_t p = 0; p < part_count; p++) {
    ex_init_params info{};
    error += ex_get_init_ext(Excn::ExodusFile(p), &info);

    local_mesh[p].title          = info.title;
    local_mesh[p].dimensionality = info.num_dim;
    local_mesh[p].nodeCount      = info.num_nodes;
    local_mesh[p].elementCount   = info.num_elem;
    local_mesh[p].blockCount     = info.num_elem_blk;
    local_mesh[p].nodesetCount   = info.num_node_sets;
    local_mesh[p].sidesetCount   = info.num_side_sets;

    if (p == 0) {
      global.title          = mytitle;
      global.dimensionality = local_mesh[p].count(Excn::DIM);
      global.blockCount     = local_mesh[p].count(Excn::EBLK);
    }
    else {
      SMART_ASSERT(global.count(Excn::DIM) == local_mesh[p].count(Excn::DIM))
      (global.count(Excn::DIM))(local_mesh[p].count(Excn::DIM))(p);
      SMART_ASSERT(global.count(Excn::EBLK) == local_mesh[p].count(Excn::EBLK))
      (global.count(Excn::EBLK))(local_mesh[p].count(Excn::EBLK))(p);
    }

    local_mesh[p].localNodeToGlobal.resize(local_mesh[p].count(Excn::NODE));
    local_mesh[p].localElementToGlobal.resize(local_mesh[p].count(Excn::ELEM));

  } // end for (p=0..part_count)

  delete[] mytitle;

  if (interface.omit_nodesets()) {
    global.nodesetCount = 0;
  }

  if (interface.omit_sidesets()) {
    global.sidesetCount = 0;
  }

  // Get database times...  Save mapping from time step to part providing that time step.
  std::vector<TimeStepMap<T>> global_times;

  T t_min = FLT_MAX;
  for (size_t p = part_count; p > 0; p--) {

    Excn::ExodusFile id(p - 1);
    bool             used = false;
    int              nts  = ex_inquire_int(id, EX_INQ_TIME);
    std::vector<T>   times(nts);
    ex_get_all_times(id, times.data());

    // A database will include all times from step 0 up to the
    // last time that is less than t_min on the following database.
    // Note that we are iterating through the databases from last to first...

    // Check delta between the last used timestep on this database
    // and the first time on the previous (later in time) database
    // (which is t_min). It must be > user-specified minimum delta.
    // Set by -interpart_minimum_time_delta {delta} command line option.

    // NOTE that this is not the delta between individual timesteps
    // within a database, it is the delta between the last timestep
    // on one database and the first on the next database.

    int i = 0;
    for (i = nts; i > 0; i--) {
      if (times[i - 1] < t_min) {
        if (used || t_min - times[i - 1] >= interface.interpart_minimum_time_delta()) {
          used = true;
          global_times.push_back(TimeStepMap<T>(p - 1, i - 1, times[i - 1]));
        }
      }
    }
    local_mesh[p - 1].timestepCount = i;
    t_min                           = t_min < times[0] ? t_min : times[0];
    if (!used) {
      std::string part = "Part " + std::to_string(p) + ": ";
      part += interface.inputFiles_[p - 1];
      std::cerr << "\nWARNING: " << part
                << " does not contain any time steps which will be used in conjoined file.\n"
                << " Current minimum time = " << t_min << ", timestep range on this part is "
                << times[0] << " to " << times[nts - 1] << "\n";
      local_mesh[p - 1].isActive = false;
    }
  }
  global.timestepCount = global_times.size();
  std::reverse(global_times.begin(), global_times.end());

  // Need these throughout run, so declare outside of this block...
  // TODO(gdsjaar): Add these to the "Mesh" class.
  std::vector<Excn::Block>              glob_blocks(global.count(Excn::EBLK));
  std::vector<std::vector<Excn::Block>> blocks(part_count);

  std::vector<Excn::SideSet<INT>>              glob_ssets;
  std::vector<std::vector<Excn::SideSet<INT>>> sidesets(part_count);

  std::vector<Excn::NodeSet<INT>>              glob_nsets;
  std::vector<std::vector<Excn::NodeSet<INT>>> nodesets(part_count);
  {
    // Now, build the reverse global node map which permits access of the
    // local id given the global id.
    std::vector<INT> global_node_map;
    if (interface.ignore_coordinates()) {
      build_reverse_node_map(local_mesh, &global, part_count, global_node_map);
    }
    else {
      GlobalMap global_node_coord_map;
      build_reverse_node_map(local_mesh, &global, part_count, global_node_coord_map);
      global_node_map.reserve(global_node_coord_map.size());
      for (const auto &ni : global_node_coord_map) {
        global_node_map.push_back(ni.id);
      }
    }

    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }

    // ****************************************************************************
    // 5. Get Block information including element attributes
    // must check for zero length blocks
    get_element_blocks(local_mesh, global, blocks, glob_blocks);

    GlobalElemMap global_element_map(global.count(Excn::ELEM));
    build_reverse_element_map(local_mesh, blocks, glob_blocks, &global, part_count,
                              global_element_map);

    //
    //    NOTE:  Node set/side set information can be different for each part
    /************************************************************************/
    // 7. Get Side sets
    if (!interface.omit_sidesets()) {
      if (debug_level & 1) {
        std::cerr << time_stamp(tsFormat);
      }
      get_sideset_metadata(local_mesh, sidesets, glob_ssets);
      if (global.count(Excn::SSET) != glob_ssets.size()) {
        global.sidesetCount = glob_ssets.size();
      }
    }

    /************************************************************************/
    // 6. Get Node sets
    if (!interface.omit_nodesets()) {
      if (debug_level & 1) {
        std::cerr << time_stamp(tsFormat);
      }
      get_nodesets(global.count(Excn::NODE), local_mesh, nodesets, glob_nsets);
      if (global.count(Excn::NSET) != glob_nsets.size()) {
        global.nodesetCount = glob_nsets.size();
      }
    }

    /************************************************************************/
    // Start writing the output file...

    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }
    Excn::CommunicationMetaData comm_data;

    // Create the output file...
    Excn::ExodusFile::create_output(interface);

    put_mesh_summary(global);

    get_put_qa(Excn::ExodusFile(0), Excn::ExodusFile::output(), global_times, interface);

    Excn::Internals exodus(Excn::ExodusFile::output(), Excn::ExodusFile::max_name_length());

    exodus.write_meta_data(global, glob_blocks, glob_nsets, glob_ssets, comm_data);

    // Output bulk mesh data....
    put_nodesets(glob_nsets);

    // c.2.  Write Global Node Number Map
    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }

    error = ex_put_id_map(Excn::ExodusFile::output(), EX_NODE_MAP, global_node_map.data());

    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }
    if (!global_element_map.empty()) {
      std::vector<INT> global_map(global.count(Excn::ELEM));
      for (size_t i = 0; i < global.count(Excn::ELEM); i++) {
        global_map[i] = global_element_map[i].first;
      }
      ex_put_id_map(Excn::ExodusFile::output(), EX_ELEM_MAP, global_map.data());
    }

    T dummy = 0;
    put_element_blocks(local_mesh, blocks, glob_blocks, dummy);
    get_put_sidesets(glob_ssets);
  }
  // ************************************************************************
  // 2. Get Coordinate Info.
  {
    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }

    error += get_put_coordinates(global, part_count, local_mesh, (T)0);

    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }
    std::cerr << "Wrote coordinate information...\n";
  }
  // ####################TRANSIENT DATA SECTION###########################
  // ***********************************************************************
  // 9. Get Variable Information and names

  if (debug_level & 1) {
    std::cerr << time_stamp(tsFormat);
  }

  //  I. read number of variables for each type. Note that exodusII does not
  //     provide for history variables
  //  NOTE: it is assumed that every part has the same global, nodal,
  //        and element lists

  bool            add_n_status = interface.nodal_status_variable() != "NONE";
  bool            add_e_status = interface.element_status_variable() != "NONE";
  Excn::Variables global_vars(Excn::GLOBAL);
  Excn::Variables nodal_vars(Excn::NODE, add_n_status);
  Excn::Variables element_vars(Excn::EBLK, add_e_status);
  Excn::Variables nodeset_vars(Excn::NSET);
  Excn::Variables sideset_vars(Excn::SSET);

  {
    Excn::ExodusFile id(0);

    get_variable_params(id, global_vars, interface.global_var_names());
    get_variable_params(id, nodal_vars, interface.node_var_names());
    get_variable_params(id, element_vars, interface.elem_var_names());
    get_variable_params(id, nodeset_vars, interface.nset_var_names());
    get_variable_params(id, sideset_vars, interface.sset_var_names());

    get_truth_table(global, blocks, glob_blocks, element_vars, 4);
    filter_truth_table(id, global, glob_blocks, element_vars, interface.elem_var_names());

    get_truth_table(global, nodesets, glob_nsets, nodeset_vars, 32);
    filter_truth_table(id, global, glob_nsets, nodeset_vars, interface.nset_var_names());

    get_truth_table(global, sidesets, glob_ssets, sideset_vars, 16);
    filter_truth_table(id, global, glob_ssets, sideset_vars, interface.sset_var_names());
  }

  // There is a slightly tricky situation here. The truthTable block order
  // is based on the ordering of the blocks on the input databases.
  // These blocks may have been reordered on output to make the 'offset'
  // variables line up correctly when the element ids are mapped back to
  // the "overall global" order.  There is not much problem since most
  // calls  outputting block-related items pass the id and don't rely
  // on ordering. However, the truth table is one of the exceptions
  // and we need to reorder the truth table to match the output block
  // order. After this call, we can use the original ordering, so just
  // need a temporary vector here...
  if (global_vars.count(Excn::OUT) + nodal_vars.count(Excn::OUT) + element_vars.count(Excn::OUT) +
          nodeset_vars.count(Excn::OUT) + sideset_vars.count(Excn::OUT) >
      0) {

    std::vector<int> elem_truth_table;
    std::vector<int> nset_truth_table;
    std::vector<int> sset_truth_table;
    create_output_truth_table(glob_blocks, element_vars, elem_truth_table);
    create_output_truth_table(glob_nsets, nodeset_vars, nset_truth_table);
    create_output_truth_table(glob_ssets, sideset_vars, sset_truth_table);

    ex_put_all_var_param(
        Excn::ExodusFile::output(), global_vars.count(Excn::OUT), nodal_vars.count(Excn::OUT),
        element_vars.count(Excn::OUT), elem_truth_table.data(), nodeset_vars.count(Excn::OUT),
        nset_truth_table.data(), sideset_vars.count(Excn::OUT), sset_truth_table.data());
  }

  // II. read/write the variable names
  int combined_status_variable_index = 0;
  {
    Excn::ExodusFile id(0);
    get_put_variable_names(id, Excn::ExodusFile::output(), global_vars, interface);
    get_put_variable_names(id, Excn::ExodusFile::output(), nodal_vars, interface);
    get_put_variable_names(id, Excn::ExodusFile::output(), element_vars, interface,
                           &combined_status_variable_index);
    get_put_variable_names(id, Excn::ExodusFile::output(), nodeset_vars, interface);
    get_put_variable_names(id, Excn::ExodusFile::output(), sideset_vars, interface);
  }
  ex_update(Excn::ExodusFile::output());

  /**********************************************************************/
  // 10. Get Transient Data
  //     This routine reads in a time dump from an EXODUSII file

  size_t num_time_steps = global.count(Excn::TIME);

  if (debug_level & 1) {
    std::cerr << time_stamp(tsFormat);
  }

  std::vector<T> global_values(global_vars.count(Excn::IN));
  std::vector<T> output_global_values(global_vars.count(Excn::OUT));

  // Determine maximum number of entities on any part
  int max_ent = find_max_entity_count(part_count, local_mesh, global, blocks, nodesets, sidesets);
  std::vector<T> values(max_ent);

  // Stage II.  Extracting transient variable data.
  //            loop over time steps

  // Determine if user wants a subset of timesteps transferred to the output file.
  // Time steps for output file
  double start_time = seacas_timer();

  // Used for output formatting at end of loop.
  int percent_width = 0;
  int field_width   = 3;
  if (num_time_steps > 100) {
    percent_width = 1;
    field_width   = 5;
  }

  int element_width  = get_width(global.count(Excn::ELEM));
  int step_width     = get_width(num_time_steps);
  int part_width     = get_width(part_count + 1);
  int loc_step_width = 0;
  for (auto &global_time : global_times) {
    if (loc_step_width < global_time.localStepNumber + 1) {
      loc_step_width = global_time.localStepNumber + 1;
    }
  }
  loc_step_width = get_width(loc_step_width);

  size_t time_step_out = 0;
  for (size_t time_step = 0; time_step < num_time_steps; time_step++) {
    time_step_out++;

    T                time_val = 0;
    size_t           p        = global_times[time_step].partNumber;
    Excn::ExodusFile id(p);

    // read in and write out the time step information
    error += ex_get_time(id, global_times[time_step].localStepNumber + 1, &time_val);
    SMART_ASSERT(time_val == global_times[time_step].timeValue)
    (time_step)(time_val)(global_times[time_step].timeValue);
    error += ex_put_time(Excn::ExodusFile::output(), time_step_out, &time_val);

    if (global_vars.count(Excn::OUT) > 0) {
      if (debug_level & 1) {
        std::cerr << time_stamp(tsFormat) << "Global Variables...\n";
      }
      error += ex_get_var(id, global_times[time_step].localStepNumber + 1, EX_GLOBAL, 0, 0,
                          global_vars.count(), (void *)global_values.data());
      // Map ...
      for (int ig = 0; ig < global_vars.count(Excn::IN); ig++) {
        if (global_vars.index_[ig] > 0) {
          SMART_ASSERT(ig < (int)global_values.size());
          output_global_values[global_vars.index_[ig] - 1] = global_values[ig];
        }
      }
      error += ex_put_var(Excn::ExodusFile::output(), time_step_out, EX_GLOBAL, 1, 0,
                          global_vars.count(Excn::OUT), output_global_values.data());
    }

    // ========================================================================
    // Nodal Values...
    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat) << "Nodal Variables...\n";
    }

    if (nodal_vars.count(Excn::OUT) > 0) {
      size_t         node_count = local_mesh[p].count(Excn::NODE);
      std::vector<T> master_nodal_values(global.count(Excn::NODE));

      int offset = nodal_vars.addStatus ? 1 : 0;
      for (int i = 0; i < nodal_vars.count(Excn::OUT) - offset;
           i++) { // Last output variable may be status
        for (int j = 0; j < nodal_vars.count(Excn::IN); j++) {
          if (nodal_vars.index_[j] - 1 == i) {
            std::fill(master_nodal_values.begin(), master_nodal_values.end(), 0.0);

            error += ex_get_var(id, global_times[time_step].localStepNumber + 1, EX_NODAL, j + 1, 0,
                                node_count, values.data());

            // copy values to master nodal value information
            for (size_t jj = 0; jj < node_count; jj++) {
              // Map local nodal value to global location...
              size_t nodal_value               = local_mesh[p].localNodeToGlobal[jj];
              master_nodal_values[nodal_value] = values[jj];
            }
            error += ex_put_var(Excn::ExodusFile::output(), time_step_out, EX_NODAL, i + 1, 0,
                                global.count(Excn::NODE), master_nodal_values.data());
            break;
          }
        }

        // Fill "node_status" variable -- 'alive' for alive; 1-alive for dead.
        // It is the last output variable...
        if (nodal_vars.addStatus) {
          SMART_ASSERT(alive == 0.0 || alive == 1.0)(alive);
          std::fill(master_nodal_values.begin(), master_nodal_values.end(), (1.0 - alive));
          for (size_t j = 0; j < node_count; j++) {
            // Map local nodal value to global location...
            size_t nodal_value               = local_mesh[p].localNodeToGlobal[j];
            master_nodal_values[nodal_value] = alive;
          }

          error += ex_put_var(Excn::ExodusFile::output(), time_step_out, EX_NODAL,
                              nodal_vars.count(Excn::OUT), 0, global.count(Excn::NODE),
                              master_nodal_values.data());
        }
      }
    }

    // ========================================================================
    // Extracting element transient variable data
    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat) << "Element Variables...\n";
    }
    if (element_vars.count(Excn::IN) > 0) {
      read_write_master_values(element_vars, global, glob_blocks, local_mesh, blocks, values, p, id,
                               global_times[time_step].localStepNumber, time_step_out);
    }

    // Add element status variable...
    // Use the output time step for writing data
    if (interface.element_status_variable() != "NONE") {
      add_status_variable(Excn::ExodusFile::output(), global, blocks[p], glob_blocks,
                          local_mesh[p].localElementToGlobal, time_step_out,
                          element_vars.index_[element_vars.count(Excn::IN)], alive,
                          combined_status_variable_index);
    }

    // ========================================================================
    // Extracting sideset transient variable data
    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat) << "Sideset Variables...\n";
    }
    if (sideset_vars.count(Excn::IN) > 0) {
      read_write_master_values(sideset_vars, global, glob_ssets, local_mesh, sidesets, values, p,
                               id, global_times[time_step].localStepNumber, time_step_out);
    }

    // ========================================================================
    // Extracting nodeset transient variable data
    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat) << "Nodeset Variables...\n";
    }
    if (nodeset_vars.count(Excn::IN) > 0) {
      read_write_master_values(nodeset_vars, global, glob_nsets, local_mesh, nodesets, values, p,
                               id, global_times[time_step].localStepNumber, time_step_out);
    }

    // ========================================================================
    ex_update(Excn::ExodusFile::output());

    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }

    std::ios::fmtflags f(std::cerr.flags());
    std::cerr << "Step " << std::setw(step_width) << time_step + 1 << "/" << num_time_steps
              << ", time " << std::scientific << std::setprecision(4) << time_val << " (Part "
              << std::setw(part_width) << p + 1 << "/" << part_count << ", step "
              << std::setw(loc_step_width) << global_times[time_step].localStepNumber + 1 << ")"
              << "  Active Elem: " << std::setw(element_width) << local_mesh[p].count(Excn::ELEM);

    double cur_time            = seacas_timer();
    double elapsed             = cur_time - start_time;
    double time_per_step       = elapsed / time_step_out;
    double percentage_done     = (time_step_out * 100.0) / global.count(Excn::TIME);
    double estimated_remaining = time_per_step * (global.count(Excn::TIME) - time_step_out);

    std::cerr << "  [" << std::fixed << std::setw(field_width) << std::setprecision(percent_width)
              << percentage_done << "%, Elapsed=" << format_time(elapsed)
              << ", ETA=" << format_time(estimated_remaining) << "]\n";
    if (debug_level & 1) {
      std::cerr << "\n";
    }
    std::cerr.flags(f);
  }

  /*************************************************************************/
  // EXIT program
  if (debug_level & 1) {
    std::cerr << time_stamp(tsFormat);
  }
  std::cerr << "******* END *******\n";
  return (error);
}

namespace {
  template <typename T>
  void get_put_qa(int id, int id_out, const std::vector<TimeStepMap<T>> &global_times,
                  Excn::SystemInterface &interface)
  {
    // NOTE: Assuming info and QA records for all parts
    int error = 0;

    int info_string_len = MAX_LINE_LENGTH;

    size_t num_info_records = ex_inquire_int(id, EX_INQ_INFO);
    //    size_t extra_info = global_times.size() + 2 + 1;
    size_t extra_info = 2 + 1;

    char **info_records = get_name_array(num_info_records + extra_info, info_string_len);

    if (num_info_records > 0) {
      error += ex_get_info(id, info_records);
    }

    // Add an info record for CONJOIN
    add_info_record(info_records[num_info_records], MAX_LINE_LENGTH);

#if 0
    // Add time/part mapping...
    for (size_t i=0; i < global_times.size(); i++) {
      std::ostringstream os;
      os << "Step " << std::setw(2) << i+1 << ", time "
	 << std::scientific << std::setprecision(4) << global_times[i].timeValue
	 << " (Part " << global_times[i].partNumber+1
	 << ", step "  << global_times[i].localStepNumber+1 << ")  File: "
	 << interface.inputFiles_[global_times[i].partNumber];
      strncpy(info_records[num_info_records+1+i], os.str().c_str(), MAX_LINE_LENGTH);
      info_records[num_info_records+1+i][MAX_LINE_LENGTH] = '\0';
    }
#endif

    error += ex_put_info(id_out, num_info_records + extra_info, info_records);

    free_name_array(info_records, num_info_records + extra_info);

    // II. Get and store QA records, if they exist
    struct qa_element
    {
      char *qa_record[1][4];
    };

    int  num_qa_records = ex_inquire_int(id, EX_INQ_QA);
    auto qaRecord       = new qa_element[num_qa_records + 1];
    for (int i = 0; i < num_qa_records + 1; i++) {
      for (int j = 0; j < 4; j++) {
        qaRecord[i].qa_record[0][j]    = new char[MAX_STR_LENGTH + 1];
        qaRecord[i].qa_record[0][j][0] = '\0';
      }
    }
    if (num_qa_records) {
      error += ex_get_qa(id, qaRecord[0].qa_record);
    }

    char buffer[MAX_STR_LENGTH + 1];

    strncpy(qaRecord[num_qa_records].qa_record[0][0], qainfo[0], MAX_STR_LENGTH); // Code
    strncpy(qaRecord[num_qa_records].qa_record[0][1], qainfo[2], MAX_STR_LENGTH); // Version

    time_t date_time = time(nullptr);
    strftime(buffer, MAX_STR_LENGTH, "%Y/%m/%d", localtime(&date_time));

    strncpy(qaRecord[num_qa_records].qa_record[0][2], buffer, MAX_STR_LENGTH);

    strftime(buffer, MAX_STR_LENGTH, "%H:%M:%S", localtime(&date_time));
    strncpy(qaRecord[num_qa_records].qa_record[0][3], buffer, MAX_STR_LENGTH);

    error += ex_put_qa(id_out, num_qa_records + 1, qaRecord[0].qa_record);

    for (int i = 0; i < num_qa_records + 1; i++) {
      for (int j = 0; j < 4; j++) {
        delete[] qaRecord[i].qa_record[0][j];
      }
    }
    delete[] qaRecord;
  }

  template <typename T, typename INT>
  int get_put_coordinates(Excn::Mesh<INT> &global, size_t part_count,
                          std::vector<Excn::Mesh<INT>> &local_mesh, T /* dummy */)
  {
    SMART_ASSERT(sizeof(T) == Excn::ExodusFile::io_word_size());
    std::vector<T> x(global.count(Excn::NODE));
    std::vector<T> y(global.count(Excn::NODE));
    std::vector<T> z(global.count(Excn::NODE));

    if (debug_level & 8) {
      std::fill(x.begin(), x.end(), FILL_VALUE);
      std::fill(y.begin(), y.end(), FILL_VALUE);
      std::fill(z.begin(), z.end(), FILL_VALUE);
    }

    int error = 0;
    for (size_t p = 0; p < part_count; p++) {
      error += get_coordinates(Excn::ExodusFile(p), global.count(Excn::DIM),
                               local_mesh[p].count(Excn::NODE), local_mesh[p].localNodeToGlobal, p,
                               x, y, z);
    } // end for p=0..part_count

    // Get Coordinate Names
    // NOTE: Assuming coordinate names should be
    // the same for all files/parts therefore, only one
    // file/part needs to be loaded
    error += get_put_coordinate_names(Excn::ExodusFile(0), Excn::ExodusFile::output(),
                                      global.count(Excn::DIM));
    // Write out coordinate information
    error += ex_put_coord(Excn::ExodusFile::output(), x.data(), y.data(), z.data());
    return error;
  }

  int get_put_coordinate_names(int in, int out, int dimensionality)
  {
    int    error            = 0;
    char **coordinate_names = get_name_array(dimensionality, Excn::ExodusFile::max_name_length());

    error += ex_get_coord_names(in, coordinate_names);
    error += ex_put_coord_names(out, coordinate_names);
    std::cerr << "Wrote coordinate names..." << '\n';

    free_name_array(coordinate_names, dimensionality);
    return error;
  }

  template <typename T, typename INT>
  int get_coordinates(int id, int dimensionality, size_t num_nodes,
                      std::vector<INT> local_node_to_global, size_t part, std::vector<T> &x,
                      std::vector<T> &y, std::vector<T> &z)
  {
    SMART_ASSERT(sizeof(T) == Excn::ExodusFile::io_word_size());
    SMART_ASSERT(local_node_to_global.size() == num_nodes);
    int            error = 0;
    std::vector<T> local_x(num_nodes);
    std::vector<T> local_y(num_nodes);
    std::vector<T> local_z(num_nodes);

    error += ex_get_coord(id, local_x.data(), local_y.data(), local_z.data());

    // Check for 2D or 3D coordinates
    if (dimensionality == 3) {
      if (debug_level & 8) {
        for (size_t i = 0; i < num_nodes; i++) {
          INT node = local_node_to_global[i];
          if (x[node] != FILL_VALUE && y[node] != FILL_VALUE && z[node] != FILL_VALUE) {
            if (!approx_equal(x[node], local_x[i]) || !approx_equal(y[node], local_y[i]) ||
                !approx_equal(z[node], local_z[i])) {
              std::ios::fmtflags f(std::cerr.flags());
              std::cerr << "\nWARNING: Node " << node + 1
                        << " has different coordinates in at least two parts.\n"
                        << "         this may indicate that this id has been reused in the current "
                           "part.\n"
                        << "         cur value = " << std::scientific << std::setprecision(6)
                        << std::setw(14) << x[node] << std::setw(14) << y[node] << std::setw(14)
                        << z[node] << "\n"
                        << "         new value = " << std::setw(14) << local_x[i] << std::setw(14)
                        << local_y[i] << std::setw(14) << local_z[i] << " from part " << part
                        << '\n';

              std::cerr.flags(f);
            }
          }
        }
      }

      for (size_t i = 0; i < num_nodes; i++) {
        // The following WILL overwrite x[node],y[node],z[node] if node is the
        // same for different parts
        INT node = local_node_to_global[i];
        x[node]  = local_x[i];
        y[node]  = local_y[i];
        z[node]  = local_z[i];
      }
    }
    else {
      if (debug_level & 8) {
        for (size_t i = 0; i < num_nodes; i++) {
          INT node = local_node_to_global[i];
          if (x[node] != FILL_VALUE && y[node] != FILL_VALUE) {
            if (!approx_equal(x[node], local_x[i]) || !approx_equal(y[node], local_y[i])) {
              std::ios::fmtflags f(std::cerr.flags());
              std::cerr << "\nWARNING: Node " << node + 1
                        << " has different coordinates in at least two parts.\n"
                        << "         this may indicate that this id has been reused in the current "
                           "part.\n"
                        << "         cur value = " << std::scientific << std::setprecision(6)
                        << std::setw(14) << x[node] << std::setw(14) << y[node] << "\n"
                        << "         new value = " << std::setw(14) << local_x[i] << std::setw(14)
                        << local_y[i] << " from part " << part << '\n';
              std::cerr.flags(f);
            }
          }
        }
      }

      for (size_t i = 0; i < num_nodes; i++) {
        // The following WILL overwrite x[node],y[node] if node is the same for
        // different parts
        INT node = local_node_to_global[i];
        x[node]  = local_x[i];
        y[node]  = local_y[i];
      }
    }
    return error;
  }

  template <typename INT>
  void get_element_blocks(const std::vector<Excn::Mesh<INT>> &   local_mesh,
                          const Excn::Mesh<INT> &                global,
                          std::vector<std::vector<Excn::Block>> &blocks,
                          std::vector<Excn::Block> &             glob_blocks)
  {
    size_t part_count = local_mesh.size();
    for (size_t ip = 0; ip < part_count; ip++) {
      blocks[ip].resize(local_mesh[ip].count(Excn::EBLK));
    }

    std::vector<INT> block_id(global.count(Excn::EBLK));

    int error = 0;
    for (size_t p = 0; p < part_count; p++) {
      Excn::ExodusFile id(p);

      error += ex_get_ids(id, EX_ELEM_BLOCK, block_id.data());

      if (p == 0) {
        for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
          glob_blocks[b].id = block_id[b];
        }
      }

      for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
        ex_block block_param{};
        block_param.id   = glob_blocks[b].id;
        block_param.type = EX_ELEM_BLOCK;

        error += ex_get_block_param(id, &block_param);

        std::vector<char> name(Excn::ExodusFile::max_name_length() + 1);
        ex_get_name(id, EX_ELEM_BLOCK, glob_blocks[b].id, name.data());

        blocks[p][b].id = glob_blocks[b].id;
        if (name[0] != '\0') {
          blocks[p][b].name_ = name.data();
          if (p == 0) {
            glob_blocks[b].name_ = name.data();
          }
        }

        // Find position of this block id in the local mesh
        for (size_t lb = 0; lb < global.count(Excn::EBLK); lb++) {
          if (block_id[lb] == glob_blocks[b].id) {
            blocks[p][b].position_ = lb;
            break;
          }
        }

        if (block_param.num_entry > 0) {
          blocks[p][b].elementCount    = block_param.num_entry;
          blocks[p][b].nodesPerElement = block_param.num_nodes_per_entry;
          blocks[p][b].attributeCount  = block_param.num_attribute;
          blocks[p][b].offset_         = block_param.num_entry;
          strncpy(blocks[p][b].elType, block_param.topology, MAX_STR_LENGTH);

          // NOTE: This is not correct, but fixed below.
          glob_blocks[b].elementCount += block_param.num_entry;

          if (glob_blocks[b].nodesPerElement == 0) {
            glob_blocks[b].nodesPerElement = block_param.num_nodes_per_entry;
          }

          if (glob_blocks[b].attributeCount == 0) {
            glob_blocks[b].attributeCount = block_param.num_attribute;
          }

          glob_blocks[b].position_ = b;
          strncpy(glob_blocks[b].elType, block_param.topology, MAX_STR_LENGTH);
        }

        if (block_param.num_attribute > 0 && glob_blocks[b].attributeNames.empty()) {
          // Get attribute names.  Assume the same on all parts
          // on which the block exists.
          char **names =
              get_name_array(block_param.num_attribute, Excn::ExodusFile::max_name_length());

          ex_get_attr_names(id, EX_ELEM_BLOCK, block_id[b], names);
          for (int i = 0; i < block_param.num_attribute; i++) {
            glob_blocks[b].attributeNames.emplace_back(names[i]);
          }
          free_name_array(names, block_param.num_attribute);
        }
      }
    } // end for p=0..part_count

    // Convert block_offset from elements/block/part to true offset
    for (size_t p = 0; p < part_count; p++) {
      if (debug_level & 4) {
        std::cerr << "\nElement block info for part " << p << "..." << '\n';
      }
      // Number of elements per block in local block order...
      std::vector<size_t> local_order_entity_count(global.count(Excn::EBLK));
      for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
        int local_order                       = blocks[p][b].position_;
        local_order_entity_count[local_order] = blocks[p][b].entity_count();
      }
      size_t sum = 0;
      for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
        size_t save                 = local_order_entity_count[b];
        local_order_entity_count[b] = sum;
        sum += save;
      }
      for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
        int local_order      = blocks[p][b].position_;
        blocks[p][b].offset_ = local_order_entity_count[local_order];

        if (debug_level & 4) {
          std::cerr << "\tBlock " << b << ", Id = " << glob_blocks[b].id;
          std::cerr << ", Name = '" << blocks[p][b].name_;
          std::cerr << "', Elements = " << std::setw(8) << blocks[p][b].entity_count();
          std::cerr << ", Nodes/element = " << blocks[p][b].nodesPerElement;
          std::cerr << ", Attributes = " << blocks[p][b].attributeCount;
          std::cerr << ", Position = " << blocks[p][b].position_;
          std::cerr << ", Offset = " << blocks[p][b].offset_ << '\n';
        }
      }
    }
  }

  template <typename T, typename INT>
  void put_element_blocks(std::vector<Excn::Mesh<INT>> &         local_mesh,
                          std::vector<std::vector<Excn::Block>> &blocks,
                          std::vector<Excn::Block> &             glob_blocks, T /* dummy */)
  {
    SMART_ASSERT(sizeof(T) == Excn::ExodusFile::io_word_size());
    int global_num_blocks = glob_blocks.size();

    if (debug_level & 1) {
      std::cerr << time_stamp(tsFormat);
    }
    std::cerr << "\nReading and Writing element connectivity & attributes\n";

    for (int b = 0; b < global_num_blocks; b++) {

      if (debug_level & 4) {
        std::cerr << "\nOutput element block info for...\n";
        std::cerr << "Block " << b << ", Id = " << glob_blocks[b].id;
        std::cerr << ", Name = '" << glob_blocks[b].name_;
        std::cerr << "', Elements = " << std::setw(8) << glob_blocks[b].entity_count();
        std::cerr << ", Nodes/element = " << glob_blocks[b].nodesPerElement;
        std::cerr << ", Attributes = " << glob_blocks[b].attributeCount << '\n';
      }

      if (debug_level & 4) {
        std::cerr << "B" << b << ":\t" << std::flush;
      }

      size_t max_nodes = glob_blocks[b].entity_count();
      max_nodes *= glob_blocks[b].nodesPerElement;
      std::vector<INT> block_linkage(max_nodes);

      // Initialize attributes list, if it exists
      std::vector<T> attributes(glob_blocks[b].attributeCount * glob_blocks[b].entity_count());

      int    error      = 0;
      size_t part_count = local_mesh.size();
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);

        size_t global_pos;
        size_t global_block_pos;

        if (blocks[p][b].entity_count() > 0) { // non-zero length block

          if (debug_level & 4) {
            std::cerr << "#" << std::flush;
          }
          size_t maximum_nodes = blocks[p][b].entity_count();
          maximum_nodes *= blocks[p][b].nodesPerElement;
          std::vector<INT> local_linkage(maximum_nodes);

          INT bid = blocks[p][b].id;
          error   = ex_get_conn(id, EX_ELEM_BLOCK, bid, local_linkage.data(), nullptr, nullptr);
          if (error < 0) {
            std::cerr << "ERROR: Cannot get element block connectivity for block " << bid
                      << " on part " << p << ".\n";
          }
          size_t pos                     = 0;
          size_t goffset                 = glob_blocks[b].offset_;
          size_t element_count           = blocks[p][b].entity_count();
          size_t boffset                 = blocks[p][b].offset_;
          size_t npe                     = blocks[p][b].nodesPerElement;
          INT *  part_loc_elem_to_global = local_mesh[p].localElementToGlobal.data();
          INT *  part_loc_node_to_global = local_mesh[p].localNodeToGlobal.data();

          for (size_t e = 0; e < element_count; e++) {
            global_block_pos = part_loc_elem_to_global[(e + boffset)] - goffset;
            global_pos       = global_block_pos * npe;

            for (size_t n = 0; n < npe; n++) {
              size_t node = part_loc_node_to_global[local_linkage[pos++] - 1];
              if (debug_level & 4) {
                SMART_ASSERT(block_linkage[global_pos] == (int)node + 1 ||
                             block_linkage[global_pos] == 0);
              }
              block_linkage[global_pos++] = node + 1;
            }
          }

          // Get attributes list,  if it exists
          if (blocks[p][b].attributeCount > 0) {
            SMART_ASSERT(blocks[p][b].attributeCount == glob_blocks[b].attributeCount)
            (p)(b)(blocks[p][b].attributeCount)(glob_blocks[b].attributeCount);

            size_t         max_attr = blocks[p][b].entity_count() * blocks[p][b].attributeCount;
            std::vector<T> local_attr(max_attr);

            error += ex_get_attr(id, EX_ELEM_BLOCK, blocks[p][b].id, local_attr.data());

            pos = 0;

            size_t att_count = blocks[p][b].attributeCount;
            for (size_t e = 0; e < element_count; e++) {
              // global_pos is global position within this element block...
              global_block_pos = local_mesh[p].localElementToGlobal[(e + boffset)] - goffset;
              global_pos       = global_block_pos * att_count;
              for (size_t n = 0; n < att_count; n++) {
                attributes[global_pos++] = local_attr[pos++];
              }
            }
          }
        } // end if blocks[p][b].entity_count() (non-zero length block)
        else if (debug_level & 4) {
          std::cerr << ".";
        }
      } // end for p=0..part_count-1

      // Verify that connectivity has been set for all elements...
      for (size_t i = 0; i < block_linkage.size(); i++) {
        SMART_ASSERT(block_linkage[i] > 0)(block_linkage[i])(i)(b)(glob_blocks[b].id);
      }

      // Write out block info
      int id_out = Excn::ExodusFile::output(); // output file identifier

      if (!block_linkage.empty()) {
        error += ex_put_conn(id_out, EX_ELEM_BLOCK, glob_blocks[b].id, block_linkage.data(),
                             nullptr, nullptr);
      }

      // Write out attributes list if it exists
      if (glob_blocks[b].attributeCount > 0) {
        error += ex_put_attr(id_out, EX_ELEM_BLOCK, glob_blocks[b].id, attributes.data());
      } // end for b=0..global_num_blocks-1
      if (debug_level & 4) {
        std::cerr << '\n';
      }
    }
    std::cerr << '\n';
  }

  template <typename INT>
  void build_reverse_element_map(std::vector<Excn::Mesh<INT>> &         local_mesh,
                                 std::vector<std::vector<Excn::Block>> &blocks,
                                 std::vector<Excn::Block> &glob_blocks, Excn::Mesh<INT> *global,
                                 size_t part_count, GlobalElemMap &global_element_map)
  {
    int error = 0;
    // Create the map that maps from a local part element to the
    // global map. This combines the mapping local part element to
    // 'global id' and then 'global id' to global position. The
    // mapping is now a direct lookup instead of a lookup followed by
    // a reverse map.
    //

    // We iterate through each element block a part at a time to
    // determine which elements are in each block. Note that an
    // element will possibly exist in multiple parts, so we need to do
    // the sort/uniquify/shrink on each element block.  The operations
    // for each element block will be:
    // 1. get elements for the block for each part into a vector.
    // 2. sort/uniquify/shrink that vector and copy into the
    // 'global_element_map' in the positions following the previous block.
    //
    // This is similar to what was done above in a global sense,
    // except it is done an element block at a time. The only need for
    // doing them both is that it lets us detect whether an element id
    // was reused and appears in both blocks.  Currently, just error
    // out if that happens...  In the future, change the id to an
    // unused value and continue....

    size_t                     tot_size = 0;
    std::vector<GlobalElemMap> global_element_numbers(part_count);
    for (size_t p = 0; p < part_count; p++) {
      Excn::ExodusFile id(p);
      global_element_numbers[p].resize(local_mesh[p].count(Excn::ELEM));
      std::vector<INT> ids(local_mesh[p].count(Excn::ELEM));
      error += ex_get_id_map(id, EX_ELEM_MAP, ids.data());
      for (size_t i = 0; i < local_mesh[p].count(Excn::ELEM); i++) {
        global_element_numbers[p][i] = std::make_pair(ids[i], size_t(0));
      }
      tot_size += local_mesh[p].count(Excn::ELEM);
    }
    global_element_map.resize(tot_size);

    size_t goffset = 0;
    for (auto &glob_block : glob_blocks) {

      size_t block_size = 0;
      for (size_t p = 0; p < part_count; p++) {
        size_t lb = 0;
        for (; lb < blocks[0].size(); lb++) {
          if (glob_block.id == blocks[p][lb].id) {
            break;
          }
        }
        SMART_ASSERT(glob_block.id == blocks[p][lb].id);
        block_size += blocks[p][lb].entity_count();
      }
      GlobalElemMap block_element_map(block_size);

      size_t poffset = 0;
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);

        size_t lb = 0;
        for (; lb < blocks[0].size(); lb++) {
          if (glob_block.id == blocks[p][lb].id) {
            break;
          }
        }
        SMART_ASSERT(glob_block.id == blocks[p][lb].id);

        // Get connectivity array for this element block...
        size_t element_count  = blocks[p][lb].entity_count();
        size_t nodes_per_elem = blocks[p][lb].nodesPerElement;
        size_t maximum_nodes  = element_count * nodes_per_elem;

        std::vector<INT> local_linkage(maximum_nodes);

        size_t bid = blocks[p][lb].id;
        error      = ex_get_conn(id, EX_ELEM_BLOCK, bid, local_linkage.data(), nullptr, nullptr);
        if (error < 0) {
          std::cerr << "ERROR: Cannot get element block connectivity for block " << bid
                    << " on part " << p << ".\n";
        }

        // Convert connectivity to global node numbers.
        for (size_t i = 0; i < element_count * nodes_per_elem; i++) {
          INT local_node   = local_linkage[i];
          INT global_node  = local_mesh[p].localNodeToGlobal[local_node - 1] + 1;
          local_linkage[i] = global_node;
        }

        // Have element global ids for all elements in this block,
        // and connectivity.  Can now create out "eleminfo" for these elements.
        size_t con_offset = 0;
        size_t boffset    = blocks[p][lb].offset_;
        for (size_t i = 0; i < element_count; i++) {
          size_t adler_crc = adler(0, &local_linkage[con_offset], nodes_per_elem * sizeof(int));
          global_element_numbers[p][boffset + i].second = adler_crc;
          block_element_map[poffset + i]                = global_element_numbers[p][boffset + i];
          con_offset += nodes_per_elem;
        }
        poffset += element_count;
      }

      // Sort, uniquify, shrink 'block_element_map' and the result
      // then contains the list of elements in this block...
      uniquify(block_element_map);

      size_t block_total_num_elements = block_element_map.size();
      glob_block.elementCount         = block_total_num_elements;
      glob_block.offset_              = goffset;

      // Copy into the global_element_map...
      std::copy(block_element_map.begin(), block_element_map.end(), &global_element_map[goffset]);
      goffset += block_total_num_elements;
    }

    global->elementCount = goffset;
    global_element_map.resize(goffset);

    size_t max_id        = global_element_map[global->elementCount - 1].first;
    bool   is_contiguous = max_id == global_element_map.size();
    // std::cerr  << "Element id map " << (is_contiguous ? "is" : "is not") << " contiguous.\n";

    // The global_element_map may or may not be globally sorted; however, each
    // block is sorted, so if we do the iteration by blocks, we can
    // use lower_bound instead of doing global searches...
    for (auto &glob_block : glob_blocks) {

      auto gm_begin = global_element_map.begin() + glob_block.offset_;
      auto gm_end   = gm_begin + glob_block.elementCount;
      auto cur_pos  = gm_begin;
      for (size_t p = 0; p < part_count; p++) {
        size_t lb = 0;
        for (; lb < blocks[0].size(); lb++) {
          if (glob_block.id == blocks[p][lb].id) {
            break;
          }
        }
        SMART_ASSERT(glob_block.id == blocks[p][lb].id);

        size_t element_count = blocks[p][lb].entity_count();
        size_t boffset       = blocks[p][lb].offset_;
        for (size_t i = 0; i < element_count; i++) {
          std::pair<int, size_t> global_element = global_element_numbers[p][boffset + i];

          if (cur_pos == gm_end || *cur_pos != global_element) {
            auto iter = std::lower_bound(gm_begin, gm_end, global_element);
            SMART_ASSERT(iter != gm_end);
            cur_pos = iter;
          }
          size_t element_value                            = cur_pos - gm_begin;
          local_mesh[p].localElementToGlobal[i + boffset] = element_value + glob_block.offset_;
          ++cur_pos;
        }
      }
    }

    // Update the element ids to give a unique, non-repeating set.  If
    // contiguous, then there is nothing to do.  If not contiguous,
    // then need to determine if there are any repeats (id reuse) and
    // if so, generate a new id for the repeated uses.  Note that
    // there is a possibility that elements in two or more element
    // blocks will have the same element id, so we generate a vector
    // containing a pair<id, position_in_global_element_map>, sort it,
    // and then use it to detect duplicates and map them to a new id.

    if (!is_contiguous) {
      std::vector<std::pair<size_t, size_t>> id_pos(global_element_map.size());
      for (size_t i = 0; i < global->elementCount; i++) {
        id_pos[i].first  = global_element_map[i].first;
        id_pos[i].second = i;
      }
      std::sort(id_pos.begin(), id_pos.end());

      max_id = id_pos.back().first;
      // Check again for contiguous ids since we now have a sorted list...
      is_contiguous = max_id == global_element_map.size();

      if (!is_contiguous) {
        int    repeat_found = 0;
        size_t id_last      = id_pos[0].first;

        for (size_t i = 1; i < global->elementCount; i++) {
          if (id_pos[i].first == id_last) {
            global_element_map[id_pos[i].second].first = ++max_id;
            repeat_found++;
          }
          else {
            id_last = id_pos[i].first;
          }
        }
        if (repeat_found > 0) {
          std::cerr << repeat_found
                    << " duplicate element ids were found. Their ids have been "
                       "renumbered to remove duplicates.\n";
        }
      }
    }
  }

  template <typename INT>
  void build_reverse_node_map(std::vector<Excn::Mesh<INT>> &local_mesh, Excn::Mesh<INT> *global,
                              size_t part_count, GlobalMap &global_node_map)
  {
    // Append all local node maps to the global node map.
    // Sort the global node map
    // Remove duplicates.
    // Position within map is now the map...
    // When building the local-part node to global id, use binary_search...

    // Global node map and count.
    std::vector<std::vector<NodeInfo>> global_nodes(part_count);

    size_t tot_size = 0;
    for (size_t p = 0; p < part_count; p++) {
      tot_size += local_mesh[p].count(Excn::NODE);
      global_nodes[p].resize(local_mesh[p].count(Excn::NODE));
    }
    global_node_map.resize(tot_size);

    size_t offset = 0;
    for (size_t p = 0; p < part_count; p++) {
      std::vector<double> x(local_mesh[p].count(Excn::NODE));
      std::vector<double> y(local_mesh[p].count(Excn::NODE));
      std::vector<double> z(local_mesh[p].count(Excn::NODE));
      std::vector<INT>    nid(local_mesh[p].count(Excn::NODE));

      Excn::ExodusFile id(p);
      ex_get_id_map(id, EX_NODE_MAP, nid.data());
      ex_get_coord(id, x.data(), y.data(), z.data());
      for (size_t i = 0; i < local_mesh[p].count(Excn::NODE); i++) {
        global_nodes[p][i] = NodeInfo(nid[i], x[i], y[i], z[i]);
      }
      std::copy(global_nodes[p].begin(), global_nodes[p].end(), &global_node_map[offset]);
      offset += local_mesh[p].count(Excn::NODE);
    }
    // Now, sort the global_node_map array and remove duplicates...
    uniquify(global_node_map);

    global->nodeCount = global_node_map.size();

    // See whether the node numbers are contiguous.  If so, we can map
    // the nodes back to their original location. Since the nodes are
    // sorted and there are no duplicates, we just need to see if the id
    // at global_node_map.size() == global_node_map.size();
    INT  max_id        = global_node_map[global->nodeCount - 1].id;
    bool is_contiguous = (int64_t)max_id == static_cast<int64_t>(global_node_map.size());
    std::cerr << "Node map " << (is_contiguous ? "is" : "is not") << " contiguous.\n";

    // Create the map that maps from a local part node to the
    // global map. This combines the mapping local part node to
    // 'global id' and then 'global id' to global position. The
    // mapping is now a direct lookup instead of a lookup followed by
    // a reverse map.
    if (is_contiguous) {
      for (size_t p = 0; p < part_count; p++) {
        size_t node_count = local_mesh[p].count(Excn::NODE);
        for (size_t i = 0; i < node_count; i++) {
          const NodeInfo &global_node        = global_nodes[p][i];
          local_mesh[p].localNodeToGlobal[i] = global_node.id - 1;
        }
      }
    }
    else {
      auto cur_pos = global_node_map.begin();
      for (size_t p = 0; p < part_count; p++) {
        size_t node_count = local_mesh[p].count(Excn::NODE);
        for (size_t i = 0; i < node_count; i++) {
          NodeInfo global_node = global_nodes[p][i];

          if (cur_pos == global_node_map.end() || *cur_pos != global_node) {
            auto iter =
                std::lower_bound(global_node_map.begin(), global_node_map.end(), global_node);
            if (iter == global_node_map.end()) {
              NodeInfo n = global_node;
              std::cerr << "Bad Node in build_reverse_node_map: " << n.id
                        << "\tat location: " << n.x << "\t" << n.y << "\t" << n.z << "\n";
              exit(EXIT_FAILURE);
            }
            cur_pos = iter;
          }
          size_t nodal_value                 = cur_pos - global_node_map.begin();
          local_mesh[p].localNodeToGlobal[i] = nodal_value;
          ++cur_pos;
        }
      }
    }

    // Update the nodal ids to give a unique, non-repeating set.  If contiguous, then
    // there is nothing to do.  If not contiguous, then need to determine if there are any
    // repeats (id reuse) and if so, generate a new id for the repeated uses.
    // A duplicate id would have the same id, but different x y z position.
    if (!is_contiguous) {
      bool   repeat_found = false;
      size_t id_last      = global_node_map[0].id;
      for (size_t i = 1; i < global->nodeCount; i++) {
        if (global_node_map[i].id == id_last) {
          global_node_map[i].id = ++max_id;
          repeat_found          = true;
        }
        else {
          id_last = global_node_map[i].id;
        }
      }
      if (repeat_found) {
        std::cerr << "Duplicate node ids were found. Their ids have been renumbered to remove "
                     "duplicates. If the part meshes should be identical, maybe use the "
                     "--ignore_coordinate option.\n";
      }
    }
  }

  template <typename INT>
  void build_reverse_node_map(std::vector<Excn::Mesh<INT>> &local_mesh, Excn::Mesh<INT> *global,
                              size_t part_count, std::vector<INT> &global_node_map)
  {
    // Append all local node maps to the global node map.
    // Sort the global node map
    // Remove duplicates.
    // Position within map is now the map...
    // When building the local-part node to global id, use binary_search...

    // Global node map and count.
    std::vector<std::vector<INT>> global_nodes(part_count);

    size_t tot_size = 0;
    for (size_t p = 0; p < part_count; p++) {
      tot_size += local_mesh[p].count(Excn::NODE);
      global_nodes[p].resize(local_mesh[p].count(Excn::NODE));
    }
    global_node_map.resize(tot_size);

    size_t offset = 0;
    for (size_t p = 0; p < part_count; p++) {
      Excn::ExodusFile id(p);
      ex_get_id_map(id, EX_NODE_MAP, global_nodes[p].data());
      std::copy(global_nodes[p].begin(), global_nodes[p].end(), &global_node_map[offset]);
      offset += local_mesh[p].count(Excn::NODE);
    }

    // Now, sort the global_node_map array and remove duplicates...
    uniquify(global_node_map);

    global->nodeCount = global_node_map.size();

    // See whether the node numbers are contiguous.  If so, we can map
    // the nodes back to their original location. Since the nodes are
    // sorted and there are no duplicates, we just need to see if the id
    // at global_node_map.size() == global_node_map.size();
    INT  max_id        = global_node_map[global->nodeCount - 1];
    bool is_contiguous = (int64_t)max_id == static_cast<int64_t>(global_node_map.size());
    std::cerr << "Node map " << (is_contiguous ? "is" : "is not") << " contiguous.\n";

    // Create the map that maps from a local part node to the
    // global map. This combines the mapping local part node to
    // 'global id' and then 'global id' to global position. The
    // mapping is now a direct lookup instead of a lookup followed by
    // a reverse map.
    if (is_contiguous) {
      for (size_t p = 0; p < part_count; p++) {
        size_t node_count = local_mesh[p].count(Excn::NODE);
        for (size_t i = 0; i < node_count; i++) {
          INT global_node                    = global_nodes[p][i];
          local_mesh[p].localNodeToGlobal[i] = global_node - 1;
        }
      }
    }
    else {
      auto cur_pos = global_node_map.begin();
      for (size_t p = 0; p < part_count; p++) {
        size_t node_count = local_mesh[p].count(Excn::NODE);
        for (size_t i = 0; i < node_count; i++) {
          INT global_node = global_nodes[p][i];

          if (cur_pos == global_node_map.end() || *cur_pos != global_node) {
            auto iter =
                std::lower_bound(global_node_map.begin(), global_node_map.end(), global_node);
            if (iter == global_node_map.end()) {
              INT n = global_node;
              std::cerr << "Bad Node in build_reverse_node_map: " << n << "\n";
              exit(EXIT_FAILURE);
            }
            cur_pos = iter;
          }
          size_t nodal_value                 = cur_pos - global_node_map.begin();
          local_mesh[p].localNodeToGlobal[i] = nodal_value;
          ++cur_pos;
        }
      }
    }
  }

  void get_put_variable_names(int id, int out, Excn::Variables &vars, Excn::SystemInterface &si,
                              int *combined_status_variable_index)
  {
    if (vars.count(Excn::OUT) > 0) {

      char **output_name_list =
          get_name_array(vars.count(Excn::OUT), Excn::ExodusFile::max_name_length());

      int num_vars       = vars.index_.size();
      int extra          = vars.addStatus ? 1 : 0;
      int num_input_vars = num_vars - extra;

      char **input_name_list = get_name_array(num_vars, Excn::ExodusFile::max_name_length());
      if (num_input_vars > 0) {
        int error = ex_get_variable_names(id, vars.type(), num_input_vars, input_name_list);
        if (error != EX_NOERR) {
          std::cerr << "ERROR: Cannot get " << vars.label() << " variable names\n";
          exit(EXIT_FAILURE);
        }
      }

      std::string status;
      if (vars.type() == EX_ELEM_BLOCK || vars.type() == EX_NODAL) {
        if (vars.type() == EX_ELEM_BLOCK) {
          status = si.element_status_variable();
          if (status != "NONE") {
            strcpy(input_name_list[num_vars - 1], status.c_str());
          }
        }
        else if (vars.type() == EX_NODAL) {
          status = si.nodal_status_variable();
          if (status != "NONE") {
            strcpy(input_name_list[num_vars - 1], status.c_str());
          }
        }
      }

      // Iterate through the 'var_index' and transfer
      // Assume that the number of pointers is limited to
      // the number of results variables
      size_t maxlen = 0;
      for (int i = 0; i < num_vars; i++) {
        if (vars.index_[i] > 0) {
          strcpy(output_name_list[vars.index_[i] - 1], input_name_list[i]);
          if (strlen(input_name_list[i]) > maxlen) {
            maxlen = strlen(input_name_list[i]);
          }
        }
      }

      // See if any of the variable names conflict with the status variable name...
      if (status != "NONE") {
        for (size_t i = 0; i < static_cast<size_t>(vars.count(Excn::OUT)) - 1; i++) {
          if (case_compare(output_name_list[i], status) == 0) {
            // Error -- duplicate element variable names on output database.
            std::cerr << "\nERROR: A " << vars.label()
                      << " variable already exists on the input database with the "
                      << "same name as the status variable '" << status
                      << "'. This is not allowed.\n\n";
            exit(EXIT_FAILURE);
          }
        }
      }

      maxlen += 2;
      // Assume 8 characters for initial tab...
      int width  = si.screen_width();
      int nfield = (width - 8) / maxlen;
      if (nfield < 1) {
        nfield = 1;
      }

      std::cerr << "Found " << vars.count(Excn::OUT) << " " << vars.label() << " variables.\n";
      {
        int i    = 0;
        int ifld = 1;
        std::cerr << "\t";
        std::ios::fmtflags f(std::cerr.flags());
        while (i < vars.count(Excn::OUT)) {
          std::cerr << std::setw(maxlen) << std::left << output_name_list[i++];
          if (++ifld > nfield && i < vars.count(Excn::OUT)) {
            std::cerr << "\n\t";
            ifld = 1;
          }
        }
        std::cerr << "\n\n";
        std::cerr.flags(f); // Reset back to what it was.
      }

      ex_put_variable_names(out, vars.type(), vars.count(Excn::OUT), output_name_list);

      // KLUGE: Handle finding combined status index variable here since it is
      // the only place we have the list of output variable names...
      if (vars.type() == EX_ELEM_BLOCK && (combined_status_variable_index != nullptr)) {
        *combined_status_variable_index = 0;
        const std::string &comb_stat    = si.combined_mesh_status_variable();
        if (!comb_stat.empty()) {
          for (int i = 0; i < vars.count(Excn::OUT); i++) {
            if (case_compare(comb_stat, output_name_list[i]) == 0) {
              *combined_status_variable_index = i + 1;
              break;
            }
          }
        }
      }
      free_name_array(output_name_list, vars.count(Excn::OUT));
      free_name_array(input_name_list, num_vars);
    }
  }

  void get_variable_params(int id, Excn::Variables &vars, const StringIdVector &variable_list)
  {
    // Determines the number of variables of type 'type()' that will
    // be written to the output database. The 'variable_list' vector
    // specifies a possibly empty list of variable names that the user
    // wants transferred to the output database. If 'variable_list' is
    // empty, then all variables of that type will be transferred; if
    // the 'variable_list' size is 1 and it contains the string 'NONE',
    // then no variables of that type will be transferred; if size is 1
    // and it contains the string 'ALL', then all variables of that type
    // will be transferred.
    //
    // Returns the number of variables which will be output Also creates
    // a 'var_index'.  The var_index is zero-based and of size
    // 'input_variable_count'. If:
    // var_index[i] ==0, variable not written to output database
    // var_index[i] > 0, variable written; is variable 'var_index[i]'

    // If 'type' is ELEMENT or NODE, then reserve space for the 'status' variable.
    int extra = vars.addStatus ? 1 : 0;
    int num_vars;
    ex_get_variable_param(id, vars.type(), &num_vars);

    vars.index_.resize(num_vars + extra);

    // Create initial index which defaults to no output...
    std::fill(vars.index_.begin(), vars.index_.end(), 0);
    // ...Except for the status variable (if any)
    if (extra == 1) {
      vars.index_[num_vars] = 1;
    }

    // If 'variable_list' is empty or specified 'ALL', then all
    // variables are to be output
    if (variable_list.empty() ||
        (variable_list.size() == 1 && case_compare(variable_list[0].first, "all") == 0)) {
      for (size_t i = 0; i < vars.index_.size(); i++) {
        vars.index_[i] = i + 1;
      }
      vars.outputCount = num_vars + extra;
      return;
    }

    // Another possibility is user specifies "NONE" for the variable
    // list so no variables will be written.  Just return 0.
    if (variable_list.size() == 1 && case_compare(variable_list[0].first, "none") == 0) {
      vars.outputCount = extra;
      return;
    }

    // At this point, the variable_list specifies at least one
    // variable to be output to the database.
    // Get the list that the user entered and the list of files
    // from the input database...
    {
      StringVector exo_names = get_exodus_variable_names(id, vars.type(), num_vars);

      // Iterate 'variable_list' and find position in 'exo_names'.  If
      // not found, there is an error -- set index to -1; otherwise set
      // index to position (1-based) in variable_list.  Others are
      // already set to '0' by default initialization.

      // The variable_list may contain multiple entries for each
      // variable if the user is specifying output only on certain
      // element blocks...
      std::string var_name;
      int         var_count = 0;
      for (auto &elem : variable_list) {
        if (var_name == elem.first) {
          continue;
        }
        var_name   = elem.first;
        bool found = false;
        for (size_t j = 0; j < exo_names.size() && !found; j++) {
          if (case_compare(exo_names[j], var_name) == 0) {
            found          = true;
            vars.index_[j] = ++var_count;
          }
        }
        if (!found) {
          std::cerr << "ERROR: Variable '" << elem.first << "' is not valid." << '\n';
          exit(EXIT_FAILURE);
        }
      }
      // Count non-zero entries in var_index;
      int nz_count = 0;
      for (auto &elem : vars.index_) {
        if (elem > 0) {
          nz_count++;
        }
      }
      SMART_ASSERT(nz_count == var_count + extra)(nz_count)(var_count);

      if (vars.addStatus) {
        vars.index_[num_vars] = nz_count; // Already counted above...
      }
      vars.outputCount = nz_count;
      return;
    }
  }

  template <typename INT> void put_mesh_summary(const Excn::Mesh<INT> &mesh)
  {
    // Write out Mesh info

    std::cout << " Title: " << mesh.title << "\n\n";
    std::cout << " Number of coordinates per node       =" << std::setw(9) << mesh.count(Excn::DIM)
              << "\n";
    std::cout << " Number of nodes                      =" << std::setw(9) << mesh.count(Excn::NODE)
              << "\n";
    std::cout << " Number of elements                   =" << std::setw(9) << mesh.count(Excn::ELEM)
              << "\n";
    std::cout << " Number of element blocks             =" << std::setw(9) << mesh.count(Excn::EBLK)
              << "\n\n";
    std::cout << " Number of nodal point sets           =" << std::setw(9) << mesh.count(Excn::NSET)
              << "\n";
    std::cout << " Number of element side sets          =" << std::setw(9) << mesh.count(Excn::SSET)
              << "\n\n";
  }

  template <typename INT>
  void get_nodesets(size_t total_node_count, std::vector<Excn::Mesh<INT>> &local_mesh,
                    std::vector<std::vector<Excn::NodeSet<INT>>> &nodesets,
                    std::vector<Excn::NodeSet<INT>> &             glob_sets)
  {
    // Find number of nodesets in the global model...
    std::set<INT>    set_ids;
    std::vector<INT> ids;

    size_t part_count = local_mesh.size();
    int    bad_ns     = 0;
    {
      int ns_count;
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);
        ns_count = ex_inquire_int(id, EX_INQ_NODE_SETS);

        // Get the ids for these
        ids.resize(ns_count);
        ex_get_ids(id, EX_NODE_SET, ids.data());

        for (int iset = 0; iset < ns_count; iset++) {
          if (ids[iset] != 0) {
            set_ids.insert(ids[iset]);
          }
          else {
            bad_ns++;
          }
        }
      }
    }

    if (bad_ns != 0) {
      std::cerr << "ERROR: There were " << bad_ns
                << " nodesets (counting all files) which had an id equal to "
                << "0 which is not allowed.\n";
    }

    if (set_ids.empty()) {
      return;
    }

    // set_ids now contains the union of all nodeset ids...
    glob_sets.resize(set_ids.size());
    size_t gnset_size = set_ids.size();
    ids.resize(gnset_size);

    {
      size_t i = 0;
      for (auto set_id : set_ids) {
        glob_sets[i].id        = set_id;
        glob_sets[i].position_ = i;
        i++;
      }
    }

    {
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);

        nodesets[p].resize(set_ids.size());

        // Get the ids again so we can map current order back to file order...
        ex_get_ids(id, EX_NODE_SET, ids.data());
        int ns_count = ex_inquire_int(id, EX_INQ_NODE_SETS);

        for (int i = 0; i < ns_count; i++) {
          nodesets[p][i].id = ids[i];

          // Find which global nodeset this id corresponds to...
          size_t gi = gnset_size;
          for (size_t j = 0; j < gnset_size; j++) {
            if (ids[i] == glob_sets[j].id) {
              gi = j;
              break;
            }
          }
          SMART_ASSERT(gi != gnset_size);
          glob_sets[gi].position_  = i;
          nodesets[p][i].position_ = gi;

          // Get the parameters for this nodeset...
          if (ex_int64_status(id) & EX_BULK_INT64_API) {
            int64_t node_count, df_count;
            ex_get_set_param(id, EX_NODE_SET, nodesets[p][i].id, &node_count, &df_count);
            nodesets[p][i].nodeCount = node_count;
            nodesets[p][i].dfCount   = df_count;
          }
          else {
            int node_count, df_count;
            ex_get_set_param(id, EX_NODE_SET, nodesets[p][i].id, &node_count, &df_count);
            nodesets[p][i].nodeCount = node_count;
            nodesets[p][i].dfCount   = df_count;
          }

          std::vector<char> name(Excn::ExodusFile::max_name_length() + 1);
          ex_get_name(id, EX_NODE_SET, nodesets[p][i].id, name.data());
          if (name[0] != '\0') {
            nodesets[p][i].name_ = name.data();
            if (glob_sets[gi].name_.empty()) {
              glob_sets[gi].name_ = name.data();
            }
          }

          if (debug_level & 32) {
            std::cerr << "Part " << p + 1 << " ";
            nodesets[p][i].dump();
          }
        }
      }
    }

    verify_set_position_mapping("nodeset", part_count, glob_sets, nodesets);

    {
      // Now get the nodeset nodes and df.
      // Currently ignore the DF.  Could add if people need it...

      // This is inefficient since the part loop is on
      // the inside...  The other ordering would use more memory...

      std::vector<INT> ns_nodes;
      for (size_t ns = 0; ns < set_ids.size(); ns++) {

        std::vector<INT> glob_ns_nodes(total_node_count + 1);
        std::fill(glob_ns_nodes.begin(), glob_ns_nodes.end(), 0);

        size_t lns = glob_sets[ns].position_;
        for (size_t p = 0; p < part_count; p++) {
          Excn::ExodusFile id(p);

          glob_sets[ns].name_ = nodesets[p][lns].name_;

          size_t size = nodesets[p][lns].entity_count();
          if (size > 0) {
            ns_nodes.resize(size);
            ex_get_set(id, EX_NODE_SET, nodesets[p][lns].id, ns_nodes.data(), nullptr);

            // The node ids are in local space -- map to global
            for (size_t iset = 0; iset < size; iset++) {
              size_t global_node         = local_mesh[p].localNodeToGlobal[ns_nodes[iset] - 1] + 1;
              glob_ns_nodes[global_node] = 1;
            }
          }
        }
        // Count number of nonzero entries and transfer to the
        // output nodeset
        // NOTE: global_node above is 1-based.
        glob_sets[ns].nodeCount = std::accumulate(glob_ns_nodes.begin(), glob_ns_nodes.end(), 0);
        glob_sets[ns].nodeSetNodes.resize(glob_sets[ns].entity_count());
        glob_sets[ns].dfCount = 0;

        size_t j = 0;
        for (size_t i = 1; i <= total_node_count; i++) {
          if (glob_ns_nodes[i] == 1) {
            glob_sets[ns].nodeSetNodes[j++] = i;
            glob_ns_nodes[i]                = j;
          }
        }
        SMART_ASSERT(j == glob_sets[ns].entity_count())(j)(glob_sets[ns].entity_count())(ns);

        // See if we need a map from local nodeset position to global nodeset position
        // Only needed if there are nodeset variables (or dist factors).
        // Assume all files have same number of variables...
        int num_vars;
        {
          Excn::ExodusFile id(0);
          ex_get_variable_param(id, EX_NODE_SET, &num_vars);
        }
        if (num_vars > 0) {
          for (size_t p = 0; p < part_count; p++) {
            Excn::ExodusFile id(p);
            // Get the nodelist, but store it in nodeOrderMap.
            // global_pos = nodeOrderMap[i]
            Excn::NodeSet<INT> &nset   = nodesets[p][ns];
            size_t              nnodes = nset.entity_count();
            nset.nodeOrderMap.resize(nnodes);
            ex_get_set(id, EX_NODE_SET, nset.id, nset.nodeOrderMap.data(), nullptr);

            for (size_t i = 0; i < nnodes; i++) {
              size_t local_node    = nset.nodeOrderMap[i];                                // 1-based
              size_t global_node   = local_mesh[p].localNodeToGlobal[local_node - 1] + 1; // 1-based
              size_t global_pos    = glob_ns_nodes[global_node];                          // 1-based
              nset.nodeOrderMap[i] = global_pos - 1;
            }
#if 0
	    if (debug_level & 32)
	      nset.dump_order();
#endif
          }
        }
      }
    }
  }

  template <typename INT> void put_nodesets(std::vector<Excn::NodeSet<INT>> &glob_sets)
  {
    int exoid = Excn::ExodusFile::output();

    if (debug_level & 32) {
      std::cerr << "\nOutput NodeSets:\n";
    }
    for (auto &glob_set : glob_sets) {
      ex_put_set(exoid, EX_NODE_SET, glob_set.id, glob_set.nodeSetNodes.data(), nullptr);
      //    ex_put_node_set_dist_fact(exoid, glob_sets[ns].id, &glob_sets[ns].distFactors[0]);
      // Done with the memory; clear out the vector containing the bulk data nodes and distFactors.
      clear(glob_set.nodeSetNodes);
      clear(glob_set.distFactors);

      if (debug_level & 32) {
        glob_set.dump();
      }
    }
  }

  template <typename INT>
  void get_sideset_metadata(std::vector<Excn::Mesh<INT>> &                local_mesh,
                            std::vector<std::vector<Excn::SideSet<INT>>> &sets,
                            std::vector<Excn::SideSet<INT>> &             glob_ssets)
  {
    // Find number of sidesets in the global model...
    std::set<int>    set_ids;
    std::vector<INT> ids;

    size_t part_count = local_mesh.size();
    int    bad_ss     = 0;
    {
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);
        int              ss_count = ex_inquire_int(id, EX_INQ_SIDE_SETS);

        // Get the ids for these
        ids.resize(ss_count);
        ex_get_ids(id, EX_SIDE_SET, ids.data());

        for (int i = 0; i < ss_count; i++) {
          if (ids[i] != 0) {
            set_ids.insert(ids[i]);
          }
          else {
            bad_ss++;
          }
        }
      }
    }

    if (bad_ss != 0) {
      std::cerr
          << "ERROR: There were " << bad_ss
          << " sidesets (counting all files) which had an id equal to 0 which is not allowed.\n";
    }

    if (set_ids.empty()) {
      return;
    }

    // set_ids now contains the union of all sideset ids...
    glob_ssets.resize(set_ids.size());
    size_t gsset_size = set_ids.size();
    ids.resize(gsset_size);

    {
      size_t i = 0;
      for (auto set_id : set_ids) {
        glob_ssets[i].id        = set_id;
        glob_ssets[i].position_ = i;
        i++;
      }
    }

    {
      std::vector<char> name(Excn::ExodusFile::max_name_length() + 1);
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);

        sets[p].resize(set_ids.size());

        // Get the ids again so we can map current order back to file order...
        ex_get_ids(id, EX_SIDE_SET, ids.data());

        int ss_count = ex_inquire_int(id, EX_INQ_SIDE_SETS);
        for (int i = 0; i < ss_count; i++) {
          sets[p][i].id = ids[i];

          // Find which global sideset this id corresponds to...
          size_t gi = gsset_size;
          for (size_t j = 0; j < gsset_size; j++) {
            if (ids[i] == glob_ssets[j].id) {
              gi = j;
              break;
            }
          }
          SMART_ASSERT(gi != gsset_size);
          glob_ssets[gi].position_ = i;
          sets[p][i].position_     = gi;

          // Get the parameters for this sideset...
          ex_set set{};
          set.type                     = EX_SIDE_SET;
          set.id                       = sets[p][i].id;
          set.entry_list               = nullptr;
          set.extra_list               = nullptr;
          set.distribution_factor_list = nullptr;
          int error                    = ex_get_sets(id, 1, &set);
          if (error != EX_NOERR) {
            std::cerr << "ERROR: Cannot get side set with id " << set.id << "\n";
            exit(EXIT_FAILURE);
          }

          sets[p][i].sideCount = set.num_entry;
          sets[p][i].dfCount   = set.num_distribution_factor;
          glob_ssets[gi].sideCount += sets[p][i].entity_count();
          glob_ssets[gi].dfCount += sets[p][i].dfCount;

          ex_get_name(id, EX_SIDE_SET, sets[p][i].id, name.data());
          if (name[0] != '\0') {
            sets[p][i].name_     = name.data();
            glob_ssets[gi].name_ = name.data();
          }
        }
      }

      verify_set_position_mapping("sideset", part_count, glob_ssets, sets);

      // See if we need a map from local sideset position to global sideset position
      // Only needed if there are sideeset variables (or dist factors which are currently ignored).
      // Assume all files have same number of variables...
      bool need_sideset_map = false;
      {
        int              num_vars;
        Excn::ExodusFile id(0);
        ex_get_variable_param(id, EX_SIDE_SET, &num_vars);
        need_sideset_map = num_vars > 0;
      }

      // Now get the sideset elements/sides. Ignoring DF for now...
      for (size_t ss = 0; ss < set_ids.size(); ss++) {
        // This is maximum possible size; will probably be reduced by duplicate elimination...
        typedef std::vector<std::pair<INT, INT>> ElemSideMap;
        ElemSideMap                              elem_side(glob_ssets[ss].sideCount);
        int                                      ss_id  = glob_ssets[ss].id;
        size_t                                   offset = 0;

        size_t lss = glob_ssets[ss].position_;
        for (size_t p = 0; p < part_count; p++) {

          Excn::ExodusFile id(p);
          sets[p][lss].elems.resize(sets[p][lss].sideCount);
          sets[p][lss].sides.resize(sets[p][lss].sideCount);
          ex_get_set(id, EX_SIDE_SET, ss_id, sets[p][lss].elems.data(), sets[p][lss].sides.data());

          // Add these to the elem_side vector...
          for (size_t i = 0; i < sets[p][lss].sideCount; i++) {
            size_t global_elem = local_mesh[p].localElementToGlobal[sets[p][lss].elems[i] - 1] + 1;
            elem_side[offset + i] = std::make_pair((INT)global_elem, (INT)sets[p][lss].sides[i]);
          }
          offset += sets[p][lss].sideCount;
        }

        uniquify(elem_side);

        // Set the output sideset definition...
        glob_ssets[ss].sideCount = elem_side.size();
        glob_ssets[ss].dfCount   = 0;
        glob_ssets[ss].elems.resize(elem_side.size());
        glob_ssets[ss].sides.resize(elem_side.size());

        // Populate the global sideset elements and sides...
        for (size_t i = 0; i < elem_side.size(); i++) {
          std::tie(glob_ssets[ss].elems[i], glob_ssets[ss].sides[i]) = elem_side[i];
        }

        if (need_sideset_map) {
          // For this sideset in each part, figure out the mapping for
          // its (elem, side, variable) position into the corresponding
          // global sideset position...

          // Try the lower_bound searching of elem_side for now.  If
          // inefficient, fix later...
          for (size_t p = 0; p < part_count; p++) {
            sets[p][lss].elemOrderMap.resize(sets[p][lss].sideCount);
            for (size_t i = 0; i < sets[p][lss].sideCount; i++) {
              size_t global_elem =
                  local_mesh[p].localElementToGlobal[sets[p][lss].elems[i] - 1] + 1;
              std::pair<INT, INT> es = std::make_pair((INT)global_elem, (INT)sets[p][lss].sides[i]);

              auto   iter = std::lower_bound(elem_side.begin(), elem_side.end(), es);
              size_t pos  = iter - elem_side.begin();
              sets[p][lss].elemOrderMap[i] = pos;
            }
          }
        }
      }

      // Calculate sideset offset
      for (size_t b = 0; b < glob_ssets.size(); b++) {
        size_t sum = 0;
        for (size_t p = 0; p < part_count; p++) {
          sets[p][b].offset_ = sum;
          sum += sets[p][b].entity_count();

          if (debug_level & 16) {
            std::cerr << "Part " << p + 1 << " ";
            sets[p][b].dump();
          }
        }
      }

      // Free some memory which is no longer needed...
      // (Could move up into sideset loop above)
      for (size_t p = 0; p < part_count; p++) {
        for (auto &elem : sets[p]) {
          clear(elem.elems);
          clear(elem.sides);
          clear(elem.distFactors);
        }
      }
    }
  }

  template <typename INT> void get_put_sidesets(std::vector<Excn::SideSet<INT>> &glob_ssets)
  {
    int exoid = Excn::ExodusFile::output(); // output file identifier
    for (auto &glob_sset : glob_ssets) {
      ex_put_set(exoid, EX_SIDE_SET, glob_sset.id, const_cast<INT *>(glob_sset.elems.data()),
                 const_cast<INT *>(&glob_sset.sides[0]));
      if (glob_sset.dfCount > 0) {
        ex_put_set_dist_fact(exoid, EX_SIDE_SET, glob_sset.id,
                             reinterpret_cast<void *>(glob_sset.distFactors.data()));
      }
    }

    for (auto &glob_sset : glob_ssets) {
      clear(glob_sset.elems);
      clear(glob_sset.sides);
      clear(glob_sset.distFactors);
    }
  }

  template <typename T, typename INT>
  void add_status_variable(int id_out, const Excn::Mesh<INT> &global,
                           const std::vector<Excn::Block> &blocks,
                           const std::vector<Excn::Block> &glob_blocks,
                           const std::vector<INT> &local_element_to_global, int step, int variable,
                           T alive, int combined_variable_index)
  {
    SMART_ASSERT(sizeof(T) == Excn::ExodusFile::io_word_size());
    std::vector<T> status;
    SMART_ASSERT(alive == 0.0 || alive == 1.0)(alive);

    for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
      status.resize(glob_blocks[b].entity_count());
      std::fill(status.begin(), status.end(), (1.0 - alive));
      size_t boffset       = blocks[b].offset_;
      size_t goffset       = glob_blocks[b].offset_;
      size_t element_count = blocks[b].entity_count();
      for (size_t e = 0; e < element_count; e++) {
        size_t global_block_pos  = local_element_to_global[(e + boffset)] - goffset;
        status[global_block_pos] = alive;
      }

      // If combining this status variable with a mesh status
      // variable, do that now...
      if (combined_variable_index > 0) {
        std::vector<T> mesh_status(glob_blocks[b].entity_count());
        ex_get_var(id_out, step, EX_ELEM_BLOCK, combined_variable_index, glob_blocks[b].id,
                   glob_blocks[b].entity_count(), mesh_status.data());

        if (alive == 1.0) {
          for (size_t i = 0; i < glob_blocks[b].entity_count(); i++) {
            status[i] = status[i] * mesh_status[i];
          }
        }
        else {
          SMART_ASSERT(alive == 0.0);
          for (size_t i = 0; i < glob_blocks[b].entity_count(); i++) {
            status[i] = 1.0 - ((1.0 - status[i]) * (1.0 - mesh_status[i]));
          }
        }
      }
      ex_put_var(id_out, step, EX_ELEM_BLOCK, variable, glob_blocks[b].id,
                 glob_blocks[b].entity_count(), status.data());
    }
  }

  StringVector get_exodus_variable_names(int id, ex_entity_type elType, size_t var_count)
  {
    // Allocate space for variable names...
    char **name_list = get_name_array(var_count, Excn::ExodusFile::max_name_length());
    int    error     = ex_get_variable_names(id, elType, var_count, name_list);
    if (error != EX_NOERR) {
      std::cerr << "ERROR: Cannot get variable names\n";
      exit(EXIT_FAILURE);
    }

    StringVector names(var_count);
    for (size_t j = 0; j < var_count; j++) {
      compress_white_space(name_list[j]);
      names[j] = std::string(name_list[j]);
    }

    free_name_array(name_list, var_count);
    return names;
  }

  template <typename T, typename INT>
  void filter_truth_table(int id, Excn::Mesh<INT> &global, std::vector<T> &glob_blocks,
                          Excn::Variables &vars, const StringIdVector &variable_names)
  {
    // This routine checks the 'variable_names' list to see if the user
    // has restricted the output of certain variables to certain element
    // blocks. If so, then the truth table is modified to match the
    // users request.
    if (variable_names.empty()) {
      return;
    }

    // Check for a non-zero id entry in the variable_names list which
    // will signify a user-specified element block.
    bool found_it = false;
    for (size_t i = 0; i < variable_names.size() && !found_it; i++) {
      if (variable_names[i].second > 0) {
        found_it = true;
      }
    }
    if (!found_it) {
      return;
    }

    // At this point, know that there is at least one block-restricted
    // variable output specification. For each one, modify the global
    // truth table to match the specification.
    StringVector exo_names = get_exodus_variable_names(id, vars.type(), vars.count());

    std::string var_name;
    int         out_position = -1;
    for (auto &variable_name : variable_names) {
      if (variable_name.second > 0) {
        if (var_name != variable_name.first) {
          var_name = variable_name.first;
          // Find which exodus variable matches this name
          out_position = -1;
          for (size_t j = 0; j < exo_names.size(); j++) {
            if (case_compare(exo_names[j], var_name) == 0) {
              out_position = vars.index_[j] - 1;
              break;
            }
          }
          if (out_position < 0) {
            std::cerr << "ERROR: Variable '" << variable_name.first
                      << "' does not exist on any block in this database.\n";
            exit(EXIT_FAILURE);
          }

          // Set all truth table entries for this variable to negative
          // of current value and then iterate over specified blocks and
          // set those positive.  This way can make sure that the
          // variable truly exists for the block that the user specified.
          found_it = false;
          for (size_t b = 0; b < global.count(vars.objectType); b++) {
            if (glob_blocks[b].id == variable_name.second) {
              if (glob_blocks[b].truthTable[out_position] == 0) {
                std::cerr << "ERROR: Variable '" << variable_name.first
                          << "' does not exist on block " << variable_name.second << ".\n";
                exit(EXIT_FAILURE);
              }
              else {
                found_it = true;
                break;
              }
            }
            else {
              // User does not want this variable output on these blocks...
              glob_blocks[b].truthTable[out_position] = 0;
            }
          }

          if (!found_it) {
            std::cerr << "ERROR: User-specified block id of " << variable_name.second
                      << " for variable '" << variable_name.first << "' does not exist.\n";
            exit(EXIT_FAILURE);
          }
        }
      }
    }
  }

  template <typename T, typename INT>
  void get_truth_table(Excn::Mesh<INT> &global, std::vector<std::vector<T>> &blocks,
                       std::vector<T> &glob_blocks, Excn::Variables &vars, int debug)
  {
    // read truth table - sum across all parts since many will
    // have values set to zero for zero length blocks the element
    // variable truth table organized as a 2D array:
    // [global.count(Excn::EBLK)][num_elem_vars]

    Excn::ObjectType object_type = vars.objectType;

    if (vars.count(Excn::OUT) > 0) {
      // For each input exodus file, get it's truth table and fill
      // in the location in the output truth table...

      size_t part_count = blocks.size();
      for (size_t p = 0; p < part_count; p++) {
        Excn::ExodusFile id(p);

        for (size_t b = 0; b < global.count(object_type); b++) {
          size_t gb = 0;
          for (; gb < global.count(object_type); gb++) {
            if (glob_blocks[gb].id == blocks[p][b].id) {
              break;
            }
          }
          SMART_ASSERT(glob_blocks[gb].id == blocks[p][b].id);

          if (p == 0) {
            glob_blocks[gb].truthTable.resize(vars.count(Excn::OUT));
          }

          if (vars.count(Excn::IN) > 0) {
            blocks[p][b].truthTable.resize(vars.count(Excn::IN));
            ex_get_object_truth_vector(id, vars.type(), blocks[p][b].id, vars.count(Excn::IN),
                                       blocks[p][b].truthTable.data());

            // Find global block corresponding to this block. (Ids match)
            for (int j = 0; j < vars.count(Excn::IN); j++) {
              if (vars.index_[j] > 0) {
                glob_blocks[gb].truthTable[vars.index_[j] - 1] += blocks[p][b].truthTable[j];
              }
            }
          }
          if (vars.addStatus) {
            glob_blocks[gb].truthTable[vars.count(Excn::OUT) - 1] = 1;
          }
        }
      }

      // reset truth table values that may be greater than 1
      for (size_t b = 0; b < global.count(object_type); b++) {
        for (int j = 0; j < vars.count(Excn::OUT); j++) {
          if (glob_blocks[b].truthTable[j] > 0) {
            glob_blocks[b].truthTable[j] = 1;
          }
        }
      }

      if (debug_level & debug) {
        std::cerr << "Truth table for " << vars.label() << "\t" << vars.count(Excn::OUT)
                  << " variables\t" << global.count(object_type) << " sets\n";
        for (size_t b = 0; b < global.count(object_type); b++) {
          for (int j = 0; j < vars.count(Excn::OUT); j++) {
            std::cerr << glob_blocks[b].truthTable[j];
          }
          std::cerr << '\n';
        }
      }
    }
  }

  int case_compare(const std::string &s1, const std::string &s2)
  {
    const char *c1 = s1.c_str();
    const char *c2 = s2.c_str();
    for (;;) {
      if (::toupper(*c1) != ::toupper(*c2)) {
        return (::toupper(*c1) - ::toupper(*c2));
      }
      if (*c1 == '\0') {
        return 0;
      }
      c1++;
      c2++;
    }
  }

  void add_info_record(char *info_record, int size)
  {
    // Add 'uname' output to the passed in character string.
    // Maximum size of string is 'size' (not including terminating nullptr)
    // This is used as information data in the concatenated results file
    // to help in tracking when/where/... the file was created
    struct utsname sys_info
    {
    };
    uname(&sys_info);

    std::string info = "CONJOIN: ";
    info += sys_info.nodename;
    info += ", OS: ";
    info += sys_info.sysname;
    info += " ";
    info += sys_info.release;
    info += ", ";
    info += sys_info.version;
    info += ", Machine: ";
    info += sys_info.machine;
    const char *sinfo = info.c_str();
    strncpy(info_record, sinfo, size);
    info_record[size] = '\0';
  }

  inline bool is_whitespace(char c)
  {
    static char white_space[] = {' ', '\t', '\n', '\r', ',', '\0'};
    return (std::strchr(white_space, c) != nullptr);
  }

  void compress_white_space(char *str)
  {
    char *ibuf = str;
    char *obuf = str;

    int i   = 0;
    int cnt = 0;

    // Don't process an empty string.
    if (str == nullptr) {
      return;
    }

    // Skip leading...
    while (*ibuf != 0 && is_whitespace(*ibuf)) {
      ++ibuf;
    }

    while (*ibuf != 0) {
      if (is_whitespace(*ibuf) && cnt > 0) {
        ibuf++;
      }
      else {
        if (!is_whitespace(*ibuf)) {
          cnt = 0;
        }
        else {
          *ibuf = ' ';
          cnt   = 1;
        }
        obuf[i++] = *ibuf++;
      }
    }
    obuf[i--] = '\0';

    // Skip trailing whitespace
    while (i > 0 && is_whitespace(obuf[i])) {
      obuf[i--] = '\0';
    }
  }

  std::string time_stamp(const std::string &format)
  {
    if (format == "") {
      return std::string("");
    }

    const int   length = 256;
    static char time_string[length];

    time_t     calendar_time = time(nullptr);
    struct tm *local_time    = localtime(&calendar_time);

    int error = strftime(time_string, length, format.c_str(), local_time);
    if (error != 0) {
      time_string[length - 1] = '\0';
      return std::string(time_string);
    }

    return std::string("[ERROR]");
  }

  std::string format_time(double seconds)
  {
    char suffix = 'u';
    if (seconds > 0.0 && seconds < 1.0) {
      return " <1s";
    }

    if (seconds > 86400) {
      suffix = 'd';
      seconds /= 86400.;
    }
    else if (seconds > 3600) {
      suffix = 'h';
      seconds /= 3600.;
    }
    else if (seconds > 60) {
      suffix = 'm';
      seconds /= 60.;
    }
    else {
      suffix = 's';
    }
    std::ostringstream os;
    os << std::showpoint << std::setprecision(2) << seconds << suffix;
    return os.str();
  }

  int get_width(size_t max_value)
  {
    // Returns the field width which will accommodate the
    // largest value.
    int width = 0;
    if (max_value >= 10) {
      width = int(log10(static_cast<double>(max_value)));
    }
    return width + 1;
  }

  template <typename T, typename INT>
  void map_element_vars(size_t loffset, size_t goffset, size_t entity_count, std::vector<T> &values,
                        std::vector<T> &global_values, INT *part_loc_elem_to_global)
  {
    // copy values to master element value information
    for (size_t j = 0; j < entity_count; j++) {
      size_t global_block_pos         = part_loc_elem_to_global[(j + loffset)] - goffset;
      global_values[global_block_pos] = values[j];
    }
  }

  template <typename T, typename U>
  void map_sideset_vars(U & /*unused*/, size_t /*unused*/, size_t /*unused*/,
                        std::vector<T> & /*unused*/, std::vector<T> & /*unused*/)
  {
    throw std::runtime_error("Internal Error.");
  }

  template <typename INT>
  void map_sideset_vars(Excn::SideSet<INT> &local_set, size_t entity_count,
                        size_t glob_entity_count, std::vector<double> &values,
                        std::vector<double> &global_values)
  {
    // copy values to master nodeset value information
    for (size_t j = 0; j < entity_count; j++) {
      size_t global_loc = local_set.elemOrderMap[j];
      SMART_ASSERT(global_loc < glob_entity_count);
      global_values[global_loc] = values[j];
    }
  }

  template <typename INT>
  void map_sideset_vars(Excn::SideSet<INT> &local_set, size_t entity_count,
                        size_t glob_entity_count, std::vector<float> &values,
                        std::vector<float> &global_values)
  {
    // copy values to master nodeset value information
    for (size_t j = 0; j < entity_count; j++) {
      size_t global_loc = local_set.elemOrderMap[j];
      SMART_ASSERT(global_loc < glob_entity_count);
      global_values[global_loc] = values[j];
    }
  }

  template <typename T, typename U>
  void map_nodeset_vars(U & /*unused*/, size_t /*unused*/, size_t /*unused*/,
                        std::vector<T> & /*unused*/, std::vector<T> & /*unused*/)
  {
    throw std::runtime_error("Internal Error.");
  }

  template <typename INT>
  void map_nodeset_vars(Excn::NodeSet<INT> &local_set, size_t entity_count,
                        size_t glob_entity_count, std::vector<double> &values,
                        std::vector<double> &global_values)
  {
    // copy values to master nodeset value information
    for (size_t j = 0; j < entity_count; j++) {
      size_t global_loc = local_set.nodeOrderMap[j];
      SMART_ASSERT(global_loc < glob_entity_count);
      global_values[global_loc] = values[j];
    }
  }

  template <typename INT>
  void map_nodeset_vars(Excn::NodeSet<INT> &local_set, size_t entity_count,
                        size_t glob_entity_count, std::vector<float> &values,
                        std::vector<float> &global_values)
  {
    // copy values to master nodeset value information
    for (size_t j = 0; j < entity_count; j++) {
      size_t global_loc = local_set.nodeOrderMap[j];
      SMART_ASSERT(global_loc < glob_entity_count);
      global_values[global_loc] = values[j];
    }
  }

  template <typename T, typename U, typename INT>
  int read_write_master_values(Excn::Variables &vars, const Excn::Mesh<INT> &global,
                               std::vector<U> &              global_sets,
                               std::vector<Excn::Mesh<INT>> &local_mesh,
                               std::vector<std::vector<U>> &local_sets, std::vector<T> &values,
                               size_t p, Excn::ExodusFile &id, int time_step, int time_step_out)
  {
    int error  = 0;
    int id_out = Excn::ExodusFile::output(); // output file identifier

    size_t max_size = 0;
    for (size_t b = 0; b < global.count(vars.objectType); b++) {
      if (max_size < global_sets[b].entity_count()) {
        max_size = global_sets[b].entity_count();
      }
    }
    std::vector<T> master_values(max_size);

    // Only needed for element, but haven't cleaned this up yet...
    INT *part_loc_elem_to_global = local_mesh[p].localElementToGlobal.data();

    for (int i = 0; i < vars.count(Excn::IN); i++) {
      if (vars.index_[i] > 0) {
        int ivar = vars.index_[i] - 1;

        for (size_t b = 0; b < global.count(vars.objectType); b++) {
          size_t lb = 0;
          for (; lb < global.count(vars.objectType); lb++) {
            if (global_sets[b].id == local_sets[p][lb].id) {
              break;
            }
          }
          SMART_ASSERT(global_sets[b].id == local_sets[p][lb].id);

          if (global_sets[b].truthTable[ivar] && local_sets[p][lb].entity_count() > 0) {

            int entity_count = local_sets[p][lb].entity_count();

            if (local_sets[p][lb].truthTable[i] > 0) {
              error += ex_get_var(id, time_step + 1, exodus_object_type(vars.objectType), i + 1,
                                  local_sets[p][lb].id, entity_count, values.data());

              switch (vars.objectType) {
              case Excn::EBLK:
                map_element_vars(local_sets[p][lb].offset_, global_sets[b].offset_, entity_count,
                                 values, master_values, part_loc_elem_to_global);
                break;

              case Excn::SSET:
                map_sideset_vars(local_sets[p][lb], entity_count, global_sets[b].entity_count(),
                                 values, master_values);
                break;

              case Excn::NSET:
                map_nodeset_vars(local_sets[p][lb], entity_count, global_sets[b].entity_count(),
                                 values, master_values);
                break;
              default: break;
              }
            }
            ex_put_var(id_out, time_step_out, exodus_object_type(vars.objectType), ivar + 1,
                       global_sets[b].id, global_sets[b].entity_count(), master_values.data());
          }
        }
      }
    }
    return error;
  }

  template <typename U>
  void create_output_truth_table(std::vector<U> &global_sets, Excn::Variables &vars,
                                 std::vector<int> &truth_table)
  {
    truth_table.resize(global_sets.size() * vars.count(Excn::OUT));
    for (auto &global_set : global_sets) {
      int bout = global_set.position_;
      SMART_ASSERT(bout >= 0);
      for (int j = 0; j < vars.count(Excn::OUT); j++) {
        int out_ttable_loc          = (bout * vars.count(Excn::OUT)) + j;
        truth_table[out_ttable_loc] = global_set.truthTable[j];
      }
    }
  }

  template <typename INT>
  size_t find_max_entity_count(size_t part_count, std::vector<Excn::Mesh<INT>> &local_mesh,
                               const Excn::Mesh<INT> &                       global,
                               std::vector<std::vector<Excn::Block>> &       blocks,
                               std::vector<std::vector<Excn::NodeSet<INT>>> &nodesets,
                               std::vector<std::vector<Excn::SideSet<INT>>> &sidesets)
  {
    size_t max_ent = local_mesh[0].count(Excn::NODE);
    for (size_t p = 1; p < part_count; p++) {
      if ((size_t)local_mesh[p].count(Excn::NODE) > max_ent) {
        max_ent = local_mesh[p].count(Excn::NODE);
      }
    }

    for (size_t p = 0; p < part_count; p++) {
      for (size_t b = 0; b < global.count(Excn::EBLK); b++) {
        if (blocks[p][b].entity_count() > max_ent) {
          max_ent = blocks[p][b].entity_count();
        }
      }
    }

    // Nodesets...
    for (size_t p = 0; p < part_count; p++) {
      for (size_t b = 0; b < global.count(Excn::NSET); b++) {
        if (nodesets[p][b].entity_count() > max_ent) {
          max_ent = nodesets[p][b].entity_count();
        }
      }
    }

    // Sidesets...
    for (size_t p = 0; p < part_count; p++) {
      for (size_t b = 0; b < global.count(Excn::SSET); b++) {
        if (sidesets[p][b].entity_count() > max_ent) {
          max_ent = sidesets[p][b].entity_count();
        }
      }
    }
    return max_ent;
  }
} // namespace
