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

   Program:   ParaView
   Module:    vtkCFSReader.cxx

   Copyright (c) 2005,2006 Sandia Corporation, Kitware Inc.
   All rights reserved.

   ParaView is a free software; you can redistribute it and/or modify it
   under the terms of the ParaView license version 1.2.

   See License_v1.2.txt for the full ParaView license.
   A copy of this license can be obtained by contacting
   Kitware Inc.
   28 Corporate Drive
   Clifton Park, NY 12065
   USA

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

/*-------------------------------------------------------------------------
   Copyright (c) 2022 Verein zur Foerderung der Software openCFS
-------------------------------------------------------------------------*/

#include <vtkBiQuadraticQuad.h>
#include <vtkCFSReader.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkDoubleArray.h>
#include <vtkErrorCode.h>
#include <vtkHexahedron.h>
#include <vtkInformation.h>
#include <vtkInformationStringKey.h>
#include <vtkInformationVector.h>
#include <vtkIntArray.h>
#include <vtkLine.h>
#include <vtkMath.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
#include <vtkPyramid.h>
#include <vtkQuad.h>
#include <vtkQuadraticEdge.h>
#include <vtkQuadraticHexahedron.h>
#include <vtkQuadraticPyramid.h>
#include <vtkQuadraticQuad.h>
#include <vtkQuadraticTetra.h>
#include <vtkQuadraticTriangle.h>
#include <vtkQuadraticWedge.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkTetra.h>
#include <vtkTriQuadraticHexahedron.h>
#include <vtkTriangle.h>
#include <vtkUnsignedIntArray.h>
#include <vtkUnstructuredGrid.h>
#include <vtkVertex.h>
#include <vtkWedge.h>
#include <vtksys/SystemTools.hxx>

#include <list>
#include <sstream>

vtkStandardNewMacro(vtkCFSReader);

// Initialize Information Keys
vtkInformationKeyMacro(vtkCFSReader, CFS_RESULT_NAME, String);
vtkInformationKeyMacro(vtkCFSReader, CFS_DOF_NAMES, StringVector);
vtkInformationKeyMacro(vtkCFSReader, CFS_DEFINED_ON, Integer);
vtkInformationKeyMacro(vtkCFSReader, CFS_ENTITY_NAME, String);
vtkInformationKeyMacro(vtkCFSReader, CFS_ENTRY_TYPE, Integer);
vtkInformationKeyMacro(vtkCFSReader, CFS_STEP_NUMS, IntegerVector);
vtkInformationKeyMacro(vtkCFSReader, CFS_STEP_VALUES, DoubleVector);
vtkInformationKeyMacro(vtkCFSReader, CFS_UNIT, String);
vtkInformationKeyMacro(vtkCFSReader, CFS_ENTITY_IDS, StringVector);
vtkInformationKeyMacro(vtkCFSReader, CFS_MULTI_SEQ_INDEX, Integer);
vtkInformationKeyMacro(vtkCFSReader, CFS_ANALYSIS_TYPE, Integer);

//-----------------------------------------------------------------------------
vtkCFSReader::vtkCFSReader()
{
  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(1);
}

//-----------------------------------------------------------------------------
vtkCFSReader::~vtkCFSReader()
{
  if (this->MBDataSet)
  {
    this->MBDataSet->Delete();
  }
  if (this->MBActiveDataSet)
  {
    this->MBActiveDataSet->Delete();
  }
}

//-----------------------------------------------------------------------------
// taken from vtkSetStringMacro and enhanced to reset history result info
void vtkCFSReader::SetFileName(const char* arg)
{
  // do nothing if nothing changes
  if (this->FileName.empty() && arg == nullptr)
  {
    return;
  }
  if (!this->FileName.empty() && arg != nullptr && strcmp(this->FileName.c_str(), arg) == 0)
  {
    return;
  }

  // something changed
  this->Reader.CloseFile();

  if (arg != nullptr)
  {
    this->FileName.assign(arg);
  }
  else
  {
    this->FileName.clear();
  }
  this->Modified();
  this->Hdf5InfoRead = false; // re-initialize info on history results
}

//-----------------------------------------------------------------------------
const char* vtkCFSReader::GetFileName()
{
  if (this->FileName.empty())
  {
    return nullptr;
  }
  else
  {
    return this->FileName.c_str();
  }
}

//-----------------------------------------------------------------------------
void vtkCFSReader::CloseFile()
{
  this->Reader.CloseFile();

  this->FileName.clear();
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetMultiSequenceStep(int step)
{
  if (step == static_cast<int>(this->MultiSequenceStep))
  {
    return;
  }
  if (step > this->MultiSequenceRange[1] || step < this->MultiSequenceRange[0])
  {
    vtkErrorMacro(<< "Please enter a valid multisequence step between "
                  << this->MultiSequenceRange[0] << " and " << this->MultiSequenceRange[1] << "!"
                  << "selected value was " << step << "\n");
  }

  this->MultiSequenceStep = step;
  this->MSStepChanged = true;

  // In addition trigger resetting the data value arrays, as the results in
  // the new msStep can be different than from the previous one.
  this->ResetDataArrays = true;

  this->Modified();
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetTimeStep(unsigned int step)
{
  vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting TimeStep to " << step);

  // security check: immediately leave, if no time / frequency steps are available
  if (this->StepVals.empty())
  {
    return;
  }

  step++;

  if (this->TimeStep != step && this->StepVals.size() >= step)
  {
    this->TimeStep = step;
    this->TimeOrFrequencyValue = this->StepVals[step - 1];
    this->TimeOrFrequencyValueStr = std::to_string(this->TimeOrFrequencyValue);

    // update pipeline
    this->Modified();
  }
}

//-----------------------------------------------------------------------------
const char* vtkCFSReader::GetTimeOrFrequencyValueStr()
{
  return (!this->HarmonicData || this->HarmonicDataAsModeShape == 0)
    ? "-"
    : this->TimeOrFrequencyValueStr.c_str();
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetAddDimensionsToArrayNames(int flag)
{
  this->AddDimensionsToArrayNames = flag;

  // In addition trigger resetting the data value arrays
  this->ResetDataArrays = true;

  // update pipeline
  this->Modified();
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetHarmonicDataAsModeShape(int flag)
{
  this->HarmonicDataAsModeShape = flag;

  if (flag == 1 && this->HarmonicData)
  {
    this->TimeStepNumberRange[1] = static_cast<int>(this->NumberOfTimeSteps - 1);
  }
  else
  {
    this->TimeStepNumberRange[1] = 0;
  }

  // In addition trigger resetting the data value arrays
  this->ResetDataArrays = true;

  // update pipeline
  this->Modified();
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetFillMissingResults(int flag)
{
  this->FillMissingResults = flag;

  // In addition trigger resetting the data value arrays
  this->ResetDataArrays = true;

  // update pipeline
  this->Modified();
}

//-----------------------------------------------------------------------------
int vtkCFSReader::GetRegionArrayStatus(const char* name)
{
  assert(name != nullptr);
  std::map<std::string, int>::const_iterator it = this->RegionSwitch.find(name);
  if (it != this->RegionSwitch.end())
  {
    return it->second;
  }
  vtkErrorMacro(<< "Region '" << name << "' not found.");
  return 0;
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetRegionArrayStatus(const char* name, int status)
{
  assert(name != nullptr);
  this->RegionSwitch[name] = status;

  // IMPORTANT: Let VTK know, that the pipeline needs an update
  // (otherwise nothing would happen)
  this->Modified();

  // In addition, set flag to false, that regions have to be updated
  this->ActiveRegionsChanged = true;
}

//-----------------------------------------------------------------------------
const char* vtkCFSReader::GetRegionArrayName(int index)
{
  // we assume the vector is properly traversed
  assert(index < this->RegionNames.size());
  return this->RegionNames[index].c_str();
}

//-----------------------------------------------------------------------------
int vtkCFSReader::GetNamedNodeArrayStatus(const char* name)
{
  assert(name != nullptr);
  std::map<std::string, int>::const_iterator it = this->NamedNodeSwitch.find(name);
  if (it != this->NamedNodeSwitch.end())
  {
    return it->second;
  }
  vtkErrorMacro(<< "Named nodes '" << name << "' not found.");
  return 0;
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetNamedNodeArrayStatus(const char* name, int status)
{
  assert(name != nullptr);
  this->NamedNodeSwitch[name] = status;
  this->Modified();
  // In addition, set flag to false, that regions have to be updated
  this->ActiveRegionsChanged = true;
}

//-----------------------------------------------------------------------------
const char* vtkCFSReader::GetNamedNodeArrayName(int index)
{
  assert(index < this->NamedNodeNames.size());
  return this->NamedNodeNames[index].c_str();
}

//-----------------------------------------------------------------------------
int vtkCFSReader::GetNamedElementArrayStatus(const char* name)
{
  assert(name != nullptr);
  std::map<std::string, int>::const_iterator it = this->NamedElementSwitch.find(name);
  if (it != this->NamedElementSwitch.end())
  {
    return it->second;
  }
  vtkErrorMacro(<< "Named elems '" << name << "' not found.");
  return 0;
}

//-----------------------------------------------------------------------------
void vtkCFSReader::SetNamedElementArrayStatus(const char* name, int status)
{
  assert(name != nullptr);
  this->NamedElementSwitch[name] = status;
  this->Modified();
  // In addition, set flag to false, that regions have to be updated
  this->ActiveRegionsChanged = true;
}

//-----------------------------------------------------------------------------
const char* vtkCFSReader::GetNamedElementArrayName(int index)
{
  assert(index < this->NamedElementNames.size());
  return this->NamedElementNames[index].c_str();
}

//-----------------------------------------------------------------------------
int vtkCFSReader::CanReadFile(const char* fname)
{
  if (fname == nullptr)
  {
    return 0;
  }

  // Try to "load" the file (which internally just opens the mesh group)
  // to see, if this is a true CFS HDF5 file.
  std::string fileName = std::string(fname);
  H5CFS::Hdf5Reader temp;
  try
  {
    temp.LoadFile(fileName);
    // If this succeeds, we assume that this is a true CFS HDF5 file.
    return 1;
  }
  catch (...)
  {
    return 0;
  }
}

//-----------------------------------------------------------------------------
void vtkCFSReader::ReadHdf5Informations()
{
  const int enableGroups = 0;

  // Load hdf file (only done first time)
  if (!this->Hdf5InfoRead)
  {
    try
    {
      this->Reader.LoadFile(this->FileName);
      this->Dimension = static_cast<int>(this->Reader.GetDimension());
      this->GridOrder = static_cast<int>(this->Reader.GetGridOrder());
      this->RegionNames = this->Reader.GetAllRegionNames();
      this->NamedNodeNames = this->Reader.GetNodeNames();
      this->NamedElementNames = this->Reader.GetElementNames();

      for (auto& name : this->RegionNames)
      {
        auto it = this->RegionSwitch.find(name);
        if (it == this->RegionSwitch.end())
        {
          this->RegionSwitch[name] = 1;
        }
      }

      for (auto& name : this->NamedNodeNames)
      {
        auto it = this->NamedNodeSwitch.find(name);
        if (it == this->NamedNodeSwitch.end())
        {
          this->NamedNodeSwitch[name] = enableGroups;
        }
      }

      for (auto& name : this->NamedElementNames)
      {
        auto it = this->NamedElementSwitch.find(name);
        if (it == this->NamedElementSwitch.end())
        {
          this->NamedElementSwitch[name] = enableGroups;
        }
      }

      // store info on multisequence steps with mesh results
      std::map<unsigned int, unsigned int> numSteps;
      this->Reader.GetNumberOfMultiSequenceSteps(this->ResAnalysisType, numSteps, false);

      // Get available mesh results
      for (auto msIt : this->ResAnalysisType)
      {
        this->Reader.GetResultTypes(msIt.first, this->MeshResultInfos[msIt.first], false);
      }

      // store info on multisequence steps for history results
      this->Reader.GetNumberOfMultiSequenceSteps(this->HistAnalysisType, numSteps, true);

      // Get available history results
      for (auto msIt : this->HistAnalysisType)
      {
        this->Reader.GetResultTypes(msIt.first, this->HistResultInfos[msIt.first], true);
      }
    }
    catch (const std::runtime_error& er)
    {
      vtkErrorMacro(<< "Caught exception when trying to read hdf5 info from file '"
                    << this->FileName << "' : '" << er.what() << "'");
      return;
    }

    this->Hdf5InfoRead = true;
    this->MSStepChanged = true;
    this->IsInitialized = false;
  } // info red
}

//-----------------------------------------------------------------------------
int vtkCFSReader::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  // Read basic information from HDF5
  this->ReadHdf5Informations();
  if (!this->Hdf5InfoRead)
  {
    return 0;
  }

  try
  {
    // Check for new multisequence step
    if (this->MSStepChanged)
    {
      // Request number of multiSequenceSteps

      // adjust range for
      if (!this->ResAnalysisType.empty())
      {
        this->MultiSequenceRange[0] = static_cast<int>(this->ResAnalysisType.begin()->first);
        this->MultiSequenceRange[1] = static_cast<int>(this->ResAnalysisType.rbegin()->first);

        // check, if any result at all are present
        // Run over all result types
        std::map<unsigned int, double> globSteps;

        if (this->ResAnalysisType.find(this->MultiSequenceStep) != this->ResAnalysisType.end())
        {
          // store current analysis type
          this->AnalysisType = this->ResAnalysisType[this->MultiSequenceStep];
          this->HarmonicData = (this->AnalysisType == H5CFS::HARMONIC ||
            this->AnalysisType == H5CFS::EIGENVALUE || this->AnalysisType == H5CFS::EIGENFREQUENCY);

          for (auto& mri : MeshResultInfos[this->MultiSequenceStep])
          {
            std::map<unsigned int, double> actStepVals;
            this->Reader.GetStepValues(this->MultiSequenceStep, mri->name, actStepVals, false);
            std::map<unsigned int, double>::const_iterator it;
            for (it = actStepVals.begin(); it != actStepVals.end(); it++)
            {
              globSteps[it->first] = it->second;
            }
          }
        }

        this->StepNumbers.resize(globSteps.size());
        this->StepVals.resize(globSteps.size());
        auto itStep = globSteps.begin();
        for (unsigned int i = 0; itStep != globSteps.end(); ++itStep, ++i)
        {
          this->StepNumbers[i] = itStep->first;
          this->StepVals[i] = itStep->second;
        }
      }
      this->NumberOfTimeSteps = static_cast<unsigned int>(this->StepNumbers.size());
      if (this->HarmonicData && this->HarmonicDataAsModeShape != 0)
      {
        this->TimeStepNumberRange[1] = static_cast<int>(this->NumberOfTimeSteps - 1);
      }
      else
      {
        this->TimeStepNumberRange[1] = 0;
      }
      if (this->NumberOfTimeSteps != 0)
      {
        this->TimeStepValuesRange[0] = this->StepVals[0];
        this->TimeStepValuesRange[1] = this->StepVals[NumberOfTimeSteps - 1];
      }

      this->MSStepChanged = false;
    } // if multisequence step has changed
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when trying to request information from file '"
                  << this->FileName << "' : '" << er.what() << "'");
    return 0;
  }

  if (this->NumberOfTimeSteps > 0)
  {
    if (this->HarmonicData && this->HarmonicDataAsModeShape != 0)
    {
      // From harmonic data we can produce a continuous time data set which ranges from 0 to one
      // All other kinds of data are discrete in time/freq.
      this->TimeStepValuesRange[0] = 0;
      this->TimeStepValuesRange[1] = 1.0;

      outputVector->GetInformationObject(0)->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    }
    else
    {
      outputVector->GetInformationObject(0)->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
        &(this->StepVals[0]), static_cast<int>(this->NumberOfTimeSteps));
    }

    outputVector->GetInformationObject(0)->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(),
      static_cast<double*>(this->TimeStepValuesRange), 2);
  }

  return 1;
}

//-----------------------------------------------------------------------------
int vtkCFSReader::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  // obtain detail for what data is actually requested
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkMultiBlockDataSet* output = vtkMultiBlockDataSet::GetData(outInfo);

  // ===================
  // Time Stepping Part
  // ===================
  // Attention: Only change value of pipeline / step values, if
  // there are any results present at all.
  if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()) != 0 &&
    !this->StepVals.empty())
  {
    RequestedTimeValue = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
    int nSteps = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    double* rsteps = outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    int timeI = 0;
    while (timeI < nSteps - 1 && rsteps[timeI] < this->RequestedTimeValue)
    {
      timeI++;
    }

    if (this->HarmonicData && this->HarmonicDataAsModeShape > 0)
    {
      outInfo->Set(vtkDataObject::DATA_TIME_STEP(), this->RequestedTimeValue);
    }
    else
    {
      outInfo->Set(vtkDataObject::DATA_TIME_STEP(), this->StepVals[timeI]);
      this->TimeStep = timeI + 1;
      this->TimeOrFrequencyValue = this->StepVals[timeI];
      this->TimeOrFrequencyValueStr = std::to_string(this->TimeOrFrequencyValue);
    }
  }

  // Read in data from file
  this->ReadFile(output);

  return 1;
}

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

  // vtkGetStringMacro
  os << indent << "File Name: " << this->FileName << "\n";

  // vtkGetMacro of interest
  os << indent << "Multi Sequence Step: " << this->GetMultiSequenceStep() << "\n";
  os << indent << "Time Step: " << this->GetTimeStep() << "\n";
  os << indent << "Complex Model Real" << this->GetComplexModeReal() << "\n";
  os << indent << "Complex Model Imaginary" << this->GetComplexModeImaginary() << "\n";
  os << indent << "Complex Model Amplitude" << this->GetComplexModeAmplitude() << "\n";
  os << indent << "Complex Model Phase" << this->GetComplexModePhase() << "\n";

  // further properties of interest
  os << indent << "Num Steps: " << this->GetNumberOfSteps() << "\n";
  os << indent << "Time Freq: " << this->GetTimeOrFrequencyValueStr() << "\n";
  os << indent << "Grid Dimensions: " << this->GetGridDimension() << "\n";
  os << indent << "Grid Order: " << this->GetGridOrder() << "\n";
  os << indent << "Number Region Arrays: " << this->GetNumberOfRegionArrays() << "\n";
  os << indent << "Number Named Node Arrays: " << this->GetNumberOfNamedNodeArrays() << "\n";
  os << indent << "Number Named Element Arrays: " << this->GetNumberOfNamedElementArrays() << "\n";
}

//----------------------------------------------------------------------------
void vtkCFSReader::ReadFile(vtkMultiBlockDataSet* output)
{
  try
  {
    // The first time, the reader is initialized, we have to
    // fill the  multi block data set
    if (!this->IsInitialized)
    {
      this->MBDataSet = vtkMultiBlockDataSet::New();
      this->MBDataSet->ShallowCopy(output);

      // Assemble number of blocks, which is
      // #regions, #named elemens, #named nodes
      unsigned int numBlocks = static_cast<unsigned int>(
        this->RegionNames.size() + this->NamedElementNames.size() + this->NamedNodeNames.size());

      for (unsigned int i = 0, n = numBlocks; i < n; i++)
      {
        vtkNew<vtkUnstructuredGrid> newGrid;
        this->MBDataSet->SetBlock(i, newGrid.GetPointer());
      }

      // initialize mapping from
      // (region,globalNodeNumber)->regionLocalNodeNumber
      this->NodeMap.resize(numBlocks);

      // read nodal connectivity
      this->ReadNodes(this->MBDataSet);

      // read element definition
      this->ReadCells(this->MBDataSet);

      // Now the grid is finalized and we store it initially also
      // to the activeSet
      this->MBActiveDataSet = vtkMultiBlockDataSet::New();
      this->MBActiveDataSet->ShallowCopy(this->MBDataSet);

      // now it should be definitely initialized
      this->IsInitialized = true;
    }

    // If region status has changed, we have to perform an update
    // of the active regions
    if (this->ActiveRegionsChanged || this->ResetDataArrays)
    {
      this->UpdateActiveRegions();
    }

    // only read data, if time/frequency steps are defined
    if (this->NumberOfTimeSteps > 0)
    {
      // read nodal data
      this->ReadNodeCellData(this->MBActiveDataSet, true);

      // read cell data
      this->ReadNodeCellData(this->MBActiveDataSet, false);
    }
    // Use the active set as base for the current set
    output->ShallowCopy(this->MBActiveDataSet);
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when trying read file '" << this->FileName << "' : '"
                  << er.what() << "'");
  }
}

//-----------------------------------------------------------------------------
void vtkCFSReader::UpdateActiveRegions()
{
  // Idea:
  // Create new multiblock data set, based on the internal stored one and copy only those
  // grid-blocks, which are active.
  vtkNew<vtkMultiBlockDataSet> newSet;

  //  Combine region and group info
  unsigned int numBlocks = static_cast<unsigned int>(
    this->RegionNames.size() + this->NamedElementNames.size() + this->NamedNodeNames.size());
  std::vector<std::string> names(numBlocks);
  std::map<std::string, int> switches;
  for (unsigned int iBlock = 0; iBlock < numBlocks; ++iBlock)
  {
    if (iBlock < this->RegionNames.size())
    {
      std::string name = this->RegionNames[iBlock];
      names[iBlock] = name;
      switches[name] = this->RegionSwitch[name];
    }
    else if (iBlock < this->RegionNames.size() + this->NamedElementNames.size())
    {
      unsigned int index = static_cast<unsigned int>(iBlock - this->RegionNames.size());
      std::string name = this->NamedElementNames[index];
      names[iBlock] = name;
      switches[name] = this->NamedElementSwitch[name];
    }
    else
    {
      unsigned int index = static_cast<unsigned int>(
        iBlock - this->RegionNames.size() - this->NamedElementNames.size());
      std::string name = this->NamedNodeNames[index];
      names[iBlock] = name;
      switches[name] = this->NamedNodeSwitch[name];
    }
  }

  //  Run over all active regions / element groups and copy grid
  unsigned int index = 0;
  for (unsigned int i = 0; i < numBlocks; i++)
  {
    // check, if region is active
    if (switches[names[i]] != 0)
    {
      // create new grid with output and make shallow copy
      vtkNew<vtkUnstructuredGrid> newGrid;
      newGrid->ShallowCopy(this->MBDataSet->GetBlock(i));

      // try to find the corresponding block of this region/group in the current
      // active set in order to copy its data to this block.
      unsigned int numActBlocks = this->MBActiveDataSet->GetNumberOfBlocks();
      unsigned int activeIndex = 0;
      bool found = false;
      for (unsigned int j = 0; j < numActBlocks; ++j)
      {
        const char* blockName =
          this->MBActiveDataSet->GetMetaData(j)->Get(vtkCompositeDataSet::NAME());
        if (blockName != nullptr)
        {
          if (std::string(blockName) == names[i])
          {
            found = true;
            activeIndex = j;
          }
        }
      }

      if (found && !this->ResetDataArrays)
      {
        vtkUnstructuredGrid* actGrid =
          vtkUnstructuredGrid::SafeDownCast(this->MBActiveDataSet->GetBlock(activeIndex));
        newGrid->GetPointData()->ShallowCopy(actGrid->GetPointData());
        newGrid->GetCellData()->ShallowCopy(actGrid->GetCellData());
      }

      // Get the array containing the actual node numbers for the current region
      newGrid->GetPointData()->SetActiveScalars("origNodeNums");
      vtkUnsignedIntArray::SafeDownCast(newGrid->GetPointData()->GetScalars());
      newGrid->GetPointData()->SetActiveScalars("");

      newSet->SetBlock(index, newGrid);

      // Tag new block with region name
      newSet->GetMetaData(index)->Set(vtkCompositeDataSet::NAME(), names[i].c_str());
      index++;
    } // if block is active
  }   // loop: all blocks (= regions, elem/node groups)

  // in the end, replace the active set by the current one
  this->MBActiveDataSet->ShallowCopy(newSet.GetPointer());

  // in the end, re-set the status of the update flag
  this->ActiveRegionsChanged = false;
  this->ResetDataArrays = false;
}

//----------------------------------------------------------------------------
void vtkCFSReader::ReadNodeCellData(vtkMultiBlockDataSet* output, bool isNode)
{
  H5CFS::EntityType entType = isNode ? H5CFS::NODE : H5CFS::ELEMENT;
  const double pi = vtkMath::Pi();
  double h180degOverPi = 180.0 / pi;

  //  Combine region and group info
  unsigned int numBlocks = static_cast<unsigned int>(
    this->RegionNames.size() + this->NamedElementNames.size() + this->NamedNodeNames.size());
  std::vector<std::string> names(numBlocks);
  std::map<std::string, int> switches;
  for (unsigned int iBlock = 0; iBlock < numBlocks; ++iBlock)
  {
    if (iBlock < this->RegionNames.size())
    {
      std::string name = this->RegionNames[iBlock];
      names[iBlock] = name;
      switches[name] = this->RegionSwitch[name];
    }
    else if (iBlock < this->RegionNames.size() + this->NamedElementNames.size())
    {
      unsigned int index = static_cast<unsigned int>(iBlock - this->RegionNames.size());
      std::string name = this->NamedElementNames[index];
      names[iBlock] = name;
      switches[name] = this->NamedElementSwitch[name];
    }
    else
    {
      unsigned int index = static_cast<unsigned int>(
        iBlock - this->RegionNames.size() - this->NamedElementNames.size());
      std::string name = this->NamedNodeNames[index];
      names[iBlock] = name;
      switches[name] = this->NamedNodeSwitch[name];
    }
  }

  try
  {
    if (MeshResultInfos.find(this->MultiSequenceStep) == this->MeshResultInfos.end())
    {
      vtkWarningMacro(
        "Multisequence step index " + std::to_string(this->MultiSequenceStep) + " is invalid.");
      return;
    }

    // get result infos for 1ms step and step 1

    // Read all result types occuring in the current multisequence step
    std::set<std::string> resultNames;
    std::map<std::string, std::shared_ptr<H5CFS::ResultInfo>> nameInfo;
    for (const auto& mri : this->MeshResultInfos[this->MultiSequenceStep])
    {
      std::string& resultName = mri->name;
      resultNames.insert(resultName);
      nameInfo[resultName] = mri;
    }

    std::shared_ptr<H5CFS::Result> result;
    //  Loop over all result names
    for (const auto& resultName : resultNames)
    {
      // Ensure, that in the nodal case we only pick nodal results.
      // In the case of element results, we accept elements and surface elements
      if ((isNode && (nameInfo[resultName]->listType != H5CFS::NODE)) ||
        (!isNode &&
          (nameInfo[resultName]->listType != H5CFS::ELEMENT &&
            nameInfo[resultName]->listType != H5CFS::SURF_ELEM)))
      {
        continue;
      }

      //  Loop over all active blocks
      unsigned int blockIndex = 0;
      for (unsigned int iBlock = 0; iBlock < numBlocks; iBlock++)
      {
        // variable used to write the result
        std::string writeName = resultName;

        // check, if block is active. If not, just leave
        if (switches[names[iBlock]] == 0)
        {
          continue;
        }

        // get grip of grid of current block
        vtkUnstructuredGrid* actGrid =
          vtkUnstructuredGrid::SafeDownCast(output->GetBlock(blockIndex++));

        //  Find resultinfo for this result and this block
        bool found = false;
        std::shared_ptr<H5CFS::ResultInfo> actInfo;
        for (auto& mri : this->MeshResultInfos[this->MultiSequenceStep])
        {
          if (mri->name == resultName && mri->listName == names[iBlock] && mri->listType == entType)
          {
            actInfo = mri;
            found = true;
            break;
          }
        }

        // if result for this block was not found and we should not generate
        // "empty" datasets for missing result, we just leave
        if (!found)
        {
          if (this->FillMissingResults != 0)
          {
            actInfo = nameInfo[resultName];
          }
          else
          {
            continue;
          }
        }

        std::vector<unsigned int> entities = this->Reader.GetEntities(entType, names[iBlock]);
        unsigned int numEntities = static_cast<unsigned int>(entities.size());
        unsigned int numDofs = static_cast<unsigned int>(actInfo->dofNames.size());
        std::vector<double>* ptRealVals = nullptr;
        std::vector<double>* ptImaginaryVals = nullptr;
        bool isComplex = false;
        std::vector<double> dummyVals;
        if (found)
        {
          result = std::make_shared<H5CFS::Result>(H5CFS::Result());
          result->resultInfo = actInfo;
          try
          {
            this->Reader.GetMeshResult(
              this->MultiSequenceStep, this->StepNumbers[this->TimeStep - 1], result.get());
          }
          catch (...)
          {
            // If we are here, the current result is not defined on the current
            // step and we should NOT fill with empty results. In this case
            // we perform the following:
            // a) if the current time/freq step is smaller than the first
            //    one the result is defined, we update to the first defined step
            // b) otherwise we simply return and keep the result from the
            //    previous steps
            std::map<unsigned int, double> steps;
            this->Reader.GetStepValues(this->MultiSequenceStep, actInfo->name, steps, false);
            if (this->StepNumbers[this->TimeStep - 1] < steps.begin()->first)
            {
              this->Reader.GetMeshResult(
                this->MultiSequenceStep, steps.begin()->first, result.get());
            }
            else
            {
              // in this case we are in the middle of the steps and just keep
              // the old value
              break;
            }
          }
          isComplex = result->isComplex;
          ptRealVals = &(result->realVals);
          ptImaginaryVals = &(result->imaginaryVals);
        }
        else
        {
          dummyVals.resize(numEntities * numDofs);
          ptRealVals = &(dummyVals);
          ptImaginaryVals = &(dummyVals);
          // NOTE: we have to find a different strategy how to
          // find out about complex valued results
          isComplex = (this->AnalysisType == H5CFS::HARMONIC);
        }

        std::vector<double>& realVals = *(ptRealVals);

        std::vector<vtkDoubleArray*> vals;

        // REAL-VALUED
        if (!isComplex)
        {
          if (!(this->HarmonicData && this->HarmonicDataAsModeShape != 0))
          {
            if (this->AddDimensionsToArrayNames != 0)
            {
              writeName += " [";
              writeName += actInfo->unit;
              writeName += "]";
            }
            vals.push_back(SaveToArray(realVals, actInfo->dofNames, numEntities, writeName));
          }
          else
          { // We obviously have a real-valued eigenvalue solution
            unsigned int numEntries = static_cast<unsigned int>(realVals.size());

            writeName = actInfo->name;
            writeName += "-mode";
            if (this->AddDimensionsToArrayNames != 0)
            {
              writeName += " [";
              writeName += actInfo->unit;
              writeName += "]";
            }

            double reRot = cos(-2 * pi * this->RequestedTimeValue);

            for (unsigned int i = 0; i < numEntries; i++)
            {
              realVals[i] *= reRot;
            }

            vals.push_back(SaveToArray(realVals, actInfo->dofNames, numEntities, writeName));
          }
        }
        else
        {
          std::vector<double>& imaginaryVals = *ptImaginaryVals;
          if (this->ComplexModeReal > 0)
          {
            std::string resName = actInfo->name + "Real";
            if (this->AddDimensionsToArrayNames != 0)
            {
              resName += " [";
              resName += actInfo->unit;
              resName += "]";
            }
            vals.push_back(SaveToArray(realVals, actInfo->dofNames, numEntities, resName));
          }
          if (this->ComplexModeImaginary > 0)
          {
            std::string resName = actInfo->name + "Imag";
            if (this->AddDimensionsToArrayNames != 0)
            {
              resName += " [";
              resName += actInfo->unit;
              resName += "]";
            }
            vals.push_back(SaveToArray(imaginaryVals, actInfo->dofNames, numEntities, resName));
          }
          if (this->ComplexModeAmplitude > 0)
          {
            unsigned int numEntries = static_cast<unsigned int>(realVals.size());
            std::string amplitudeName = actInfo->name + "Ampl";
            if (this->AddDimensionsToArrayNames != 0)
            {
              amplitudeName += " [";
              amplitudeName += actInfo->unit;
              amplitudeName += "]";
            }
            std::vector<double> amplitude(numEntries);
            for (unsigned int i = 0; i < numEntries; i++)
            {
              amplitude[i] = hypot(realVals[i], imaginaryVals[i]);
            }
            vals.push_back(SaveToArray(amplitude, actInfo->dofNames, numEntities, amplitudeName));
          }
          if (this->ComplexModePhase > 0)
          {
            unsigned int numEntries = static_cast<unsigned int>(realVals.size());
            std::string phaseName = actInfo->name + "Phase";
            if (this->AddDimensionsToArrayNames != 0)
            {
              phaseName += " [deg]";
            }
            std::vector<double> phase(numEntries);
            for (unsigned int i = 0; i < numEntries; i++)
            {
              phase[i] = (std::abs(imaginaryVals[i]) > 1e-16)
                ? std::atan2(imaginaryVals[i], realVals[i]) * h180degOverPi
                : (realVals[i] < 0.0) ? 180 : 0;
            }
            vals.push_back(SaveToArray(phase, actInfo->dofNames, numEntities, phaseName));
          }
          // Compute time continuous field from harmonic data by multiplying complex
          // result field with rotating phasor
          if (this->HarmonicDataAsModeShape != 0)
          {
            unsigned int numEntries = static_cast<unsigned int>(realVals.size());
            writeName = actInfo->name;
            if (AddDimensionsToArrayNames != 0)
            {
              writeName += " [";
              writeName += actInfo->unit;
              writeName += "]";
            }
            std::vector<double> modeShape(numEntries);
            // Compute the rotation pointer from the RequestedTimeValue which is in the
            // interval [0,1]
            double reRot = cos(2 * pi * this->RequestedTimeValue);
            double imRot = sin(2 * pi * this->RequestedTimeValue);
            for (unsigned int i = 0; i < numEntries; i++)
            {
              // Just compute the real part of the field
              modeShape[i] = realVals[i] * reRot - imaginaryVals[i] * imRot;
            }
            // vals5 = SaveToArray( modeShape, actInfo->dofNames, numEntities, writeName );
            vals.push_back(SaveToArray(modeShape, actInfo->dofNames, numEntities, writeName));
          } // harmonic data as mode shape
        }   // if isComplex

        // save values to grid
        if (entType == H5CFS::NODE)
        {
          for (auto& val : vals)
          {
            actGrid->GetPointData()->AddArray(val);
          }
        }
        else
        {
          for (auto& val : vals)
          {
            actGrid->GetCellData()->AddArray(val);
          }
        }
        // delete values (de-crement reference counter)
        for (auto& val : vals)
        {
          val->Delete();
        }
      } // loop over active blocks
    }   // loop over result types
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when reading results: '" << er.what() << "'");
  }
}

//-----------------------------------------------------------------------------
vtkDoubleArray* vtkCFSReader::SaveToArray(const std::vector<double>& vals,
  const std::vector<std::string>& dofNames, unsigned int numEntities, const std::string& name)
{
  vtkDoubleArray* ret = vtkDoubleArray::New();
  unsigned int numDofs = static_cast<unsigned int>(dofNames.size());
  // Case 1: numDofs is 1 (scalar) or >= 3 (real vector/tensor)
  if (numDofs == 1 || numDofs >= 3)
  {
    ret->SetNumberOfComponents(static_cast<int>(numDofs));
    for (unsigned int i = 0; i < numDofs; ++i)
    {
      ret->SetComponentName(i, dofNames[i].c_str());
    }
    ret->SetNumberOfTuples(numEntities);
    ret->SetName(name.c_str());
    double* retPtr = ret->GetPointer(0);
    unsigned int numEntries = numDofs * numEntities;
    for (unsigned int j = 0; j < numEntries; j++)
    {
      retPtr[j] = vals[j];
    }
  }
  else
  {
    // Case 2: numDofs (vector field in 2D)
    // -> add artificial 3rd component, which is 0
    ret->SetNumberOfComponents(3);
    for (unsigned int i = 0; i < numDofs; ++i)
    {
      ret->SetComponentName(i, dofNames[i].c_str());
    }
    ret->SetComponentName(2, "");
    ret->SetNumberOfTuples(numEntities);
    ret->SetName(name.c_str());
    double* retPtr = ret->GetPointer(0);
    unsigned int index = 0;
    for (unsigned int iEnt = 0; iEnt < numEntities; iEnt++)
    {
      retPtr[index++] = vals[iEnt * 2 + 0];
      retPtr[index++] = vals[iEnt * 2 + 1];
      retPtr[index++] = 0.0;
    }
  }
  return ret;
}

//-----------------------------------------------------------------------------
void vtkCFSReader::ReadCells(vtkMultiBlockDataSet* output)
{
  std::vector<H5CFS::ElemType> types;
  std::vector<std::vector<unsigned int>> connect;

  try
  {
    this->Reader.GetElements(types, connect);
    vtkDebugMacro(<< "Read in the element definitions");

    // loop over all regions
    for (unsigned int iRegion = 0; iRegion < this->RegionNames.size(); iRegion++)
    {
      const std::vector<unsigned int> regionElements =
        this->Reader.GetElementsOfRegion(this->RegionNames[iRegion]);

      vtkUnstructuredGrid* actGrid = vtkUnstructuredGrid::SafeDownCast(output->GetBlock(iRegion));

      // Add a result vector containing the original element numbers
      // to the current unstructured grid.
      vtkNew<vtkUnsignedIntArray> origElementNumbers;
      origElementNumbers->SetNumberOfValues(regionElements.size());
      origElementNumbers->SetName("origElementNums");
      actGrid->GetCellData()->AddArray(origElementNumbers.GetPointer());

      // loop over all elements
      for (unsigned int i = 0; i < regionElements.size(); i++)
      {
        origElementNumbers->SetValue(i, regionElements[i]);
      }

      this->AddElements(actGrid, iRegion, regionElements, types, connect);
    }

    unsigned int offset = static_cast<unsigned int>(this->RegionNames.size());

    // loop over all element groups
    for (unsigned int iGroup = 0; iGroup < this->NamedElementNames.size(); iGroup++)
    {
      std::string groupName = this->NamedElementNames[iGroup];
      std::vector<unsigned int> groupElements = this->Reader.GetNamedElements(groupName);

      vtkUnstructuredGrid* actGrid =
        vtkUnstructuredGrid::SafeDownCast(output->GetBlock(iGroup + offset));
      this->AddElements(actGrid, iGroup + offset, groupElements, types, connect);
    }

    offset += static_cast<unsigned int>(this->NamedElementNames.size());

    //  loop over all nodal groups
    for (unsigned int iGroup = 0; iGroup < static_cast<unsigned int>(this->NamedNodeNames.size());
         iGroup++)
    {
      std::string groupName = this->NamedNodeNames[iGroup];
      std::vector<unsigned int> groupNodes = this->Reader.GetNamedNodes(groupName);
      unsigned int numGroupNodes = static_cast<unsigned int>(groupNodes.size());

      // setup array with virtual element numbers and virtual connqectivity
      std::vector<unsigned int> vElementNumbers(numGroupNodes);
      std::vector<std::vector<unsigned int>> vElemConnect(numGroupNodes);
      std::vector<H5CFS::ElemType> vElemType(numGroupNodes);
      for (unsigned int iNode = 0; iNode < numGroupNodes; ++iNode)
      {
        vElementNumbers[iNode] = iNode + 1;
        vElemConnect[iNode].resize(1);
        vElemConnect[iNode][0] = groupNodes[iNode];
        vElemType[iNode] = H5CFS::ET_POINT;
      }

      vtkUnstructuredGrid* actGrid =
        vtkUnstructuredGrid::SafeDownCast(output->GetBlock(iGroup + offset));
      AddElements(actGrid, iGroup + offset, vElementNumbers, vElemType, vElemConnect);
    }
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when reading cells: '" << er.what() << "'");
  }
}

//-----------------------------------------------------------------------------
void vtkCFSReader::AddElements(vtkUnstructuredGrid* actGrid, unsigned int blockIndex,
  const std::vector<unsigned int>& elems, std::vector<H5CFS::ElemType>& types,
  std::vector<std::vector<unsigned int>>& connect)
{
  VTKCellType elemType = VTK_EMPTY_CELL;
  try
  {
    unsigned int numElements = static_cast<unsigned int>(elems.size());
    actGrid->Allocate(numElements);

    // loop over all elements
    for (unsigned int i = 0; i < numElements; i++)
    {
      unsigned int actElem = elems[i];

      // create correct vtk element
      elemType = GetCellIdType(types[actElem - 1]);
      if (elemType == VTK_EMPTY_CELL)
      {
        vtkErrorMacro(<< "FE type " << types[actElem - 1] << " not implemented yet");
      }

      // set connectivity
      std::array<vtkIdType, 27> nodeIds;
      for (unsigned int j = 0, n = static_cast<unsigned int>(connect[actElem - 1].size()); j < n;
           j++)
      {
        vtkDebugMacro(<< "addding nodeNum" << connect[actElem - 1][j]);
        nodeIds[j] = this->NodeMap[blockIndex][connect[actElem - 1][j] - 1] - 1;
      }
      if (elemType == VTK_TRIQUADRATIC_HEXAHEDRON)
      {
        // top
        //   7--14--6
        //   |      |
        //  15  25  13
        //   |      |
        //   4--12--5

        //   middle
        //  19--23--18
        //   |      |
        //  20  26  21
        //   |      |
        //  16--22--17

        //  bottom
        //   3--10--2
        //   |      |
        //  11  24  9
        //   |      |
        //   0-- 8--1

        nodeIds[20] = this->NodeMap[blockIndex][connect[actElem - 1][23] - 1] - 1;
        nodeIds[21] = this->NodeMap[blockIndex][connect[actElem - 1][21] - 1] - 1;
        nodeIds[22] = this->NodeMap[blockIndex][connect[actElem - 1][20] - 1] - 1;
        nodeIds[23] = this->NodeMap[blockIndex][connect[actElem - 1][22] - 1] - 1;
      }

      actGrid->InsertNextCell(elemType, connect[actElem - 1].size(), nodeIds.data());
    } // loop: elements
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when adding elements: '" << er.what() << "'");
  }
}

//-----------------------------------------------------------------------------
void vtkCFSReader::ReadNodes(vtkMultiBlockDataSet* output)
{
  try
  {
    // get all nodal coordinates
    std::vector<std::vector<double>> coords;
    this->Reader.GetNodeCoordinates(coords);

    // loop over all regions
    for (unsigned int iRegion = 0; iRegion < this->RegionNames.size(); iRegion++)
    {
      NodeMap[iRegion].resize(coords.size());
      const std::vector<unsigned int> regionNodes =
        this->Reader.GetNodesOfRegion(this->RegionNames[iRegion]);
      unsigned int numRegionNodes = static_cast<unsigned int>(regionNodes.size());

      vtkUnstructuredGrid* actGrid = vtkUnstructuredGrid::SafeDownCast(output->GetBlock(iRegion));

      vtkNew<vtkPoints> points;
      points->SetNumberOfPoints(numRegionNodes);

      // Add a result vector containing the original node numbers
      // to the current unstructured grid.
      vtkNew<vtkUnsignedIntArray> origNodeNumbers;
      origNodeNumbers->SetNumberOfValues(numRegionNodes);
      origNodeNumbers->SetName("origNodeNums");
      actGrid->GetPointData()->AddArray(origNodeNumbers.GetPointer());

      for (unsigned int i = 0; i < numRegionNodes; i++)
      {
        const std::vector<double>& p = coords[regionNodes[i] - 1];
        origNodeNumbers->SetValue(i, regionNodes[i]);
        NodeMap[iRegion][regionNodes[i] - 1] = i + 1;
        points->SetPoint(i, p[0], p[1], p[2]);
      }
      actGrid->SetPoints(points.GetPointer());
    }

    // loop over all element groups and set up nodal mapping
    unsigned int pos = static_cast<unsigned int>(RegionNames.size());
    for (unsigned int iGroup = 0; iGroup < this->NamedElementNames.size(); iGroup++, pos++)
    {
      std::string groupName = this->NamedElementNames[iGroup];
      std::vector<unsigned int> groupNodes = this->Reader.GetNamedNodes(groupName);

      this->NodeMap[pos].resize(coords.size());
      unsigned int numGroupNodes = static_cast<unsigned int>(groupNodes.size());

      vtkUnstructuredGrid* actGrid = vtkUnstructuredGrid::SafeDownCast(output->GetBlock(pos));

      vtkNew<vtkPoints> points;
      points->SetNumberOfPoints(numGroupNodes);

      for (unsigned int i = 0; i < numGroupNodes; i++)
      {
        const std::vector<double>& p = coords[groupNodes[i] - 1];
        this->NodeMap[pos][groupNodes[i] - 1] = i + 1;
        points->SetPoint(i, p[0], p[1], p[2]);
      }
      actGrid->SetPoints(points.GetPointer());
    }

    // loop over all nodal groups and set up nodal mapping
    for (unsigned int iGroup = 0; iGroup < this->NamedNodeNames.size(); iGroup++, pos++)
    {
      std::string groupName = this->NamedNodeNames[iGroup];
      std::vector<unsigned int> groupNodes = this->Reader.GetNamedNodes(groupName);
      NodeMap[pos].resize(coords.size());
      unsigned int numGroupNodes = static_cast<unsigned int>(groupNodes.size());

      vtkUnstructuredGrid* actGrid = vtkUnstructuredGrid::SafeDownCast(output->GetBlock(pos));

      vtkNew<vtkPoints> points;
      points->SetNumberOfPoints(numGroupNodes);

      for (unsigned int i = 0; i < numGroupNodes; i++)
      {
        const std::vector<double>& p = coords[groupNodes[i] - 1];
        this->NodeMap[pos][groupNodes[i] - 1] = i + 1;
        points->SetPoint(i, p[0], p[1], p[2]);
      }
      actGrid->SetPoints(points.GetPointer());
    }
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when reading nodes: '" << er.what() << "'");
  }

  vtkDebugMacro(<< "Finished reading nodes");
}

//-----------------------------------------------------------------------------
VTKCellType vtkCFSReader::GetCellIdType(H5CFS::ElemType cfsType)
{
  VTKCellType type = VTK_EMPTY_CELL;
  switch (cfsType)
  {
    case H5CFS::ET_UNDEF:
      break;
    case H5CFS::ET_POINT:
      type = VTK_VERTEX;
      break;
    case H5CFS::ET_LINE2:
      type = VTK_LINE;
      break;
    case H5CFS::ET_LINE3:
      type = VTK_QUADRATIC_EDGE;
      break;
    case H5CFS::ET_TRIA3:
      type = VTK_TRIANGLE;
      break;
    case H5CFS::ET_TRIA6:
      type = VTK_QUADRATIC_TRIANGLE;
      break;
    case H5CFS::ET_QUAD4:
      type = VTK_QUAD;
      break;
    case H5CFS::ET_QUAD8:
      type = VTK_QUADRATIC_QUAD;
      break;
    case H5CFS::ET_QUAD9:
      type = VTK_BIQUADRATIC_QUAD;
      break;
    case H5CFS::ET_TET4:
      type = VTK_TETRA;
      break;
    case H5CFS::ET_TET10:
      type = VTK_QUADRATIC_TETRA;
      break;
    case H5CFS::ET_HEXA8:
      type = VTK_HEXAHEDRON;
      break;
    case H5CFS::ET_HEXA20:
      type = VTK_QUADRATIC_HEXAHEDRON;
      break;
    case H5CFS::ET_HEXA27:
      type = VTK_TRIQUADRATIC_HEXAHEDRON;
      break;
    case H5CFS::ET_PYRA5:
      type = VTK_PYRAMID;
      break;
    case H5CFS::ET_PYRA13:
    case H5CFS::ET_PYRA14:
      type = VTK_QUADRATIC_PYRAMID;
      break;
    case H5CFS::ET_WEDGE6:
      type = VTK_WEDGE;
      break;
    case H5CFS::ET_WEDGE15:
      type = VTK_QUADRATIC_WEDGE;
      break;
    case H5CFS::ET_WEDGE18:
      type = VTK_BIQUADRATIC_QUADRATIC_WEDGE;
      break;
  }
  return type;
}

//-----------------------------------------------------------------------------
void vtkCFSReader::GetNodeCoordinates(vtkDoubleArray* vtkCoords)
{

  // Read basic information from HDF5
  ReadHdf5Informations();
  if (!this->Hdf5InfoRead)
  {
    return;
  }

  std::vector<std::vector<double>> coords;
  try
  {
    this->Reader.GetNodeCoordinates(coords);
  }
  catch (const std::runtime_error& er)
  {
    vtkErrorMacro(<< "Caught exception when obtaining node coords: '" << er.what() << "'");
    return;
  }

  unsigned int numNodes = static_cast<unsigned int>(coords.size());

  vtkCoords->SetNumberOfComponents(3);
  vtkCoords->SetNumberOfTuples(numNodes);

  if (vtkCoords->HasStandardMemoryLayout())
  {
    double* vtkPtr = vtkCoords->GetPointer(0);
    unsigned int idx = 0;

    for (unsigned int i = 0; i < numNodes; ++i)
    {
      vtkPtr[idx] = coords[i][0];
      vtkPtr[idx + 1] = coords[i][1];
      vtkPtr[idx + 2] = coords[i][2];
      idx += 3;
    }
  }
  else
  {
    vtkErrorMacro(<< "Contiguous memory layout is required.");
  }
}
