// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause

#include "vtkLSDynaPartCollection.h"
#include "LSDynaMetaData.h"
#include "vtkLSDynaPart.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkIdTypeArray.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkStringArray.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"

#include <algorithm>
#include <list>
#include <vector>

//------------------------------------------------------------------------------
VTK_ABI_NAMESPACE_BEGIN
class vtkLSDynaPartCollection::LSDynaPartStorage
{
protected:
  //---------------------------------------------------------------------------
  // stores the number of cells for a given part
  // this struct is meant to resemble a run length encoding style of storage
  // of mapping cell ids to the part that holds those cells
  struct PartInfo
  {
    PartInfo(vtkLSDynaPart* p, const int& type, const vtkIdType& pId, const vtkIdType& start,
      const vtkIdType& npts)
      : numCells(1)
      , // we are inserting the first cell when we create this so start with 1
      startId(start)
      , cellStructureSize(npts)
      , // start with number of points in first cell
      partId(pId)
    {
      // we store the part id our selves because we can have null parts
      // if the user has disabled reading that part
      this->part = p;
      if (this->part)
      {
        this->part->SetPartType(type);
      }
    }

    vtkIdType numCells;          // number of cells in this continuous block
    vtkIdType startId;           // the global index to start of this block
    vtkIdType cellStructureSize; // stores the size of the cell array for this section
    vtkIdType partId; // id of the part this block represents, because the part can be nullptr
    vtkLSDynaPart* part;
  };

  //---------------------------------------------------------------------------
  struct PartInsertion
  {
    PartInsertion()
      : numCellsInserted(0)
    {
    }

    PartInsertion(std::vector<PartInfo>* pInfo)
      : numCellsInserted(0)
    {
      this->pIt = pInfo->begin();
    }

    // increments the numCells, and when needed increments to the next part
    void inc()
    {
      ++numCellsInserted;
      if ((*pIt).numCells == numCellsInserted)
      {
        ++pIt;
        numCellsInserted = 0;
      }
    }

    std::vector<PartInfo>::iterator pIt;
    vtkIdType numCellsInserted;
  };
  //---------------------------------------------------------------------------

public:
  LSDynaPartStorage(const vtkIdType& numMaterials)
    : NumParts(numMaterials)
    , PartIteratorLoc(0)
  {
    // a part represents a single material. A part type is
    this->Info = new std::vector<PartInfo>[LSDynaMetaData::NUM_CELL_TYPES];
    this->CellInsertionIterators = new PartInsertion[LSDynaMetaData::NUM_CELL_TYPES];
    this->Parts = new vtkLSDynaPart*[numMaterials];
    for (vtkIdType i = 0; i < numMaterials; ++i)
    {
      this->Parts[i] = nullptr;
    }
  }
  ~LSDynaPartStorage()
  {
    for (vtkIdType i = 0; i < this->NumParts; ++i)
    {
      if (this->Parts[i])
      {
        this->Parts[i]->Delete();
        this->Parts[i] = nullptr;
      }
    }
    delete[] this->Parts;
    delete[] this->CellInsertionIterators;
    delete[] this->Info;
  }

  //---------------------------------------------------------------------------
  vtkIdType GetNumParts() const { return NumParts; }

  //---------------------------------------------------------------------------
  void RegisterCell(const int& partType, const vtkIdType& matId, const vtkIdType& npts)
  {
    if (!this->Info[partType].empty())
    {
      PartInfo* info = &this->Info[partType].back();
      if (info->partId == matId)
      {
        // append to this item
        ++info->numCells;
        info->cellStructureSize += npts;
      }
      else
      {
        // add a new item
        // PartInfo sets the part type!
        PartInfo newInfo(
          this->Parts[matId], partType, matId, (info->startId + info->numCells), npts);
        this->Info[partType].push_back(newInfo);
      }
    }
    else
    {
      // PartInfo sets the part type!
      PartInfo newInfo(this->Parts[matId], partType, matId, 0, npts);
      this->Info[partType].push_back(newInfo);
    }
  }

  //---------------------------------------------------------------------------
  void ConstructPart(const vtkIdType& index, const std::string& name, const int& materialId,
    const int& numGlobalNodes, const int& wordSize)
  {
    vtkLSDynaPart* p = vtkLSDynaPart::New();
    p->InitPart(name, index, materialId, numGlobalNodes, wordSize);
    this->Parts[index] = p;
  }

  //---------------------------------------------------------------------------
  void InitCellInsertion()
  {
    // we build up an array of cell insertion iterators
    // that point to the first element of each part type info
    for (int i = 0; i < LSDynaMetaData::NUM_CELL_TYPES; ++i)
    {
      if (!this->Info[i].empty())
      {
        PartInsertion partIt(&this->Info[i]);
        this->CellInsertionIterators[i] = partIt;
      }
    }
  }

  //---------------------------------------------------------------------------
  void InsertCell(
    const int& partType, const int& cellType, const vtkIdType& npts, vtkIdType conn[8])
  {
    // get the correct iterator from the array of iterations
    if (this->CellInsertionIterators[partType].pIt->part)
    {
      // only insert the cell if the part is turned on
      this->CellInsertionIterators[partType].pIt->part->AddCell(cellType, npts, conn);
    }
    this->CellInsertionIterators[partType].inc();
  }

  //---------------------------------------------------------------------------
  bool PartExists(const vtkIdType& index) const
  {
    if (index < 0 || index > this->NumParts)
    {
      return false;
    }
    return (this->Parts[index] != nullptr && this->Parts[index]->HasCells());
  }

  //---------------------------------------------------------------------------
  vtkLSDynaPart* GetPart(const vtkIdType& index) { return this->Parts[index]; }

  //---------------------------------------------------------------------------
  vtkUnstructuredGrid* GetPartGrid(const vtkIdType& index)
  {
    return this->Parts[index]->GenerateGrid();
  }

  //---------------------------------------------------------------------------
  void InitPartIteration(const int& partType)
  {
    for (vtkIdType i = 0; i < this->NumParts; ++i)
    {
      if (this->Parts[i] && this->Parts[i]->PartType() == partType)
      {
        PartIteratorLoc = i;
        this->PartIterator = this->Parts[i];
        return;
      }
    }
    // failed to find a part that matches the type
    PartIteratorLoc = -1;
    this->PartIterator = nullptr;
  }

  //---------------------------------------------------------------------------
  bool GetNextPart(vtkLSDynaPart*& part)
  {
    if (!this->PartIterator)
    {
      part = nullptr;
      return false;
    }
    part = this->PartIterator;

    // clear iterator before we search for the next part
    vtkIdType pos = this->PartIteratorLoc + 1;
    this->PartIterator = nullptr;
    this->PartIteratorLoc = -1;

    // find the next part
    for (vtkIdType i = pos; i < this->NumParts; i++)
    {
      if (this->Parts[i] && this->Parts[i]->PartType() == part->PartType())
      {
        this->PartIteratorLoc = i;
        this->PartIterator = this->Parts[i];
        break;
      }
    }
    return true;
  }

  //---------------------------------------------------------------------------
  void AllocateParts()
  {
    vtkIdType numCells = 0, cellLength = 0;
    for (vtkIdType i = 0; i < this->NumParts; ++i)
    {
      vtkLSDynaPart* part = this->Parts[i];

      if (part)
      {
        bool canBeAllocated = this->GetInfoForPart(part, numCells, cellLength);
        if (canBeAllocated)
        {
          part->AllocateCellMemory(numCells, cellLength);
        }
        else
        {
          // this part has no cells allocated to it, so remove it now.
          part->Delete();
          this->Parts[i] = nullptr;
        }
      }
    }
    // Only needed when debugging
    // this->DumpPartInfo();
  }

  //---------------------------------------------------------------------------
  bool GetInfoForPart(vtkLSDynaPart* part, vtkIdType& numCells, vtkIdType& cellArrayLength) const

  {
    // verify that the part is valid
    numCells = 0;
    cellArrayLength = 0;
    bool validPart = part->hasValidType();
    if (!validPart)
    {
      // we return early because an invalid type would
      // cause the Info array to be accessed out of bounds
      return validPart;
    }

    // give a part type and a material id
    // walk the run length encoding to determe the total size
    std::vector<PartInfo>::const_iterator it;
    for (it = this->Info[part->PartType()].begin(); it != this->Info[part->PartType()].end(); ++it)
    {
      const PartInfo* info = &(*it);
      if (info->partId == part->GetPartId())
      {
        validPart = true;
        numCells += info->numCells;
        cellArrayLength += info->cellStructureSize;
      }
    }

    return validPart;
  }

  //---------------------------------------------------------------------------
  void DumpPartInfo()
  {
    for (int i = 0; i < LSDynaMetaData::NUM_CELL_TYPES; ++i)
    {
      // now lets dump all the part info

      std::cout << "For Info index: " << i << std::endl;
      std::cout << "We have " << this->Info[i].size() << " info entries" << std::endl;

      std::vector<PartInfo>::const_iterator it;
      for (it = this->Info[i].begin(); it != this->Info[i].end(); ++it)
      {
        const PartInfo* info = &(*it);
        if (info->part != nullptr)
        {
          std::cout << "The material id is: " << info->partId << std::endl;
          std::cout << "The numCells is:    " << info->numCells << std::endl;
          std::cout << std::endl;
          std::cout << "The Part is :" << std::endl;
          info->part->PrintSelf(cout, vtkIndent().GetNextIndent());
          std::cout << std::endl;
          std::cout << std::endl;
        }
      }
    }
  }

  //---------------------------------------------------------------------------
  void InitCellIteration(const int& partType, int pos = 0)
  {
    this->CellIteratorEnd = this->Info[partType].end();
    if (!this->Info[partType].empty())
    {
      this->CellIterator = this->Info[partType].begin();
    }
    else
    {
      this->CellIterator = this->Info[partType].end();
    }

    while (pos > 0 && this->CellIterator != this->CellIteratorEnd)
    {
      pos -= (*this->CellIterator).numCells;
      if (pos > 0)
      {
        ++this->CellIterator;
      }
    }
  }

  //---------------------------------------------------------------------------
  bool GetNextCellPart(vtkIdType& startId, vtkIdType& numCells, vtkLSDynaPart*& part)
  {
    if (this->CellIterator == this->CellIteratorEnd)
    {
      return false;
    }

    startId = (*this->CellIterator).startId;
    numCells = (*this->CellIterator).numCells;
    part = (*this->CellIterator).part;
    ++this->CellIterator;
    return true;
  }

  //---------------------------------------------------------------------------
  void FinalizeTopology()
  {
    for (vtkIdType i = 0; i < this->NumParts; ++i)
    {
      vtkLSDynaPart* part = this->Parts[i];
      if (part && part->HasCells())
      {
        part->BuildToplogy();
      }
      else if (part)
      {
        part->Delete();
        this->Parts[i] = nullptr;
      }
    }
  }

  //---------------------------------------------------------------------------
  void DisableDeadCells()
  {
    for (vtkIdType i = 0; i < this->NumParts; ++i)
    {
      vtkLSDynaPart* part = this->Parts[i];
      if (part && part->HasCells())
      {
        part->DisableDeadCells();
      }
    }
  }

  //---------------------------------------------------------------------------
  void PrintSelf(ostream& os, vtkIndent indent)
  {
    for (vtkIdType i = 0; i < this->NumParts; ++i)
    {
      os << indent << "Part Number " << i << std::endl;
      if (this->PartExists(i))
      {
        vtkLSDynaPart* part = this->Parts[i];
        part->PrintSelf(os, indent.GetNextIndent());
      }
      else
      {
        os << indent.GetNextIndent() << "Does not exist." << std::endl;
      }
    }
  }

protected:
  vtkIdType NumParts;

  // stores all the parts for this collection.
  vtkLSDynaPart** Parts;

  // maps cell indexes which are tracked by output type to the part
  // Since cells are ordered the same between the cell connectivity data block
  // and the state block in the d3plot format we only need to know which part
  // the cell is part of.
  // This info is constant for each time step
  std::vector<PartInfo>* Info;
  PartInsertion* CellInsertionIterators;

  std::vector<PartInfo>::const_iterator CellIterator, CellIteratorEnd;
  vtkLSDynaPart* PartIterator;
  vtkIdType PartIteratorLoc;
};

vtkStandardNewMacro(vtkLSDynaPartCollection);
//------------------------------------------------------------------------------
vtkLSDynaPartCollection::vtkLSDynaPartCollection()
{
  this->MetaData = nullptr;
  this->Storage = nullptr;
  this->MinIds = nullptr;
  this->MaxIds = nullptr;
}

//------------------------------------------------------------------------------
vtkLSDynaPartCollection::~vtkLSDynaPartCollection()
{
  delete this->Storage;
  delete[] this->MinIds;
  delete[] this->MaxIds;
  this->MetaData = nullptr;
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::PrintSelf(ostream& os, vtkIndent indent)
{
  // just needs to print all public accessible ivars
  this->Superclass::PrintSelf(os, indent);

  // number of parts
  os << indent << "Number of Parts: " << this->GetNumberOfParts() << std::endl;

  // print self for each part
  this->Storage->PrintSelf(os, indent.GetNextIndent());
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::InitCollection(
  LSDynaMetaData* metaData, vtkIdType* mins, vtkIdType* maxs)
{
  delete this->Storage;
  delete[] this->MinIds;
  delete[] this->MaxIds;

  // reserve enough space for the grids. Each node
  // will have a part allocated, since we don't know yet
  // how the cells map to parts.
  this->Storage = new LSDynaPartStorage(static_cast<vtkIdType>(metaData->PartIds.size()));

  this->MinIds = new vtkIdType[LSDynaMetaData::NUM_CELL_TYPES];
  this->MaxIds = new vtkIdType[LSDynaMetaData::NUM_CELL_TYPES];

  // We only have to map the cell ids between min and max, so we
  // skip into the proper place
  for (int i = 0; i < LSDynaMetaData::NUM_CELL_TYPES; ++i)
  {
    this->MinIds[i] = (mins != nullptr) ? mins[i] : 0;
    this->MaxIds[i] = (maxs != nullptr) ? maxs[i] : metaData->NumberOfCells[i];
  }

  if (metaData)
  {
    this->MetaData = metaData;
    this->BuildPartInfo();
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::BuildPartInfo()
{
  // we iterate on part materials as those are those are from 1 to num Parts.
  // the part ids are the user material ids

  std::vector<int>::const_iterator partMIt;
  std::vector<int>::const_iterator materialIdIt = this->MetaData->PartIds.begin();
  std::vector<int>::const_iterator statusIt = this->MetaData->PartStatus.begin();
  std::vector<std::string>::const_iterator nameIt = this->MetaData->PartNames.begin();

  for (partMIt = this->MetaData->PartMaterials.begin();
       partMIt != this->MetaData->PartMaterials.end();
       ++partMIt, ++statusIt, ++nameIt, ++materialIdIt)
  {
    if (*statusIt)
    {
      // make the index contain a part
      this->Storage->ConstructPart((*partMIt) - 1, *nameIt, *materialIdIt,
        this->MetaData->NumberOfNodes, this->MetaData->Fam.GetWordSize());
    }
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::RegisterCellIndexToPart(
  const int& partType, const vtkIdType& matId, const vtkIdType&, const vtkIdType& npts)

{
  this->Storage->RegisterCell(partType, matId - 1, npts);
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::AllocateParts()
{
  this->Storage->AllocateParts();
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::InitCellInsertion()
{
  this->Storage->InitCellInsertion();
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::InsertCell(const int& partType, const vtkIdType&, const int& cellType,
  const vtkIdType& npts, vtkIdType conn[8])
{

  this->Storage->InsertCell(partType, cellType, npts, conn);
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::SetCellDeadFlags(
  const int& partType, vtkUnsignedCharArray* death, const int& deadCellsAsGhostArray)
{
  // go through and flag each part cell as deleted or not.
  // this means breaking up this array into an array for each part
  if (!death)
  {
    return;
  }

  // The array that passed in from the reader only contains the subset
  // of the full data that we are interested in so we don't have to adjust
  // any indices
  this->Storage->InitCellIteration(partType);
  vtkIdType numCells, startId;
  vtkLSDynaPart* part;
  unsigned char* dead = static_cast<unsigned char*>(death->GetVoidPointer(0));
  while (this->Storage->GetNextCellPart(startId, numCells, part))
  {
    // perfectly valid to have a nullptr part being returned
    // just skip it as the user doesn't want it loaded.
    if (part)
    {
      part->EnableDeadCells(deadCellsAsGhostArray);
      part->SetCellsDeadState(dead, numCells);
    }
    dead += numCells;
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::AddProperty(const LSDynaMetaData::LSDYNA_TYPES& type,
  const char* name, const int& offset, const int& numComps)
{

  vtkLSDynaPart* part = nullptr;
  this->Storage->InitPartIteration(type);
  while (this->Storage->GetNextPart(part))
  {
    if (part)
    {
      part->AddCellProperty(name, offset, numComps);
    }
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::FillCellProperties(float* buffer,
  const LSDynaMetaData::LSDYNA_TYPES& type, const vtkIdType& startId, const vtkIdType& numCells,
  const int& numPropertiesInCell)
{
  this->FillCellArray(buffer, type, startId, numCells, numPropertiesInCell);
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::FillCellProperties(double* buffer,
  const LSDynaMetaData::LSDYNA_TYPES& type, const vtkIdType& startId, const vtkIdType& numCells,
  const int& numPropertiesInCell)
{
  this->FillCellArray(buffer, type, startId, numCells, numPropertiesInCell);
}

//------------------------------------------------------------------------------
template <typename T>
void vtkLSDynaPartCollection::FillCellArray(T* buffer, const LSDynaMetaData::LSDYNA_TYPES& type,
  const vtkIdType& startId, vtkIdType numCells, const int& numPropertiesInCell)
{
  // we only need to iterate the array for the subsection we need
  T* loc = buffer;
  vtkIdType size, globalStartId;
  vtkLSDynaPart* part;
  this->Storage->InitCellIteration(type, startId);
  while (this->Storage->GetNextCellPart(globalStartId, size, part))
  {
    vtkIdType start = std::max(globalStartId, startId);
    vtkIdType end = std::min(globalStartId + size, startId + numCells);
    if (end < start)
    {
      break;
    }
    vtkIdType is = end - start;
    if (part)
    {
      part->ReadCellProperties(loc, is, numPropertiesInCell);
    }
    loc += is * numPropertiesInCell;
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::ReadCellUserIds(
  const LSDynaMetaData::LSDYNA_TYPES& type, const int& status)
{

  vtkIdType numCells, numSkipStart, numSkipEnd;
  this->GetPartReadInfo(type, numCells, numSkipStart, numSkipEnd);

  if (!status)
  {
    // skip this part type
    this->MetaData->Fam.SkipWords(numSkipStart + numCells + numSkipEnd);
    return;
  }

  this->MetaData->Fam.SkipWords(numSkipStart);
  vtkIdType numChunks = this->MetaData->Fam.InitPartialChunkBuffering(numCells, 1);
  vtkIdType startId = 0;
  if (this->MetaData->Fam.GetWordSize() == 8 && numCells > 0)
  {
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      vtkIdType chunkSize = this->MetaData->Fam.GetNextChunk(LSDynaFamily::Float);
      vtkIdType numCellsInChunk = chunkSize;
      vtkIdType* buf = this->MetaData->Fam.GetBufferAs<vtkIdType>();
      this->FillCellUserId(buf, type, startId, numCellsInChunk);
      startId += numCellsInChunk;
    }
  }
  else if (numCells > 0)
  {
    for (vtkIdType i = 0; i < numChunks; ++i)
    {
      vtkIdType chunkSize = this->MetaData->Fam.GetNextChunk(LSDynaFamily::Float);
      vtkIdType numCellsInChunk = chunkSize;
      int* buf = this->MetaData->Fam.GetBufferAs<int>();
      this->FillCellUserId(buf, type, startId, numCellsInChunk);
      startId += numCellsInChunk;
    }
  }
  this->MetaData->Fam.SkipWords(numSkipEnd);

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

//------------------------------------------------------------------------------
template <typename T>
void vtkLSDynaPartCollection::FillCellUserIdArray(
  T* buffer, const LSDynaMetaData::LSDYNA_TYPES& type, const vtkIdType& startId, vtkIdType numCells)
{
  const int numWordsPerIdType(this->MetaData->Fam.GetWordSize() / sizeof(T));
  // we only need to iterate the array for the subsection we need
  T* loc = buffer;
  vtkIdType size, globalStartId;
  vtkLSDynaPart* part;
  this->Storage->InitCellIteration(type, startId);
  while (this->Storage->GetNextCellPart(globalStartId, size, part))
  {
    vtkIdType start = std::max(globalStartId, startId);
    vtkIdType end = std::min(globalStartId + size, startId + numCells);
    if (end < start)
    {
      break;
    }
    vtkIdType is = (end - start) * numWordsPerIdType;
    if (part)
    {
      part->EnableCellUserIds();
      for (vtkIdType i = 0; i < is; i += numWordsPerIdType)
      {
        part->SetNextCellUserIds((vtkIdType)loc[i]);
      }
    }
    // perfectly valid to have a nullptr part being returned
    // just skip it as the user doesn't want it loaded.
    loc += is;
  }
}

//------------------------------------------------------------------------------
bool vtkLSDynaPartCollection::IsActivePart(const int& id) const
{
  return this->Storage->PartExists(id);
}

//------------------------------------------------------------------------------
vtkUnstructuredGrid* vtkLSDynaPartCollection::GetGridForPart(const int& index) const
{
  return this->Storage->GetPartGrid(index);
}

//------------------------------------------------------------------------------
int vtkLSDynaPartCollection::GetNumberOfParts() const
{
  return static_cast<int>(this->Storage->GetNumParts());
}
//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::DisbleDeadCells()
{
  this->Storage->DisableDeadCells();
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::GetPartReadInfo(const int& partType, vtkIdType& numberOfCells,
  vtkIdType& numCellsToSkipStart, vtkIdType& numCellsToSkipEnd) const
{
  vtkIdType size = this->MaxIds[partType] - this->MinIds[partType];
  if (size <= 0)
  {
    numberOfCells = 0;
    // skip everything
    numCellsToSkipStart = this->MetaData->NumberOfCells[partType];
    numCellsToSkipEnd = 0; // no reason to skip anything else
  }
  else
  {
    numberOfCells = size;
    numCellsToSkipStart = this->MinIds[partType];
    numCellsToSkipEnd =
      this->MetaData->NumberOfCells[partType] - (numberOfCells + numCellsToSkipStart);
  }
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::FinalizeTopology()
{
  this->Storage->FinalizeTopology();
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::ReadPointUserIds(const vtkIdType& numTuples, const char* name)
{
  this->SetupPointPropertyForReading(numTuples, 1, name, true, true, false, false);
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::ReadPointProperty(const vtkIdType& numTuples,
  const vtkIdType& numComps, const char* name, const bool& isProperty, const bool& isGeometryPoints,
  const bool& isRoadPoints)
{
  this->SetupPointPropertyForReading(
    numTuples, numComps, name, false, isProperty, isGeometryPoints, isRoadPoints);
}

//------------------------------------------------------------------------------
void vtkLSDynaPartCollection::SetupPointPropertyForReading(const vtkIdType& numTuples,
  const vtkIdType& numComps, const char* name, const bool& isIdType, const bool& isProperty,
  const bool& isGeometryPoints, const bool& isRoadPoints)
{
  if (!isProperty && !isGeometryPoints && !isRoadPoints)
  {
    // don't read arrays the user didn't request, just skip them
    this->MetaData->Fam.SkipWords(numTuples * numComps);
    return;
  }

  // If this is a geometeric point property it needs to apply
  // to the following
  // BEAM,SHELL,THICK_SHELL,SOLID,Particles
  // if it is road surface it only applies to RigidSurfaceData part types
  vtkLSDynaPart* part = nullptr;
  vtkLSDynaPart** validParts = new vtkLSDynaPart*[this->Storage->GetNumParts()];
  vtkIdType idx = 0;
  if (!isRoadPoints)
  {
    enum LSDynaMetaData::LSDYNA_TYPES validCellTypes[5] = { LSDynaMetaData::PARTICLE,
      LSDynaMetaData::BEAM, LSDynaMetaData::SHELL, LSDynaMetaData::THICK_SHELL,
      LSDynaMetaData::SOLID };
    for (int i = 0; i < 5; ++i)
    {
      this->Storage->InitPartIteration(validCellTypes[i]);
      while (this->Storage->GetNextPart(part))
      {
        part->AddPointProperty(name, numComps, isIdType, isProperty, isGeometryPoints);
        validParts[idx++] = part;
      }
    }
  }
  else
  {
    // is a road point
    this->Storage->InitPartIteration(LSDynaMetaData::ROAD_SURFACE);
    while (this->Storage->GetNextPart(part))
    {
      part->AddPointProperty(name, numComps, isIdType, isProperty, isGeometryPoints);
      validParts[idx++] = part;
    }
  }

  if (idx <= 0)
  {
    // don't do anything as we have no valid parts
  }
  else if (this->MetaData->Fam.GetWordSize() == 8)
  {
    this->FillPointProperty<double>(numTuples, numComps, validParts, idx);
  }
  else
  {
    this->FillPointProperty<float>(numTuples, numComps, validParts, idx);
  }

  delete[] validParts;
}

namespace
{
// this function is used to sort a collection of parts
// based on the max and min global point ids that the part
// we use both to enforce better weak ordering
bool sortPartsOnGlobalIds(const vtkLSDynaPart* p1, const vtkLSDynaPart* p2)
{
  return p1->GetMaxGlobalPointId() < p2->GetMaxGlobalPointId();
}
}

//------------------------------------------------------------------------------
template <typename T>
void vtkLSDynaPartCollection::FillPointProperty(
  const vtkIdType& numTuples, const vtkIdType& numComps, vtkLSDynaPart** parts, vtkIdType numParts)
{
  LSDynaMetaData* p = this->MetaData;

  // construct the sorted array of parts so we only
  // have to iterate a subset that are interested in the points we have
  // are reading in.
  std::list<vtkLSDynaPart*> sortedParts(parts, parts + numParts);
  std::list<vtkLSDynaPart*>::iterator partIt;

  sortedParts.sort(sortPartsOnGlobalIds);

  // find the max as the subset of points
  const vtkIdType maxGlobalPoint(sortedParts.back()->GetMaxGlobalPointId());
  vtkIdType minGlobalPoint = maxGlobalPoint;
  for (partIt = sortedParts.begin(); partIt != sortedParts.end(); ++partIt)
  {
    minGlobalPoint = std::min((*partIt)->GetMinGlobalPointId(), minGlobalPoint);
  }

  const vtkIdType realNumberOfTuples(maxGlobalPoint - minGlobalPoint);
  const vtkIdType numPointsToSkipStart(minGlobalPoint);
  const vtkIdType numPointsToSkipEnd(numTuples - (realNumberOfTuples + minGlobalPoint));

  vtkIdType offset = numPointsToSkipStart;
  const vtkIdType numPointsToRead(1048576);
  const vtkIdType loopTimes(realNumberOfTuples / numPointsToRead);
  const vtkIdType leftOver(realNumberOfTuples % numPointsToRead);
  const vtkIdType bufferChunkSize(numPointsToRead * numComps);

  T* buf = nullptr;
  p->Fam.SkipWords(numPointsToSkipStart * numComps);
  for (vtkIdType j = 0; j < loopTimes; ++j, offset += numPointsToRead)
  {
    p->Fam.BufferChunk(LSDynaFamily::Float, bufferChunkSize);
    buf = p->Fam.GetBufferAs<T>();

    partIt = sortedParts.begin();
    while (partIt != sortedParts.end() && (*partIt)->GetMaxGlobalPointId() < offset)
    {
      // remove all parts from the list that have already been
      // filled by previous loops
      sortedParts.pop_front();
      partIt = sortedParts.begin();
    }

    while (partIt != sortedParts.end())
    {
      // only read the points which have a point that lies within this section
      // so we stop once the min is larger than our max id
      (*partIt)->ReadPointBasedProperty(buf, numPointsToRead, numComps, offset);
      ++partIt;
    }
  }
  if (leftOver > 0 && !sortedParts.empty())
  {
    p->Fam.BufferChunk(LSDynaFamily::Float, leftOver * numComps);
    buf = p->Fam.GetBufferAs<T>();
    for (partIt = sortedParts.begin(); partIt != sortedParts.end(); ++partIt)
    {
      (*partIt)->ReadPointBasedProperty(buf, leftOver, numComps, offset);
    }
  }
  p->Fam.SkipWords(numPointsToSkipEnd * numComps);
}
VTK_ABI_NAMESPACE_END
