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

  Program:   Visualization Toolkit
  Module:    vtkLSDynaReader.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.
----------------------------------------------------------------------------*/

// NOTE TO DEVELOPERS: ========================================================
//
// This is a really big reader.
// It uses several classes defined in Utilities/LSDyna:
// - LSDynaFamily:
//    A class to abstract away I/O from families of output files.
//    This performs the actual reads and writes plus any required byte swapping.
//    Also contains a subclass, LSDynaFamilyAdaptLevel, used to store
//    file+offset
//    information for each mesh adaptation's state info.
// - LSDynaMetaData:
//    A class to hold metadata about a particular file (such as time steps,
//    the start of state information for each time step, the number of
//    adaptive remeshes, and the large collection of constants that determine
//    the available attributes). It contains an LSDynaFamily instance.

// It also uses a helper vtk class
// - vtkLSDynaSummaryParser:
//    A class to parse XML summary files containing part names and their IDs.
//    This class is used by vtkLSDynaReader::ReadInputDeckXML().

// This class is preceded by some file-static constants and utility routines.

#include "vtkLSDynaReader.h"
#include "LSDynaFamily.h"
#include "LSDynaMetaData.h"
#include "vtkLSDynaPartCollection.h"
#include "vtkLSDynaSummaryParser.h"

#include "vtksys/FStream.hxx"
#include "vtksys/SystemTools.hxx"

#include <algorithm>
#include <cassert>
#include <map>
#include <string>
#include <vector>

#include "vtkCellType.h"
#include "vtkDataObject.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationDoubleVectorKey.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"
#include "vtkVectorOperators.h"

vtkStandardNewMacro(vtkLSDynaReader);

// Names of vtkDataArrays provided with grid:
#define LS_ARRAYNAME_DEATH "Death"
#define LS_ARRAYNAME_USERID "UserID"
#define LS_ARRAYNAME_SPECIES_BLNK "SpeciesXX"
#define LS_ARRAYNAME_SPECIES_FMT "Species%02d"
#define LS_ARRAYNAME_SPECIES_01 "Species01"
#define LS_ARRAYNAME_SPECIES_02 "Species02"
#define LS_ARRAYNAME_SPECIES_03 "Species03"
#define LS_ARRAYNAME_SPECIES_04 "Species04"
#define LS_ARRAYNAME_SPECIES_05 "Species05"
#define LS_ARRAYNAME_SPECIES_06 "Species06"
#define LS_ARRAYNAME_SPECIES_07 "Species07"
#define LS_ARRAYNAME_SPECIES_08 "Species08"
#define LS_ARRAYNAME_SPECIES_09 "Species09"
#define LS_ARRAYNAME_SPECIES_10 "Species10"
#define LS_ARRAYNAME_TEMPERATURE "Temperature"
#define LS_ARRAYNAME_DEFLECTION "Deflected Coordinates"
#define LS_ARRAYNAME_VELOCITY "Velocity"
#define LS_ARRAYNAME_ACCELERATION "Acceleration"
#define LS_ARRAYNAME_PRESSURE "Pressure"
#define LS_ARRAYNAME_VORTICITY "Vorticity"
#define LS_ARRAYNAME_VORTICITY_X "Vorticity_X"
#define LS_ARRAYNAME_VORTICITY_Y "Vorticity_Y"
#define LS_ARRAYNAME_VORTICITY_Z "Vorticity_Z"
#define LS_ARRAYNAME_RESULTANTVORTICITY "ResVorticity"
#define LS_ARRAYNAME_ENSTROPHY "Enstrophy"
#define LS_ARRAYNAME_HELICITY "Helicity"
#define LS_ARRAYNAME_STREAMFUNCTION "StreamFunc"
#define LS_ARRAYNAME_ENTHALPY "Enthalpy"
#define LS_ARRAYNAME_DENSITY "Density"
#define LS_ARRAYNAME_TURBULENTKE "TurbulentKE"
#define LS_ARRAYNAME_DISSIPATION "Dissipation"
#define LS_ARRAYNAME_EDDYVISCOSITY "EddyVisc"
#define LS_ARRAYNAME_RADIUSOFINFLUENCE "InfluenceRadius"
#define LS_ARRAYNAME_NUMNEIGHBORS "NumberOfNeighbors"
#define LS_ARRAYNAME_SEGMENTID "SegmentID"
#define LS_ARRAYNAME_STRAIN "Strain"
#define LS_ARRAYNAME_STRESS "Stress"
#define LS_ARRAYNAME_EPSTRAIN "EffPlastStrn"
#define LS_ARRAYNAME_INTEGRATIONPOINT "IntPtData"
#define LS_ARRAYNAME_RESULTANTS "Resultants"
#define LS_ARRAYNAME_ELEMENTMISC "ElementMisc"
#define LS_ARRAYNAME_INTERNALENERGY "InternalEnergy"
#define LS_ARRAYNAME_AXIALFORCE "AxialForce"
#define LS_ARRAYNAME_SHEARRESULTANT "ShearResultant"
#define LS_ARRAYNAME_BENDINGRESULTANT "BendingResultant"
#define LS_ARRAYNAME_TORSIONRESULTANT "TorsionResultant"
#define LS_ARRAYNAME_NORMALRESULTANT "NormalResultant"
#define LS_ARRAYNAME_AXIALSTRAIN "AxialStrain"
#define LS_ARRAYNAME_AXIALSTRESS "AxialStress"
#define LS_ARRAYNAME_SHEARSTRAIN "ShearStrain"
#define LS_ARRAYNAME_SHEARSTRESS "ShearStress"
#define LS_ARRAYNAME_PLASTICSTRAIN "PlasticStrain"
#define LS_ARRAYNAME_THICKNESS "Thickness"
#define LS_ARRAYNAME_MASS "Mass"
#define LS_ARRAYNAME_VOLUME_FRACTION_FMT "VolumeFraction%02d"
#define LS_ARRAYNAME_DOMINANT_GROUP "DominantGroup"
#define LS_ARRAYNAME_SPECIES_MASS_FMT "SpeciesMass%02d"
#define LS_ARRAYNAME_MATERIAL "Material"

// Possible material options
#define LS_MDLOPT_NONE 0
#define LS_MDLOPT_POINT 1
#define LS_MDLOPT_CELL 2

#ifdef VTK_LSDYNA_DBG_MULTIBLOCK
static void vtkDebugMultiBlockStructure(vtkIndent indent, vtkMultiGroupDataSet* mbds);
#endif // VTK_LSDYNA_DBG_MULTIBLOCK

namespace
{
const char* vtkLSDynaCellTypes[] = { "Point", "Beam", "Shell", "Thick Shell", "Solid", "Rigid Body",
  "Road Surface" };

// Read in lines until one that's
// - not empty, and
// - not a comment
// is encountered. Return with that text stored in \a line.
// If an error or EOF is hit, return 0. Otherwise, return 1.
int vtkLSNextSignificantLine(istream& deck, std::string& line)
{
  while (deck.good())
  {
    std::getline(deck, line, '\n');
    if (!line.empty() && line[0] != '$')
    {
      return 1;
    }
  }
  return 0;
}

void vtkLSTrimWhitespace(std::string& line)
{
  std::string::size_type llen = line.length();
  while (llen &&
    (line[llen - 1] == ' ' || line[llen - 1] == '\t' || line[llen - 1] == '\r' ||
      line[llen - 1] == '\n'))
  {
    --llen;
  }

  std::string::size_type nameStart = 0;
  while (nameStart < llen && (line[nameStart] == ' ' || line[nameStart] == '\t'))
  {
    ++nameStart;
  }

  line = line.substr(nameStart, llen - nameStart);
}

void vtkLSDowncaseFirstWord(std::string& downcased, const std::string& line)
{
  std::string::size_type i;
  std::string::value_type chr;
  int leadingSpace = 0;
  downcased = "";
  for (i = 0; i < line.length(); ++i)
  {
    chr = tolower(line[i]);
    if (chr == ' ' || chr == '\t')
    {
      if (leadingSpace)
      { // We've trimmed leading whitespace already, so we're done with the word.
        return;
      }
    }
    else
    {
      leadingSpace = 1;
      if (chr == ',')
      { // We're at a separator (other than whitespace). No need to continue.
        return;
      }
    }
    downcased += chr;
  }
}

void vtkLSSplitString(std::string& input, std::vector<std::string>& splits, const char* separators)
{
  std::string::size_type posBeg = 0;
  std::string::size_type posEnd;
  do
  {
    posEnd = input.find_first_of(separators, posBeg);
    if (posEnd > posBeg)
    {
      // don't include empty entries in splits.
      // NOTE: This means ",comp,1, ,3" with separators ", " yields "comp","1","3", not
      // "","comp","1","","","3".
      splits.push_back(input.substr(posBeg, posEnd - posBeg));
    }
    posBeg = input.find_first_not_of(separators, posEnd);
  } while (posBeg != std::string::npos);
}

template <int hostBitSize, int fileBitSize, int cellLength>
struct Converter
{
  // general use case that the host
  // bit size and file bit size are the same
  vtkIdType* convert(vtkIdType* buff, const vtkIdType&) { return buff; }
};

template <int cellLength>
struct Converter<8, 4, cellLength>
{
  // specilization of 64bit machine and 32bit file
  // so we have to copy each item individually
  vtkIdType* convert(int* buff, const vtkIdType& size)
  {
    for (vtkIdType i = 0; i < size; ++i)
    {
      this->Conn[i] = static_cast<vtkIdType>(buff[i]);
    }
    return Conn;
  }
  vtkIdType Conn[cellLength];
};

template <int cellLength>
struct Converter<4, 8, cellLength>
{
  // specilization for reading 64 bit files on a 32 bit machine
  // which means reading the bottom half of the long long
  vtkIdType* convert(int* buff, const vtkIdType& size)
  {
    vtkIdType idx = 0;
    for (vtkIdType i = 0; i < size; i += 2, ++idx)
    {
      this->Conn[idx] = static_cast<vtkIdType>(buff[i]);
    }
    return Conn;
  }
  vtkIdType Conn[cellLength];
};

template <int type, int wordSize, int cellLength>
struct FillBlock
{
  Converter<sizeof(vtkIdType), wordSize, cellLength> BC;

  template <typename T>
  FillBlock(T* buff, vtkLSDynaPartCollection* parts, LSDynaMetaData* p,
    const vtkIdType& numWordsPerCell, const int& cellType)
  {
    // determine the relationship between the file bit size and
    // the host machine bit size. This allows us to read 64 bit files on a
    // 32 bit machine. The Converter allows us to easily convert 32bit
    // arrays to 64bit arrays
    const int numWordsPerIdType(p->Fam.GetWordSize() / sizeof(T));
    const vtkIdType numFileWordsPerCell(numWordsPerCell * numWordsPerIdType);
    const vtkIdType offsetToMatId(numWordsPerIdType * (numWordsPerCell - 1));
    vtkIdType* conn;

    vtkIdType nc = 0, j = 0, matlId = 0;
    vtkIdType numCellsToSkip = 0, numCellsToSkipEnd = 0, chunkSize = 0;

    // get from the part the read information for this lsdyna block type
    parts->GetPartReadInfo(type, nc, numCellsToSkip, numCellsToSkipEnd);

    p->Fam.SkipWords(numFileWordsPerCell * numCellsToSkip); // skip to the right start id

    // buffer the amount in small chunks so we don't create a massive buffer
    vtkIdType numChunks = p->Fam.InitPartialChunkBuffering(nc, numWordsPerCell);
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      chunkSize = p->Fam.GetNextChunk(LSDynaFamily::Int);
      buff = p->Fam.GetBufferAs<T>();
      for (j = 0; j < chunkSize; j += numWordsPerCell)
      {
        conn = BC.convert(buff, offsetToMatId);
        buff += offsetToMatId;
        matlId = static_cast<vtkIdType>(*buff);
        buff += numWordsPerIdType;
        parts->InsertCell(type, matlId, cellType, cellLength, conn);
      }
    }
    p->Fam.SkipWords(numFileWordsPerCell * numCellsToSkipEnd);
  }
};

template <int wordSize, int cellLength>
struct FillBlock<LSDynaMetaData::SOLID, wordSize, cellLength>
{
  Converter<sizeof(vtkIdType), wordSize, cellLength> BC;

  template <typename T>
  FillBlock(T* buff, vtkLSDynaPartCollection* parts, LSDynaMetaData* p,
    const vtkIdType& numWordsPerCell, const int& vtkNotUsed(cellType))
  {
    // determine the relationship between the file bit size and
    // the host machine bit size. This allows us to read 64 bit files on a
    // 32 bit machine. The Converter allows us to easily convert 32bit
    // arrays to 64bit arrays
    const int numWordsPerIdType(p->Fam.GetWordSize() / sizeof(T));
    const vtkIdType numFileWordsPerCell(numWordsPerCell * numWordsPerIdType);
    const vtkIdType offsetToMatId(numWordsPerIdType * cellLength);
    vtkIdType* conn;

    // This is a read solids template specialization since it has a special use
    // case for cell length based on the connectivity mapping
    vtkIdType nc = 0, j = 0, matlId = 0;
    vtkIdType numCellsToSkip = 0, numCellsToSkipEnd = 0, chunkSize = 0;

    // get from the part the read information for this lsdyna block type
    parts->GetPartReadInfo(LSDynaMetaData::SOLID, nc, numCellsToSkip, numCellsToSkipEnd);

    p->Fam.SkipWords(numFileWordsPerCell * numCellsToSkip); // skip to the right start id

    // buffer the amount in small chunks so we don't create a massive buffer
    vtkIdType numChunks = p->Fam.InitPartialChunkBuffering(nc, numWordsPerCell);
    vtkIdType npts = 0;
    int ctype = 0;
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      chunkSize = p->Fam.GetNextChunk(LSDynaFamily::Int);
      buff = p->Fam.GetBufferAs<T>();
      for (j = 0; j < chunkSize; j += numWordsPerCell)
      {
        conn = BC.convert(buff, offsetToMatId);
        buff += offsetToMatId;
        matlId = static_cast<vtkIdType>(*buff);
        buff += numWordsPerIdType;

        // Detect repeated connectivity entries to determine element type
        if (conn[3] == conn[7])
        {
          ctype = VTK_TETRA;
          npts = 4;
        }
        else if (conn[4] == conn[7])
        {
          ctype = VTK_PYRAMID;
          npts = 5;
        }
        else if (conn[5] == conn[7])
        {
          ctype = VTK_WEDGE;
          npts = 6;
        }
        else
        {
          ctype = VTK_HEXAHEDRON;
          npts = 8;
        }

        // push this cell back into the unstructured grid for this part(if the part is active)
        parts->InsertCell(LSDynaMetaData::SOLID, matlId, ctype, npts, conn);
      }
    }
    p->Fam.SkipWords(numFileWordsPerCell * numCellsToSkipEnd);
  }
};

template <int wordSize, int cellLength>
struct FillBlock<LSDynaMetaData::SHELL, wordSize, cellLength>
{
  Converter<sizeof(vtkIdType), wordSize, cellLength> BC;

  template <typename T>
  FillBlock(T* buff, vtkLSDynaPartCollection* parts, LSDynaMetaData* p,
    const vtkIdType& numWordsPerCell, const int& cellType)
  {
    // determine the relationship between the file bit size and
    // the host machine bit size. This allows us to read 64 bit files on a
    // 32 bit machine. The Converter allows us to easily convert 32bit
    // arrays to 64bit arrays
    const int numWordsPerIdType(p->Fam.GetWordSize() / sizeof(T));
    const vtkIdType numFileWordsPerCell(numWordsPerCell * numWordsPerIdType);
    const vtkIdType offsetToMatId(numWordsPerIdType * cellLength);
    vtkIdType* conn;

    // This is a read RIGID_BODY and SHELL template specialization since it
    // has a weird weaving of cell types
    bool haveRigidMaterials = (p->Dict["MATTYP"] != 0) && !p->RigidMaterials.empty();

    vtkIdType nc = 0, j = 0, matlId = 0;
    vtkIdType numCellsToSkip = 0, numCellsToSkipEnd = 0, chunkSize = 0;

    // get from the part the read information for this lsdyna block type
    parts->GetPartReadInfo(LSDynaMetaData::SHELL, nc, numCellsToSkip, numCellsToSkipEnd);

    p->Fam.SkipWords(numFileWordsPerCell * numCellsToSkip); // skip to the right start id

    // buffer the amount in small chunks so we don't create a massive buffer
    vtkIdType numChunks = p->Fam.InitPartialChunkBuffering(nc, numWordsPerCell);
    int pType = 0;
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      chunkSize = p->Fam.GetNextChunk(LSDynaFamily::Int);
      buff = p->Fam.GetBufferAs<T>();
      for (j = 0; j < chunkSize; j += numWordsPerCell)
      {
        conn = BC.convert(buff, offsetToMatId);
        buff += offsetToMatId;
        matlId = static_cast<vtkIdType>(*buff);
        buff += numWordsPerIdType;

        if (haveRigidMaterials && p->RigidMaterials.find(matlId) == p->RigidMaterials.end())
        {
          pType = LSDynaMetaData::RIGID_BODY;
        }
        else
        {
          pType = LSDynaMetaData::SHELL;
        }
        parts->InsertCell(pType, matlId, cellType, cellLength, conn);
      }
    }
    p->Fam.SkipWords(numFileWordsPerCell * numCellsToSkipEnd);
  }
};

template <int wordSize, int cellLength>
struct FillBlock<LSDynaMetaData::ROAD_SURFACE, wordSize, cellLength>
{
  template <typename T>
  FillBlock(
    T*, vtkLSDynaPartCollection* parts, LSDynaMetaData* p, const vtkIdType&, const int& cellType)
  {
    // This is a ROAD_SURFACE specialization
    // has a weird weaving of cell types
    vtkIdType nc = 0, segId = 0, segSz = 0;
    vtkIdType numCellsToSkip = 0, numCellsToSkipEnd = 0;

    // get from the part the read information for this lsdyna block type
    parts->GetPartReadInfo(LSDynaMetaData::SHELL, nc, numCellsToSkip, numCellsToSkipEnd);

    // the road surface format is horrible for parallel reading.
    // we don't know the number of cells in each road surface.
    // only the total number of cells. So we have to do some fun stuff to correctly skip
    // this is unoptimized since I don't have any road surface data
    vtkIdType currentCell = 0;
    vtkIdType conn[4];
    for (vtkIdType i = 0; i < p->Dict["NSURF"]; ++i)
    {
      p->Fam.BufferChunk(LSDynaFamily::Int, 2);
      segId = p->Fam.GetNextWordAsInt();
      segSz = p->Fam.GetNextWordAsInt();
      p->Fam.BufferChunk(LSDynaFamily::Int, 4 * segSz);
      for (vtkIdType t = 0; t < segSz; ++t, ++currentCell)
      {
        if (currentCell >= numCellsToSkip)
        {
          for (int j = 0; j < 4; ++j)
          {
            conn[j] = p->Fam.GetNextWordAsInt() - 1;
          }
          parts->InsertCell(LSDynaMetaData::ROAD_SURFACE, segId, cellType, 4, conn);
        }
      }
    }
  }
};

}
// =================================================== Start of public interface
vtkLSDynaReader::vtkLSDynaReader()
{
  this->P = new LSDynaMetaData;
  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(1);
  this->TimeStepRange[0] = 0;
  this->TimeStepRange[1] = 0;
  this->DeformedMesh = 1;
  this->RemoveDeletedCells = 1;
  this->DeletedCellsAsGhostArray = 0;
  this->InputDeck = nullptr;
  this->Parts = nullptr;
}

vtkLSDynaReader::~vtkLSDynaReader()
{
  this->ResetPartsCache();
  this->SetInputDeck(nullptr);
  delete this->P;
  this->P = nullptr;
}

void vtkLSDynaReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "Title: \"" << this->GetTitle() << "\"" << endl;
  os << indent << "InputDeck: " << (this->InputDeck ? this->InputDeck : "(null)") << endl;
  os << indent << "DeformedMesh: " << (this->DeformedMesh ? "On" : "Off") << endl;
  os << indent << "RemoveDeletedCells: " << (this->RemoveDeletedCells ? "On" : "Off") << endl;
  os << indent << "TimeStepRange: " << this->TimeStepRange[0] << ", " << this->TimeStepRange[1]
     << endl;

  if (this->P)
  {
    os << indent << "PrivateData: " << this->P << endl;
  }
  else
  {
    os << indent << "PrivateData: (none)" << endl;
  }
  os << indent
     << "Show Deleted Cells as Ghost Cells: " << (this->DeletedCellsAsGhostArray ? "On" : "Off")
     << endl;

  os << indent << "Dimensionality: " << this->GetDimensionality() << endl;
  os << indent << "Nodes: " << this->GetNumberOfNodes() << endl;
  os << indent << "Cells: " << this->GetNumberOfCells() << endl;

  os << indent << "PointArrays: ";
  for (int i = 0; i < this->GetNumberOfPointArrays(); ++i)
  {
    os << this->GetPointArrayName(i) << " ";
  }
  os << endl;

  os << "CellArrays: " << endl;
  for (int ct = 0; ct < LSDynaMetaData::NUM_CELL_TYPES; ++ct)
  {
    os << vtkLSDynaCellTypes[ct] << ":" << endl;
    for (int i = 0; i < this->GetNumberOfCellArrays(ct); ++i)
    {
      os << this->GetCellArrayName(ct, i) << " ";
    }
    os << endl;
  }
  os << endl;

  os << indent << "Time Steps: " << this->GetNumberOfTimeSteps() << endl;
  for (int j = 0; j < this->GetNumberOfTimeSteps(); ++j)
  {
    os.precision(5);
    os.width(12);
    os << this->GetTimeValue(j);
    if ((j + 1) % 8 == 0 && j != this->GetNumberOfTimeSteps() - 1)
    {
      os << endl << indent;
    }
    else
    {
      os << " ";
    }
  }
  os << endl;
}

void vtkLSDynaReader::Dump(ostream& os)
{
  vtkIndent indent;
  os << indent << "Title: \"" << this->GetTitle() << "\"" << endl
     << indent << "DeformedMesh: " << (this->DeformedMesh ? "On" : "Off") << endl
     << indent << "RemoveDeletedCells: " << (this->RemoveDeletedCells ? "On" : "Off") << endl
     << indent << "TimeStepRange: " << this->TimeStepRange[0] << ", " << this->TimeStepRange[1]
     << endl
     << indent << "PrivateData: " << this->P << endl
     << indent << "Dimensionality: " << this->GetDimensionality() << endl
     << indent << "Nodes: " << this->GetNumberOfNodes() << endl
     << indent << "Cells: " << this->GetNumberOfCells() << endl
     << indent << "PointArrays:    ";
  for (int i = 0; i < this->GetNumberOfPointArrays(); ++i)
  {
    os << this->GetPointArrayName(i) << " ";
  }
  os << endl << "CellArrays:" << endl;
  for (int ct = 0; ct < LSDynaMetaData::NUM_CELL_TYPES; ++ct)
  {
    os << vtkLSDynaCellTypes[ct] << ":" << endl;
    for (int i = 0; i < this->GetNumberOfCellArrays(ct); ++i)
    {
      os << this->GetCellArrayName(ct, i) << " ";
    }
    os << endl;
  }
  os << endl;

  os << indent << "Time Steps:       " << this->GetNumberOfTimeSteps() << endl;
  for (int j = 0; j < this->GetNumberOfTimeSteps(); ++j)
  {
    os.precision(5);
    os.width(12);
    os << this->GetTimeValue(j);
    if ((j + 1) % 8 == 0 && j != this->GetNumberOfTimeSteps() - 1)
    {
      os << endl << indent;
    }
    else
    {
      os << " ";
    }
  }
  os << endl;
}

void vtkLSDynaReader::DebugDump()
{
  this->Dump(cout);
}

int vtkLSDynaReader::CanReadFile(const char* fname)
{
  if (!fname)
    return 0;

  std::string dbDir = vtksys::SystemTools::GetFilenamePath(fname);
  std::string dbName = vtksys::SystemTools::GetFilenameName(fname);
  std::string dbExt;
  std::string::size_type dot;
  LSDynaMetaData* p = new LSDynaMetaData;
  int result = 0;

  // GetFilenameExtension doesn't look for the rightmost "." ... do it ourselves.
  dot = dbName.rfind('.');
  if (dot != std::string::npos)
  {
    dbExt = dbName.substr(dot);
  }
  else
  {
    dbExt = "";
  }

  p->Fam.SetDatabaseDirectory(dbDir);

  if (dbExt == ".k" || dbExt == ".lsdyna")
  {
    p->Fam.SetDatabaseBaseName("/d3plot");
  }
  else
  {
    vtksys::SystemTools::Stat_t st;
    if (vtksys::SystemTools::Stat(fname, &st) == 0)
    {
      dbName.insert(0, "/");
      p->Fam.SetDatabaseBaseName(dbName);
    }
    else
    {
      p->Fam.SetDatabaseBaseName("/d3plot");
    }
  }
  // If the time step is set before RequestInformation is called, we must
  // read the header information immediately in order to determine whether
  // the timestep that's been passed is valid. If it's not, we ignore it.
  if (!p->FileIsValid)
  {
    if (p->Fam.GetDatabaseDirectory().empty())
    {
      result = -1;
    }
    else
    {
      if (p->Fam.GetDatabaseBaseName().empty())
      {
        p->Fam.SetDatabaseBaseName("/d3plot"); // not a bad assumption.
      }
      p->Fam.ScanDatabaseDirectory();
      if (p->Fam.GetNumberOfFiles() < 1)
      {
        result = -1;
      }
      else
      {
        if (p->Fam.DetermineStorageModel() != 0)
          result = 0;
        else
          result = 1;
      }
    }
  }
  delete p;

  return result > 0; // -1 and 0 are both problems, 1 indicates success.
}

void vtkLSDynaReader::SetDatabaseDirectory(const std::string& f)
{
  this->SetDatabaseDirectory(f.c_str());
}

void vtkLSDynaReader::SetDatabaseDirectory(const char* f)
{
  vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting DatabaseDirectory to " << f);
  if (!f)
  {
    if (!this->P->Fam.GetDatabaseDirectory().empty())
    { // no string => no database directory
      this->P->Reset();
      this->SetInputDeck(nullptr);
      this->ResetPartsCache();
      this->Modified();
    }
    return;
  }
  if (strcmp(this->P->Fam.GetDatabaseDirectory().c_str(), f) != 0)
  {
    this->P->Reset();
    this->SetInputDeck(nullptr);
    this->P->Fam.SetDatabaseDirectory(std::string(f));
    this->ResetPartsCache();
    this->Modified();
  }
}

#if defined(VTK_LEGACY_REMOVE)
std::string vtkLSDynaReader::GetDatabaseDirectory()
{
  return this->P->Fam.GetDatabaseDirectory();
}
#else
const char* vtkLSDynaReader::GetDatabaseDirectory()
{
  static VTK_THREAD_LOCAL std::string surrogate;
  surrogate = this->P->Fam.GetDatabaseDirectory();
  return surrogate.c_str();
}
#endif

int vtkLSDynaReader::IsDatabaseValid()
{
  return this->P->FileIsValid;
}

void vtkLSDynaReader::SetFileName(const std::string& f)
{
  this->SetFileName(f.c_str());
}

void vtkLSDynaReader::SetFileName(const char* f)
{
  std::string dbDir = vtksys::SystemTools::GetFilenamePath(f);
  std::string dbName = vtksys::SystemTools::GetFilenameName(f);
  std::string dbExt;
  std::string::size_type dot;

  // GetFilenameExtension doesn't look for the rightmost "." ... do it ourselves.
  dot = dbName.rfind('.');
  if (dot != std::string::npos)
  {
    dbExt = dbName.substr(dot);
  }
  else
  {
    dbExt = "";
  }

  this->SetDatabaseDirectory(dbDir.c_str());

  if (dbExt == ".k" || dbExt == ".lsdyna")
  {
    this->SetInputDeck(f);
    this->P->Fam.SetDatabaseBaseName("/d3plot");
  }
  else
  {
    vtksys::SystemTools::Stat_t st;
    if (vtksys::SystemTools::Stat(f, &st) == 0)
    {
      dbName.insert(0, "/");
      this->P->Fam.SetDatabaseBaseName(dbName);
    }
    else
    {
      this->P->Fam.SetDatabaseBaseName("/d3plot");
    }
  }
}

#if defined(VTK_LEGACY_REMOVE)
std::string vtkLSDynaReader::GetFileName()
{
  std::string filename = this->P->Fam.GetDatabaseDirectory() + "/d3plot";
  return filename;
}
#else
const char* vtkLSDynaReader::GetFileName()
{
  static VTK_THREAD_LOCAL std::string surrogate;
  surrogate = this->P->Fam.GetDatabaseDirectory() + "/d3plot";
  return surrogate.c_str();
}
#endif

char* vtkLSDynaReader::GetTitle()
{
  return this->P->Title;
}

int vtkLSDynaReader::GetDimensionality()
{
  return this->P->Dimensionality;
}

void vtkLSDynaReader::SetTimeStep(vtkIdType t)
{
  LSDynaMetaData* p = this->P;
  if (p->CurrentState == t)
  {
    return;
  }

  // If the time step is set before RequestInformation is called, we must
  // read the header information immediately in order to determine whether
  // the timestep that's been passed is valid. If it's not, we ignore it.
  if (!p->FileIsValid)
  {
    if (p->Fam.GetDatabaseDirectory().empty())
    {
      vtkErrorMacro("You haven't set the LS-Dyna database directory!");
      return;
    }

    p->Fam.SetDatabaseBaseName("/d3plot"); // force this for now.
    p->Fam.ScanDatabaseDirectory();
    if (p->Fam.GetNumberOfFiles() < 1)
    {
      p->FileIsValid = 0;
      return;
    }
    p->Fam.DetermineStorageModel();
    p->MaxFileLength = p->FileSizeFactor * 512 * 512 * p->Fam.GetWordSize();
    p->FileIsValid = 1;

    // OK, now we have a list of files. Next, determine the length of the
    // state vector (#bytes of data stored per time step):
    this->ReadHeaderInformation(0);

    // Finally, we can loop through and determine where all the state
    // vectors start for each time step.
    this->ScanDatabaseTimeSteps();
  }

  // Now, make sure we update the dictionary to contain information
  // relevant to the adaptation level that matches the requested timestep.
  if (t >= 0 && t < (int)p->TimeValues.size())
  {
    if (p->Fam.GetCurrentAdaptLevel() != p->Fam.TimeAdaptLevel(t))
    {
      if (this->ReadHeaderInformation(p->Fam.TimeAdaptLevel(t)) == 0)
      {
        // unable to read the header information for the adaptation level corresponding
        // to the requested time step
        return;
      }
    }
  }

  p->CurrentState = t;
  this->Modified();
}

vtkIdType vtkLSDynaReader::GetTimeStep()
{
  return this->P->CurrentState;
}

vtkIdType vtkLSDynaReader::GetNumberOfTimeSteps()
{
  return (vtkIdType)this->P->TimeValues.size();
}

double vtkLSDynaReader::GetTimeValue(vtkIdType s)
{
  if (s < 0 || s >= (vtkIdType)this->P->TimeValues.size())
  {
    return -1.0;
  }

  return this->P->TimeValues[s];
}

vtkIdType vtkLSDynaReader::GetNumberOfNodes()
{
  return this->P->NumberOfNodes;
}

vtkIdType vtkLSDynaReader::GetNumberOfCells()
{
  vtkIdType tmp = 0;
  for (int c = 0; c < LSDynaMetaData::NUM_CELL_TYPES; ++c)
  {
    tmp += this->P->NumberOfCells[c];
  }
  return tmp;
}

vtkIdType vtkLSDynaReader::GetNumberOfSolidCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::SOLID];
}

vtkIdType vtkLSDynaReader::GetNumberOfThickShellCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::THICK_SHELL];
}

vtkIdType vtkLSDynaReader::GetNumberOfShellCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::SHELL];
}

vtkIdType vtkLSDynaReader::GetNumberOfRigidBodyCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::RIGID_BODY];
}

vtkIdType vtkLSDynaReader::GetNumberOfRoadSurfaceCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::ROAD_SURFACE];
}

vtkIdType vtkLSDynaReader::GetNumberOfBeamCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::BEAM];
}

vtkIdType vtkLSDynaReader::GetNumberOfParticleCells()
{
  return this->P->NumberOfCells[LSDynaMetaData::PARTICLE];
}

vtkIdType vtkLSDynaReader::GetNumberOfContinuumCells()
{
  vtkIdType tmp = 0;
  for (int c = LSDynaMetaData::PARTICLE + 1; c < LSDynaMetaData::NUM_CELL_TYPES; ++c)
  {
    tmp += this->P->NumberOfCells[c];
  }
  return tmp;
}

// =================================== Point array queries
int vtkLSDynaReader::GetNumberOfPointArrays()
{
  return (int)this->P->PointArrayNames.size();
}

const char* vtkLSDynaReader::GetPointArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->PointArrayNames.size())
    return nullptr;

  return this->P->PointArrayNames[a].c_str();
}

int vtkLSDynaReader::GetPointArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->PointArrayStatus.size())
    return 0;

  return this->P->PointArrayStatus[a];
}

void vtkLSDynaReader::SetPointArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->PointArrayStatus.size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->PointArrayStatus[a])
    return;

  this->P->PointArrayStatus[a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

int vtkLSDynaReader::GetNumberOfComponentsInPointArray(int a)
{
  if (a < 0 || a >= (int)this->P->PointArrayStatus.size())
    return 0;

  return this->P->PointArrayComponents[a];
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfCellArrays(int ct)
{
  return (int)this->P->CellArrayNames[ct].size();
}

const char* vtkLSDynaReader::GetCellArrayName(int ct, int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[ct].size())
    return nullptr;

  return this->P->CellArrayNames[ct][a].c_str();
}

int vtkLSDynaReader::GetCellArrayStatus(int ct, int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[ct].size())
    return 0;

  return this->P->CellArrayStatus[ct][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInCellArray(int ct, int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[ct].size())
    return 0;

  return this->P->CellArrayComponents[ct][a];
}

void vtkLSDynaReader::SetCellArrayStatus(int ct, int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[ct].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[ct][a])
    return;

  this->P->CellArrayStatus[ct][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfSolidArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::SOLID].size();
}

const char* vtkLSDynaReader::GetSolidArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::SOLID].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::SOLID][a].c_str();
}

int vtkLSDynaReader::GetSolidArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::SOLID].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::SOLID][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInSolidArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::SOLID].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::SOLID][a];
}

void vtkLSDynaReader::SetSolidArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::SOLID].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::SOLID][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::SOLID][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfThickShellArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::THICK_SHELL].size();
}

const char* vtkLSDynaReader::GetThickShellArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::THICK_SHELL].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::THICK_SHELL][a].c_str();
}

int vtkLSDynaReader::GetThickShellArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::THICK_SHELL].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::THICK_SHELL][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInThickShellArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::THICK_SHELL].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::THICK_SHELL][a];
}

void vtkLSDynaReader::SetThickShellArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::THICK_SHELL].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::THICK_SHELL][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::THICK_SHELL][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfShellArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::SHELL].size();
}

const char* vtkLSDynaReader::GetShellArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::SHELL].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::SHELL][a].c_str();
}

int vtkLSDynaReader::GetShellArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::SHELL].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::SHELL][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInShellArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::SHELL].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::SHELL][a];
}

void vtkLSDynaReader::SetShellArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::SHELL].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::SHELL][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::SHELL][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfRigidBodyArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::RIGID_BODY].size();
}

const char* vtkLSDynaReader::GetRigidBodyArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::RIGID_BODY].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::RIGID_BODY][a].c_str();
}

int vtkLSDynaReader::GetRigidBodyArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::RIGID_BODY].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::RIGID_BODY][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInRigidBodyArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::RIGID_BODY].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::RIGID_BODY][a];
}

void vtkLSDynaReader::SetRigidBodyArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::RIGID_BODY].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::RIGID_BODY][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::RIGID_BODY][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfRoadSurfaceArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::ROAD_SURFACE].size();
}

const char* vtkLSDynaReader::GetRoadSurfaceArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::ROAD_SURFACE].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::ROAD_SURFACE][a].c_str();
}

int vtkLSDynaReader::GetRoadSurfaceArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::ROAD_SURFACE].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::ROAD_SURFACE][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInRoadSurfaceArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::ROAD_SURFACE].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::ROAD_SURFACE][a];
}

void vtkLSDynaReader::SetRoadSurfaceArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::ROAD_SURFACE].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::ROAD_SURFACE][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::ROAD_SURFACE][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfBeamArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::BEAM].size();
}

const char* vtkLSDynaReader::GetBeamArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::BEAM].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::BEAM][a].c_str();
}

int vtkLSDynaReader::GetBeamArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::BEAM].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::BEAM][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInBeamArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::BEAM].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::BEAM][a];
}

void vtkLSDynaReader::SetBeamArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::BEAM].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::BEAM][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::BEAM][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Cell array queries
int vtkLSDynaReader::GetNumberOfParticleArrays()
{
  return (int)this->P->CellArrayNames[LSDynaMetaData::PARTICLE].size();
}

const char* vtkLSDynaReader::GetParticleArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayNames[LSDynaMetaData::PARTICLE].size())
    return nullptr;

  return this->P->CellArrayNames[LSDynaMetaData::PARTICLE][a].c_str();
}

int vtkLSDynaReader::GetParticleArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::PARTICLE].size())
    return 0;

  return this->P->CellArrayStatus[LSDynaMetaData::PARTICLE][a];
}

int vtkLSDynaReader::GetNumberOfComponentsInParticleArray(int a)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::PARTICLE].size())
    return 0;

  return this->P->CellArrayComponents[LSDynaMetaData::PARTICLE][a];
}

void vtkLSDynaReader::SetParticleArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->CellArrayStatus[LSDynaMetaData::PARTICLE].size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->CellArrayStatus[LSDynaMetaData::PARTICLE][a])
    return;

  this->P->CellArrayStatus[LSDynaMetaData::PARTICLE][a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// =================================== Part queries
int vtkLSDynaReader::GetNumberOfPartArrays()
{
  return (int)this->P->PartNames.size();
}

const char* vtkLSDynaReader::GetPartArrayName(int a)
{
  if (a < 0 || a >= (int)this->P->PartNames.size())
    return nullptr;

  return this->P->PartNames[a].c_str();
}

int vtkLSDynaReader::GetPartArrayStatus(int a)
{
  if (a < 0 || a >= (int)this->P->PartStatus.size())
    return 0;

  return this->P->PartStatus[a];
}

void vtkLSDynaReader::SetPartArrayStatus(int a, int stat)
{
  if (a < 0 || a >= (int)this->P->PartStatus.size())
  {
    vtkWarningMacro("Cannot set status of non-existent point array " << a);
    return;
  }

  if (stat == this->P->PartStatus[a])
    return;

  this->P->PartStatus[a] = stat;
  this->ResetPartsCache();
  this->Modified();
}

// ===================================
void vtkLSDynaReader::ResetPartsCache()
{
  if (this->Parts)
  {
    this->Parts->Delete();
    this->Parts = nullptr;
  }
}

// =================================== Read the control word header for the current adaptation level
int vtkLSDynaReader::ReadHeaderInformation(int curAdapt)
{
  LSDynaMetaData* p = this->P;

  // =================================== Control Word Section
  p->Fam.SkipToWord(LSDynaFamily::ControlSection, curAdapt /*timestep*/, 0);
  p->Fam.BufferChunk(LSDynaFamily::Char, 10);
  memcpy(p->Title, p->Fam.GetNextWordAsChars(), 40 * sizeof(char));
  p->Title[40] = '\0';

  p->Fam.SkipToWord(LSDynaFamily::ControlSection, curAdapt /*timestep*/, 13);

  // Release number: Release number in character*4 form: 50 for R5.0, 511c for R5.1.1c.
  p->Fam.BufferChunk(LSDynaFamily::Char, 1);
  memcpy(p->ReleaseNumber, p->Fam.GetNextWordAsChars(), 4 * sizeof(char));
  p->ReleaseNumber[4] = '\0';

  // Version: Code version, floating number, eg 960.0 it is used to distinguish the
  // floating point format, like cray, ieee, and dpieee.
  p->Fam.BufferChunk(LSDynaFamily::Float, 1);
  p->CodeVersion = p->Fam.GetNextWordAsFloat();

  p->Fam.BufferChunk(LSDynaFamily::Int, 49);
  p->Dict["NDIM"] = p->Fam.GetNextWordAsInt();
  p->Dict["NUMNP"] = p->Fam.GetNextWordAsInt();
  p->Dict["ICODE"] = p->Fam.GetNextWordAsInt();
  p->Dict["NGLBV"] = p->Fam.GetNextWordAsInt();
  p->Dict["IT"] = p->Fam.GetNextWordAsInt();
  p->Dict["IU"] = p->Fam.GetNextWordAsInt();
  p->Dict["IV"] = p->Fam.GetNextWordAsInt();
  p->Dict["IA"] = p->Fam.GetNextWordAsInt();
  p->Dict["NEL8"] = p->Fam.GetNextWordAsInt();
  p->Dict["NUMMAT8"] = p->Fam.GetNextWordAsInt();
  p->Fam.GetNextWordAsInt(); // BLANK
  p->Fam.GetNextWordAsInt(); // BLANK
  p->Dict["NV3D"] = p->Fam.GetNextWordAsInt();
  p->Dict["NEL2"] = p->Fam.GetNextWordAsInt();
  p->Dict["NUMMAT2"] = p->Fam.GetNextWordAsInt();
  p->Dict["NV1D"] = p->Fam.GetNextWordAsInt();
  p->Dict["NEL4"] = p->Fam.GetNextWordAsInt();
  p->Dict["NUMMAT4"] = p->Fam.GetNextWordAsInt();
  p->Dict["NV2D"] = p->Fam.GetNextWordAsInt();
  p->Dict["NEIPH"] = p->Fam.GetNextWordAsInt();
  p->Dict["NEIPS"] = p->Fam.GetNextWordAsInt();
  p->Dict["MAXINT"] = p->Fam.GetNextWordAsInt();
  // do MDLOPT here?
  p->Dict["NMSPH"] = p->Fam.GetNextWordAsInt();
  p->Dict["EDLOPT"] = p->Dict["NMSPH"]; // EDLOPT is not standard
  p->Dict["NGPSPH"] = p->Fam.GetNextWordAsInt();
  p->Dict["NARBS"] = p->Fam.GetNextWordAsInt();
  p->Dict["NELT"] = p->Fam.GetNextWordAsInt();
  p->Dict["NUMMATT"] = p->Fam.GetNextWordAsInt();
  p->Dict["NV3DT"] = p->Fam.GetNextWordAsInt();
  p->Dict["IOSHL(1)"] = p->Fam.GetNextWordAsInt() == 1000 ? 1 : 0;
  p->Dict["IOSHL(2)"] = p->Fam.GetNextWordAsInt() == 1000 ? 1 : 0;
  p->Dict["IOSHL(3)"] = p->Fam.GetNextWordAsInt() == 1000 ? 1 : 0;
  p->Dict["IOSHL(4)"] = p->Fam.GetNextWordAsInt() == 1000 ? 1 : 0;
  p->Dict["IALEMAT"] = p->Fam.GetNextWordAsInt();
  p->Dict["NCFDV1"] = p->Fam.GetNextWordAsInt();
  p->Dict["NCFDV2"] = p->Fam.GetNextWordAsInt();
  p->Dict["NADAPT"] = p->Fam.GetNextWordAsInt();
  p->Dict["NMMAT"] = p->Fam.GetNextWordAsInt();

  // NUMFLUID: Total number of ALE fluid groups. Fluid density and
  // volume fractions output as history variables, and a flag
  // for the dominant group. If negative multi-material
  // species mass for each group is also output. Order is: rho,
  // vf1, ... vfn, dvf flag, m1, ... mn. Density is at position 8
  // after the location for plastic strain. Any element material
  // history variables are written before the Ale variables, and
  // the six element strains components after these if
  // ISTRN=1.
  p->Dict["NUMFLUID"] = p->Fam.GetNextWordAsInt();

  // Compute the derived values in this->P
  // =========================================== Control Word Section Processing
  int itmp;
  vtkIdType iddtmp;
  char ctmp[128]; // temp space for generating keywords (i.e. isphfg) and attribute names (i.e.,
                  // StressIntPt3)

  // --- Initialize some values
  p->ReadRigidRoadMvmt = 0;
  p->PreStateSize = 64 * p->Fam.GetWordSize();
  p->StateSize = p->Fam.GetWordSize(); // Account for "time word"
  p->Dimensionality = p->Dict["NDIM"];
  switch (p->Dimensionality)
  {
    case 2:
    case 3:
      p->Dict["MATTYP"] = 0;
      p->ConnectivityUnpacked = 0;
      break;
    case 7:
      p->ReadRigidRoadMvmt = 1;
      VTK_FALLTHROUGH;
    case 5:
      p->Dict["MATTYP"] = 1;
      p->ConnectivityUnpacked = 1;
      p->Dimensionality = 3;
      break;
    case 4:
      p->ConnectivityUnpacked = 1;
      p->Dict["MATTYP"] = 0;
      p->Dimensionality = 3;
      break;
    default:
      vtkErrorMacro("Unknown Dimensionality " << p->Dimensionality << " encountered");
      p->FileIsValid = 0;
      return 0;
  }

  // FIXME Are these marks valid since we are marking the word past the end of the chunk?
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::StaticSection);
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::MaterialTypeData);
  if (p->Dict["MATTYP"] != 0)
  {
    p->Fam.BufferChunk(LSDynaFamily::Int, 2);
    p->Dict["NUMRBE"] = p->Fam.GetNextWordAsInt();
    p->Dict["NUMMAT"] = p->Fam.GetNextWordAsInt();
  }
  else
  {
    p->Dict["NUMRBE"] = 0;
    p->Dict["NUMMAT"] = 0;
  }
  p->NumberOfNodes = p->Dict["NUMNP"];

  p->NumberOfCells[LSDynaMetaData::RIGID_BODY] = p->Dict["NUMRBE"];
  p->NumberOfCells[LSDynaMetaData::SOLID] = p->Dict["NEL8"];
  p->NumberOfCells[LSDynaMetaData::THICK_SHELL] = p->Dict["NELT"];
  p->NumberOfCells[LSDynaMetaData::SHELL] = p->Dict["NEL4"];
  p->NumberOfCells[LSDynaMetaData::BEAM] = p->Dict["NEL2"];
  p->NumberOfCells[LSDynaMetaData::PARTICLE] = p->Dict["NMSPH"];

  p->StateSize += p->Dict["NGLBV"] * p->Fam.GetWordSize();

  if (p->Dict["IT"] != 0)
  {
    p->AddPointArray(LS_ARRAYNAME_TEMPERATURE, 1, 1);
    p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
  }
  if (p->Dict["IU"] != 0)
  {
    p->AddPointArray(LS_ARRAYNAME_DEFLECTION, p->Dimensionality, 1);
    p->StateSize += p->NumberOfNodes * p->Dimensionality * p->Fam.GetWordSize();
  }
  if (p->Dict["IV"] != 0)
  {
    p->AddPointArray(LS_ARRAYNAME_VELOCITY, p->Dimensionality, 1);
    p->StateSize += p->NumberOfNodes * p->Dimensionality * p->Fam.GetWordSize();
  }
  if (p->Dict["IA"] != 0)
  {
    p->AddPointArray(LS_ARRAYNAME_ACCELERATION, p->Dimensionality, 1);
    p->StateSize += p->NumberOfNodes * p->Dimensionality * p->Fam.GetWordSize();
  }
  p->Dict["cfdPressure"] = 0;
  p->Dict["cfdVort"] = 0;
  p->Dict["cfdXVort"] = 0;
  p->Dict["cfdYVort"] = 0;
  p->Dict["cfdZVort"] = 0;
  p->Dict["cfdRVort"] = 0;
  p->Dict["cfdEnstrophy"] = 0;
  p->Dict["cfdHelicity"] = 0;
  p->Dict["cfdStream"] = 0;
  p->Dict["cfdEnthalpy"] = 0;
  p->Dict["cfdDensity"] = 0;
  p->Dict["cfdTurbKE"] = 0;
  p->Dict["cfdDiss"] = 0;
  p->Dict["cfdEddyVisc"] = 0;
  itmp = p->Dict["NCFDV1"];
  if (itmp & 2)
  {
    p->AddPointArray(LS_ARRAYNAME_PRESSURE, 1, 1);
    p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
    p->Dict["cfdPressure"] = 1;
  }
  if ((itmp & 28) == 28)
  {
    p->AddPointArray(LS_ARRAYNAME_VORTICITY, 3, 1);
    p->StateSize += p->NumberOfNodes * 3 * p->Fam.GetWordSize();
    p->Dict["cfdVort"] = 1;
    p->Dict["cfdXVort"] = 1;
    p->Dict["cfdYVort"] = 1;
    p->Dict["cfdZVort"] = 1;
  }
  else
  { // OK, we don't have all the vector components... maybe we have some of them?
    if (itmp & 4)
    {
      p->AddPointArray(LS_ARRAYNAME_VORTICITY_X, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdXVort"] = 1;
    }
    if (itmp & 8)
    {
      p->AddPointArray(LS_ARRAYNAME_VORTICITY_Y, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdYVort"] = 1;
    }
    if (itmp & 16)
    {
      p->AddPointArray(LS_ARRAYNAME_VORTICITY_Z, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdZVort"] = 1;
    }
  }
  if (itmp & 32)
  {
    p->AddPointArray(LS_ARRAYNAME_RESULTANTVORTICITY, 1, 1);
    p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
    p->Dict["cfdRVort"] = 1;
    if (itmp & 64)
    {
      p->AddPointArray(LS_ARRAYNAME_ENSTROPHY, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdEnstrophy"] = 1;
    }
    if (itmp & 128)
    {
      p->AddPointArray(LS_ARRAYNAME_HELICITY, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdHelicity"] = 1;
    }
    if (itmp & 256)
    {
      p->AddPointArray(LS_ARRAYNAME_STREAMFUNCTION, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdStream"] = 1;
    }
    if (itmp & 512)
    {
      p->AddPointArray(LS_ARRAYNAME_ENTHALPY, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdEnthalpy"] = 1;
    }
    if (itmp & 1024)
    {
      p->AddPointArray(LS_ARRAYNAME_DENSITY, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdDensity"] = 1;
    }
    if (itmp & 2048)
    {
      p->AddPointArray(LS_ARRAYNAME_TURBULENTKE, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdTurbKE"] = 1;
    }
    if (itmp & 4096)
    {
      p->AddPointArray(LS_ARRAYNAME_DISSIPATION, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdDiss"] = 1;
    }
    if (itmp & 1040384)
    {
      p->AddPointArray(LS_ARRAYNAME_EDDYVISCOSITY, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      p->Dict["cfdEddyVisc"] = 1;
    }
  }

  char sname[] = LS_ARRAYNAME_SPECIES_BLNK;
  iddtmp = p->Dict["NCFDV2"];
  for (itmp = 1; itmp < 11; ++itmp)
  {
    if (iddtmp & (static_cast<vtkIdType>(1) << itmp))
    {
      snprintf(sname, sizeof(sname), LS_ARRAYNAME_SPECIES_FMT, itmp);
      p->AddPointArray(sname, 1, 1);
      p->StateSize += p->NumberOfNodes * p->Fam.GetWordSize();
      snprintf(sname, sizeof(sname), "cfdSpec%02d", itmp);
      p->Dict[sname] = 1;
    }
    else
    {
      snprintf(sname, sizeof(sname), "cfdSpec%02d", itmp);
      p->Dict[sname] = 0;
    }
  }

  // Solid element state size   FIXME: 7 + NEIPH should really be NV3D (in case things change)
  p->StateSize +=
    (7 + p->Dict["NEIPH"]) * p->NumberOfCells[LSDynaMetaData::SOLID] * p->Fam.GetWordSize();
  // Thick shell state size
  p->StateSize +=
    p->Dict["NV3DT"] * p->NumberOfCells[LSDynaMetaData::THICK_SHELL] * p->Fam.GetWordSize();
  // (Thin) shell state size (we remove rigid body shell element state below)
  p->StateSize += p->Dict["NV2D"] * p->NumberOfCells[LSDynaMetaData::SHELL] * p->Fam.GetWordSize();
  // Beam state size
  p->StateSize += p->Dict["NV1D"] * p->NumberOfCells[LSDynaMetaData::BEAM] * p->Fam.GetWordSize();

  // ================================================ Static Information Section
  // This is marked above so we can read NUMRBE in time to do StateSize calculations
  // ================================================ Material Type Data Section
  // This is marked above so we can read NUMRBE in time to do StateSize calculations
  if (p->Dict["MATTYP"] != 0)
  {
    // If there are rigid body elements, they won't have state data and we must
    // reduce the state size
    p->StateSize -= p->Dict["NV2D"] * p->NumberOfCells[LSDynaMetaData::RIGID_BODY];

    p->Fam.BufferChunk(LSDynaFamily::Int, p->Dict["NUMMAT"]);
    for (itmp = 0; itmp < p->Dict["NUMMAT"]; ++itmp)
    {
      p->RigidMaterials.insert(p->Fam.GetNextWordAsInt());
    }
    p->PreStateSize += (2 + p->Dict["NUMMAT"]) * p->Fam.GetWordSize();
  }

  // process the deletion array
  // save the position we currently have as it is the offset to jump to
  // when reading deletion
  p->ElementDeletionOffset = p->StateSize / p->Fam.GetWordSize();

  int mdlopt;
  int intpts2;
  mdlopt = p->Dict["MAXINT"];
  if (mdlopt >= 0 && mdlopt <= 10000)
  {
    intpts2 = mdlopt;
    mdlopt = LS_MDLOPT_NONE;
  }
  else if (mdlopt < -10000)
  {
    intpts2 = -mdlopt - 10000;
    mdlopt = LS_MDLOPT_CELL;

    // WARNING: This needs NumberOfCells[LSDynaMetaData::RIGID_BODY] set, which relies on NUMRBE
    p->StateSize += this->GetNumberOfContinuumCells() * p->Fam.GetWordSize();
  }
  else if (mdlopt > 10000)
  {
    intpts2 = mdlopt - 10000;
    mdlopt = LS_MDLOPT_CELL;

    // WARNING: This needs NumberOfCells[LSDynaMetaData::RIGID_BODY] set, which relies on NUMRBE
    p->StateSize += this->GetNumberOfContinuumCells() * p->Fam.GetWordSize();
  }
  else
  {
    intpts2 = -mdlopt;
    mdlopt = LS_MDLOPT_POINT;
    // p->AddPointArray( LS_ARRAYNAME_DEATH, 1, 1 );
    p->StateSize += this->GetNumberOfNodes() * p->Fam.GetWordSize();
  }
  p->Dict["MDLOPT"] = mdlopt;
  p->Dict["_MAXINT_"] = intpts2;
  if (p->Dict["NV2D"] > 0)
  {
    if (p->Dict["NV2D"] -
        (p->Dict["_MAXINT_"] * (6 * p->Dict["IOSHL(1)"] + p->Dict["IOSHL(2)"] + p->Dict["NEIPS"]) +
          8 * p->Dict["IOSHL(3)"] + 4 * p->Dict["IOSHL(4)"]) >
      1)
    {
      p->Dict["ISTRN"] = 1;
    }
    else
    {
      p->Dict["ISTRN"] = 0;
    }
  }
  else if (p->Dict["NELT"] > 0)
  {
    if (p->Dict["NV3D"] -
        p->Dict["_MAXINT_"] * (6 * p->Dict["IOSHL(1)"] + p->Dict["IOSHL(2)"] + p->Dict["NEIPS"]) >
      1)
    {
      p->Dict["ISTRN"] = 1;
    }
    else
    {
      p->Dict["ISTRN"] = 0;
    }
  }
  else
  {
    p->Dict["ISTRN"] = 0;
  }

  // OK, we are done processing the header (control) section.

  p->SPHStateOffset = p->StateSize / p->Fam.GetWordSize();
  // ============================================ Fluid Material ID Data Section
  // IALEMAT offset
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::FluidMaterialIdData);
  p->PreStateSize += p->Dict["IALEMAT"];
  p->Fam.BufferChunk(LSDynaFamily::Int, p->Dict["IALEMAT"]);
  for (itmp = 0; itmp < p->Dict["IALEMAT"]; ++itmp)
  {
    p->FluidMaterials.insert(p->Fam.GetNextWordAsInt());
  }
  // p->Fam.SkipToWord( LSDynaFamily::FluidMaterialIdData, curAdapt, p->Dict["IALEMAT"] );

  // =======================  Smooth Particle Hydrodynamics Element Data Section
  // Only when NMSPH > 0
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::SPHElementData);
  if (p->NumberOfCells[LSDynaMetaData::PARTICLE] > 0)
  {
    p->Fam.BufferChunk(LSDynaFamily::Int, 1);
    vtkIdType sphAttributes = p->Fam.GetNextWordAsInt();
    p->Dict["isphfg(1)"] = sphAttributes;
    if (sphAttributes >= 9)
    {
      p->Fam.BufferChunk(LSDynaFamily::Int, sphAttributes - 1); // should be 9 or 10
      // Dyna docs call statePerParticle "NUM_SPH_DATA":
      int statePerParticle = 1; // start at 1 because we always have material ID of particle.
      for (itmp = 2; itmp <= sphAttributes; ++itmp)
      {
        int numComponents = p->Fam.GetNextWordAsInt();
        snprintf(ctmp, sizeof(ctmp), "isphfg(%d)", itmp);
        p->Dict[ctmp] = numComponents;
        statePerParticle += numComponents;
      }
      p->Dict["NUM_SPH_DATA"] = statePerParticle;
      p->StateSize +=
        statePerParticle * p->NumberOfCells[LSDynaMetaData::PARTICLE] * p->Fam.GetWordSize();
    }
    else
    {
      p->FileIsValid = 0;
      return 0;
    }
    p->Fam.SkipToWord(LSDynaFamily::SPHElementData, curAdapt, p->Dict["isphfg(1)"]);
    p->PreStateSize += p->Dict["isphfg(1)"] * p->Fam.GetWordSize();
  }

  // ===================================================== Geometry Data Section
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::GeometryData);
  iddtmp =
    this->GetNumberOfNodes() * p->Dimensionality * p->Fam.GetWordSize(); // Size of nodes on disk
  iddtmp +=
    p->NumberOfCells[LSDynaMetaData::SOLID] * 9 * p->Fam.GetWordSize(); // Size of hexes on disk
  iddtmp += p->NumberOfCells[LSDynaMetaData::THICK_SHELL] * 9 *
    p->Fam.GetWordSize(); // Size of thick shells on disk
  iddtmp +=
    p->NumberOfCells[LSDynaMetaData::SHELL] * 5 * p->Fam.GetWordSize(); // Size of quads on disk
  iddtmp +=
    p->NumberOfCells[LSDynaMetaData::BEAM] * 6 * p->Fam.GetWordSize(); // Size of beams on disk
  p->PreStateSize += iddtmp;
  p->Fam.SkipToWord(
    LSDynaFamily::GeometryData, curAdapt, iddtmp / p->Fam.GetWordSize()); // Skip to end of geometry

  // =========== User Material, Node, And Element Identification Numbers Section
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::UserIdData);
  if (p->Dict["NARBS"] != 0)
  {
    // For sequential material numbering
    // NARBS = 10+NUMNP+NEL8+NEL2+NEL4+NELT+3*NMMAT   (even though the 3*NMMAT numbers are unused)
    // For arbitrary material numbering (NSORT<0)
    // NARBS = 16+NUMNP+NEL8+NEL2+NEL4+NELT+3*NMMAT   (Note 6 extra params)

    p->Fam.BufferChunk(LSDynaFamily::Int, 10);
    p->Dict["NSORT"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRH"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRB"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRS"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRT"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSORTD"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRHD"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRBD"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRSD"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSRTD"] = p->Fam.GetNextWordAsInt();

    if (p->Dict["NSORT"] < 0)
    {
      p->Fam.BufferChunk(LSDynaFamily::Int, 6);
      p->Dict["NSRMA"] = p->Fam.GetNextWordAsInt();
      p->Dict["NSRMU"] = p->Fam.GetNextWordAsInt();
      p->Dict["NSRMP"] = p->Fam.GetNextWordAsInt();
      p->Dict["NSRTM"] = p->Fam.GetNextWordAsInt();
      p->Dict["NUMRBS"] = p->Fam.GetNextWordAsInt();
      p->Dict["NMMAT"] = p->Fam.GetNextWordAsInt();
    }

    iddtmp = p->Dict["NARBS"] * p->Fam.GetWordSize();
    p->PreStateSize += iddtmp;
    p->Fam.SkipToWord(LSDynaFamily::UserIdData, curAdapt, p->Dict["NARBS"]);
  }
  else
  {
    p->Dict["NSORT"] = 0;
  }
  // Break from convention and read in actual array values (as opposed to just summary information)
  // about the material IDs. This is required because the reader must present part names after
  // RequestInformation is called and that cannot be done without a list of material IDs.
  this->ReadUserMaterialIds();

  // ======================================= Adapted Element Parent List Section
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::AdaptedParentData);
  p->Fam.SkipToWord(LSDynaFamily::AdaptedParentData, curAdapt, 2 * p->Dict["NADAPT"]);
  iddtmp = 2 * p->Dict["NADAPT"] * p->Fam.GetWordSize();
  p->PreStateSize += iddtmp;

  // ============== Smooth Particle Hydrodynamics Node And Material List Section
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::SPHNodeData);
  iddtmp = 2 * p->NumberOfCells[LSDynaMetaData::PARTICLE] * p->Fam.GetWordSize();
  p->PreStateSize += iddtmp;
  p->Fam.SkipToWord(
    LSDynaFamily::SPHNodeData, curAdapt, 2 * p->NumberOfCells[LSDynaMetaData::PARTICLE]);

  // =========================================== Rigid Road Surface Data Section
  p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::RigidSurfaceData);
  if (p->Dict["NDIM"] > 5)
  {
    p->Fam.BufferChunk(LSDynaFamily::Int, 4);
    p->PreStateSize += 4 * p->Fam.GetWordSize();
    p->Dict["NNODE"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSEG"] = p->Fam.GetNextWordAsInt();
    p->Dict["NSURF"] = p->Fam.GetNextWordAsInt();
    p->Dict["MOTION"] = p->Fam.GetNextWordAsInt();
    iddtmp = 4 * p->Dict["NNODE"] * p->Fam.GetWordSize();
    p->PreStateSize += iddtmp;
    p->Fam.SkipWords(4 * p->Dict["NNODE"]);
    p->NumberOfCells[LSDynaMetaData::ROAD_SURFACE] = p->Dict["NSEG"];

    for (itmp = 0; itmp < p->Dict["NSURF"]; ++itmp)
    {
      p->Fam.BufferChunk(LSDynaFamily::Int, 2);
      p->Fam.GetNextWordAsInt();          // Skip SURFID
      iddtmp = p->Fam.GetNextWordAsInt(); // Read SURFNSEG[SURFID]
      p->RigidSurfaceSegmentSizes.push_back(iddtmp);
      p->PreStateSize += (2 + 4 * iddtmp) * p->Fam.GetWordSize();
      p->Fam.SkipWords(4 * iddtmp);
    }

    if (p->Dict["NSEG"] > 0)
    {
      p->AddCellArray(LSDynaMetaData::ROAD_SURFACE, LS_ARRAYNAME_SEGMENTID, 1, 1);
      // FIXME: p->AddRoadPointArray( LSDynaMetaData::ROAD_SURFACE, LS_ARRAYNAME_USERID, 1, 1 );
    }

    if (p->Dict["MOTION"])
    {
      p->StateSize += 6 * p->Dict["NSURF"] * p->Fam.GetWordSize();
    }
  }

  // if ( curAdapt == 0 )
  {
    p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::EndOfStaticSection);
    p->Fam.MarkSectionStart(curAdapt, LSDynaFamily::TimeStepSection);
  }
  p->Fam.SetStateSize(p->StateSize / p->Fam.GetWordSize());

  // ==========================================================================
  // Now that we've read the header, create a list of the cell arrays for
  // each output mesh.

  // Currently, the LS-Dyna reader only gives users a few knobs to control which cell variables are
  // loaded. It is a difficult problem since many attributes only occur on some element types, there
  // are many dyna flags that conditionally omit results, and some quantities are repeated over
  // differing numbers of points for different types of cells. Given the complexity, we punt by
  // defining some knobs for "types" of data and define related fields. In a perfect world,
  // finer-grained control would be available.
  //
  // As an example: if there are any
  // - 3-D cells, OR
  // - non-rigid 2-D cells with IOSHL(1)==1, OR
  // - beam cells with NV1D > 6, OR
  // - SPH cells with isphfg(4)==6
  // then there will be stresses present

  // Every cell always has a material type
  // FIXME: Is this true? Rigid bodies may be an exception, in which
  // case we need to check that the number of cells in the other 5 meshes sum to >0

  // Total number of ALE fluid groups. Fluid density and
  // volume fractions output as history variables, and a flag
  // for the dominant group. If negative multi-material
  // species mass for each group is also output. Order is: rho,
  // vf1, … vfn, dvf flag, m1, … mn. Density is at position 8
  // after the location for plastic strain. Any element material
  // history variables are written before the Ale variables, and
  // the six element strains components after these if
  // ISTRN=1.
  int numGroups = std::abs(static_cast<int>(p->Dict["NUMFLUID"]));
  bool hasMass = (p->Dict["NUMFLUID"] < 0);
  int numALEvalues = (numGroups > 0) ? 1 + numGroups + 1 + (hasMass ? numGroups : 0) : 0;

  if (p->Dict["NARBS"])
  {
    p->AddPointArray(LS_ARRAYNAME_USERID, 1, 1);
  }

  if (p->NumberOfCells[LSDynaMetaData::PARTICLE])
  {
    // One value is always output which is the material number as a floating point number for each
    // particle. If this value is negative then the particle has been deleted from the model.
    p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_MATERIAL, 1, 1);
    // p->AddCellArray( LSDynaMetaData::PARTICLE, LS_ARRAYNAME_DEATH, 1, 1 );  // There is no
    // mention of this in manual... Was this the -ve Material case?
    if (p->Dict["isphfg(2)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_RADIUSOFINFLUENCE, 1, 1);
    }
    if (p->Dict["isphfg(3)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_PRESSURE, 1, 1);
    }
    if (p->Dict["isphfg(4)"] == 6)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_STRESS, 6, 1);
    }
    if (p->Dict["isphfg(5)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_EPSTRAIN, 1, 1);
    }
    if (p->Dict["isphfg(6)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_DENSITY, 1, 1);
    }
    if (p->Dict["isphfg(7)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_INTERNALENERGY, 1, 1);
    }
    if (p->Dict["isphfg(8)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_NUMNEIGHBORS, 1, 1);
    }
    if (p->Dict["isphfg(9)"] == 6)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_STRAIN, 6, 1);
    }
    if (p->Dict["isphfg(10)"] == 1)
    {
      p->AddCellArray(LSDynaMetaData::PARTICLE, LS_ARRAYNAME_MASS, 1, 1);
    }
  }

  if (p->NumberOfCells[LSDynaMetaData::BEAM])
  {
    //    p->AddCellArray( LSDynaMetaData::BEAM, LS_ARRAYNAME_MATERIAL, 1, 1 );
    //    if ( p->Dict["MDLOPT"] == LS_MDLOPT_CELL )
    //      {
    //      p->AddCellArray( LSDynaMetaData::BEAM, LS_ARRAYNAME_DEATH, 1, 1 );
    //      }

    if (p->Dict["NARBS"] != 0)
    {
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_USERID, 1, 1);
    }
    if (p->Dict["NV1D"] >= 6)
    {
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_AXIALFORCE, 1, 1);
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_SHEARRESULTANT, 2, 1);
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_BENDINGRESULTANT, 2, 1);
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_TORSIONRESULTANT, 1, 1);
    }
    if (p->Dict["NV1D"] > 6)
    {
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_SHEARSTRESS, 2, 1);
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_AXIALSTRESS, 1, 1);
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_PLASTICSTRAIN, 1, 1);
      p->AddCellArray(LSDynaMetaData::BEAM, LS_ARRAYNAME_AXIALSTRAIN, 1, 1);
    }
  }

  if (p->NumberOfCells[LSDynaMetaData::SHELL])
  {
    //    p->AddCellArray( LSDynaMetaData::SHELL, LS_ARRAYNAME_MATERIAL, 1, 1 );
    //    if ( p->Dict["MDLOPT"] == LS_MDLOPT_CELL )
    //      {
    //      p->AddCellArray( LSDynaMetaData::SHELL, LS_ARRAYNAME_DEATH, 1, 1 );
    //      }
    if (p->Dict["NARBS"] != 0)
    {
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_USERID, 1, 1);
    }
    if (p->Dict["IOSHL(1)"])
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_STRESS, 6, 1);
    if (p->Dict["IOSHL(2)"])
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_EPSTRAIN, 1, 1);

    int neips = p->Dict["NEIPS"];
    int extraValues = neips;
    if (extraValues > 0)
    {
      // Any element material history variables are written before the Ale variables, and the six
      // element strains components after these if ISTRN=1
      int materialValues = extraValues - numALEvalues;
      if (materialValues > 0)
      {
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_INTEGRATIONPOINT, materialValues, 1);
        extraValues -= materialValues;
      }

      if ((numALEvalues > 0) && (extraValues >= numALEvalues))
      {
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_DENSITY, 1, 1);
        extraValues--;

        for (int g = 0; g < numGroups; ++g)
        {
          snprintf(ctmp, sizeof(ctmp), LS_ARRAYNAME_VOLUME_FRACTION_FMT, static_cast<int>(g + 1));
          p->AddCellArray(LSDynaMetaData::SHELL, ctmp, 1, 1);
          extraValues--;
        }

        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_DOMINANT_GROUP, 1, 1);
        extraValues--;

        for (int g = 0; hasMass && (g < numGroups); ++g)
        {
          snprintf(ctmp, sizeof(ctmp), LS_ARRAYNAME_SPECIES_MASS_FMT, static_cast<int>(g + 1));
          p->AddCellArray(LSDynaMetaData::SHELL, ctmp, 1, 1);
          extraValues--;
        }
      }
    }

    if (p->Dict["_MAXINT_"] >= 2)
    {
      if (p->Dict["IOSHL(1)"])
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_STRESS "InnerSurf", 6, 1);
      if (p->Dict["IOSHL(2)"])
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_EPSTRAIN "InnerSurf", 1, 1);
      if (neips)
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_INTEGRATIONPOINT "InnerSurf", neips, 1);
    }
    if (p->Dict["_MAXINT_"] >= 3)
    {
      if (p->Dict["IOSHL(1)"])
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_STRESS "OuterSurf", 6, 1);
      if (p->Dict["IOSHL(2)"])
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_EPSTRAIN "OuterSurf", 1, 1);
      if (neips)
        p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_INTEGRATIONPOINT "OuterSurf", neips, 1);
    }

    for (itmp = 3; itmp < p->Dict["_MAXINT_"]; ++itmp)
    {
      if (p->Dict["IOSHL(1)"])
      {
        snprintf(ctmp, sizeof(ctmp), "%sIntPt%d", LS_ARRAYNAME_STRESS, itmp + 1);
        p->AddCellArray(LSDynaMetaData::SHELL, ctmp, 6, 1);
      }
      if (p->Dict["IOSHL(2)"])
      {
        snprintf(ctmp, sizeof(ctmp), "%sIntPt%d", LS_ARRAYNAME_EPSTRAIN, itmp + 1);
        p->AddCellArray(LSDynaMetaData::SHELL, ctmp, 1, 1);
      }
      if (neips)
      {
        snprintf(ctmp, sizeof(ctmp), "%sIntPt%d", LS_ARRAYNAME_INTEGRATIONPOINT, itmp + 1);
        p->AddCellArray(LSDynaMetaData::SHELL, ctmp, neips, 1);
      }
    }
    if (p->Dict["IOSHL(3)"])
    {
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_NORMALRESULTANT, 3, 1);
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_SHEARRESULTANT, 2, 1);
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_BENDINGRESULTANT, 3, 1);
    }
    if (p->Dict["IOSHL(4)"])
    {
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_THICKNESS, 1, 1);
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_ELEMENTMISC, 2, 1);
    }

    if (p->Dict["ISTRN"])
    {
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_STRAIN "InnerSurf", 6, 1);
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_STRAIN "OuterSurf", 6, 1);
    }
    if (!p->Dict["ISTRN"] || (p->Dict["NV2D"] >= 45))
    {
      p->AddCellArray(LSDynaMetaData::SHELL, LS_ARRAYNAME_INTERNALENERGY, 1, 1);
    }
  }

  if (p->NumberOfCells[LSDynaMetaData::THICK_SHELL])
  {
    //    p->AddCellArray( LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_MATERIAL, 1, 1 );
    //    if ( p->Dict["MDLOPT"] == LS_MDLOPT_CELL )
    //      {
    //      p->AddCellArray( LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_DEATH, 1, 1 );
    //      }
    if (p->Dict["NARBS"] != 0)
    {
      p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_USERID, 1, 1);
    }
    if (p->Dict["IOSHL(1)"])
    {
      if (p->Dict["_MAXINT_"] >= 3)
      {
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_STRESS, 6, 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_STRESS "InnerSurf", 6, 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_STRESS "OuterSurf", 6, 1);
      }
      for (itmp = 3; itmp < p->Dict["_MAXINT_"]; ++itmp)
      {
        snprintf(ctmp, sizeof(ctmp), "%sIntPt%d", LS_ARRAYNAME_STRESS, itmp + 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, ctmp, 6, 1);
      }
    }
    if (p->Dict["IOSHL(2)"])
    {
      if (p->Dict["_MAXINT_"] >= 3)
      {
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_EPSTRAIN, 1, 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_EPSTRAIN "InnerSurf", 1, 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_EPSTRAIN "OuterSurf", 1, 1);
      }
      for (itmp = 3; itmp < p->Dict["_MAXINT_"]; ++itmp)
      {
        snprintf(ctmp, sizeof(ctmp), "%sIntPt%d", LS_ARRAYNAME_EPSTRAIN, itmp + 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, ctmp, 1, 1);
      }
    }
    if (p->Dict["NEIPS"])
    {
      int neips = p->Dict["NEIPS"];
      if (p->Dict["_MAXINT_"] >= 3)
      {
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_INTEGRATIONPOINT, neips, 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_INTEGRATIONPOINT, neips, 1);
        p->AddCellArray(
          LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_INTEGRATIONPOINT "InnerSurf", neips, 1);
        p->AddCellArray(
          LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_INTEGRATIONPOINT "OuterSurf", neips, 1);
      }
      for (itmp = 3; itmp < p->Dict["_MAXINT_"]; ++itmp)
      {
        snprintf(ctmp, sizeof(ctmp), "%sIntPt%d", LS_ARRAYNAME_INTEGRATIONPOINT, itmp + 1);
        p->AddCellArray(LSDynaMetaData::THICK_SHELL, ctmp, 6, 1);
      }
    }
    if (p->Dict["ISTRN"])
    {
      p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_STRAIN "InnerSurf", 6, 1);
      p->AddCellArray(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_STRAIN "OuterSurf", 6, 1);
    }
  }

  if (p->NumberOfCells[LSDynaMetaData::SOLID])
  {
    //    p->AddCellArray( LSDynaMetaData::SOLID, LS_ARRAYNAME_MATERIAL, 1, 1 );
    //    if ( p->Dict["MDLOPT"] == LS_MDLOPT_CELL )
    //      {
    //      p->AddCellArray( LSDynaMetaData::SOLID, LS_ARRAYNAME_DEATH, 1, 1 );
    //      }
    if (p->Dict["NARBS"] != 0)
    {
      p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_USERID, 1, 1);
    }
    p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_STRESS, 6, 1);
    p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_EPSTRAIN, 1, 1);

    int extraValues = p->Dict["NEIPH"];
    if (extraValues > 0)
    {
      int strainValues = (p->Dict["ISTRN"] == 1) ? 6 : 0; // last six values are strain.

      // Any element material history variables are written before the Ale variables, and the six
      // element strains components after these if ISTRN=1
      int materialValues = extraValues - (numALEvalues + strainValues);
      if (materialValues > 0)
      {
        p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_INTEGRATIONPOINT, materialValues, 1);
        extraValues -= materialValues;
      }

      if ((numALEvalues > 0) && (extraValues >= numALEvalues))
      {
        p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_DENSITY, 1, 1);
        extraValues--;

        for (int g = 0; g < numGroups; ++g)
        {
          snprintf(ctmp, sizeof(ctmp), LS_ARRAYNAME_VOLUME_FRACTION_FMT, static_cast<int>(g + 1));
          p->AddCellArray(LSDynaMetaData::SOLID, ctmp, 1, 1);
          extraValues--;
        }

        p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_DOMINANT_GROUP, 1, 1);
        extraValues--;

        for (int g = 0; hasMass && (g < numGroups); ++g)
        {
          snprintf(ctmp, sizeof(ctmp), LS_ARRAYNAME_SPECIES_MASS_FMT, static_cast<int>(g + 1));
          p->AddCellArray(LSDynaMetaData::SOLID, ctmp, 1, 1);
          extraValues--;
        }
      }
      if ((strainValues > 0) && (extraValues >= strainValues))
      {
        p->AddCellArray(LSDynaMetaData::SOLID, LS_ARRAYNAME_STRAIN, strainValues, 1);
      }
    }
  }

  // Only try reading the keyword file if we don't have part names.
  if (curAdapt == 0 && p->PartNames.empty())
  {
    this->ResetPartInfo();

    int result = this->ReadInputDeck();

    if (result == 0)
    {
      // we failed to read the input deck so we are going to read the first binary file for part
      // names
      this->ReadPartTitlesFromRootFile();
    }
  }

  return -1;
}

int vtkLSDynaReader::ScanDatabaseTimeSteps()
{
  LSDynaMetaData* p = this->P;

  // ======================================================= State Data Sections
  // The 2 lines below are now in ReadHeaderInformation:
  // p->Fam.MarkSectionStart( curAdapt, LSDynaFamily::TimeStepSection );
  // p->Fam.SetStateSize( p->StateSize / p->Fam.GetWordSize() );
  // It may be useful to call
  // p->JumpToMark( LSDynaFamily::TimeStepSection );
  // here.
  if (p->Fam.GetStateSize() <= 0)
  {
    vtkErrorMacro("Database has bad state size (" << p->Fam.GetStateSize() << ").");
    return 1;
  }

  // Discover the number of states and record the time value for each.
  int ntimesteps = 0;
  double time;
  int itmp = 1;
  int lastAdapt = 0;
  do
  {
    if (p->Fam.BufferChunk(LSDynaFamily::Float, 1) == 0)
    {
      time = p->Fam.GetNextWordAsFloat();
      if (time != LSDynaFamily::EOFMarker)
      {
        p->Fam.MarkTimeStep();
        p->TimeValues.push_back(time);
        // fprintf( stderr, "%d %f\n", (int) p->TimeValues.size() - 1, time ); fflush(stderr);
        if (p->Fam.SkipToWord(LSDynaFamily::TimeStepSection, ntimesteps++, p->Fam.GetStateSize()))
        {
          itmp = 0;
        }
      }
      else
      {
        if (p->Fam.AdvanceFile())
        {
          itmp = 0;
        }
        else
        {
          if (ntimesteps == 0)
          {
            // First time step was an EOF marker... move the marker to the
            // beginning of the first real time step...
            p->Fam.MarkSectionStart(lastAdapt, LSDynaFamily::TimeStepSection);
          }
        }
        int nextAdapt = p->Fam.GetCurrentAdaptLevel();
        if (nextAdapt != lastAdapt)
        {
          // Read the next static header section... state size has changed.
          p->Fam.MarkSectionStart(nextAdapt, LSDynaFamily::ControlSection);
          this->ReadHeaderInformation(nextAdapt);
          // p->Fam.SkipToWord( LSDynaFamily::EndOfStaticSection, nextAdapt, 0 );
          lastAdapt = nextAdapt;
        }
      }
    }
    else
    {
      itmp = 0;
    }
  } while (itmp);

  this->TimeStepRange[0] = 0;
  this->TimeStepRange[1] = ntimesteps ? ntimesteps - 1 : 0;

  return -1;
}

// =================================== Provide information about the database to the pipeline
int vtkLSDynaReader::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(iinfo), vtkInformationVector* oinfo)
{
  LSDynaMetaData* p = this->P;
  // If the time step is set before RequestInformation is called, we must
  // read the header information immediately in order to determine whether
  // the timestep that's been passed is valid. If it's not, we ignore it.
  if (!p->FileIsValid)
  {
    if (p->Fam.GetDatabaseDirectory().empty())
    {
      // fail silently for CanReadFile()'s sake.
      // vtkErrorMacro( "You haven't set the LS-Dyna database directory!" );
      return 1;
    }

    if (p->Fam.GetDatabaseBaseName().empty())
    {
      p->Fam.SetDatabaseBaseName("/d3plot"); // not a bad assumption.
    }
    p->Fam.ScanDatabaseDirectory();
    if (p->Fam.GetNumberOfFiles() < 1)
    {
      p->FileIsValid = 0;
      return 1;
    }
    p->Fam.DetermineStorageModel();
    p->MaxFileLength = p->FileSizeFactor * 512 * 512 * p->Fam.GetWordSize();
    p->FileIsValid = 1;

    // OK, now we have a list of files. Next, determine the length of the
    // state vector (#bytes of data stored per time step):
    this->ReadHeaderInformation(0);

    // Finally, we can loop through and determine where all the state
    // vectors start for each time step.
    // This will call ReadHeaderInformation when it encounters any
    // mesh adaptations (so that it can get the new state vector size).
    this->ScanDatabaseTimeSteps();
  }

  if (p->TimeValues.empty())
  {
    vtkErrorMacro("No valid time steps in the LS-Dyna database");
    return 0;
  }

  // Clamp timestep to be valid here.
  if (p->CurrentState < 0)
  {
    p->CurrentState = 0;
  }
  else if (p->CurrentState >= static_cast<vtkIdType>(p->TimeValues.size()))
  {
    p->CurrentState = static_cast<vtkIdType>(p->TimeValues.size() - 1);
  }

  int newAdaptLevel = p->Fam.TimeAdaptLevel(p->CurrentState);
  if (p->Fam.GetCurrentAdaptLevel() != newAdaptLevel)
  {
    // Requested time step has a different mesh adaptation than current one.
    // Update the header information so that queries like GetNumberOfCells() work properly.
    int result;
    result = this->ReadHeaderInformation(newAdaptLevel);
    if (result >= 0)
    {
      this->ResetPartsCache();
      return result;
    }
  }

  // Every output object has all the time steps.
  vtkInformation* outInfo = oinfo->GetInformationObject(0);
  outInfo->Set(
    vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &p->TimeValues[0], (int)p->TimeValues.size());
  double timeRange[2];
  timeRange[0] = p->TimeValues[0];
  timeRange[1] = p->TimeValues[p->TimeValues.size() - 1];
  outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);

  return 1;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadTopology()
{
  bool readTopology = false;
  if (!this->Parts)
  {
    readTopology = true;
    this->Parts = vtkLSDynaPartCollection::New();
    this->Parts->InitCollection(this->P, nullptr, nullptr);
  }
  if (!readTopology)
  {
    return 0;
  }

  if (this->ReadPartSizes())
  {
    vtkErrorMacro("Could not read cell sizes.");
    return 1;
  }

  if (this->ReadConnectivityAndMaterial())
  {
    vtkErrorMacro("Could not read connectivity.");
    return 1;
  }

  this->Parts->FinalizeTopology();

  if (this->ReadNodes())
  {
    vtkErrorMacro("Could not read static node values.");
    return 1;
  }

  // we need to read the user ids after we have read the topology
  // so we know how many cells are in each part
  if (this->ReadUserIds())
  {
    vtkErrorMacro("Could not read user node/element IDs.");
    return 1;
  }

  return 0;
}

// ============================================================= Section parsing
int vtkLSDynaReader::ReadNodes()
{
  LSDynaMetaData* p = this->P;

  // Always read geometry, even if the mesh will be deformed using deflected coordinates
  // (LS_ARRAYNAME_DEFLECTION). That way, we can compute deflection array later on.
  p->Fam.SkipToWord(LSDynaFamily::GeometryData, p->Fam.GetCurrentAdaptLevel(), 0);
  this->Parts->ReadPointProperty(p->NumberOfNodes, p->Dimensionality, nullptr, false, true, false);

  // Note that in any event we still have to read the rigid road coordinates.
  // If the mesh is deformed each state will have the points so see ReadState
  if (p->ReadRigidRoadMvmt)
  {
    vtkIdType nnode = p->Dict["NNODE"];
    p->Fam.SkipToWord(LSDynaFamily::RigidSurfaceData, p->Fam.GetCurrentAdaptLevel(), 4 + nnode);
    this->Parts->ReadPointProperty(nnode, 3, nullptr, false, false, true);
  }

  return 0;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadUserIds()
{
  // Below here is code that runs when user node or element numbers are present
  int arbitraryMaterials = this->P->Dict["NSORT"] < 0 ? 1 : 0;
  vtkIdType isz = this->GetNumberOfNodes();
  if (arbitraryMaterials)
  {
    this->P->Fam.SkipToWord(LSDynaFamily::UserIdData, this->P->Fam.GetCurrentAdaptLevel(), 16);
  }
  else
  {
    this->P->Fam.SkipToWord(LSDynaFamily::UserIdData, this->P->Fam.GetCurrentAdaptLevel(), 10);
  }

  // Node numbers
  bool nodeIdStatus = this->GetPointArrayStatus(LS_ARRAYNAME_USERID) == 1;
  if (nodeIdStatus)
  {
    this->Parts->ReadPointUserIds(isz, LS_ARRAYNAME_USERID);
  }

  // FIXME: This won't work if Rigid Body and Shell elements are interleaved (which I now believe
  // they are)
  this->Parts->ReadCellUserIds(
    LSDynaMetaData::SOLID, this->GetCellArrayStatus(LSDynaMetaData::SOLID, LS_ARRAYNAME_USERID));
  this->Parts->ReadCellUserIds(
    LSDynaMetaData::BEAM, this->GetCellArrayStatus(LSDynaMetaData::BEAM, LS_ARRAYNAME_USERID));
  this->Parts->ReadCellUserIds(
    LSDynaMetaData::SHELL, this->GetCellArrayStatus(LSDynaMetaData::SHELL, LS_ARRAYNAME_USERID));
  this->Parts->ReadCellUserIds(LSDynaMetaData::THICK_SHELL,
    this->GetCellArrayStatus(LSDynaMetaData::THICK_SHELL, LS_ARRAYNAME_USERID));
  this->Parts->ReadCellUserIds(LSDynaMetaData::RIGID_BODY,
    this->GetCellArrayStatus(LSDynaMetaData::RIGID_BODY, LS_ARRAYNAME_USERID));
  return 0;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadDeletion()
{
  enum LSDynaMetaData::LSDYNA_TYPES validCellTypes[4] = { LSDynaMetaData::SOLID,
    LSDynaMetaData::THICK_SHELL, LSDynaMetaData::SHELL, LSDynaMetaData::BEAM };

  if (!this->RemoveDeletedCells)
  {
    // this functions doesn't have to lead the reader to a certain
    // position in the files
    this->Parts->DisbleDeadCells();
    return 0;
  }

  LSDynaMetaData* p = this->P;
  vtkUnsignedCharArray* death;
  switch (p->Dict["MDLOPT"])
  {
    case LS_MDLOPT_POINT:
      vtkErrorMacro("We currently only support cell death");
      break;
    case LS_MDLOPT_CELL:
      for (int i = 0; i < 4; ++i)
      {
        const LSDynaMetaData::LSDYNA_TYPES type = validCellTypes[i];
        vtkIdType numCells, numSkipStart, numSkipEnd;
        this->Parts->GetPartReadInfo(type, numCells, numSkipStart, numSkipEnd);

        death = vtkUnsignedCharArray::New();
        death->SetName(LS_ARRAYNAME_DEATH);
        death->SetNumberOfComponents(1);
        death->SetNumberOfTuples(numCells);

        p->Fam.SkipWords(numSkipStart);
        this->ReadDeletionArray(death, 0, 1);
        p->Fam.SkipWords(numSkipEnd);
        this->Parts->SetCellDeadFlags(type, death, this->DeletedCellsAsGhostArray);
        death->Delete();
      }

      // we are now at the position to read the SPH deletion info from the sph state info
      if (p->NumberOfCells[LSDynaMetaData::PARTICLE] > 0)
      {
        const LSDynaMetaData::LSDYNA_TYPES type = LSDynaMetaData::PARTICLE;
        vtkIdType numCells, numSkipStart, numSkipEnd;
        this->Parts->GetPartReadInfo(type, numCells, numSkipStart, numSkipEnd);

        death = vtkUnsignedCharArray::New();
        death->SetName(LS_ARRAYNAME_DEATH);
        death->SetNumberOfComponents(1);
        death->SetNumberOfTuples(numCells);

        p->Fam.SkipWords(numSkipStart);
        // we are really reading the material id as the death flag
        // and each particle has twenty words of info, so we have to skip 19
        // since luckily material id is first
        this->ReadDeletionArray(death, 0, 20);
        p->Fam.SkipWords(numSkipEnd);
        this->Parts->SetCellDeadFlags(type, death, this->DeletedCellsAsGhostArray);
        death->Delete();
      }
      break;
    case LS_MDLOPT_NONE:
    default:
      // do nothing.
      break;
  }
  return 0;
}

//------------------------------------------------------------------------------
void vtkLSDynaReader::ReadDeletionArray(vtkUnsignedCharArray* arr, const int& pos, const int& size)
{
  // setup to do a block read, way faster than converting each
  // float/double individually
  LSDynaMetaData* p = this->P;
  vtkIdType startId = 0;
  vtkIdType numChunks = p->Fam.InitPartialChunkBuffering(arr->GetNumberOfTuples(), size);
  if (p->Fam.GetWordSize() == 8)
  {
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      vtkIdType chunkSize = p->Fam.GetNextChunk(LSDynaFamily::Float);
      vtkIdType numCellsInChunk = chunkSize / size;
      double* dbuf = p->Fam.GetBufferAs<double>();
      this->FillDeletionArray(dbuf, arr, startId, numCellsInChunk, pos, size);
      startId += numCellsInChunk;
    }
  }
  else
  {
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      vtkIdType chunkSize = p->Fam.GetNextChunk(LSDynaFamily::Float);
      vtkIdType numCellsInChunk = chunkSize / size;
      float* fbuf = p->Fam.GetBufferAs<float>();
      this->FillDeletionArray(fbuf, arr, startId, numCellsInChunk, pos, size);
      startId += numCellsInChunk;
    }
  }
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadState(vtkIdType step)
{
  // remember C style return so zero is pass
  if (this->ReadNodeStateInfo(step))
  {
    vtkErrorMacro("Problem reading state point information.");
    return 1;
  }
  if (this->ReadCellStateInfo(step))
  {
    vtkErrorMacro("Problem reading state cell information.");
    return 1;
  }
  if (this->ReadDeletion())
  {
    vtkErrorMacro("Problem reading state deletion information.");
    return 1;
  }
  return 0;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadNodeStateInfo(vtkIdType step)
{
  LSDynaMetaData* p = this->P;

  // Skip global variables for now
  p->Fam.SkipToWord(LSDynaFamily::TimeStepSection, step, 1 + p->Dict["NGLBV"]);

  // Read nodal data ===========================================================
  std::vector<std::string> names;
  std::vector<int> cmps;
  // Important: push_back in the order these are interleaved on disk
  // Note that temperature and deflection are swapped relative to the order they
  // are specified in the header section.
  const char* aNames[] = {
    LS_ARRAYNAME_DEFLECTION,
    LS_ARRAYNAME_TEMPERATURE,
    LS_ARRAYNAME_VELOCITY,
    LS_ARRAYNAME_ACCELERATION,
    LS_ARRAYNAME_PRESSURE,
    LS_ARRAYNAME_VORTICITY_X,
    LS_ARRAYNAME_VORTICITY_Y,
    LS_ARRAYNAME_VORTICITY_Z,
    LS_ARRAYNAME_RESULTANTVORTICITY,
    LS_ARRAYNAME_ENSTROPHY,
    LS_ARRAYNAME_HELICITY,
    LS_ARRAYNAME_STREAMFUNCTION,
    LS_ARRAYNAME_ENTHALPY,
    LS_ARRAYNAME_DENSITY,
    LS_ARRAYNAME_TURBULENTKE,
    LS_ARRAYNAME_DISSIPATION,
    LS_ARRAYNAME_EDDYVISCOSITY,
    LS_ARRAYNAME_SPECIES_01,
    LS_ARRAYNAME_SPECIES_02,
    LS_ARRAYNAME_SPECIES_03,
    LS_ARRAYNAME_SPECIES_04,
    LS_ARRAYNAME_SPECIES_05,
    LS_ARRAYNAME_SPECIES_06,
    LS_ARRAYNAME_SPECIES_07,
    LS_ARRAYNAME_SPECIES_08,
    LS_ARRAYNAME_SPECIES_09,
    LS_ARRAYNAME_SPECIES_10,
  };

  const char* aDictNames[] = {
    "IU",
    "IT",
    "IV",
    "IA",
    "cfdPressure",
    "cfdXVort",
    "cfdYVort",
    "cfdZVort",
    "cfdRVort",
    "cfdEnstrophy",
    "cfdHelicity",
    "cfdStream",
    "cfdEnthalpy",
    "cfdDensity",
    "cfdTurbKE",
    "cfdDiss",
    "cfdEddyVisc",
    "cfdSpec01",
    "cfdSpec02",
    "cfdSpec03",
    "cfdSpec04",
    "cfdSpec05",
    "cfdSpec06",
    "cfdSpec07",
    "cfdSpec08",
    "cfdSpec09",
    "cfdSpec10",
  };
  int aComponents[] = { -1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1 };
  int vppt = 0; // values per point
  int allVortPresent = p->Dict["cfdXVort"] && p->Dict["cfdYVort"] && p->Dict["cfdZVort"];

  for (int nvnum = 0; nvnum < (int)(sizeof(aComponents) / sizeof(aComponents[0])); ++nvnum)
  {
    if (p->Dict[aDictNames[nvnum]])
    {
      if (allVortPresent &&
        !strncmp(LS_ARRAYNAME_VORTICITY, aNames[nvnum], sizeof(LS_ARRAYNAME_VORTICITY)))
      {
        // turn the vorticity components from individual scalars into one vector (with a hack)
        if (nvnum < 7)
          continue;
        aComponents[nvnum] = 3;
        aNames[nvnum] = LS_ARRAYNAME_VORTICITY;
      }
      names.emplace_back(aNames[nvnum]);
      cmps.push_back(aComponents[nvnum] == -1 ? p->Dimensionality : aComponents[nvnum]);
      vppt += cmps.back();
    }
  }

  if (vppt != 0)
  {
    for (size_t i = 0; i < cmps.size(); i++)
    {
      // Note, we don't do anything special when reading deflected coordinates here.
      // See `ComputeDeflectionAndUpdateGeometry` for computing of deflection and
      // updating for geometry if requested to be deformed.
      bool valid = this->GetPointArrayStatus(names[i].c_str()) != 0;
      this->Parts->ReadPointProperty(p->NumberOfNodes, cmps[i], names[i].c_str(), valid);
    }

    // clear the buffer as it will be very large and not needed
    p->Fam.ClearBuffer();
  }
  return 0;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadCellStateInfo(vtkIdType vtkNotUsed(step))
{

  LSDynaMetaData* p = this->P;

  // ENN (Total element data for state) = Sum of
  //      NEL8  (# 8 node solid elems)           * NV3D
  //      NELT  (# 8 node thick shell elems)     * NV3DT
  //      NEL2  (# 2 node one-dimensional elems) * NV1D
  //      NEL4  (# 4 node shells elems)          * NV2D
  //      NMSPH (#of SPH Nodes)                  * NUM_SPH_VARS
  // where:
  //      NV3D = 7 + NEIPH  ( NEIPH = # ALE Data + # Strain Data)
  //      NV2D = MAXINT* (6*IOSHL(1) + 1*IOSHL(2) + NEIPS) +8*IOSHL(3) + 4*IOSHL(4) + 12*ISTRN
  //

  // Instead of repeating the specification of what array values are read for each cell type here
  // we can use the Property lists which we generated in ReadHeaderInformation().

  LSDynaMetaData::LSDYNA_TYPES celltypes[] = { LSDynaMetaData::SOLID, LSDynaMetaData::THICK_SHELL,
    LSDynaMetaData::BEAM, LSDynaMetaData::SHELL };
  vtkIdType cells[] = { p->Dict["NEL8"], p->Dict["NELT"], p->Dict["NEL2"], p->Dict["NEL4"] };
  vtkIdType cellVals[] = { p->Dict["NV3D"], p->Dict["NV3DT"], p->Dict["NV1D"], p->Dict["NV2D"] };

  // Be careful to exclude arrays which are not part State data
  unsigned int firstStateArrayNdx =
    (p->Dict["NARBS"] > 0) ? 1 : 0; // Skip first Array if it is UserIds

  for (int i = 0; i < 4; i++)
  {
    if (cells[i] > 0)
    {
      LSDynaMetaData::LSDYNA_TYPES celltype = celltypes[i];
      int startPos = 0;
      for (unsigned int a = firstStateArrayNdx; a < p->CellArrayNames[celltype].size(); a++)
      {
        int numComps = this->GetNumberOfComponentsInCellArray(celltype, a);
        // std::cout << setw(3) << numComps << " " << this->GetCellArrayName(celltype,a) <<
        // std::endl;
        if (this->GetCellArrayStatus(celltype, a))
          this->Parts->AddProperty(
            celltype, this->GetCellArrayName(celltype, a), startPos, numComps);
        startPos += numComps;
      }
      this->ReadCellProperties(celltype, cellVals[i]);
    }
  }

#undef VTK_LS_CELLARRAY
  return 0;
}

//------------------------------------------------------------------------------
void vtkLSDynaReader::ReadCellProperties(const int& type, const int& numTuples)
{
  LSDynaMetaData::LSDYNA_TYPES t = static_cast<LSDynaMetaData::LSDYNA_TYPES>(type);
  vtkIdType numCells, numSkipStart, numSkipEnd;
  this->Parts->GetPartReadInfo(type, numCells, numSkipStart, numSkipEnd);

  this->P->Fam.SkipWords(numSkipStart * numTuples);
  vtkIdType numChunks = this->P->Fam.InitPartialChunkBuffering(numCells, numTuples);
  vtkIdType startId = 0;
  if (this->P->Fam.GetWordSize() == 8 && numCells > 0)
  {
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      // we need offsets!
      vtkIdType chunkSize = this->P->Fam.GetNextChunk(LSDynaFamily::Float);
      vtkIdType numCellsInChunk = chunkSize / numTuples;
      double* dbuf = this->P->Fam.GetBufferAs<double>();
      this->Parts->FillCellProperties(dbuf, t, startId, numCellsInChunk, numTuples);
      startId += numCellsInChunk;
    }
  }
  else if (numCells > 0)
  {
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      vtkIdType chunkSize = this->P->Fam.GetNextChunk(LSDynaFamily::Float);
      vtkIdType numCellsInChunk = chunkSize / numTuples;
      float* fbuf = this->P->Fam.GetBufferAs<float>();
      this->Parts->FillCellProperties(fbuf, t, startId, numCellsInChunk, numTuples);
      startId += numCellsInChunk;
    }
  }
  this->P->Fam.SkipWords(numSkipEnd * numTuples);

  // clear the buffer as it will be very large and not needed
  this->P->Fam.ClearBuffer();
}

int vtkLSDynaReader::ReadSPHState(vtkIdType vtkNotUsed(step))
{
  LSDynaMetaData* p = this->P;

  // Make sure we are at the start of the SPH state data
  p->Fam.SkipToWord(LSDynaFamily::TimeStepSection, p->CurrentState, 0);
  p->Fam.SkipWords(p->SPHStateOffset);

#define VTK_LS_SPHARRAY(cond, celltype, arrayname, numComps)                                       \
  if ((cond) && this->GetCellArrayStatus(celltype, arrayname))                                     \
  {                                                                                                \
    this->Parts->AddProperty(celltype, arrayname, startPos, numComps);                             \
  }                                                                                                \
  if (cond)                                                                                        \
    startPos += (numComps);

  // Smooth Particle ========================================================

  // currently have a bug when reading SPH properties disabling for now
  int startPos = 0; // used to keep track of the startpos between calls to VTK_LS_CELLARRAY
  VTK_LS_SPHARRAY(
    true, LSDynaMetaData::PARTICLE, LS_ARRAYNAME_MATERIAL, 1); // always keep material id / death
  VTK_LS_SPHARRAY(
    p->Dict["isphfg(2)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_RADIUSOFINFLUENCE, 1);
  VTK_LS_SPHARRAY(p->Dict["isphfg(3)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_PRESSURE, 1);
  VTK_LS_SPHARRAY(p->Dict["isphfg(4)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_STRESS, 6);
  VTK_LS_SPHARRAY(p->Dict["isphfg(5)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_EPSTRAIN, 1);
  VTK_LS_SPHARRAY(p->Dict["isphfg(6)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_DENSITY, 1);
  VTK_LS_SPHARRAY(p->Dict["isphfg(7)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_INTERNALENERGY, 1);
  VTK_LS_SPHARRAY(p->Dict["isphfg(8)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_NUMNEIGHBORS, 1);
  VTK_LS_SPHARRAY(p->Dict["isphfg(9)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_STRAIN, 6);
  VTK_LS_SPHARRAY(p->Dict["isphfg(10)"], LSDynaMetaData::PARTICLE, LS_ARRAYNAME_MASS, 1);

  //  std::cout << "NUM_SPH_DATA: " << p->Dict["NUM_SPH_DATA"] << "start Pos is " << startPos <<
  //  std::endl;
  this->ReadCellProperties(LSDynaMetaData::PARTICLE, p->Dict["NUM_SPH_DATA"]);

#undef VTK_LS_SPHARRAY
  return 0;
}

int vtkLSDynaReader::ReadUserMaterialIds()
{
  LSDynaMetaData* p = this->P;
  vtkIdType m, numMats;

  p->MaterialsOrdered.clear();
  p->MaterialsUnordered.clear();
  p->MaterialsLookup.clear();
  // Does the file contain arbitrary material IDs?

  if ((p->Dict["NARBS"] > 0) && (p->Dict["NSORT"] < 0))
  { // Yes, it does. Read them.

    // Skip over arbitrary node and element IDs:
    vtkIdType skipIds =
      p->Dict["NUMNP"] + p->Dict["NEL8"] + p->Dict["NEL2"] + p->Dict["NEL4"] + p->Dict["NELT"];
    p->Fam.SkipToWord(LSDynaFamily::UserIdData, p->Fam.GetCurrentAdaptLevel(), 16 + skipIds);

    // in some cases the number of materials in NMMAT is incorrect since we are loading
    // SPH materials.
    numMats = p->Dict["NMMAT"];

    // Read in material ID lists:
    p->Fam.BufferChunk(LSDynaFamily::Int, numMats * 3);
    for (m = 0; m < numMats; ++m)
    {
      p->MaterialsOrdered.push_back(p->Fam.GetNextWordAsInt());
    }
    for (m = 0; m < numMats; ++m)
    {
      p->MaterialsUnordered.push_back(p->Fam.GetNextWordAsInt());
    }
    for (m = 0; m < numMats; ++m)
    {
      p->MaterialsLookup.push_back(p->Fam.GetNextWordAsInt());
    }
  }
  else
  {
    numMats = p->Dict["NUMMAT8"] + p->Dict["NUMMATT"] + p->Dict["NUMMAT4"] + p->Dict["NUMMAT2"] +
      p->Dict["NGPSPH"];
    // No, it doesn't. Fabricate a list of sequential IDs
    // construct the (trivial) material lookup tables
    for (m = 1; m <= numMats; ++m)
    {
      p->MaterialsOrdered.push_back(m);
      p->MaterialsUnordered.push_back(m);
      p->MaterialsLookup.push_back(m);
    }
  }
  return 0;
}

int vtkLSDynaReader::ReadPartTitlesFromRootFile()
{
  /*
  The extra data is written at the end of the following files:
  d3plot, d3part and intfor files, and the header and part titles are written directly after the
  EOF (= -999999.0) marker.
  Value Length Description
  -------------------------------

  NTYPE 1 entity type = 90001
  NUMPROP 1 number of parts

  For NUMPROP parts:
  IDP 1 part id
  PTITLE 18 Part title (72 characters)

  For NUMPROP parts:
  NTYPE 1 entity type = 90000
  HEAD 18 Header title (72 characters)
  */

  LSDynaMetaData* p = this->P;
  if (p->PreStateSize <= 0)
  {
    vtkErrorMacro("Database has bad pre state size(" << p->PreStateSize << ").");
    return 1;
  }

  // when called this method is at the right spot to read the part names
  vtkIdType currentFileLoc = p->Fam.GetCurrentFWord();
  vtkIdType currentAdaptLevel = p->Fam.GetCurrentAdaptLevel();

  p->Fam.BufferChunk(LSDynaFamily::Float, 1);
  double eofM = p->Fam.GetNextWordAsFloat();
  if (eofM != LSDynaFamily::EOFMarker)
  {
    // we failed to find a marker stop on the part names
    p->Fam.SkipToWord(LSDynaFamily::ControlSection, currentAdaptLevel, currentFileLoc);
    return 1;
  }

  // make sure that the root files has room left for the amount of data we are going to request
  // if it doesn't we know it can't have part names
  vtkIdType numParts = static_cast<vtkIdType>(p->PartIds.size());
  vtkIdType partTitlesByteSize =
    p->Fam.GetWordSize() * (2 + numParts); // NType + NUMPRop + (header part ids)
  partTitlesByteSize +=
    (numParts * 72); // names are constant at 72 bytes each independent of word size

  vtkIdType fileSize = p->Fam.GetFileSize(0);
  if (fileSize < partTitlesByteSize + p->Fam.GetCurrentFWord())
  {
    // this root file doesn't part names
    p->Fam.SkipToWord(LSDynaFamily::ControlSection, currentAdaptLevel, currentFileLoc);
    return 1;
  }

  // we can now safely read the part titles
  p->Fam.SkipWords(2); // skip types and num of parts
  vtkIdType nameWordSize = 72 / p->Fam.GetWordSize();
  for (vtkIdType i = 0; i < numParts; ++i)
  {
    p->Fam.BufferChunk(LSDynaFamily::Int, 1);
    p->Fam.GetNextWordAsInt(); // vtkIdType partId

    p->Fam.BufferChunk(LSDynaFamily::Char, nameWordSize);
    std::string name(p->Fam.GetNextWordAsChars(), 72);
    if (!name.empty() && name[0] != ' ')
    {
      // strip the name to the subset that
      size_t found = name.find_last_not_of(' ');
      if (found != std::string::npos)
      {
        name = name.substr(0, found + 1);
      }
      // get the right part id
      p->PartNames[i] = name;
    }
  }
  p->Fam.SkipToWord(LSDynaFamily::ControlSection, currentAdaptLevel, currentFileLoc);
  return 0;
}

void vtkLSDynaReader::ResetPartInfo()
{
  LSDynaMetaData* p = this->P;
  p->PartNames.clear();
  p->PartIds.clear();
  p->PartMaterials.clear();
  p->PartStatus.clear();

  // Create simple part names as place holders
  int mat = 1, realMat = 1;
  int i;
  int N;
  char partLabel[64];
  int arbitraryMaterials = p->Dict["NSORT"] < 0 ? 1 : 0;

#define VTK_LSDYNA_PARTLABEL(dict, fmt)                                                            \
  N = p->Dict[dict];                                                                               \
  for (i = 0; i < N; ++i, ++mat)                                                                   \
  {                                                                                                \
    if (arbitraryMaterials)                                                                        \
    {                                                                                              \
      if (mat < static_cast<int>(p->MaterialsOrdered.size()))                                      \
      {                                                                                            \
        realMat = p->MaterialsOrdered[mat - 1];                                                    \
      }                                                                                            \
      else                                                                                         \
      {                                                                                            \
        realMat = mat;                                                                             \
      }                                                                                            \
      snprintf(partLabel, sizeof(partLabel), fmt " (Matl%d)", mat, realMat);                       \
    }                                                                                              \
    else                                                                                           \
    {                                                                                              \
      realMat = mat;                                                                               \
      snprintf(partLabel, sizeof(partLabel), fmt, mat);                                            \
    }                                                                                              \
    p->PartNames.emplace_back(partLabel);                                                          \
    p->PartIds.emplace_back(realMat);                                                              \
    p->PartMaterials.emplace_back(mat);                                                            \
    p->PartStatus.emplace_back(1);                                                                 \
  }

  VTK_LSDYNA_PARTLABEL("NUMMAT8", "Part%d"); // was "PartSolid%d
  VTK_LSDYNA_PARTLABEL("NUMMATT", "Part%d"); // was "PartThickShell%d
  VTK_LSDYNA_PARTLABEL("NUMMAT4", "Part%d"); // was "PartShell%d
  VTK_LSDYNA_PARTLABEL("NUMMAT2", "Part%d"); // was "PartBeam%d
  VTK_LSDYNA_PARTLABEL("NGPSPH", "Part%d");  // was "PartParticle%d
  VTK_LSDYNA_PARTLABEL("NSURF", "Part%d");   // was "PartRoadSurface%d
  VTK_LSDYNA_PARTLABEL("NUMMAT", "Part%d");  // was "PartRigidBody%d

#undef VTK_LSDYNA_PARTLABEL
}

int vtkLSDynaReader::ReadInputDeck()
{
  if (!this->InputDeck)
  {
    // nothing more we can do
    return 0;
  }

  vtksys::ifstream deck(this->InputDeck, ios::in);
  if (!deck.good())
  {
    return 0;
  }

  std::string header;
  std::getline(deck, header, '\n');
  deck.seekg(0, ios::beg);
  int retval;
  if (vtksys::SystemTools::StringStartsWith(header.c_str(), "<?xml"))
  {
    retval = this->ReadInputDeckXML(deck);
  }
  else
  {
    retval = this->ReadInputDeckKeywords(deck);
  }

  return retval;
}

int vtkLSDynaReader::ReadInputDeckXML(istream& deck)
{
  vtkLSDynaSummaryParser* parser = vtkLSDynaSummaryParser::New();
  parser->MetaData = this->P;
  parser->SetStream(&deck);
  // We must be able to parse the file and end up with 1 part per material ID
  if (!parser->Parse() || this->P->GetTotalMaterialCount() != (vtkIdType)this->P->PartNames.size())
  {
    // We had a problem identifying a part, give up and start over by resetting the parts
    this->ResetPartInfo();
  }
  parser->Delete();

  return 0;
}

int vtkLSDynaReader::ReadInputDeckKeywords(istream& deck)
{
  int success = 1;
  std::map<std::string, int> parameters;
  std::string line;
  std::string lineLowercase;
  std::string partName;
  int partMaterial;
  int partId;
  int curPart = 0;

  while (
    deck.good() && vtkLSNextSignificantLine(deck, line) && curPart < (int)this->P->PartNames.size())
  {
    if (line[0] == '*')
    {
      vtkLSDowncaseFirstWord(lineLowercase, line.substr(1));
      if (vtksys::SystemTools::StringStartsWith(lineLowercase.c_str(), "part"))
      {
        // found a part
        // ... read the next non-comment line as the part name
        if (vtkLSNextSignificantLine(deck, line))
        {
          // Get rid of leading and trailing newlines, whitespace, etc.
          vtkLSTrimWhitespace(line);
          partName = line;
        }
        else
        {
          partName = "";
        }
        // ... read the next non-comment line as the part id or a reference to it.
        if (vtkLSNextSignificantLine(deck, line))
        {
          std::vector<std::string> splits;
          vtkLSSplitString(line, splits, "& ,\t\n\r");
          if (line[0] == '&')
          {
            // found a reference. look it up.
            partId = !splits.empty() ? parameters[splits[0]] : -1;
          }
          else
          {
            if (splits.empty() || sscanf(splits[0].c_str(), "%d", &partId) <= 0)
            {
              partId = -1;
            }
          }
          if (splits.size() < 3)
          {
            partMaterial = -1;
          }
          else
          {
            if (splits[2][0] == '&')
            {
              partMaterial = parameters[splits[2]];
            }
            else
            {
              if (sscanf(splits[2].c_str(), "%d", &partMaterial) <= 0)
              {
                partMaterial = -1;
              }
            }
          }
        } // read the part id or reference
        else
        {
          partId = -1;
          partMaterial = -1;
        }
        // Comment on next line: partId values need not be consecutive. FIXME: ... or even positive?
        if (!partName.empty() && partId >= 0)
        {
          this->P->PartNames[curPart] = partName;
          this->P->PartIds[curPart] = partId;
          this->P->PartMaterials[curPart] = partMaterial;
          this->P->PartStatus[curPart] = 1;
          fprintf(stderr, "%2d: Part: \"%s\" Id: %d\n", curPart, partName.c_str(), partId);
          ++curPart;
        }
        else
        {
          success = 0;
        }
      }
      else if (vtksys::SystemTools::StringStartsWith(lineLowercase.c_str(), "parameter"))
      {
        // found a reference
        // ... read the next non-comment line to decode the reference
        if (vtkLSNextSignificantLine(deck, line))
        {
          std::string paramName;
          int paramIntVal;
          // Look for "^[IiRr]\s*(\w+)\s+([\w\.-]+)" and set parameters[\2]=\1
          if (line[0] == 'I' || line[0] == 'i')
          { // We found an integer parameter. Those are the only ones we care about.
            line = line.substr(1);
            std::string::size_type paramStart = line.find_first_not_of(" \t,");
            if (paramStart == std::string::npos)
            { // ignore a bad parameter line
              continue;
            }
            std::string::size_type paramEnd = line.find_first_of(" \t,", paramStart);
            if (paramEnd == std::string::npos)
            { // found the parameter name, but no value after it
              continue;
            }
            paramName = line.substr(paramStart, paramEnd - paramStart);
            if (sscanf(line.substr(paramEnd + 1).c_str(), "%d", &paramIntVal) <= 0)
            { // unable to read id
              continue;
            }
            parameters[paramName] = paramIntVal;
          }
        }
        else
        {
          // no valid line after "*parameter" keyword. Silently ignore it.
        }
      } // "parameter line"
    }   // line starts with "*"
  }     // while ( deck.good() )

  if (success)
  {
    // Save a summary file if possible. The user can open the summary file next
    // time and not be forced to parse the entire input deck to get part IDs.
    std::string deckDir = vtksys::SystemTools::GetFilenamePath(this->InputDeck);
    std::string deckName = vtksys::SystemTools::GetFilenameName(this->InputDeck);
    std::string deckExt;
    std::string::size_type dot;
    std::string xmlSummary;

    // GetFilenameExtension doesn't look for the rightmost "." ... do it ourselves.
    dot = deckName.rfind('.');
    if (dot != std::string::npos)
    {
      deckExt = deckName.substr(dot);
      deckName = deckName.substr(0, dot);
    }
    else
    {
      deckExt = "";
    }
#ifndef _WIN32
    xmlSummary = deckDir + "/" + deckName + ".lsdyna";
#else
    xmlSummary = deckDir + "\\" + deckName + ".lsdyna";
#endif // _WIN32
    // As long as we don't kill the input deck, write the summary XML:
    if (xmlSummary != this->InputDeck)
    {
      this->WriteInputDeckSummary(xmlSummary.c_str());
    }
  }
  else
  {
    // We had a problem identifying a part, give up and reset part info
    this->ResetPartInfo();
  }

  return !success;
}

int vtkLSDynaReader::WriteInputDeckSummary(const char* fname)
{
  vtksys::ofstream xmlSummary(fname, ios::out | ios::trunc);
  if (!xmlSummary.good())
  {
    return 1;
  }

  xmlSummary << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl << "<lsdyna>" << endl;

  std::string dbDir = this->P->Fam.GetDatabaseDirectory();
  std::string dbName = this->P->Fam.GetDatabaseBaseName();
  if (this->IsDatabaseValid() && !dbDir.empty() && !dbName.empty())
  {
#ifndef _WIN32
    if (dbDir[0] == '/')
#else
    if (dbDir[0] == '\\')
#endif // _WIN32
    {
      // OK, we have an absolute path, so it should be safe to write it out.
      xmlSummary << "  <database path=\"" << dbDir.c_str() << "\" name=\"" << dbName.c_str()
                 << "\"/>" << endl;
    }
  }

  for (unsigned p = 0; p < this->P->PartNames.size(); ++p)
  {
    xmlSummary << "  <part id=\"" << this->P->PartIds[p] << "\" material_id=\""
               << this->P->PartMaterials[p] << "\" status=\"" << this->P->PartStatus[p]
               << "\"><name>" << this->P->PartNames[p].c_str() << "</name></part>" << endl;
  }

  xmlSummary << "</lsdyna>" << endl;

  return 0;
}

// ================================================== OK Already! Read the file!
int vtkLSDynaReader::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(iinfo), vtkInformationVector* oinfo)
{
  LSDynaMetaData* p = this->P;

  if (!p->FileIsValid)
  {
    // This should have been set in RequestInformation()
    return 0;
  }
  p->Fam.ClearBuffer();
  p->Fam.OpenFileHandles();

  vtkMultiBlockDataSet* mbds = nullptr;
  vtkInformation* oi = oinfo->GetInformationObject(0);
  if (!oi)
  {
    return 0;
  }

  if (oi->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
  {
    // Only return single time steps for now.
    double requestedTimeStep = oi->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
    int timeStepLen = oi->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    double* timeSteps = oi->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());

    int cnt = 0;
    while (cnt < timeStepLen - 1 && timeSteps[cnt] < requestedTimeStep)
    {
      ++cnt;
    }
    this->SetTimeStep(cnt);

    oi->Set(vtkDataObject::DATA_TIME_STEP(), p->TimeValues[p->CurrentState]);
  }

  mbds = vtkMultiBlockDataSet::SafeDownCast(oi->Get(vtkDataObject::DATA_OBJECT()));
  if (!mbds)
  {
    return 0;
  }
  this->UpdateProgress(0.01);

  if (p->Dict["MATTYP"])
  {
    // Do something with material type data
  }
  this->UpdateProgress(0.05);

  if (p->Dict["IALEMAT"])
  {
    // Do something with fluid material ID data
  }
  this->UpdateProgress(0.10);

  if (p->Dict["NMSPH"])
  {
    // Do something with smooth particle hydrodynamics element data
  }
  this->UpdateProgress(0.15);

  // Read in the topology information for caching
  this->ReadTopology();

  // Adapted element parent list
  // This isn't even implemented by LS-Dyna yet

  // Smooth Particle Hydrodynamics Node and Material List are handled in
  // ReadConnectivityAndMaterial()

  // Start of state data ===================
  // I. Node and Cell State
  this->UpdateProgress(0.6);
  if (this->ReadState(p->CurrentState))
  {
    vtkErrorMacro("Problem reading state data for time step " << p->CurrentState);
    return 1;
  }

  // III. SPH Node State
  this->UpdateProgress(0.7);
  if (this->GetNumberOfParticleCells())
  {
    if (this->ReadSPHState(p->CurrentState))
    {
      vtkErrorMacro("Problem reading smooth particle hydrodynamics state.");
      return 1;
    }
  }

  this->UpdateProgress(0.8);
  // add all the parts as child blocks to the output
  int size = this->Parts->GetNumberOfParts();
  for (int i = 0; i < size; ++i)
  {
    if (this->Parts->IsActivePart(i))
    {
      vtkUnstructuredGrid* ug = this->Parts->GetGridForPart(i);
      this->ComputeDeflectionAndUpdateGeometry(ug);

      mbds->SetBlock(i, ug);
      mbds->GetMetaData(i)->Set(vtkCompositeDataSet::NAME(), this->P->PartNames[i].c_str());
    }
    else
    {
      mbds->SetBlock(i, nullptr);
    }
  }

  this->P->Fam.ClearBuffer();
  this->UpdateProgress(1.0);
  return 1;
}

//------------------------------------------------------------------------------
template <typename T>
void vtkLSDynaReader::FillDeletionArray(T* buffer, vtkUnsignedCharArray* arr,
  const vtkIdType& start, const vtkIdType& numCells, const int& deathPos, const int& cellSize)
{
  unsigned char val;
  for (vtkIdType i = 0; i < numCells; ++i)
  {
    // Quote from LSDyna Manual:
    //"each value is set to the element material number or =0,
    // if the element is deleted"
    val = (buffer[deathPos] == 0.0) ? 1 : 0;
    buffer += cellSize;
    arr->SetTuple1(start + i, val);
  }
}

//------------------------------------------------------------------------------
template <int wordSize, typename T>
int vtkLSDynaReader::FillTopology(T* buff)
{
  // the passed in buffer is null, and only used to specialze the method
  // as pure method specialization isn't support by some compilers
  // READ PARTICLES
  this->P->Fam.SkipToWord(LSDynaFamily::SPHNodeData, this->P->Fam.GetCurrentAdaptLevel(), 0);
  FillBlock<LSDynaMetaData::PARTICLE, wordSize, 1>(buff, this->Parts, this->P, 2, VTK_VERTEX);

  // READ SOLIDS
  this->P->Fam.SkipToWord(LSDynaFamily::GeometryData, this->P->Fam.GetCurrentAdaptLevel(),
    this->P->NumberOfNodes * this->P->Dimensionality);

  // other than buff, these parameters are changed by the template specialization
  // as SOLID is a unique case
  FillBlock<LSDynaMetaData::SOLID, wordSize, 8>(buff, this->Parts, this->P, 9, VTK_HEXAHEDRON);

  // READ THICK_SHELL
  FillBlock<LSDynaMetaData::THICK_SHELL, wordSize, 8>(
    buff, this->Parts, this->P, 9, VTK_QUADRATIC_QUAD);

  // READ BEAM
  FillBlock<LSDynaMetaData::BEAM, wordSize, 2>(buff, this->Parts, this->P, 6, VTK_LINE);

  // READ SHELL and RIGID_BODY
  // uses a specialization to weave SHELL and RIGID BODY cells together
  FillBlock<LSDynaMetaData::SHELL, wordSize, 4>(buff, this->Parts, this->P, 5, VTK_QUAD);

  // Read Road Surface
  if (this->P->ReadRigidRoadMvmt)
  {
    this->P->Fam.SkipToWord(LSDynaFamily::RigidSurfaceData, this->P->Fam.GetCurrentAdaptLevel(),
      4 + 4 * this->P->Dict["NNODE"]);
    FillBlock<LSDynaMetaData::ROAD_SURFACE, wordSize, 4>(buff, this->Parts, this->P, 5, VTK_QUAD);
  }
  return 0;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadConnectivityAndMaterial()
{
  LSDynaMetaData* p = this->P;
  if (p->ConnectivityUnpacked == 0)
  {
    // FIXME
    vtkErrorMacro("Packed connectivity isn't supported yet.");
    return 1;
  }

  this->Parts->InitCellInsertion();
  if (p->Fam.GetWordSize() == 8)
  {
    vtkIdType* buf = nullptr;
    return this->FillTopology<8>(buf);
  }
  else
  {
    int* buf = nullptr;
    return this->FillTopology<4>(buf);
  }
}

//------------------------------------------------------------------------------
template <typename T, int blockType, vtkIdType numWordsPerCell, vtkIdType cellLength>
void vtkLSDynaReader::ReadBlockCellSizes()
{
  // determine the relationship between the file bit size and
  // the host machine bit size. This allows us to read 64 bit files on a
  // 32 bit machine
  const int numWordsPerIdType(this->P->Fam.GetWordSize() / sizeof(T));

  vtkIdType nc = 0, t = 0, j = 0, matlId = 0;
  vtkIdType numCellsToSkip = 0, numCellsToSkipEnd = 0, chunkSize = 0;
  const T fileNumWordsPerCell(numWordsPerCell * numWordsPerIdType);
  const T offsetToMatId(numWordsPerIdType * (numWordsPerCell - 1));
  T* buff = nullptr;

  // get from the part the read information for this lsdyna block type
  this->Parts->GetPartReadInfo(blockType, nc, numCellsToSkip, numCellsToSkipEnd);

  this->P->Fam.SkipWords(fileNumWordsPerCell * numCellsToSkip); // skip to the right start id

  // buffer the amount in small chunks so we don't create a massive buffer
  vtkIdType numChunks = this->P->Fam.InitPartialChunkBuffering(nc, numWordsPerCell);
  for (vtkIdType i = 0; i < numChunks; ++i)
  {
    chunkSize = this->P->Fam.GetNextChunk(LSDynaFamily::Int);
    buff = this->P->Fam.GetBufferAs<T>();

    for (j = 0; j < chunkSize; j += numWordsPerCell)
    {
      buff += offsetToMatId;
      matlId = static_cast<vtkIdType>(*buff);
      buff += numWordsPerIdType;
      this->Parts->RegisterCellIndexToPart(blockType, matlId, t++, cellLength);
    }
  }
  this->P->Fam.SkipWords(fileNumWordsPerCell * numCellsToSkipEnd);
}

//------------------------------------------------------------------------------
template <typename T>
int vtkLSDynaReader::FillPartSizes()
{
  // READ PARTICLES
  this->P->Fam.SkipToWord(LSDynaFamily::SPHNodeData, this->P->Fam.GetCurrentAdaptLevel(), 0);
  this->ReadBlockCellSizes<T, LSDynaMetaData::PARTICLE, 2, 1>();

  // READ SOLIDS
  this->P->Fam.SkipToWord(LSDynaFamily::GeometryData, this->P->Fam.GetCurrentAdaptLevel(),
    this->P->NumberOfNodes * this->P->Dimensionality);

  this->ReadBlockCellSizes<T, LSDynaMetaData::SOLID, 9, 8>();

  // READ THICK_SHELL
  this->ReadBlockCellSizes<T, LSDynaMetaData::THICK_SHELL, 9, 8>();

  // READ BEAM
  this->ReadBlockCellSizes<T, LSDynaMetaData::BEAM, 6, 2>();

  // READ SHELL and RIGID_BODY
  // uses a specialization to weave SHELL and RIGID BODY cells together
  this->ReadBlockCellSizes<T, LSDynaMetaData::SHELL, 5, 4>();

  // Read Road Surface
  if (this->P->ReadRigidRoadMvmt)
  {
    this->P->Fam.SkipToWord(LSDynaFamily::RigidSurfaceData, this->P->Fam.GetCurrentAdaptLevel(),
      4 + 4 * this->P->Dict["NNODE"]);
    this->ReadBlockCellSizes<T, LSDynaMetaData::ROAD_SURFACE, 5, 4>();
  }

  // now that all the registering is done tell the collection
  // it can allocate the necessary space for each part
  this->Parts->AllocateParts();
  return 0;
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ReadPartSizes()
{
  LSDynaMetaData* p = this->P;
  if (p->ConnectivityUnpacked == 0)
  {
    // FIXME
    vtkErrorMacro("Packed connectivity isn't supported yet.");
    return 1;
  }

  if (p->Fam.GetWordSize() == 8)
  {
    return this->FillPartSizes<vtkIdType>();
  }
  else
  {
    return this->FillPartSizes<int>();
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaReader::SetDeformedMesh(vtkTypeBool deformed)
{
  if (this->DeformedMesh != deformed)
  {
    this->DeformedMesh = deformed;
    this->ResetPartsCache();
    this->Modified();
  }
}

namespace
{
template <class T, int NumberOfComponents>
vtkSmartPointer<T> vtkComputeDifference(T* aArray, T* bArray)
{
  if (aArray == nullptr || bArray == nullptr)
  {
    return nullptr;
  }

  const vtkIdType numTuples = aArray->GetNumberOfTuples();
  const int numComps = aArray->GetNumberOfComponents();
  if (bArray->GetNumberOfTuples() != numTuples || bArray->GetNumberOfComponents() != numComps ||
    numComps != NumberOfComponents)
  {
    return nullptr;
  }

  vtkVector<typename T::ValueType, NumberOfComponents> tupleA, tupleB;
  vtkSmartPointer<T> result = vtkSmartPointer<T>::New();
  result->SetNumberOfComponents(NumberOfComponents);
  result->SetNumberOfTuples(numTuples);
  for (vtkIdType cc = 0, max = numTuples; cc < max; ++cc)
  {
    aArray->GetTypedTuple(cc, tupleA.GetData());
    bArray->GetTypedTuple(cc, tupleB.GetData());
    result->SetTypedTuple(cc, (tupleA - tupleB).GetData());
  }
  return result;
}
}

//------------------------------------------------------------------------------
int vtkLSDynaReader::ComputeDeflectionAndUpdateGeometry(vtkUnstructuredGrid* ug)
{
  // If LS_ARRAYNAME_DEFLECTION is preset then this computes the deflection.
  // If this->DeformedMesh is true, this additionally swaps the output geometry points
  // to be the LS_ARRAYNAME_DEFLECTION (which is the deflected coordinates).
  LSDynaMetaData* p = this->P;
  vtkDataArray* deflectedCoords =
    ug ? ug->GetPointData()->GetArray(LS_ARRAYNAME_DEFLECTION) : nullptr;
  if (deflectedCoords)
  {
    vtkSmartPointer<vtkDataArray> deflection;
    if (p->Fam.GetWordSize() == 8)
    {
      deflection =
        vtkComputeDifference<vtkDoubleArray, 3>(vtkDoubleArray::SafeDownCast(deflectedCoords),
          vtkDoubleArray::SafeDownCast(ug->GetPoints()->GetData()));
    }
    else
    {
      deflection =
        vtkComputeDifference<vtkFloatArray, 3>(vtkFloatArray::SafeDownCast(deflectedCoords),
          vtkFloatArray::SafeDownCast(ug->GetPoints()->GetData()));
    }

    if (deflection)
    {
      deflection->SetName("Deflection");
      ug->GetPointData()->AddArray(deflection);
    }

    if (this->DeformedMesh)
    {
      ug->GetPoints()->SetData(deflectedCoords);
    }
  }
  return EXIT_SUCCESS;
}