/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkExodusIIWriter.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkArrayIteratorIncludes.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCompositeDataIterator.h"
#include "vtkCompositeDataSet.h"
#include "vtkDataObject.h"
#include "vtkDataObjectTreeIterator.h"
#include "vtkDoubleArray.h"
#include "vtkExodusIIWriter.h"
#include "vtkFieldData.h"
#include "vtkFloatArray.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkModelMetadata.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkStdString.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStringArray.h"
#include "vtkThreshold.h"
#include "vtkToolkits.h" // for VTK_USE_PARALLEL
#include "vtkUnstructuredGrid.h"

#include "vtk_exodusII.h"
#include <cctype>
#include <map>
#include <ctime>
#include <sstream>

vtkObjectFactoryNewMacro (vtkExodusIIWriter);
vtkCxxSetObjectMacro (vtkExodusIIWriter, ModelMetadata, vtkModelMetadata);

namespace
{
unsigned int GetNumberOfDigits(unsigned int i)
{
  if (i<10)
  {
    return 1;
  }
  return GetNumberOfDigits(i/10)+1;
}
}

//----------------------------------------------------------------------------

vtkExodusIIWriter::vtkExodusIIWriter ()
{
  this->fid = -1;
  this->FileName = nullptr;

  this->StoreDoubles = -1;
  this->GhostLevel = 0;
  this->WriteOutBlockIdArray = 0;
  this->WriteOutGlobalElementIdArray = 0;
  this->WriteOutGlobalNodeIdArray = 0;
  this->WriteAllTimeSteps = 0;
  this->BlockIdArrayName = nullptr;
  this->ModelMetadata = nullptr;

  this->NumberOfTimeSteps = 0;
  this->CurrentTimeIndex = 0;
  this->FileTimeOffset = 0;

  this->AtLeastOneGlobalElementIdList = 0;
  this->AtLeastOneGlobalNodeIdList = 0;

  this->NumPoints = 0;
  this->NumCells = 0;

  this->PassDoubles = 1;

  this->BlockElementVariableTruthTable = nullptr;

  this->LocalNodeIdMap = nullptr;
  this->LocalElementIdMap = nullptr;
  this->TopologyChanged = false;
  this->IgnoreMetaDataWarning = 0;
}

vtkExodusIIWriter::~vtkExodusIIWriter ()
{
  this->SetModelMetadata(nullptr); // kill the reference if its there

  delete [] this->FileName;
  delete [] this->BlockIdArrayName;


  delete [] this->BlockElementVariableTruthTable;

  for (size_t i = 0; i < this->BlockIdList.size (); i ++)
  {
    this->BlockIdList[i]->UnRegister (this);
  }
}

void vtkExodusIIWriter::PrintSelf (ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  os << indent << "FileName " << (this->FileName ? this->FileName : "(none)") << endl;
  os << indent << "StoreDoubles " << this->StoreDoubles << endl;
  os << indent << "GhostLevel " << this->GhostLevel << endl;
  os << indent << "WriteOutBlockIdArray " << this->WriteOutBlockIdArray << endl;
  os << indent << "WriteOutGlobalNodeIdArray " << this->WriteOutGlobalNodeIdArray << endl;
  os << indent << "WriteOutGlobalElementIdArray " << this->WriteOutGlobalElementIdArray << endl;
  os << indent << "WriteAllTimeSteps " << this->WriteAllTimeSteps << endl;
  os << indent << "BlockIdArrayName " << (this->BlockIdArrayName ? this->BlockIdArrayName : "(none)") << endl;
  os << indent << "ModelMetadata " << (this->ModelMetadata ? "" : "(none)") << endl;
  if (this->ModelMetadata)
  {
    this->ModelMetadata->PrintSelf (os, indent.GetNextIndent ());
  }
  os << indent << "IgnoreMetaDataWarning " << this->IgnoreMetaDataWarning << endl;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::ProcessRequest (
  vtkInformation* request,
  vtkInformationVector** inputVector,
  vtkInformationVector* outputVector)
{
  if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
  {
    return this->RequestInformation(request, inputVector, outputVector);
  }
  else if(request->Has(
      vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
  {
    return this->RequestUpdateExtent(request, inputVector, outputVector);
  }
  // generate the data
  else if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
  {
    return this->RequestData(request, inputVector, outputVector);
  }

  return this->Superclass::ProcessRequest(request, inputVector, outputVector);
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::RequestInformation (
  vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector,
  vtkInformationVector* vtkNotUsed(outputVector))
{
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  if ( inInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()) )
  {
    this->NumberOfTimeSteps =
      inInfo->Length( vtkStreamingDemandDrivenPipeline::TIME_STEPS() );
  }
  else
  {
    this->NumberOfTimeSteps = 0;
  }

  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::RequestUpdateExtent (
  vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector,
  vtkInformationVector* vtkNotUsed(outputVector))
{
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  if ( this->WriteAllTimeSteps &&
       inInfo->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS()) )
  {
    double* timeSteps =
      inInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS() );
    double timeReq= timeSteps[this->CurrentTimeIndex];
    inputVector[0]->GetInformationObject(0)->Set
      ( vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP(),  timeReq);
  }
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::FillInputPortInformation (
  int vtkNotUsed (port),
  vtkInformation* info)
{
  info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(),
    "vtkCompositeDataSet");
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::RequestData (
  vtkInformation* request,
  vtkInformationVector** inputVector,
  vtkInformationVector* vtkNotUsed(outputVector))
{
  if (!this->FileName)
  {
    return 1;
  }

  vtkInformation* inInfo = inputVector[0]->GetInformationObject (0);
  this->OriginalInput = vtkDataObject::SafeDownCast (
    inInfo->Get(vtkDataObject::DATA_OBJECT ()));

  // is this the first request
  if (this->CurrentTimeIndex == 0 && this->WriteAllTimeSteps)
  {
    // Tell the pipeline to start looping.
    request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 1);
  }

  this->WriteData ();

  this->CurrentTimeIndex ++;
  if (this->CurrentTimeIndex >= this->NumberOfTimeSteps ||
      this->TopologyChanged)
  {
    this->CloseExodusFile ();
    this->CurrentTimeIndex = 0;
    if (this->WriteAllTimeSteps)
    {
      // Tell the pipeline to stop looping.
      request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 0);
    }
  }
  // still close out the file after each step written.
  if (!this->WriteAllTimeSteps)
  {
    this->CloseExodusFile ();
  }

  int localContinue = request->Get(
    vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING());
  if (this->GlobalContinueExecuting(localContinue) != localContinue)
  {
    // Some other node decided to stop the execution.
    assert (localContinue == 1);
    request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 0);
  }
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::GlobalContinueExecuting(int localContinueExecution)
{
  return localContinueExecution;
}

//----------------------------------------------------------------------------
void vtkExodusIIWriter::WriteData ()
{
  this->NewFlattenedInput.clear ();
  this->NewFlattenedNames.clear ();
  // Is it safe to assume this is the same?
  bool newHierarchy = false;
  if (!this->FlattenHierarchy (this->OriginalInput, "", newHierarchy))
  {
    vtkErrorMacro (
      "vtkExodusIIWriter::WriteData Unable to flatten hierarchy");
    return;
  }
  if (this->FlattenedInput.size () != this->NewFlattenedInput.size ())
  {
    newHierarchy = true;
  }

  // For now, we don't support changing topology even if the writer supports it.
  // The reader needs to be updated to support changing topology.
  if (newHierarchy && !this->FlattenedInput.empty())
  {
    this->TopologyChanged = true;
    vtkErrorMacro (
      "Temporal data with changing topology is not yet supported.");
    return;
  }

  // Copies over the new results data in the new objects
  this->FlattenedInput = this->NewFlattenedInput;
  this->FlattenedNames = this->NewFlattenedNames;

  this->RemoveGhostCells ();

  // move check parameters up here and then if there's a change, new file.
  if (this->WriteAllTimeSteps && this->CurrentTimeIndex != 0 && !newHierarchy)
  {
    if (!this->WriteNextTimeStep ())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData results");
    }
    return;
  }
  else
  {
    // Close out the old file, if we have one
    if (this->CurrentTimeIndex > 0)
    {
      this->CloseExodusFile ();
    }


    // The file has changed, initialize new file
    if (!this->CheckParameters ())
    {
      return;
    }

    // TODO this should increment a counter
    if (!this->CreateNewExodusFile ())
    {
      vtkErrorMacro (
        "vtkExodusIIWriter: WriteData can't create exodus file");
      return;
    }

    if (!this->WriteInitializationParameters ())
    {
      vtkErrorMacro(<< "vtkExodusIIWriter::WriteData init params");
      return;
    }

    if (!this->WriteInformationRecords())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData information records");
      return;
    }

    if (!this->WritePoints())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData points");
      return;
    }

    if (!this->WriteCoordinateNames())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData coordinate names");
      return;
    }

    if (!this->WriteGlobalPointIds())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData global point IDs");
      return;
    }

    if (!this->WriteBlockInformation())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData block information");
      return;
    }

    if (!this->WriteGlobalElementIds())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData global element IDs");
      return;
    }

    if (!this->WriteVariableArrayNames())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData variable array names");
      return;
    }

    if (!this->WriteNodeSetInformation())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData can't node sets");
      return;
    }

    if (!this->WriteSideSetInformation())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData can't side sets");
      return;
    }

    if (!this->WriteProperties())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData can't properties");
      return;
    }

    if (!this->WriteNextTimeStep ())
    {
      vtkErrorMacro("vtkExodusIIWriter::WriteData results");
      return;
    }
  }
}

//----------------------------------------------------------------------------
char *vtkExodusIIWriter::StrDupWithNew(const char *s)
{
  char *newstr = nullptr;

  if (s)
  {
    size_t len = strlen(s);
    newstr = new char [len + 1];
    strcpy(newstr, s);
  }

  return newstr;
}

//----------------------------------------------------------------------------
void vtkExodusIIWriter::StringUppercase(std::string& str)
{
  for (size_t i=0; i<str.size (); i++)
  {
    str[i] = toupper(str[i]);
  }
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::FlattenHierarchy (vtkDataObject* input, const char *name, bool& changed)
{
  if (input->IsA ("vtkMultiBlockDataSet"))
  {
    vtkMultiBlockDataSet* castObj = vtkMultiBlockDataSet::SafeDownCast(input);
    vtkSmartPointer<vtkDataObjectTreeIterator> iter;
    iter.TakeReference (castObj->NewTreeIterator ());
    iter->VisitOnlyLeavesOff ();
    iter->TraverseSubTreeOff ();
    iter->SkipEmptyNodesOff ();
    for (iter->InitTraversal ();
         !iter->IsDoneWithTraversal ();
         iter->GoToNextItem ())
    {
      name = iter->GetCurrentMetaData()->Get (vtkCompositeDataSet::NAME());
      if (name != nullptr && strstr (name, "Sets") != nullptr)
      {
        continue;
      }
      if (name == nullptr)
      {
        // avoid null references in StdString
        name = "";
      }
      if (iter->GetCurrentDataObject () && !this->FlattenHierarchy (iter->GetCurrentDataObject (), name, changed))
      {
        return 0;
      }
    }
  }
  else if (input->IsA ("vtkCompositeDataSet"))
  {
    vtkCompositeDataSet* castObj = vtkCompositeDataSet::SafeDownCast(input);
    vtkSmartPointer<vtkCompositeDataIterator> iter;
    iter.TakeReference (castObj->NewIterator ());
    vtkDataObjectTreeIterator* treeIter = vtkDataObjectTreeIterator::SafeDownCast (castObj);
    if (treeIter)
    {
      treeIter->VisitOnlyLeavesOff ();
      treeIter->TraverseSubTreeOff ();
      treeIter->SkipEmptyNodesOff ();
    }
    for (iter->InitTraversal ();
         !iter->IsDoneWithTraversal ();
         iter->GoToNextItem ())
    {
      if (iter->GetCurrentDataObject () && !this->FlattenHierarchy (iter->GetCurrentDataObject (), name, changed))
      {
        return 0;
      }
    }
  }
  else if (input->IsA ("vtkDataSet"))
  {
    vtkSmartPointer<vtkUnstructuredGrid> output =
        vtkSmartPointer<vtkUnstructuredGrid>::New ();
    if (input->IsA ("vtkUnstructuredGrid"))
    {
      output->ShallowCopy (input);
    }
    else
    {
      vtkDataSet* castObj = vtkDataSet::SafeDownCast (input);

      output->GetFieldData ()->ShallowCopy (castObj->GetFieldData ());
      output->GetPointData ()->ShallowCopy (castObj->GetPointData ());
      output->GetCellData ()->ShallowCopy (castObj->GetCellData ());

      vtkIdType numPoints = castObj->GetNumberOfPoints ();
      vtkSmartPointer<vtkPoints> outPoints = vtkSmartPointer<vtkPoints>::New ();
      outPoints->SetNumberOfPoints (numPoints);
      for (vtkIdType i = 0; i < numPoints; i ++)
      {
        outPoints->SetPoint (i, castObj->GetPoint (i));
      }
      output->SetPoints (outPoints);

      int numCells = castObj->GetNumberOfCells ();
      output->Allocate (numCells);
      vtkIdList* ptIds = vtkIdList::New ();
      for (int i = 0; i < numCells; i ++)
      {
        castObj->GetCellPoints (i, ptIds);
        output->InsertNextCell (castObj->GetCellType (i), ptIds);
      }
      ptIds->Delete ();
    }
    // check to see if we need a new exodus file
    // because the element or node count needs to be changed
    size_t checkIndex = this->NewFlattenedInput.size ();
    if (this->FlattenedInput.size () > checkIndex)
    {
      int numPoints = this->FlattenedInput[checkIndex]->GetNumberOfPoints ();
      int numCells = this->FlattenedInput[checkIndex]->GetNumberOfCells ();
      if (numPoints != output->GetNumberOfPoints () ||
          numCells != output->GetNumberOfCells ())
      {
        changed = true;
      }
    }
    else
    {
      changed = true;
    }
    this->NewFlattenedInput.push_back (output);
    if (!name)
    {
      // Setting an arbitrary name for datasets that have not been assigned one.
      name = "block";
    }
    this->NewFlattenedNames.push_back (name);
  }
  else
  {
    vtkErrorMacro (<< "Incorrect class type " << input->GetClassName () << " on input");
    return 0;
  }
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::CreateNewExodusFile()
{
  int compWordSize= (this->PassDoubles ? sizeof(double) : sizeof(float));
  int IOWordSize = (this->StoreDoubles ? sizeof(double) : sizeof(float));


  if (this->NumberOfProcesses == 1)
  {
    if (this->WriteAllTimeSteps == false || this->CurrentTimeIndex == 0)
    {
      this->fid = ex_create(this->FileName, EX_CLOBBER, &compWordSize, &IOWordSize);
      if (fid <= 0)
      {
        vtkErrorMacro (
          << "vtkExodusIIWriter: CreateNewExodusFile can't create "
          << this->FileName);
      }
    }
    else
    {
      char *myFileName = new char [VTK_MAXPATH];
      snprintf(myFileName, VTK_MAXPATH, "%s-s.%06d", this->FileName, this->CurrentTimeIndex);
      this->fid = ex_create(myFileName, EX_CLOBBER, &compWordSize, &IOWordSize);
      if (fid <= 0)
      {
        vtkErrorMacro (
          << "vtkExodusIIWriter: CreateNewExodusFile can't create "
          << myFileName);
      }
      delete [] myFileName;
    }
  }
  else
  {
    std::ostringstream myFileName;
    myFileName << this->FileName;
    if (this->WriteAllTimeSteps == false || this->CurrentTimeIndex == 0)
    {
      myFileName << ".";
    }
    else
    {
      myFileName << "-s." << std::setfill('0') << std::setw(6) <<
        this->CurrentTimeIndex << std::setw(0) << ".";
    }
    unsigned int numDigits = GetNumberOfDigits(
      static_cast<unsigned int>(this->NumberOfProcesses-1));
    myFileName << this->NumberOfProcesses << "." << std::setfill('0')
               << std::setw(numDigits) << this->MyRank;
    this->fid = ex_create(myFileName.str().c_str(), EX_CLOBBER,
                          &compWordSize, &IOWordSize);
    if (this->fid <= 0)
    {
      vtkErrorMacro (
        << "vtkExodusIIWriter: CreateNewExodusFile can't create "
        << myFileName.str());
    }
  }
  ex_set_max_name_length(this->fid, static_cast<int>(this->GetMaxNameLength()));

  // FileTimeOffset makes the time in the file relative
  // e.g., if the CurrentTimeIndex for this file is 4 and goes through 6, the
  // file will write them as 0 1 2 instead of 4 5 6
  this->FileTimeOffset = this->CurrentTimeIndex;
  return this->fid > 0;
}
//
//------------------------------------------------------------------
void vtkExodusIIWriter::CloseExodusFile()
{
  if (this->fid >= 0)
  {
    ex_close(this->fid);
    this->fid = -1;
    return;
  }
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::IsDouble ()
{
  // Determine whether we should pass single or double precision
  // floats to the Exodus Library.  We'll look through the arrays
  // and points in the input and pick the precision of the
  // first float we see.

  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {

    vtkCellData *cd = this->FlattenedInput[i]->GetCellData();
    if (cd)
    {
      int numCellArrays = cd->GetNumberOfArrays();
      for (int j=0; j<numCellArrays; j++)
      {
        vtkDataArray *a = cd->GetArray(j);
        int type = a->GetDataType();

        if (type == VTK_DOUBLE)
        {
          return 1;
        }
        else if (type == VTK_FLOAT)
        {
          return 0;
        }
      }
    }

    vtkPointData *pd = this->FlattenedInput[i]->GetPointData();
    if (pd)
    {
      int numPtArrays = pd->GetNumberOfArrays();
      for (int j=0; j<numPtArrays; j++)
      {
        vtkDataArray *a = pd->GetArray(j);
        int type = a->GetDataType();

        if (type == VTK_DOUBLE)
        {
          return 1;
        }
        else if (type == VTK_FLOAT)
        {
          return 0;
        }
      }
    }

  }
  return -1;
}

//----------------------------------------------------------------------------
void vtkExodusIIWriter::RemoveGhostCells()
{
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    vtkUnsignedCharArray *da = this->FlattenedInput[i]->GetCellGhostArray();

    if (da)
    {
      vtkThreshold *t = vtkThreshold::New();
      t->SetInputData(this->FlattenedInput[i]);
      t->ThresholdByLower(0);
      t->SetInputArrayToProcess(
        0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, vtkDataSetAttributes::GhostArrayName());

      t->Update();

      this->FlattenedInput[i] = vtkSmartPointer<vtkUnstructuredGrid>(t->GetOutput());
      t->Delete();

      this->FlattenedInput[i]->GetCellData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());
      this->FlattenedInput[i]->GetPointData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());

      this->GhostLevel = 1;
    }
    else
    {
      this->GhostLevel = 0;
    }

  }
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::CheckParametersInternal (int _NumberOfProcesses, int _MyRank)
{
   if (!this->FileName)
   {
     vtkErrorMacro("No filename specified.");
     return 0;
   }

   this->PassDoubles = this->IsDouble ();
   if (this->PassDoubles < 0)
   {
     // Can't find float types in input, assume doubles
     this->PassDoubles = 1;
   }

   if (this->StoreDoubles < 0)
   {
     // The default is to store in the
     // same precision that appears in the input.

     this->StoreDoubles = this->PassDoubles;
   }

   this->NumberOfProcesses = _NumberOfProcesses;
   this->MyRank = _MyRank;

  if (!this->CheckInputArrays ())
  {
    return 0;
  }

  if (!this->ConstructBlockInfoMap ())
  {
    return 0;
  }

  if (!this->ConstructVariableInfoMaps ())
  {
    return 0;
  }

  // Data (and possibly the topology) is different for every time step so
  // always create a new metadata.
  if (!this->CreateDefaultMetadata ())
  {
    return 0;
  }

  if (!this->ParseMetadata ())
  {
    return 0;
  }

  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::CheckParameters ()
{
  return this->CheckParametersInternal(1, 0);
}

int vtkExodusIIWriter::CheckInputArrays ()
{
  this->BlockIdList.resize (this->FlattenedInput.size ());
  this->GlobalElementIdList.resize (this->FlattenedInput.size ());
  this->AtLeastOneGlobalElementIdList = 0;
  this->GlobalNodeIdList.resize (this->FlattenedInput.size ());
  this->AtLeastOneGlobalNodeIdList = 0;

  this->NumPoints = 0;
  this->NumCells = 0;
  this->MaxId = 0;

  size_t i;
  for (i = 0; i < this->FlattenedInput.size (); i ++)
  {
    this->NumPoints += this->FlattenedInput[i]->GetNumberOfPoints ();
    int numCells = this->FlattenedInput[i]->GetNumberOfCells ();
    this->NumCells += numCells;

    vtkCellData *cd = this->FlattenedInput[i]->GetCellData();
    vtkPointData *pd = this->FlattenedInput[i]->GetPointData();

    vtkIntArray *bia = this->GetBlockIdArray (
      this->BlockIdArrayName, this->FlattenedInput[i]);
    if (bia)
    {
      this->BlockIdList[i] = bia;
      this->BlockIdList[i]->Register (this);
      // computing the max known id in order to create unique fill in values below
      for (int j = 0; j < numCells; j ++)
      {
        if (this->BlockIdList[i]->GetValue (j) > this->MaxId)
        {
          this->MaxId = this->BlockIdList[i]->GetValue (j);
        }
      }
    }
    else
    {
      // Will fill in below
      this->BlockIdList[i] = nullptr;
    }

    // Trying to find global element id
    vtkDataArray *da = cd->GetGlobalIds();
    if (!da)
    { // try finding the array explicitly named
      da = cd->GetArray ("GlobalElementId");
    }
    if (da)
    {
      vtkIdTypeArray *ia = vtkArrayDownCast<vtkIdTypeArray>(da);
      if (!ia)
      {
        vtkWarningMacro(<<
          "vtkExodusIIWriter, element ID array is not an Id array, ignoring it");
        this->GlobalElementIdList[i] = nullptr;
      }
      else
      {
        this->GlobalElementIdList[i] = ia->GetPointer(0);
        this->AtLeastOneGlobalElementIdList = 1;
      }
    }


    // Trying to find global node id
    da = pd->GetGlobalIds();
    if (!da)
    { // try finding the array explicitly named
      da = pd->GetArray ("GlobalNodeId");
    }
    if (da)
    {
      vtkIdTypeArray *ia = vtkArrayDownCast<vtkIdTypeArray>(da);
      if (!ia)
      {
        vtkWarningMacro(<<
          "vtkExodusIIWriter, node ID array is not an Id array, ignoring it");
        this->GlobalNodeIdList[i] = nullptr;
      }
      else
      {
        this->GlobalNodeIdList[i] = ia->GetPointer(0);
        this->AtLeastOneGlobalNodeIdList = 1;
      }
    }
    else
    {
      this->GlobalNodeIdList[i] = nullptr;
    }
  }

  return 1;
}

int vtkExodusIIWriter::ConstructBlockInfoMap ()
{
  // The elements in the input may not be in order by block, but we must
  // write element IDs and element variables out to the Exodus file in
  // order by block.  Create a mapping if necessary, for an ordering by
  // block to the ordering found in the input unstructured grid.
  this->CellToElementOffset.resize (this->FlattenedInput.size ());
  this->BlockInfoMap.clear ();
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    // If we weren't explicitly given the block ids, try to extract them from the
    // block id array embedded in the cell data.
    vtkIdType ncells = this->FlattenedInput[i]->GetNumberOfCells();
    if (!this->BlockIdList[i])
    {
      vtkIntArray *ia = vtkIntArray::New ();
      ia->SetNumberOfValues (ncells);
      for (int j = 0; j < ncells; j ++)
      {
        ia->SetValue (j, this->FlattenedInput[i]->GetCellType(j) + this->MaxId);
      }

      // Pretend we had it in the metadata
      this->BlockIdList[i] = ia;
      this->BlockIdList[i]->Register (this);
      ia->Delete ();

      // Also increment the MaxId so we can keep it unique
      this->MaxId += VTK_NUMBER_OF_CELL_TYPES;
    }

    // Compute all the block information mappings.
    this->CellToElementOffset[i].resize (ncells);
    for(int j=0; j<ncells; j++)
    {
      std::map<int, Block>::iterator iter =
        this->BlockInfoMap.find (this->BlockIdList[i]->GetValue (j));
                                 // shift by 1 in case there's a 0  This was removed because it changes the user supplied block ids in the metadata.
      if (iter == this->BlockInfoMap.end ())
      {
        this->CellToElementOffset[i][j] = 0;
        Block& b = this->BlockInfoMap[this->BlockIdList[i]->GetValue (j)];//CLM
        b.Name = this->FlattenedNames[i].c_str();
        b.Type = this->FlattenedInput[i]->GetCellType (j);
        b.NumElements = 1;
        b.ElementStartIndex = 0;
        switch (b.Type)
        {
          case VTK_POLY_LINE:
          case VTK_POLYGON:
          case VTK_POLYHEDRON:
            // this block contains variable numbers of nodes per element
            b.NodesPerElement = 0;
            b.EntityCounts = std::vector<int>(ncells);
            b.EntityCounts[0] = this->FlattenedInput[i]->GetCell (j)->GetNumberOfPoints ();
            b.EntityNodeOffsets = std::vector<int>(ncells);
            b.EntityNodeOffsets[0] = 0;
            break;
          default:
            b.NodesPerElement = this->FlattenedInput[i]->GetCell (j)->GetNumberOfPoints ();
        };

        // TODO this could be a push if i is different.
        b.GridIndex = i;

        // This may get pulled from the meta data below,
        // but if not, default reasonably to 0
        b.NumAttributes = 0;
        b.BlockAttributes = nullptr;
      }
      else
      {
        // TODO we should be able to deal with this, not just error out
        if (iter->second.GridIndex != i)
        {
          vtkWarningMacro ("Block ids are not unique across the hierarchy ");
        }

        this->CellToElementOffset[i][j] = iter->second.NumElements;
        int index = iter->second.NumElements;
        if (iter->second.NodesPerElement == 0)
        {
          iter->second.EntityCounts[index] =
                  this->FlattenedInput[i]->GetCell (j)->GetNumberOfPoints ();
          iter->second.EntityNodeOffsets[index] =
                  iter->second.EntityNodeOffsets[index - 1] +
                  iter->second.EntityCounts[index - 1];
        }
        iter->second.NumElements ++;
      }
    }
  }

  this->CheckBlockInfoMap();

  // Find the ElementStartIndex and the output order
  std::map<int, Block>::iterator iter;
  int runningCount = 0;
  int index = 0;
  for (iter = this->BlockInfoMap.begin ();
       iter != this->BlockInfoMap.end ();
       ++iter)
  {
    iter->second.ElementStartIndex = runningCount;
    runningCount += iter->second.NumElements;

    iter->second.OutputIndex = index;
    index ++;
  }

  return 1;
}

int vtkExodusIIWriter::ConstructVariableInfoMaps ()
{
  // Create the variable info maps.
  if (this->GlobalVariableMap.empty())
  {
    this->NumberOfScalarGlobalArrays = 0;
  }
  this->NumberOfScalarElementArrays = 0;
  this->NumberOfScalarNodeArrays = 0;
  this->BlockVariableMap.clear ();
  this->NodeVariableMap.clear ();
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    vtkFieldData *fd = this->FlattenedInput[i]->GetFieldData ();
    for (int j = 0; j < fd->GetNumberOfArrays (); j ++)
    {
      char *name = nullptr;
      if (fd->GetAbstractArray(j))
      {
        name = fd->GetAbstractArray(j)->GetName();
      }
      if (name == nullptr)
      {
        vtkWarningMacro ("Array in input field data has Null name, cannot output it");
        continue;
      }
      std::string upper (name);
      this->StringUppercase(upper);
      if (strncmp(upper.c_str (),"TITLE",5) == 0)
      {
        continue;
      }
      if (strncmp(upper.c_str (),"QA_RECORD",9) == 0)
      {
        continue;
      }
      if (strncmp(upper.c_str (),"INFO_RECORD",11) == 0)
      {
        continue;
      }
      if (strncmp(upper.c_str (),"ELEMENTBLOCKIDS",15) == 0)
      {
        continue;
      }
      int numComp = fd->GetAbstractArray(j)->GetNumberOfComponents ();
      std::map<std::string, VariableInfo>::const_iterator iter =
        this->GlobalVariableMap.find (name);
      if (iter == this->GlobalVariableMap.end ())
      {
        VariableInfo &info = this->GlobalVariableMap[name];
        info.NumComponents = numComp;
        info.OutNames.resize (info.NumComponents);
        info.ScalarOutOffset = this->NumberOfScalarGlobalArrays;
        this->NumberOfScalarGlobalArrays += numComp;
      }
      else if (iter->second.NumComponents != numComp)
      {
        vtkErrorMacro (
          "Disagreement in the hierarchy for the number of components in " <<
          name);
        return 0;
      }
    }

    vtkCellData *cd = this->FlattenedInput[i]->GetCellData();
    for (int j = 0; j < cd->GetNumberOfArrays(); j ++)
    {
      char *name = nullptr;
      if (cd->GetArray(j))
      {
        name = cd->GetArray(j)->GetName();
      }
      if (name == nullptr)
      {
        vtkWarningMacro ("Array in input cell data has Null name, cannot output it");
        continue;
      }
      std::string upper (name);
      this->StringUppercase(upper);

      if (!this->WriteOutGlobalElementIdArray &&
          cd->IsArrayAnAttribute(j) == vtkDataSetAttributes::GLOBALIDS)
      {
        continue;
      }
      if (!this->WriteOutBlockIdArray &&
          this->BlockIdArrayName && strcmp (name, this->BlockIdArrayName) == 0)
      {
        continue;
      }
      if (strncmp(upper.c_str (),"PEDIGREE",8) == 0)
      {
        continue;
      }

      int numComp = cd->GetArray(j)->GetNumberOfComponents ();
      std::map<std::string, VariableInfo>::const_iterator iter =
        this->BlockVariableMap.find (name);
      if (iter == this->BlockVariableMap.end ())
      {
        VariableInfo &info = this->BlockVariableMap[name];
        info.NumComponents = numComp;
        info.OutNames.resize (info.NumComponents);
        info.ScalarOutOffset = this->NumberOfScalarElementArrays;
        this->NumberOfScalarElementArrays += numComp;
      }
      else if (iter->second.NumComponents != numComp)
      {
        vtkErrorMacro (
          "Disagreement in the hierarchy for the number of components in " <<
          name);
        return 0;
      }
    }

    vtkPointData *pd = this->FlattenedInput[i]->GetPointData();
    for (int j = 0; j < pd->GetNumberOfArrays(); j ++)
    {
      char *name = nullptr;
      if (pd->GetArray(j))
      {
        name = pd->GetArray(j)->GetName();
      }
      if (name == nullptr)
      {
        vtkWarningMacro ("Array in input point data has Null name, cannot output it");
        continue;
      }
      std::string upper (name);
      this->StringUppercase(upper);

      // Is this array displacement?
      // If it is and we are not writing all the timesteps,
      // do not write out. It would mess up the geometry the
      // next time the file was read in.
      if (!this->WriteOutGlobalNodeIdArray &&
          pd->IsArrayAnAttribute(j) == vtkDataSetAttributes::GLOBALIDS)
      {
        continue;
      }
      if (strncmp(upper.c_str (),"PEDIGREE",8) == 0)
      {
        continue;
      }
      // Is this array displacement?
      // If it is and we are not writing all the timesteps,
      // do not write out. It would mess up the geometry the
      // next time the file was read in.
      if (!this->WriteAllTimeSteps && strncmp(upper.c_str (),"DIS",3) == 0)
      {
        continue;
      }

      int numComp = pd->GetArray(j)->GetNumberOfComponents ();
      std::map<std::string, VariableInfo>::const_iterator iter =
        this->NodeVariableMap.find (name);
      if (iter == this->NodeVariableMap.end ())
      {
        VariableInfo &info = this->NodeVariableMap[name];
        info.NumComponents = numComp;
        info.OutNames.resize (info.NumComponents);
        info.ScalarOutOffset = this->NumberOfScalarNodeArrays;
        this->NumberOfScalarNodeArrays += info.NumComponents;
      }
      else if (iter->second.NumComponents != numComp)
      {
        vtkErrorMacro (
          "Disagreement in the hierarchy for the number of components in " <<
          name);
        return 0;
      }
    }
  }

  // BLOCK/ELEMENT TRUTH TABLE
  size_t ttsize = this->BlockInfoMap.size () * this->NumberOfScalarElementArrays;
  if (ttsize > 0)
  {
    this->BlockElementVariableTruthTable = new int [ttsize];
    int index = 0;
    std::map<int, Block>::const_iterator blockIter;
    for (blockIter = this->BlockInfoMap.begin ();
         blockIter != this->BlockInfoMap.end ();
         ++blockIter)
    {
      std::map<std::string, VariableInfo>::const_iterator varIter;
      for (varIter = this->BlockVariableMap.begin ();
           varIter != this->BlockVariableMap.end ();
           ++varIter)
      {
        vtkCellData *cd =
          this->FlattenedInput[blockIter->second.GridIndex]->GetCellData();
        int truth = 0;
        if (cd)
        {
          truth = 1;
        }
        for (int c = 0; c < varIter->second.NumComponents; c ++)
        {
            this->BlockElementVariableTruthTable[index++] = truth;
        }
      }
    }
  }

  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::CreateDefaultMetadata ()
{
  // There is no metadata associated with this input.  If we have enough
  // information, we create reasonable defaults.

  vtkModelMetadata *em = vtkModelMetadata::New();

  char *title = new char [MAX_LINE_LENGTH + 1];
  time_t currentTime = time(nullptr);
  char *stime = ctime(&currentTime);

  snprintf(title, MAX_LINE_LENGTH + 1, "Created by vtkExodusIIWriter, %s", stime);

  em->SetTitle(title);

  delete [] title;

  char **dimNames = new char * [3];
  dimNames[0] = vtkExodusIIWriter::StrDupWithNew("X");
  dimNames[1] = vtkExodusIIWriter::StrDupWithNew("Y");
  dimNames[2] = vtkExodusIIWriter::StrDupWithNew("Z");
  em->SetCoordinateNames(3, dimNames);

  if (!this->CreateBlockIdMetadata(em))
  {
    return 0;
  }

  if (!this->CreateBlockVariableMetadata(em))
  {
    return 0;
  }

  this->CreateSetsMetadata(em);
  this->SetModelMetadata(em);
  em->Delete();

  return 1;
}

//----------------------------------------------------------------------------
char *vtkExodusIIWriter::GetCellTypeName(int t)
{
  if (MAX_STR_LENGTH < 32) return nullptr;
  char *nm = new char [MAX_STR_LENGTH + 1];

  switch (t)
  {
    case VTK_EMPTY_CELL:
      strcpy(nm, "empty cell");
      break;
    case VTK_VERTEX:
      strcpy(nm, "sphere");
      break;
    case VTK_POLY_VERTEX:
      strcpy(nm, "sup");
      break;
    case VTK_LINE:
      strcpy(nm, "edge");
      break;
    case VTK_POLY_LINE:
      strcpy(nm, "NSIDED");
      break;
    case VTK_TRIANGLE:
      strcpy(nm, "TRIANGLE");
      break;
    case VTK_TRIANGLE_STRIP:
      strcpy(nm, "TRIANGLE");
      break;
    case VTK_POLYGON:
      strcpy(nm, "NSIDED");
      break;
    case VTK_POLYHEDRON:
      strcpy(nm, "NFACED");
      break;
    case VTK_PIXEL:
      strcpy(nm, "sphere");
      break;
    case VTK_QUAD:
      strcpy(nm, "quad");
      break;
    case VTK_TETRA:
      strcpy(nm, "TETRA");
      break;
    case VTK_VOXEL :
      strcpy(nm, "HEX");
      break;
    case VTK_HEXAHEDRON:
      strcpy(nm, "HEX");
      break;
    case VTK_WEDGE:
      strcpy(nm, "wedge");
      break;
    case VTK_PYRAMID:
      strcpy(nm, "pyramid");
      break;
    case VTK_PENTAGONAL_PRISM:
      strcpy(nm, "pentagonal prism");
      break;
    case VTK_HEXAGONAL_PRISM:
      strcpy(nm, "hexagonal prism");
      break;
    case VTK_QUADRATIC_EDGE:
      strcpy(nm, "edge");
      break;
    case VTK_QUADRATIC_TRIANGLE:
      strcpy(nm, "triangle");
      break;
    case VTK_QUADRATIC_QUAD:
      strcpy(nm, "quad");
      break;
    case VTK_QUADRATIC_TETRA:
      strcpy(nm, "tetra");
      break;
    case VTK_QUADRATIC_HEXAHEDRON:
      strcpy(nm, "hexahedron");
      break;
    case VTK_QUADRATIC_WEDGE:
      strcpy(nm, "wedge");
      break;
    case VTK_QUADRATIC_PYRAMID:
      strcpy(nm, "pyramid");
      break;
    case VTK_CONVEX_POINT_SET:
      strcpy(nm, "convex point set");
      break;
    case VTK_PARAMETRIC_CURVE:
      strcpy(nm, "parametric curve");
      break;
    case VTK_PARAMETRIC_SURFACE:
      strcpy(nm, "parametric surface");
      break;
    case VTK_PARAMETRIC_TRI_SURFACE:
      strcpy(nm, "parametric tri surface");
      break;
    case VTK_PARAMETRIC_QUAD_SURFACE:
      strcpy(nm, "parametric quad surface");
      break;
    case VTK_PARAMETRIC_TETRA_REGION:
      strcpy(nm, "parametric tetra region");
      break;
    case VTK_PARAMETRIC_HEX_REGION:
      strcpy(nm, "paramertric hex region");
      break;
    default:
      strcpy(nm, "unknown cell type");
      break;
  }

  return nm;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::CreateBlockIdMetadata(vtkModelMetadata *em)
{
  // vtkModelMetadata frees the memory when its done so we need to create a copy
  size_t nblocks = this->BlockInfoMap.size ();
  if (nblocks < 1) return 1;
  em->SetNumberOfBlocks(static_cast<int>(nblocks));

  int *blockIds = new int [nblocks];
  char **blockNames = new char * [nblocks];
  int *numElements = new int [nblocks];
  int *numNodesPerElement = new int [nblocks];
  int *numAttributes = new int [nblocks];

  std::map<int, Block>::const_iterator iter;
  for(iter = this->BlockInfoMap.begin();
      iter != this->BlockInfoMap.end();
      ++iter)
  {
    int index = iter->second.OutputIndex;
    blockIds[index] = iter->first;
    blockNames[index] = vtkExodusIIWriter::GetCellTypeName (iter->second.Type);
    numElements[index] = iter->second.NumElements;
    numNodesPerElement[index] = iter->second.NodesPerElement;
    numAttributes[index] = 0;
  }
  em->SetBlockIds(blockIds);
  em->SetBlockElementType(blockNames);
  em->SetBlockNumberOfElements(numElements);
  em->SetBlockNodesPerElement(numNodesPerElement);
  em->SetBlockNumberOfAttributesPerElement(numAttributes);
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::CreateBlockVariableMetadata (vtkModelMetadata *em)
{
  size_t narrays = this->GlobalVariableMap.size ();
  char **flattenedNames = nullptr;
  if (narrays > 0)
  {
    flattenedNames = vtkExodusIIWriter::FlattenOutVariableNames(
                   this->NumberOfScalarGlobalArrays, this->GlobalVariableMap);
    em->SetGlobalVariableNames (this->NumberOfScalarGlobalArrays, flattenedNames);
  }

  narrays = this->BlockVariableMap.size ();
  char **nms = nullptr;
  if (narrays > 0)
  {
    nms = new char * [narrays];
    int *numComponents = new int [narrays];
    int *scalarIndex = new int [narrays];

    int index = 0;
    std::map<std::string, VariableInfo>::const_iterator var_iter;
    for (var_iter = this->BlockVariableMap.begin ();
         var_iter != this->BlockVariableMap.end ();
         ++var_iter)
    {
      nms[index] =
        vtkExodusIIWriter::StrDupWithNew (var_iter->first.c_str ());
      numComponents[index] = var_iter->second.NumComponents;
      scalarIndex[index] = var_iter->second.ScalarOutOffset;
      index ++;
    }

    flattenedNames = vtkExodusIIWriter::FlattenOutVariableNames(
                   this->NumberOfScalarElementArrays, this->BlockVariableMap);

    // these variables are now owned by em.  No deletion
    em->SetElementVariableInfo(this->NumberOfScalarElementArrays,
                  flattenedNames, static_cast<int>(narrays), nms,
                  numComponents, scalarIndex);
  }

  narrays = this->NodeVariableMap.size ();
  if (narrays > 0)
  {
    nms = new char * [narrays];
    int *numComponents = new int [narrays];
    int *scalarOutOffset = new int [narrays];

    int index = 0;
    std::map<std::string, VariableInfo>::const_iterator iter;
    for (iter = this->NodeVariableMap.begin ();
         iter != this->NodeVariableMap.end ();
         ++iter)
    {
      nms[index] =
            vtkExodusIIWriter::StrDupWithNew (iter->first.c_str ());
      numComponents[index] = iter->second.NumComponents;
      scalarOutOffset[index] = iter->second.ScalarOutOffset;
      index ++;
    }

      flattenedNames = vtkExodusIIWriter::FlattenOutVariableNames(
                    this->NumberOfScalarNodeArrays, this->NodeVariableMap);

      em->SetNodeVariableInfo(this->NumberOfScalarNodeArrays, flattenedNames,
          static_cast<int>(narrays), nms, numComponents, scalarOutOffset);
  }
  return 1;
}

int vtkExodusIIWriter::CreateSetsMetadata (vtkModelMetadata* em)
{
  bool isASideSet = false;
  bool isANodeSet = false;

  if (this->OriginalInput->IsA ("vtkMultiBlockDataSet"))
  {
    int numNodeSets = 0;
    int sumNodes = 0;
    vtkNew<vtkStringArray> nodeSetNames;
    vtkSmartPointer<vtkIntArray> nodeSetIds = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> nodeSetSizes = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> nodeSetNumDF = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> nodeIds = vtkSmartPointer<vtkIntArray>::New ();
    int numSideSets = 0;
    int sumSides = 0;
    vtkNew<vtkStringArray> sideSetNames;
    vtkSmartPointer<vtkIntArray> sideSetIds = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> sideSetSizes = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> sideSetNumDF = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> sideSetElementList = vtkSmartPointer<vtkIntArray>::New ();
    vtkSmartPointer<vtkIntArray> sideSetSideList = vtkSmartPointer<vtkIntArray>::New ();

    vtkMultiBlockDataSet* castObj = vtkMultiBlockDataSet::SafeDownCast(this->OriginalInput);
    vtkSmartPointer<vtkDataObjectTreeIterator> iter;
    iter.TakeReference (castObj->NewTreeIterator ());
    iter->VisitOnlyLeavesOff ();
    iter->TraverseSubTreeOn ();
    int node_id = 0;
    int side_id = 0;
    for (iter->InitTraversal ();
         !iter->IsDoneWithTraversal ();
         iter->GoToNextItem ())
    {
      const char *name = iter->GetCurrentMetaData()->Get (vtkCompositeDataSet::NAME());
      if (iter->GetCurrentDataObject ()->IsA ("vtkMultiBlockDataSet"))
      {
        isASideSet = (name != nullptr && strncmp (name, "Side Sets", 9) == 0);
        isANodeSet = (name != nullptr && strncmp (name, "Node Sets", 9) == 0);
      }
      else if (isANodeSet)
      {
        numNodeSets ++;
        const char* id_str = name != nullptr ? strstr (name, "ID:") : nullptr;
        if (id_str != nullptr)
        {
          id_str += 3;
          node_id = atoi (id_str);
        }
        nodeSetIds->InsertNextTuple1 (node_id);

        if (nodeSetNames->GetNumberOfValues() <= node_id)
        {
          nodeSetNames->SetNumberOfValues ((node_id + 1) * 2);
        }
        nodeSetNames->SetValue (node_id, name);

        node_id ++; // Make sure the node_id is unique if id_str is invalid
        vtkUnstructuredGrid* grid = vtkUnstructuredGrid::SafeDownCast (iter->GetCurrentDataObject ());
        vtkFieldData* field = grid->GetPointData ();
        vtkIdTypeArray* globalIds = vtkArrayDownCast<vtkIdTypeArray>(field ? field->GetArray ("GlobalNodeId") : nullptr);
        if (globalIds)
        {
          int size = 0;
          for (int c = 0; c < grid->GetNumberOfCells (); c ++)
          {
            vtkCell* cell = grid->GetCell (c);
            for (int p = 0; p < cell->GetNumberOfPoints (); p ++)
            {
              size ++;
              vtkIdType pointId = cell->GetPointId (p);
              vtkIdType gid = globalIds->GetValue (pointId);
              nodeIds->InsertNextTuple1 (gid);
              // TODO implement distribution factors
              nodeSetNumDF->InsertNextTuple1 (0);
            }
          }
          nodeSetSizes->InsertNextTuple1 (size);
          sumNodes += size;
        }
        // else
          // {
          // vtkErrorMacro ("We have a node set, but it doesn't have GlobalNodeIDs");
          // }
      }
      else if (isASideSet)
      {
        int hexSides = 0;
        int wedgeSides = 0;
        int otherSides = 0;
        int badSides = 0;
        numSideSets ++;
        const char* id_str = name != nullptr ? strstr (name, "ID:") : nullptr;
        if (id_str != nullptr)
        {
          id_str += 3;
          side_id = atoi (id_str);
        }
        sideSetIds->InsertNextTuple1 (side_id);
        if (sideSetNames->GetNumberOfValues () <= side_id)
        {
          sideSetNames->SetNumberOfValues((side_id + 1) * 2);
        }
        sideSetNames->SetValue (side_id, name);
        side_id ++; // Make sure the side_id is unique if id_str is invalid
        vtkUnstructuredGrid* grid = vtkUnstructuredGrid::SafeDownCast (iter->GetCurrentDataObject ());
        vtkFieldData* field = grid->GetCellData ();
        int cells = grid->GetNumberOfCells ();
        vtkIdTypeArray* sourceElement = vtkArrayDownCast<vtkIdTypeArray>(field ? field->GetArray ("SourceElementId") : nullptr);
        vtkIntArray* sourceSide = vtkArrayDownCast<vtkIntArray>(field ? field->GetArray ("SourceElementSide") : nullptr);
        if (sourceElement && sourceSide)
        {
          for (int c = 0; c < cells; c ++)
          {
            sideSetElementList->InsertNextTuple1 (sourceElement->GetValue (c) + 1);
            switch (GetElementType (sourceElement->GetValue (c) + 1))
            {
              case -1:
              {
                badSides ++;
                break;
              }
              case VTK_WEDGE:
              {
                int wedgeMapping[5] = {3, 4, 0, 1, 2};
                int side = wedgeMapping[sourceSide->GetValue (c)] + 1;
                sideSetSideList->InsertNextTuple1 (side);
                wedgeSides ++;
                break;
              }
              case VTK_HEXAHEDRON:
              {
                int hexMapping[6] = {3, 1, 0, 2, 4, 5};
                sideSetSideList->InsertNextTuple1 (hexMapping[sourceSide->GetValue (c)] + 1);
                hexSides ++;
                break;
              }
              default:
              {
                sideSetSideList->InsertNextTuple1 (sourceSide->GetValue (c) + 1);
                otherSides ++;
                break;
              }
            }
            // TODO implement distribution factors
            sideSetNumDF->InsertNextTuple1 (0);
          }
          sideSetSizes->InsertNextTuple1 (cells);
          sumSides += cells;
        }
        // else
          // {
          // vtkErrorMacro ("We have a side set, but it doesn't have SourceElementId or SourceElementSide");
          // }
      }
    }


    em->SetNumberOfNodeSets (numNodeSets);
    em->SetSumNodesPerNodeSet (sumNodes);
    em->SetNodeSetNames (nodeSetNames);

    int *nodeSetIds_a = new int[nodeSetIds->GetNumberOfTuples ()];
    memcpy (nodeSetIds_a, nodeSetIds->GetPointer (0), nodeSetIds->GetNumberOfTuples () * sizeof(int));
    em->SetNodeSetIds (nodeSetIds_a);

    int *nodeSetSizes_a = new int[nodeSetSizes->GetNumberOfTuples ()];
    memcpy (nodeSetSizes_a, nodeSetSizes->GetPointer (0), nodeSetSizes->GetNumberOfTuples () * sizeof(int));
    em->SetNodeSetSize (nodeSetSizes_a);

    int *nodeSetNumDF_a = new int[nodeSetNumDF->GetNumberOfTuples ()];
    memcpy (nodeSetNumDF_a, nodeSetNumDF->GetPointer (0), nodeSetNumDF->GetNumberOfTuples () * sizeof(int));
    em->SetNodeSetNumberOfDistributionFactors (nodeSetNumDF_a);

    int *nodeIds_a = new int[nodeIds->GetNumberOfTuples ()];
    memcpy (nodeIds_a, nodeIds->GetPointer (0), nodeIds->GetNumberOfTuples() * sizeof(int));
    em->SetNodeSetNodeIdList (nodeIds_a);

    em->SetNumberOfSideSets (numSideSets);
    em->SetSumSidesPerSideSet (sumSides);
    em->SetSideSetNames (sideSetNames);

    int *sideSetIds_a = new int[sideSetIds->GetNumberOfTuples ()];
    memcpy (sideSetIds_a, sideSetIds->GetPointer (0), sideSetIds->GetNumberOfTuples() * sizeof(int));
    em->SetSideSetIds (sideSetIds_a);

    int *sideSetSizes_a = new int[sideSetSizes->GetNumberOfTuples ()];
    memcpy (sideSetSizes_a, sideSetSizes->GetPointer (0), sideSetSizes->GetNumberOfTuples() * sizeof(int));
    em->SetSideSetSize (sideSetSizes_a);

    int *sideSetNumDF_a = new int[sideSetNumDF->GetNumberOfTuples ()];
    memcpy (sideSetNumDF_a, sideSetNumDF->GetPointer (0), sideSetNumDF->GetNumberOfTuples() * sizeof(int));
    em->SetSideSetNumDFPerSide (sideSetNumDF_a);

    int *sideSetElementList_a = new int[sideSetElementList->GetNumberOfTuples ()];
    memcpy (sideSetElementList_a, sideSetElementList->GetPointer (0), sideSetElementList->GetNumberOfTuples() * sizeof(int));
    em->SetSideSetElementList (sideSetElementList_a);

    int *sideSetSideList_a = new int[sideSetSideList->GetNumberOfTuples ()];
    memcpy (sideSetSideList_a, sideSetSideList->GetPointer (0), sideSetSideList->GetNumberOfTuples() * sizeof(int));
    em->SetSideSetSideList (sideSetSideList_a);
  }
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::ParseMetadata ()
{
  vtkModelMetadata* em = this->GetModelMetadata ();
  int nblocks = em->GetNumberOfBlocks ();
  int *ids = em->GetBlockIds();
  int *numAttributes = em->GetBlockNumberOfAttributesPerElement();
  float *att = em->GetBlockAttributes();
  int *attIdx = em->GetBlockAttributesIndex();
  // Extract the attribute data from the meta model.
  for (int n = 0; n < nblocks; n ++)
  {
    std::map<int, Block>::iterator iter = this->BlockInfoMap.find (ids[n]);
    if (iter == this->BlockInfoMap.end ())
    {
      vtkErrorMacro (<< "Unknown id " << ids[n] << " found in meta data");
      return 0;
    }
    iter->second.NumAttributes = numAttributes[n];
    iter->second.BlockAttributes = att + attIdx[n];
  }

  this->ConvertVariableNames (this->GlobalVariableMap);
  this->ConvertVariableNames (this->BlockVariableMap);
  this->ConvertVariableNames (this->NodeVariableMap);
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteInitializationParameters()
{
  vtkModelMetadata *em = this->GetModelMetadata();

  int dim = em->GetDimension();
  int nnsets = em->GetNumberOfNodeSets();
  int nssets = em->GetNumberOfSideSets();
  const char *title = em->GetTitle();
  int numBlocks = em->GetNumberOfBlocks();
  int rc = ex_put_init(this->fid, title, dim, this->NumPoints, this->NumCells,
                       numBlocks, nnsets, nssets);
  return rc >= 0;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteInformationRecords()
{

  vtkModelMetadata *em = this->GetModelMetadata();

  int nlines = em->GetNumberOfInformationLines();

  if (nlines > 0)
  {
    char **lines = nullptr;

    em->GetInformationLines(&lines);

    ex_put_info(this->fid, nlines, lines);
  }

  return 1;
}

template<typename T>
int vtkExodusIIWriterWritePoints (
    std::vector< vtkSmartPointer<vtkUnstructuredGrid> > input,
    int numPoints, int fid)
{
    T *px = new T [numPoints];
    T *py = new T [numPoints];
    T *pz = new T [numPoints];

    int array_index = 0;
    for (size_t i = 0; i < input.size (); i ++)
    {
      vtkPoints *pts = input[i]->GetPoints ();
      if (pts)
      {
        int npts = pts->GetNumberOfPoints ();
        vtkDataArray *da = pts->GetData ();
        for (int j = 0; j < npts; j ++)
        {
          px[array_index] = da->GetComponent(j, 0);
          py[array_index] = da->GetComponent(j, 1);
          pz[array_index] = da->GetComponent(j, 2);
          array_index ++;
        }
      }
    }

    int rc = ex_put_coord(fid, px, py, pz);

    delete [] px;
    delete [] py;
    delete [] pz;

    return rc >= 0;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WritePoints()
{
  if (this->PassDoubles)
  {
    return vtkExodusIIWriterWritePoints<double>(
                    this->FlattenedInput, this->NumPoints, this->fid);
  }
  else
  {
    return vtkExodusIIWriterWritePoints<float>(
                    this->FlattenedInput, this->NumPoints, this->fid);
  }
}

//---------------------------------------------------------
// Points and point IDs, element IDs
//---------------------------------------------------------
int vtkExodusIIWriter::WriteCoordinateNames()
{
  vtkModelMetadata *em = this->GetModelMetadata();

  int rc = ex_put_coord_names(this->fid, em->GetCoordinateNames());

  return rc >= 0;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteGlobalPointIds()
{
  if (!this->AtLeastOneGlobalNodeIdList)
  {
    return 1;
  }
  int *copyOfIds = new int [this->NumPoints];
  int index = 0;
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    vtkIdType npoints = this->FlattenedInput[i]->GetNumberOfPoints();

    vtkIdType *ids = this->GlobalNodeIdList[i];
    if (ids)
    {
      for (int j=0; j<npoints; j++)
      {
        copyOfIds[index++] = static_cast<int>(ids[j]);
      }
    }
    else
    {
      for (int j=0; j<npoints; j++)
      {
        copyOfIds[index++] = 0;
      }
    }
  }
  int rc = ex_put_node_num_map(this->fid, copyOfIds);
  delete [] copyOfIds;

  return rc >= 0;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteBlockInformation()
{
  int rc;

  std::map<int, Block>::const_iterator blockIter;
  size_t nblocks = this->BlockInfoMap.size ();

  std::vector<int*> connectivity (nblocks);

  // Use this to copy the attributes into if we need it.
  std::vector<double*> attributesD (nblocks);

  // For each block, a map from element global ID to it's location
  // within it's block in the ExodusModel object.
  for (blockIter = this->BlockInfoMap.begin ();
       blockIter != this->BlockInfoMap.end ();
       ++blockIter)
  {
    int outputIndex = blockIter->second.OutputIndex;
    int numElts = blockIter->second.NumElements;
    int numAtts = blockIter->second.NumAttributes;
    int numNodes = blockIter->second.NodesPerElement;

    int numPoints;
    if (numNodes == 0 &&
        // Number of elements is 0 if all data is on one process
        numElts > 0)
    {
      numPoints = blockIter->second.EntityNodeOffsets[numElts - 1]
                + blockIter->second.EntityCounts[numElts - 1];
    }
    else
    {
      numPoints = numElts * numNodes;
    }

    if (numElts > 0)
    {
      connectivity[outputIndex] = new int [numPoints];

      if (numAtts > 0)
      {
        if (this->PassDoubles)
        {
          attributesD[outputIndex] = new double [numElts * numAtts];
        }
      }
    }
    else
    {
      connectivity[outputIndex] = nullptr;
      attributesD[outputIndex] = nullptr;
    }
  }

  // Prepare the pointers as flat arrays for Exodus
  // connectivity and attributes, if we need doubles.
  int pointOffset = 0;
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    vtkCellArray *ca = this->FlattenedInput[i]->GetCells();
    vtkIdType *ptIds = nullptr;
    if (ca)
    {
      ptIds = ca->GetPointer ();
    }
    vtkIdTypeArray *loca = this->FlattenedInput[i]->GetCellLocationsArray();
    vtkIdType *loc = nullptr;
    if (loca)
    {
      loc = loca->GetPointer(0);
    }

    int ncells = this->FlattenedInput[i]->GetNumberOfCells();
    for (int j = 0; j < ncells; j++)
    {
      int blockId = this->BlockIdList[i]->GetValue(j);//CLM
      int blockOutIndex = this->BlockInfoMap[blockId].OutputIndex;

      int nodesPerElement = this->BlockInfoMap[blockId].NodesPerElement;
      vtkIdType elementOffset = this->CellToElementOffset[i][j];
      int offset;
      if (nodesPerElement == 0)
      {
        offset = this->BlockInfoMap[blockId].EntityNodeOffsets[j];
      }
      else
      {
        offset = elementOffset * nodesPerElement;
      }

      // the block connectivity array
      vtkIdType ptListIdx = loc[j];
      vtkIdType npts = ptIds[ptListIdx++];

      switch (this->FlattenedInput[i]->GetCellType (j))
      {
        case VTK_VOXEL: // reorder to exodus HEX type
          connectivity[blockOutIndex][offset + 0] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 1] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 3] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 2] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 4] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 5] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 7] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          connectivity[blockOutIndex][offset + 6] = pointOffset + (int) ptIds[ptListIdx++] + 1;
          break;
        default:
          for (vtkIdType p=0; p<npts; p++)
          {
            int ExodusPointId = pointOffset + (int) ptIds[ptListIdx++] + 1;
            connectivity[blockOutIndex][offset + p] = ExodusPointId;
          }
      }

      // the block element attributes
      float *att = this->BlockInfoMap[blockId].BlockAttributes;

      int numAtts = this->BlockInfoMap[blockId].NumAttributes;

      if ((numAtts == 0) || (att == nullptr)) continue;

      int attOff = (elementOffset * numAtts); // location for the element in the block

      if (this->PassDoubles)
      {
        for (int k = 0; k <numAtts; k++)
        {
          int off = attOff + k;
          // TODO verify the assumption that ModelMetadata stores
          // elements in the same order we do
          // Could probably use the global node id? but how global is global.
          attributesD[blockOutIndex][off] = static_cast<double>(att[off]);
        }
      }
    }
    pointOffset += this->FlattenedInput[i]->GetNumberOfPoints ();
  }

  // Now, finally, write out the block information
  int fail = 0;
  for (blockIter = this->BlockInfoMap.begin ();
       blockIter != this->BlockInfoMap.end ();
       ++blockIter)
  {
    char *type_name = vtkExodusIIWriter::GetCellTypeName (blockIter->second.Type);
    if (blockIter->second.NodesPerElement == 0 &&
        blockIter->second.NumElements > 0)
    {
      int numElts = blockIter->second.NumElements;
      int numPoints = blockIter->second.EntityNodeOffsets[numElts - 1]
                    + blockIter->second.EntityCounts[numElts - 1];
      rc = ex_put_elem_block(this->fid, blockIter->first,
                  type_name,
                  blockIter->second.NumElements,
                  numPoints,
                  blockIter->second.NumAttributes);
    }
    else
    {
      rc = ex_put_elem_block(this->fid, blockIter->first,
                  type_name,
                  blockIter->second.NumElements,
                  blockIter->second.NodesPerElement,
                  blockIter->second.NumAttributes);
    }
    delete [] type_name;

    if (rc < 0)
    {
      vtkErrorMacro (<< "Problem adding block with id " << blockIter->first);
      continue;
    }

    ex_put_name (this->fid, EX_ELEM_BLOCK, blockIter->first, blockIter->second.Name);

    if (blockIter->second.NumElements > 0)
    {
      rc = ex_put_elem_conn(this->fid, blockIter->first,
                      connectivity[blockIter->second.OutputIndex]);

      if (rc < 0)
      {
        vtkErrorMacro (<< "Problem writing connectivity " << blockIter->first);
        continue;
      }

      if (blockIter->second.NumAttributes != 0)
      {
        if (this->PassDoubles)
        {
          rc = ex_put_elem_attr(this->fid, blockIter->first,
                          attributesD[blockIter->second.OutputIndex]);
        }
        else
        {
          // TODO verify the assumption that ModelMetadata stores
          // elements in the same order we do
          rc = ex_put_elem_attr(this->fid, blockIter->first,
                            blockIter->second.BlockAttributes);
        }

        if (rc < 0)
        {
          continue;
        }
      }

      if (blockIter->second.NodesPerElement == 0)
      {
        rc = ex_put_entity_count_per_polyhedra (this->fid, EX_ELEM_BLOCK,
                     blockIter->first, &(blockIter->second.EntityCounts[0]));
      }
    }
  }

  for (size_t n = 0; n < nblocks; n ++)
  {
    delete [] connectivity[n];
    if (this->PassDoubles)
    {
      delete [] attributesD[n];
    }
  }
  return !fail;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteGlobalElementIds()
{
  int rc = 0;

  //if (sizeof(vtkIdType) != sizeof(int))
  //  {
  //  vtkErrorMacro(<<"vtkExodusIIWriter::WriteGlobalElementIds cannot convert vtkIdType to int.");
  //  return -1;
  //  }

  if (this->AtLeastOneGlobalElementIdList)
  {
    int *copyOfIds = new int [this->NumCells];
    memset (copyOfIds, 0, sizeof (int) * this->NumCells);
    for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
    {
      vtkIdType *ids = this->GlobalElementIdList[i];
      if (ids)
      {
        int ncells = this->FlattenedInput[i]->GetNumberOfCells();
        for (int j=0; j<ncells; j++)
        {
          int blockId = this->BlockIdList[i]->GetValue (j);
          int start = this->BlockInfoMap[blockId].ElementStartIndex;
          vtkIdType offset = this->CellToElementOffset[i][j];
          copyOfIds[start + offset] = static_cast<int>(ids[j]);
        }
      }
    }
    rc = ex_put_elem_num_map(this->fid, copyOfIds);
    delete [] copyOfIds;
  }

  return rc >= 0;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteVariableArrayNames()
{
  int rc = 0;

  //  1. We convert vector arrays to individual scalar arrays, using
  //     their original names if we have those.
  //  2. For the element variables, create the element/block truth table.

  // GLOBAL VARIABLES
  const char **outputArrayNames;
  if (this->NumberOfScalarGlobalArrays > 0)
  {
    outputArrayNames = new const char * [this->NumberOfScalarGlobalArrays];
    std::map<std::string, VariableInfo>::const_iterator iter;
    for (iter = this->GlobalVariableMap.begin ();
         iter != this->GlobalVariableMap.end ();
         ++iter)
    {
      int off = iter->second.ScalarOutOffset;
      for (int j=0; j< iter->second.NumComponents; j++)
      {
        outputArrayNames[off + j] = iter->second.OutNames[j].c_str ();
      }
    }

    rc = ex_put_var_param(this->fid, "G", this->NumberOfScalarGlobalArrays);
    if (rc < 0)
    {
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames cell variables");
      delete [] outputArrayNames;
      return 0;
    }

    rc = ex_put_var_names(this->fid, "G", this->NumberOfScalarGlobalArrays,
                          const_cast<char **>(outputArrayNames));
                          // This should be treating this read only... hopefully
    if (rc < 0)
    {
      delete [] outputArrayNames;
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames cell variables");
      return 0;
    }

    delete [] outputArrayNames;
  }

  // CELL (ELEMENT) VARIABLES
  if (this->NumberOfScalarElementArrays > 0 && this->NumCells > 0)
  {
    outputArrayNames = new const char * [this->NumberOfScalarElementArrays];

    std::map<std::string, VariableInfo>::const_iterator iter;
    for (iter = this->BlockVariableMap.begin ();
         iter != this->BlockVariableMap.end ();
         ++iter)
    {
      int off = iter->second.ScalarOutOffset;
      for (int j=0; j< iter->second.NumComponents; j++)
      {
        outputArrayNames[off + j] = iter->second.OutNames[j].c_str ();
      }
    }

    rc = ex_put_var_param(this->fid, "E", this->NumberOfScalarElementArrays);
    if (rc < 0)
    {
      delete [] outputArrayNames;
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames cell variables");
      return 0;
    }

    rc = ex_put_var_names(this->fid, "E", this->NumberOfScalarElementArrays,
                          const_cast<char **>(outputArrayNames));
                          // This should be treating this read only... hopefully
    if (rc < 0)
    {
      delete [] outputArrayNames;
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames cell variables");
      return 0;
    }

    rc = ex_put_elem_var_tab(this->fid,
                             static_cast<int>(this->BlockInfoMap.size ()),
                             this->NumberOfScalarElementArrays,
                             this->BlockElementVariableTruthTable);
    if (rc < 0)
    {
      delete [] outputArrayNames;
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames cell variables");
      return 0;
    }

    delete [] outputArrayNames;
  }

  // POINT (NODE) VARIABLES
  if (this->NumberOfScalarNodeArrays > 0 && this->NumPoints > 0)
  {
    outputArrayNames = new const char * [this->NumberOfScalarNodeArrays];

    std::map<std::string, VariableInfo>::const_iterator iter;
    for (iter = this->NodeVariableMap.begin ();
         iter != this->NodeVariableMap.end ();
         ++iter)
    {
      int off = iter->second.ScalarOutOffset;
      for (int j=0; j<iter->second.NumComponents; j++)
      {
        outputArrayNames[off + j] = iter->second.OutNames[j].c_str ();
      }
    }

    rc = ex_put_var_param(this->fid, "N", this->NumberOfScalarNodeArrays);

    if (rc < 0)
    {
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames " <<
        "failure to write " << this->NumberOfScalarNodeArrays << " arrays");
      delete [] outputArrayNames;
      return 0;
    }

    rc = ex_put_var_names(this->fid, "N", this->NumberOfScalarNodeArrays,
                          const_cast<char **>(outputArrayNames));
                          // This should not save references... hopefully
    if (rc < 0)
    {
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames " <<
        "failure to write the array names");
      delete [] outputArrayNames;
      return 0;
    }

    delete [] outputArrayNames;
  }

  // GLOBAL VARIABLES
/*
  int ngvars = mmd->GetNumberOfGlobalVariables();

  if (ngvars > 0)
    {
    char **names = mmd->GetGlobalVariableNames();

    rc = ex_put_var_param(this->fid, "G", ngvars);

    if (rc == 0)
      {
      rc = ex_put_var_names(this->fid, "G", ngvars, names);
      }

    if (rc < 0)
      {
      vtkErrorMacro(<<
        "vtkExodusIIWriter::WriteVariableArrayNames global variables");
      return 0;
      }
    }
*/
  return 1;
}

//----------------------------------------------------------------------------
void vtkExodusIIWriter::ConvertVariableNames(
             std::map<std::string, VariableInfo>& variableMap)
{
  std::map<std::string, VariableInfo>::iterator varIter;
  // Global output variable names
  for (varIter = variableMap.begin ();
       varIter != variableMap.end ();
       ++varIter)
  {
    int numComp = varIter->second.NumComponents;
    if (numComp == 1)
    {
      varIter->second.OutNames[0] = std::string (varIter->first);
    }
    else
    {
      for (int component = 0; component < numComp; component ++)
      {
        varIter->second.OutNames[component] =
                CreateNameForScalarArray (varIter->first.c_str (),
                                          component,
                                          numComp);
      }
    }
  }
}

char **vtkExodusIIWriter::FlattenOutVariableNames(
             int nScalarArrays,
             const std::map<std::string, VariableInfo>& variableMap)
{
  char **newNames = new char * [nScalarArrays];

  std::map<std::string, VariableInfo>::const_iterator iter;
  for (iter = variableMap.begin ();
       iter != variableMap.end ();
       ++iter)
  {
    for (int component = 0; component < iter->second.NumComponents; component ++)
    {
      int index = iter->second.ScalarOutOffset + component;
      newNames[index] =
          StrDupWithNew (
              CreateNameForScalarArray(iter->first.c_str (),
                                       component,
                                       iter->second.NumComponents).c_str ());
    }
  }

  return newNames;
}

//----------------------------------------------------------------------------
std::string vtkExodusIIWriter::CreateNameForScalarArray(
  const char *root, int component, int numComponents)
{
  // Naming conventions chosen to match ExodusIIReader expectations
  if (component >= numComponents)
  {
    vtkErrorMacro ("CreateNameForScalarArray: Component out of range");
    return std::string ();
  }
  if (numComponents == 1)
  {
    return std::string (root);
  }
  else if (numComponents <= 2)
  {
    std::string s (root);
    switch (component)
    {
      case 0:
        s.append ("_R");
        break;
      case 1:
        s.append ("_Z");
        break;
    }
    return s;
  }
  else if (numComponents <= 3)
  {
    std::string s (root);
    switch (component)
    {
      case 0:
        s.append ("X");
        break;
      case 1:
        s.append ("Y");
        break;
      case 2:
        s.append ("Z");
        break;
    }
    return s;
  }
  else if (numComponents <= 6)
  {
    std::string s (root);
    switch (component)
    {
      case 0:
        s.append ("XX");
        break;
      case 1:
        s.append ("XY");
        break;
      case 2:
        s.append ("XZ");
        break;
      case 3:
        s.append ("YY");
        break;
      case 4:
        s.append ("YZ");
        break;
      case 5:
        s.append ("ZZ");
        break;
    }
    return s;
  }
  else
  {
    std::string s (root);
    // assume largest for 32 bit decimal representation
    char n[11];
    snprintf (n, sizeof(n), "%10d", component);
    s.append (n);
    return s;
  }
}

//----------------------------------------------------------------------------
vtkIdType vtkExodusIIWriter::GetNodeLocalId(vtkIdType id)
{
  if (!this->LocalNodeIdMap)
  {
    this->LocalNodeIdMap = new std::map<vtkIdType, vtkIdType>;
    vtkIdType index = 0;
    for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
    {
      vtkIdType npoints = this->FlattenedInput[i]->GetNumberOfPoints();
      vtkIdType *ids = this->GlobalNodeIdList[i];
      if (ids)
      {
        for (int j=0; j<npoints; j++)
        {
          this->LocalNodeIdMap->insert(
            std::map<vtkIdType,vtkIdType>::value_type(ids[j], index));
          index ++;
        }
      }
      else
      {
        index += npoints;
      }
    }
  }

  std::map<vtkIdType,vtkIdType>::iterator mapit =
    this->LocalNodeIdMap->find(id);

  if (mapit == this->LocalNodeIdMap->end())
  {
    return -1;
  }
  else
  {
    return mapit->second;
  }
}

//-----------------------------------------------------------------------
// Side sets and node sets
//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteNodeSetInformation()
{
  int rc = 0;
  int i, j;

  vtkModelMetadata *em = this->GetModelMetadata();

  int nnsets = em->GetNumberOfNodeSets();

  if (nnsets < 1) return 1;

  int nids = em->GetSumNodesPerNodeSet();

  if (nids < 1 || !this->AtLeastOneGlobalNodeIdList)
  {
    int *buf = new int [nnsets];

    memset(buf, 0, sizeof(int) * nnsets);

    rc = ex_put_concat_node_sets(this->fid, em->GetNodeSetIds(),
              buf, buf, buf, buf, nullptr, nullptr);

    delete [] buf;

    return (rc >= 0);
  }

  int *nsSize = new int [nnsets];
  int *nsNumDF = new int [nnsets];
  int *nsIdIdx = new int [nnsets];
  int *nsDFIdx = new int [nnsets];

  int ndf = em->GetSumDistFactPerNodeSet();

  int *idBuf = new int [nids];
  float *dfBuf = nullptr;
  double *dfBufD = nullptr;

  if (ndf)
  {
    if (this->PassDoubles)
    {
      dfBufD = new double [ndf];
    }
    else
    {
      dfBuf = new float [ndf];
    }
  }

  int *emNsSize = em->GetNodeSetSize();
  int *emNumDF = em->GetNodeSetNumberOfDistributionFactors();

  int nextId = 0;
  int nextDF = 0;

  int *ids = em->GetNodeSetNodeIdList();
  float *df = em->GetNodeSetDistributionFactors();

  for (i=0; i<nnsets; i++)
  {
    nsSize[i] = 0;
    nsNumDF[i] = 0;

    nsIdIdx[i] = nextId;
    nsDFIdx[i] = nextDF;

    for (j=0; j< emNsSize[i]; j++)
    {
      int lid = this->GetNodeLocalId(*ids);
      ids ++;

      if (lid < 0) continue;

      nsSize[i]++;
      idBuf[nextId++] = lid + 1;

      if (*emNumDF > 0)
      {
        nsNumDF[i]++;

        if (this->PassDoubles)
        {
          dfBufD[nextDF++] = (double)*df;
        }
        else
        {
          dfBuf[nextDF++] = *df;
        }
      }
    }
      df ++;
      emNumDF ++;
  }

  int *node_ids = em->GetNodeSetIds();

  if (this->PassDoubles)
  {
    rc = ex_put_concat_node_sets(this->fid, node_ids,
              nsSize, nsNumDF, nsIdIdx, nsDFIdx, idBuf, dfBufD);
  }
  else
  {
    rc = ex_put_concat_node_sets(this->fid, node_ids,
              nsSize, nsNumDF, nsIdIdx, nsDFIdx, idBuf, dfBuf);
  }

  for (i = 0; i < nnsets; i ++) {
    vtkStdString name = em->GetNodeSetNames ()->GetValue (node_ids[i]);
    ex_put_name (this->fid, EX_NODE_SET, node_ids[i], name.c_str());
  }


  delete [] nsSize;
  delete [] nsNumDF;
  delete [] nsIdIdx;
  delete [] nsDFIdx;
  delete [] idBuf;
  delete [] dfBuf;
  delete [] dfBufD;

  return (rc >= 0);
}

//----------------------------------------------------------------------------
vtkIdType vtkExodusIIWriter::GetElementLocalId(vtkIdType id)
{
  if (!this->LocalElementIdMap)
  {
    this->LocalElementIdMap = new std::map<vtkIdType, vtkIdType>;
    for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
    {
      if (this->GlobalElementIdList[i])
      {
        vtkIdType ncells = this->FlattenedInput[i]->GetNumberOfCells();
        for (vtkIdType j=0; j<ncells; j++)
        {
          vtkIdType gid = this->GlobalElementIdList[i][j];
          std::map<vtkIdType,vtkIdType>::iterator mapit = this->LocalElementIdMap->find(gid);
          if (mapit != this->LocalElementIdMap->end ())
          {
            vtkErrorMacro ("Overlapping gids in the dataset");
            continue;
          }

          std::map<int, Block>::const_iterator blockIter =
                  this->BlockInfoMap.find (this->BlockIdList[i]->GetValue (j));
          if (blockIter == this->BlockInfoMap.end ())
          {
            vtkWarningMacro ("vtkExodusIIWriter: The block id map has come out of sync");
            continue;
          }

          int index = blockIter->second.ElementStartIndex + CellToElementOffset[i][j];
          this->LocalElementIdMap->insert(
            std::map<vtkIdType, vtkIdType>::value_type(gid, index));
        }
      }
    }
  }

  std::map<vtkIdType,vtkIdType>::iterator mapit = this->LocalElementIdMap->find(id);
  if (mapit == this->LocalElementIdMap->end())
  {
      return -1;
  }
  else
  {
      return mapit->second;
  }
}

int vtkExodusIIWriter::GetElementType(vtkIdType id)
{
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    if (this->GlobalElementIdList[i])
    {
      vtkIdType ncells = this->FlattenedInput[i]->GetNumberOfCells();
      for (vtkIdType j=0; j<ncells; j++)
      {
        vtkIdType gid = this->GlobalElementIdList[i][j];
        if (gid == id)
        {
          return this->FlattenedInput[i]->GetCellType (j);
        }
      }
    }
  }
  return -1;
}

//-----------------------------------------------------------------------
int vtkExodusIIWriter::WriteSideSetInformation()
{
  int i, j, k;
  int rc= 0;

  vtkModelMetadata *em = this->GetModelMetadata();

  int nssets = em->GetNumberOfSideSets();

  if (nssets < 1) return 1;

  // Cells are written out to file in a different order than
  // they appear in the input. We need a mapping from their internal
  // id in the input to their internal id in the output.

  int nids = em->GetSumSidesPerSideSet();

  if (nids < 1)
  {
    int *buf = new int [nssets];

    memset(buf, 0, sizeof(int) * nssets);

    rc = ex_put_concat_side_sets(this->fid, em->GetSideSetIds(),
              buf, buf, buf, buf, nullptr, nullptr, nullptr);

    delete [] buf;

    return (rc >= 0);
  }

  int *ssSize = new int [nssets];
  int *ssNumDF = new int [nssets];
  int *ssIdIdx = new int [nssets];
  int *ssDFIdx = new int [nssets];

  int ndf = em->GetSumDistFactPerSideSet();

  int *idBuf = new int [nids];
  int *sideBuf = new int [nids];
  float *dfBuf = nullptr;
  double *dfBufD = nullptr;

  if (ndf)
  {
    if (this->PassDoubles)
    {
      dfBufD = new double [ndf];
    }
    else
    {
      dfBuf = new float [ndf];
    }
  }

  int *emSsSize = em->GetSideSetSize();
  int *emDFIdx = em->GetSideSetDistributionFactorIndex();

  int nextId = 0;
  int nextDF = 0;

  int *ids = em->GetSideSetElementList();
  int *sides = em->GetSideSetSideList();
  int *numDFPerSide = em->GetSideSetNumDFPerSide();

  for (i=0; i<nssets; i++)
  {
    ssSize[i] = 0;
    ssNumDF[i] = 0;

    ssIdIdx[i] = nextId;
    ssDFIdx[i] = nextDF;

    if (emSsSize[i] == 0) continue;

    float *df = nullptr;

    if (ndf > 0)
    {
      df = em->GetSideSetDistributionFactors() + emDFIdx[i];
    }

    for (j=0; j< emSsSize[i]; j++)
    {
      // Have to check if this element is still in the ugrid.
      // It may have been deleted since the ExodusModel was created.

      int lid = this->GetElementLocalId(*ids);
      ids ++;

      if (lid >= 0)
      {
        ssSize[i]++;

        idBuf[nextId] = lid + 1;
        sideBuf[nextId] = *sides;
        sides ++;

        nextId++;

        if (*numDFPerSide > 0)
        {
          ssNumDF[i] += *numDFPerSide;

          if (this->PassDoubles)
          {
            for (k=0; k < *numDFPerSide; k++)
            {
              dfBufD[nextDF++] = (double)df[k];
            }
          }
          else
          {
            for (k=0; k < *numDFPerSide; k++)
            {
              dfBuf[nextDF++] = df[k];
            }
          }
        }
      }

      if (df) df += *numDFPerSide;
      numDFPerSide ++;
    }
  }

  int *sids = em->GetSideSetIds ();

  if (this->PassDoubles)
  {
    rc = ex_put_concat_side_sets(this->fid, sids,
            ssSize, ssNumDF, ssIdIdx, ssDFIdx, idBuf, sideBuf, dfBufD);
  }
  else
  {
    rc = ex_put_concat_side_sets(this->fid, sids,
            ssSize, ssNumDF, ssIdIdx, ssDFIdx, idBuf, sideBuf, dfBuf);
  }

  for (i = 0; i < nssets; i ++) {
    vtkStdString name = em->GetSideSetNames ()->GetValue (sids[i]);
    ex_put_name (this->fid, EX_SIDE_SET, sids[i], name.c_str());
  }

  delete [] ssSize;
  delete [] ssNumDF;
  delete [] ssIdIdx;
  delete [] ssDFIdx;
  delete [] idBuf;
  delete [] sideBuf;
  delete [] dfBuf;
  delete [] dfBufD;

  return rc >= 0;
}


//----------------------------------------------------------------------------
int vtkExodusIIWriter::BlockVariableTruthValue(int blockIdx, int varIdx)
{
  int tt=0;
  int nvars = this->NumberOfScalarElementArrays;
  int nblocks = static_cast<int>(this->BlockInfoMap.size ());

  if ( (blockIdx >= 0) && (blockIdx < nblocks) &&
       (varIdx >= 0) && (varIdx < nvars))
  {
    tt = this->BlockElementVariableTruthTable[(blockIdx * nvars) + varIdx];
  }
  else
  {
    vtkWarningMacro(<< "vtkExodusIIWriter::BlockVariableTruthValue invalid index");
  }

  return tt;
}

//-----------------------------------------------------------------------
// Properties
//-----------------------------------------------------------------------
int vtkExodusIIWriter::WriteProperties()
{
  int rc = 0;
  int i;

  vtkModelMetadata *em = this->GetModelMetadata();

  int nbprop = em->GetNumberOfBlockProperties();
  int nnsprop = em->GetNumberOfNodeSetProperties();
  int nssprop = em->GetNumberOfSideSetProperties();

  if (nbprop)
  {
    char **names = em->GetBlockPropertyNames();

    // Exodus library "feature".  By convention there is a property
    // array called "ID", the value of which is the ID of the block,
    // node set or side set.  This property is special.  For example,
    // if you change the property value for a block, that block's
    // block ID is changed.  I had no idea *how* special this property
    // was, however.  If you use ex_put_prop_names to tell the library
    // what your property names are, and "ID" happens to be one of those
    // names, then the library fills out the whole property array for
    // you.  Then if you follow this call with ex_put_prop_array for
    // each property array, including "ID", you get *two* arrays named
    // "ID".  This is not documented, and totally unexpected.
    //
    // ex_put_prop_names is not required, it's just more efficient to
    // call it before all the ex_put_prop_array calls.  So we are
    // not going to call it.
    //
    // rc = ex_put_prop_names(this->fid, EX_ELEM_BLOCK, nbprop, names);

    if (rc >= 0)
    {
      int *values = em->GetBlockPropertyValue();

      for (i=0; i<nbprop; i++)
      {
        rc = ex_put_prop_array(this->fid, EX_ELEM_BLOCK, names[i], values);
        if (rc) break;
        // TODO Handle the addition of Blocks not known by the metadata
        values += this->BlockInfoMap.size ();
      }
    }
  }

  if (!rc && nnsprop)
  {
    char **names = em->GetNodeSetPropertyNames();
    int nnsets = em->GetNumberOfNodeSets();

    // rc = ex_put_prop_names(this->fid, EX_NODE_SET, nnsprop, names);

    if (rc >= 0)
    {
      int *values = em->GetNodeSetPropertyValue();

      for (i=0; i<nnsprop; i++)
      {
        rc = ex_put_prop_array(this->fid, EX_NODE_SET, names[i], values);
        if (rc) break;
        values += nnsets;
      }
    }
  }

  if (!rc && nssprop)
  {
    char **names = em->GetSideSetPropertyNames();
    int nssets = em->GetNumberOfSideSets();

    // rc = ex_put_prop_names(this->fid, EX_SIDE_SET, nssprop, names);

    if (rc >= 0)
    {
      int *values = em->GetSideSetPropertyValue();

      for (i=0; i<nssprop; i++)
      {
        rc = ex_put_prop_array(this->fid, EX_SIDE_SET, names[i], values);
        if (rc) break;
        values += nssets;
      }
    }
  }

  return (rc >= 0);
}

//========================================================================
//   VARIABLE ARRAYS:
//========================================================================
template <typename iterT>
double vtkExodusIIWriterGetComponent (iterT* it, vtkIdType ind)
{
  vtkVariant v(it->GetValue (ind));
  return v.ToDouble ();
}

//----------------------------------------------------------------------------
double vtkExodusIIWriter::ExtractGlobalData (const char *name, int comp, int ts)
{
  double ret = 0.0;
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    // find the first block that matches this global data.  Assumes it's global.
    vtkDataArray *da = this->FlattenedInput[i]->GetFieldData ()->GetArray (name);
    if (da)
    {
      int numTup = da->GetNumberOfTuples ();
      if (numTup == 1)
      {
        ret = da->GetComponent (0, comp);
      }
      // Exodus doesn't support multiple tuples on the global values.
      // But the ExodusIIReader reads all timesteps into the field array
      // at every time step.  This will assume that if we have multiple tuples
      // in the array they are from an exodus file so we'll output them
      // back as expected on another read.  Not perfect...
      else if (ts < numTup)
      {
        ret = da->GetComponent (ts, comp);
      }
    }
  }
  return ret;
}

//----------------------------------------------------------------------------
void vtkExodusIIWriter::ExtractCellData (const char *name, int comp,
                                         vtkDataArray *buffer)
{
  buffer->SetNumberOfTuples (this->NumCells);
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    vtkDataArray *da = this->FlattenedInput[i]->GetCellData ()->GetArray (name);
    int ncells = this->FlattenedInput[i]->GetNumberOfCells ();
    if (da)
    {
      vtkArrayIterator *arrayIter = da->NewIterator ();
      vtkIdType ncomp = da->GetNumberOfComponents ();
      for (vtkIdType j = 0; j < ncells; j ++)
      {
        std::map<int, Block>::const_iterator blockIter =
                this->BlockInfoMap.find (this->BlockIdList[i]->GetValue (j));
        if (blockIter == this->BlockInfoMap.end ())
        {
          vtkWarningMacro ("vtkExodusIIWriter: The block id map has come out of sync");
          continue;
        }
        int index = blockIter->second.ElementStartIndex +
                CellToElementOffset[i][j];
        switch (da->GetDataType ())
        {
          vtkArrayIteratorTemplateMacro(
            buffer->SetTuple1 (index,
                               vtkExodusIIWriterGetComponent (
                                    static_cast<VTK_TT*>(arrayIter),
                                    j * ncomp + comp)));
        }
      }
      arrayIter->Delete ();
    }
    else
    {
      for (vtkIdType j = 0; j < ncells; j ++)
      {
        std::map<int, Block>::const_iterator blockIter =
                this->BlockInfoMap.find (this->BlockIdList[i]->GetValue (j));
        if (blockIter == this->BlockInfoMap.end ())
        {
          vtkWarningMacro ("vtkExodusIIWriter: The block id map has come out of sync");
          continue;
        }
        int index = blockIter->second.ElementStartIndex +
                CellToElementOffset[i][j];
        buffer->SetTuple1 (index, 0);
      }
    }
  }
}

//----------------------------------------------------------------------------
void vtkExodusIIWriter::ExtractPointData (const char *name, int comp,
                                          vtkDataArray* buffer)
{
  buffer->SetNumberOfTuples (this->NumPoints);
  int index = 0;
  for (size_t i = 0; i < this->FlattenedInput.size (); i ++)
  {
    vtkDataArray *da = this->FlattenedInput[i]->GetPointData ()->GetArray (name);
    if (da)
    {
      vtkArrayIterator *iter = da->NewIterator ();
      vtkIdType ncomp = da->GetNumberOfComponents ();
      vtkIdType nvals = ncomp * da->GetNumberOfTuples ();
      for (vtkIdType j = comp; j < nvals; j += ncomp)
      {
        switch (da->GetDataType ())
        {
          vtkArrayIteratorTemplateMacro(
            buffer->SetTuple1 (index++,
                               vtkExodusIIWriterGetComponent (
                                        static_cast<VTK_TT*>(iter), j)));
        }
      }
      iter->Delete ();
    }
    else
    {
      vtkIdType nvals = this->FlattenedInput[i]->GetNumberOfPoints ();
      for (vtkIdType j = 0; j < nvals; j ++)
      {
        buffer->SetTuple1 (index ++, 0);
      }
    }
  }
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteGlobalData (int timestep, vtkDataArray *buffer)
{
  std::map<std::string, VariableInfo>::const_iterator varIter;
  buffer->Initialize ();
  buffer->SetNumberOfComponents (1);
  buffer->SetNumberOfTuples (this->NumberOfScalarGlobalArrays);
  for (varIter = this->GlobalVariableMap.begin ();
       varIter != this->GlobalVariableMap.end ();
       ++varIter)
  {
    const char *nameIn = varIter->first.c_str ();
    int numComp = varIter->second.NumComponents;
    for (int component = 0; component < numComp; component ++)
    {
      double val = this->ExtractGlobalData (nameIn, component, timestep);
      int varOutIndex = varIter->second.ScalarOutOffset + component;
      buffer->SetComponent (varOutIndex, 0, val);
    }
  }
  int rc;
  if (buffer->IsA ("vtkDoubleArray"))
  {
    vtkDoubleArray *da = vtkArrayDownCast<vtkDoubleArray>(buffer);
    rc = ex_put_glob_vars (this->fid, timestep + 1,
                  this->NumberOfScalarGlobalArrays, da->GetPointer (0));
  }
  else /* (buffer->IsA ("vtkFloatArray")) */
  {
    vtkFloatArray *fa = vtkArrayDownCast<vtkFloatArray>(buffer);
    rc = ex_put_glob_vars (this->fid, timestep + 1,
                  this->NumberOfScalarGlobalArrays, fa->GetPointer (0));
  }
  if (rc < 0)
  {
    vtkErrorMacro(<< "vtkExodusIIWriter::WriteNextTimeStep glob vars");
    return 0;
  }
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteCellData (int timestep, vtkDataArray *buffer)
{
  std::map<std::string, VariableInfo>::const_iterator varIter;
  for (varIter = this->BlockVariableMap.begin ();
       varIter != this->BlockVariableMap.end ();
       ++varIter)
  {
    const char *nameIn = varIter->first.c_str ();
    int numComp = varIter->second.NumComponents;

    for (int component = 0; component < numComp; component ++)
    {
      buffer->Initialize ();
      this->ExtractCellData(nameIn, component, buffer);
      int varOutIndex = varIter->second.ScalarOutOffset + component;

      std::map<int, Block>::const_iterator blockIter;
      for (blockIter = this->BlockInfoMap.begin ();
           blockIter != this->BlockInfoMap.end ();
           ++blockIter)
      {
        int numElts = blockIter->second.NumElements;
        if (numElts < 1) continue;   // no cells in this block

        int defined = this->BlockVariableTruthValue(
                        blockIter->second.OutputIndex, varOutIndex);
        if (!defined) continue;    // var undefined in this block

        int id = blockIter->first;
        int start = blockIter->second.ElementStartIndex;

        int rc;
        if (buffer->IsA ("vtkDoubleArray"))
        {
          vtkDoubleArray *da = vtkArrayDownCast<vtkDoubleArray>(buffer);
          rc = ex_put_elem_var(this->fid, timestep + 1, varOutIndex + 1, id,
                                   numElts, da->GetPointer(start));
        }
        else /* (buffer->IsA ("vtkFloatArray")) */
        {
          vtkFloatArray *fa = vtkArrayDownCast<vtkFloatArray>(buffer);
          rc = ex_put_elem_var(this->fid, timestep + 1, varOutIndex + 1, id,
                                   numElts, fa->GetPointer(start));
        }

        if (rc < 0)
        {
          vtkErrorMacro(<< "vtkExodusIIWriter::WriteNextTimeStep ex_put_elem_var");
          return 0;
        }
      }
    }
  }
  return 1;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WritePointData (int timestep, vtkDataArray *buffer)
{
  if (this->NumPoints == 0)
  {
    return 1;
  }
  std::map<std::string, VariableInfo>::const_iterator varIter;
  for (varIter = this->NodeVariableMap.begin ();
       varIter != this->NodeVariableMap.end ();
       ++varIter)
  {
    const char *nameIn = varIter->first.c_str ();
    int numComp = varIter->second.NumComponents;
    for (int component = 0; component < numComp; component ++)
    {
      buffer->Initialize ();
      this->ExtractPointData(nameIn, component, buffer);
      int varOutIndex = varIter->second.ScalarOutOffset + component;
      int rc;
      if (buffer->IsA ("vtkDoubleArray"))
      {
        vtkDoubleArray *da = vtkArrayDownCast<vtkDoubleArray>(buffer);
        rc = ex_put_nodal_var(this->fid, timestep + 1, varOutIndex + 1,
                                  this->NumPoints, da->GetPointer (0));
      }
      else /* (buffer->IsA ("vtkFloatArray")) */
      {
        vtkFloatArray *fa = vtkArrayDownCast<vtkFloatArray>(buffer);
        rc = ex_put_nodal_var(this->fid, timestep + 1, varOutIndex + 1,
                                  this->NumPoints, fa->GetPointer (0));
      }

      if (rc < 0)
      {
        vtkErrorMacro(<< "vtkExodusIIWriter::WriteNextTimeStep ex_put_nodal_var");
        return 0;
      }
    }
  }
  return 1;
}

namespace
{
unsigned int GetLongestFieldDataName(vtkFieldData* fd)
{
  unsigned int maxName = 0;
  for (int i=0;i<fd->GetNumberOfArrays();i++)
  {
    unsigned int length = static_cast<unsigned int>(strlen(fd->GetArrayName(i)));
    if (length > maxName)
    {
      maxName = length;
    }
  }
  return maxName;
}

unsigned int GetLongestDataSetName(vtkDataSet* ds)
{
  unsigned int maxName = 32;
  unsigned int maxDataSetName = GetLongestFieldDataName(ds->GetPointData());
  if (maxDataSetName > maxName)
  {
    maxName = maxDataSetName;
  }
  maxDataSetName = GetLongestFieldDataName(ds->GetCellData());
  if (maxDataSetName > maxName)
  {
    maxName = maxDataSetName;
  }
  maxDataSetName = GetLongestFieldDataName(ds->GetFieldData());
  if (maxDataSetName > maxName)
  {
    maxName = maxDataSetName;
  }
  return maxName;
}
}

//----------------------------------------------------------------------------
unsigned int vtkExodusIIWriter::GetMaxNameLength()
{
  unsigned int maxName = 32;
  if (vtkMultiBlockDataSet* mb = vtkMultiBlockDataSet::SafeDownCast(this->OriginalInput))
  {
    vtkCompositeDataIterator* iter = mb->NewIterator();
    iter->SkipEmptyNodesOn();
    for (iter->InitTraversal();!iter->IsDoneWithTraversal();iter->GoToNextItem())
    {
      if (vtkDataSet* dataSet = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject()))
      {
        unsigned int maxDataSetName = GetLongestDataSetName(dataSet);
        if (maxDataSetName > maxName)
        {
          maxName = maxDataSetName;
        }
        if (vtkInformation* info = iter->GetCurrentMetaData())
        {
          if (const char* objectName = info->Get(vtkCompositeDataSet::NAME()))
          {
            maxDataSetName = static_cast<unsigned int>(strlen(objectName));
            if (maxDataSetName > maxName)
            {
              maxName = maxDataSetName;
            }
          }
        }
      }
    }
    iter->Delete();
  }
  else if (vtkDataSet* dataSet = vtkDataSet::SafeDownCast(this->OriginalInput))
  {
    maxName = GetLongestDataSetName(dataSet);
  }
  return maxName;
}

//----------------------------------------------------------------------------
int vtkExodusIIWriter::WriteNextTimeStep()
{
  int rc = 0;
  int ts = this->CurrentTimeIndex - this->FileTimeOffset;
  float tsv = 0.;
  if (this->GetInput()->GetInformation()->Has(vtkDataObject::DATA_TIME_STEP()) )
  {
    tsv = this->GetInput()->GetInformation()->Get(vtkDataObject::DATA_TIME_STEP());
  }

  vtkSmartPointer<vtkDataArray> buffer;
  if (this->PassDoubles)
  {
    double dtsv = (double)tsv;
    rc = ex_put_time(this->fid, ts + 1, &dtsv);
    if (rc < 0)
    {
      vtkErrorMacro(<< "vtkExodusIIWriter::WriteNextTimeStep time step values"
                    << " fid " << this->fid << " ts " << ts + 1 << " tsv " << tsv);
      return 0;
    }
    buffer = vtkSmartPointer<vtkDoubleArray>::New ();
  }
  else
  {
    rc = ex_put_time(this->fid, ts + 1, &tsv);
    if (rc < 0)
    {
      vtkErrorMacro(<< "vtkExodusIIWriter::WriteNextTimeStep time step values"
                    << " fid " << this->fid << " ts " << ts + 1 << " tsv " << tsv);
      return 0;
    }
    buffer = vtkSmartPointer<vtkFloatArray>::New ();
  }

  // Buffer is used to help these determine the type of the data to write out
  if (!this->WriteGlobalData (ts, buffer))
  {
    return 0;
  }
  if (!this->WriteCellData (ts, buffer))
  {
    return 0;
  }
  if (!this->WritePointData (ts, buffer))
  {
    return 0;
  }

  return 1;
}

void vtkExodusIIWriter::CheckBlockInfoMap()
{
  // no op for serial version
}

bool vtkExodusIIWriter::SameTypeOfCells (vtkIntArray* cellToBlockId,
                                         vtkUnstructuredGrid* input)
{
  if(cellToBlockId->GetNumberOfComponents () != 1 ||
      cellToBlockId->GetNumberOfTuples () != input->GetNumberOfCells())
  {
    return false;
  }
  std::map<int, int> blockIdToCellType;
  for (vtkIdType cellId = 0; cellId < cellToBlockId->GetNumberOfTuples ();
       ++cellId)
  {
    vtkIdType blockId = cellToBlockId->GetValue(cellId);
    std::map<int, int>::iterator it = blockIdToCellType.find(blockId);
    if (it == blockIdToCellType.end())
    {
      blockIdToCellType[blockId] = input->GetCellType(cellId);
    }
    else
    {
      if (it->second != input->GetCellType(cellId))
      {
        return false;
      }
    }
  }
  return true;
}

vtkIntArray* vtkExodusIIWriter::GetBlockIdArray (
  const char* name, vtkUnstructuredGrid* input)
{
  vtkDataArray *da = nullptr;
  vtkCellData *cd = input->GetCellData();
  if (name)
  {
    da = cd->GetArray(name);
  }
  if (!da)
  {
    name = "ObjectId";
    da = cd->GetArray(name);
    if (!da)
    {
      name = "ElementBlockIds";
      da = cd->GetArray(name);
    }
  }
  if (da)
  {
    vtkIntArray *ia = vtkArrayDownCast<vtkIntArray>(da);
    if (ia != nullptr && vtkExodusIIWriter::SameTypeOfCells (ia, input))
    {
      this->SetBlockIdArrayName(name);
      return ia;
    }
  }
  this->SetBlockIdArrayName(nullptr);
  if ((this->NumberOfProcesses > 1) &&
      // you don't have metadata but you have some tuples.
      cd->GetNumberOfTuples() > 0 &&
      // depending on what we're trying to write out we may not care about missing metadata
      this->IgnoreMetaDataWarning == 0)
  {
    // Parallel apps must have a global list of all block IDs, plus a
    // list of block IDs for each cell.
    vtkWarningMacro(<< "Attempting to proceed without metadata");
  }
  return nullptr;
}