vtkreader.hh

Go to the documentation of this file.

#pragma once
 
#include <iosfwd>
#include <map>
#include <memory>
#include <vector>
 
#include <dune/common/shared_ptr.hh>
#include <dune/common/typelist.hh>
#include <dune/common/typeutilities.hh>
 
#include <dune/vtk/filereader.hh>
#include <dune/vtk/types.hh>
#include <dune/vtk/utility/errors.hh>
 
// default GridCreator
#include <dune/vtk/gridcreators/continuousgridcreator.hh>
#include <dune/vtk/gridfunctions/common.hh>
 
namespace Dune::Vtk
{
 
  template <class G, class GC = Vtk::ContinuousGridCreator<G>, class FT = double>
  class VtkReader
      : public Vtk::FileReader<G, VtkReader<G, GC, FT>>
  {
    // Sections visited during the xml parsing
    enum Sections {
      NO_SECTION = 0, VTK_FILE, UNSTRUCTURED_GRID, PIECE, POINT_DATA, PD_DATA_ARRAY, CELL_DATA, CD_DATA_ARRAY,
      POINTS, POINTS_DATA_ARRAY, CELLS, CELLS_DATA_ARRAY, APPENDED_DATA, XML_NAME, XML_NAME_ASSIGN, XML_VALUE
    };
 
  public:
    // Type storing information about read data
    struct DataArrayAttributes
    {
      std::string name;
      Vtk::DataTypes type = Vtk::DataTypes::UNKNOWN;
      unsigned int components = 1;
      std::uint64_t offset = 0;
      Sections section = NO_SECTION;
    };
 
  public:
    using Grid = G;
 
    using GridCreator = GC;
 
    using FieldType = FT;
 
    using GlobalCoordinate = typename GridCreator::GlobalCoordinate;
 
    // Template representing a grid-function that is created in getPointData() and getCellData()
    // with Context either Vtk::PointContext or Vtk::CellContext, respectively.
    // To each GridCreator a GridFunction is associated, see, e.g. Vtk::ContinuousGridFunction
    // or Vtk::LagrangeGridFunction.
    template <class Context, class Range>
    using GridFunction = typename Vtk::AssociatedGridFunction<GridCreator, Range, Context>::type;
 
  public:
 
    template <class... Args,
      std::enable_if_t<std::is_constructible<GridCreator, Args...>::value,int> = 0>
    explicit VtkReader (Args&&... args)
      : VtkReader(std::make_shared<GridCreator>(std::forward<Args>(args)...))
    {}
 
    explicit VtkReader (GridCreator& creator)
      : VtkReader(stackobject_to_shared_ptr(creator))
    {}
 
    explicit VtkReader (GridCreator&& creator)
      : VtkReader(std::make_shared<GridCreator>(std::move(creator)))
    {}
 
    explicit VtkReader (std::shared_ptr<GridCreator> creator)
      : creator_(std::move(creator))
    {}
 
 
    void read (std::string const& filename, bool fillCreator = true);
 
    GridCreator& gridCreator ()
    {
      return *creator_;
    }
 
    std::unique_ptr<Grid> createGrid () const
    {
      return creator_->createGrid();
    }
 
    template <class Range = DynamicVector<FieldType>>
    GridFunction<Vtk::PointContext,Range> getPointData (std::string const& name) const
    {
      auto data_it = dataArray_.find("PointData." + name);
      auto point_it = pointData_.find("PointData." + name);
      VTK_ASSERT_MSG(data_it != dataArray_.end() && point_it != pointData_.end(),
        "The data to extract is not found in point-data. Try `getCellData()` instead!");
      VTK_ASSERT(data_it->second.section == POINT_DATA);
 
      return {*creator_, point_it->second, data_it->second.name, data_it->second.components,
              data_it->second.type, vec_types, vec_offsets, vec_connectivity};
    }
 
    std::vector<DataArrayAttributes> getPointDataAttributes () const
    {
      std::vector<DataArrayAttributes> attributes;
      attributes.reserve(pointData_.size());
      for (auto const& da : dataArray_) {
        if (da.second.section == POINT_DATA)
          attributes.push_back(da.second);
      }
      return attributes;
    }
 
    template <class Range = DynamicVector<FieldType>>
    GridFunction<Vtk::CellContext,Range> getCellData (std::string const& name) const
    {
      auto data_it = dataArray_.find("CellData." + name);
      auto cell_it = cellData_.find("CellData." + name);
      VTK_ASSERT_MSG(data_it != dataArray_.end() && cell_it != cellData_.end(),
        "The data to extract is not found in cell-data. Try `getPointData()` instead!");
      VTK_ASSERT(data_it->second.section == CELL_DATA);
 
      return {*creator_, cell_it->second, data_it->second.name, data_it->second.components,
              data_it->second.type, vec_types, vec_offsets, vec_connectivity};
    }
 
    std::vector<DataArrayAttributes> getCellDataAttributes () const
    {
      std::vector<DataArrayAttributes> attributes;
      attributes.reserve(cellData_.size());
      for (auto const& da : dataArray_) {
        if (da.second.section == CELL_DATA)
          attributes.push_back(da.second);
      }
      return attributes;
    }
 
 
 
    void readSerialFileFromStream (std::ifstream& input, bool create = true);
 
 
    void readParallelFileFromStream (std::ifstream& input, int rank, int size, bool create = true);
 
    void fillGridCreator (bool insertPieces = true);
 
 
    std::vector<std::string> const& pieces () const
    {
      return pieces_;
    }
 
#ifndef DOXYGEN
    // Implementation of the FileReader interface
    static void fillFactoryImpl (GridFactory<Grid>& factory, std::string const& filename)
    {
      VtkReader reader{factory};
      reader.read(filename);
    }
#endif
 
  private:
    // Read values stored on the cells with ID `id`
    Sections readCellData (std::ifstream& input, std::string id);
 
    template <class F, class H>
    void readCellDataAppended (MetaType<F>, MetaType<H>, std::ifstream& input, std::string id);
 
    // Read values stored on the points with ID `id`
    Sections readPointData (std::ifstream& input, std::string id);
 
    template <class F, class H>
    void readPointDataAppended (MetaType<F>, MetaType<H>, std::ifstream& input, std::string id);
 
 
    // Read vertex coordinates from `input` stream and store in into `factory`
    Sections readPoints (std::ifstream& input, std::string id);
 
    template <class F, class H>
    void readPointsAppended (MetaType<F>, MetaType<H>, std::ifstream& input);
 
 
    // Read cell type, cell offsets and connectivity from `input` stream
    Sections readCells (std::ifstream& input, std::string id);
 
    template <class H>
    void readCellsAppended (MetaType<H>, std::ifstream& input);
 
    // Read data from appended section in vtk file, starting from `offset`
    template <class FloatType, class HeaderType>
    void readAppended (std::ifstream& input, std::vector<FloatType>& values, HeaderType offset);
 
    // Test whether line belongs to section
    bool isSection (std::string line,
                    std::string key,
                    Sections current,
                    Sections parent = NO_SECTION) const
    {
      bool result = line.substr(1, key.length()) == key;
      if (result && current != parent)
        DUNE_THROW(Exception , "<" << key << "> in wrong section." );
      return result;
    }
 
    // Convert a section into a string
    std::string toString (Sections s) const;
 
    // Find beginning of appended binary data
    std::uint64_t findAppendedDataPosition (std::ifstream& input) const;
 
    // Read attributes from current xml tag
    std::map<std::string, std::string> parseXml (std::string const& line, bool& closed);
 
    // clear all vectors
    void clear ();
 
    auto comm () const
    {
      return MPIHelper::getCommunication();
    }
 
  private:
    std::shared_ptr<GridCreator> creator_;
 
    Vtk::FormatTypes format_ = Vtk::FormatTypes::UNKNOWN;
 
    Vtk::CompressorTypes compressor_ = Vtk::CompressorTypes::NONE;
 
    // Temporary data to construct the grid elements
    std::vector<GlobalCoordinate> vec_points;
    std::vector<std::uint64_t> vec_point_ids;   //< Global unique vertex ID
    std::vector<std::uint8_t> vec_types;        //< VTK cell type ID
    std::vector<std::int64_t> vec_offsets;      //< offset of vertices of cell
    std::vector<std::int64_t> vec_connectivity; //< vertex indices of cell
 
    std::size_t numberOfCells_ = 0;   //< Number of cells in the grid
    std::size_t numberOfPoints_ = 0;  //< Number of vertices in the grid
 
    // offset information for appended data
    // map: id -> {Name,DataType,NumberOfComponents,Offset}
    std::map<std::string, DataArrayAttributes> dataArray_;
 
    // storage for internal point and cell data
    // map: id -> vector of data entries per point/cell
    std::map<std::string, std::vector<FieldType>> pointData_;
    std::map<std::string, std::vector<FieldType>> cellData_;
 
    // vector of filenames of parallel pieces
    std::vector<std::string> pieces_;
 
    std::uint64_t offset0_ = 0;
 
    bool read_ = false;
  };
 
  // deduction guides
  template <class Grid>
  VtkReader (GridFactory<Grid>&)
    -> VtkReader<Grid, ContinuousGridCreator<Grid>>;
 
  template <class GridCreator,
    class = std::void_t<typename GridCreator::Grid>>
  VtkReader (GridCreator&)
    -> VtkReader<typename GridCreator::Grid, GridCreator>;
 
  template <class GridCreator,
    class = std::void_t<typename GridCreator::Grid>>
  VtkReader (GridCreator&&)
    -> VtkReader<typename GridCreator::Grid, GridCreator>;
 
  template <class GridCreator,
    class = std::void_t<typename GridCreator::Grid>>
  VtkReader (std::shared_ptr<GridCreator>)
    -> VtkReader<typename GridCreator::Grid, GridCreator>;
 
} // end namespace Dune::Vtk
 
namespace Dune
{
  template <class Grid, class GC = Vtk::ContinuousGridCreator<Grid>, class FieldType = double>
  class [[deprecated("Use Vtk::VtkReader instead.")]] VtkReader
    : public Vtk::VtkReader<Grid,GC,FieldType>
  {
    using Base = Vtk::VtkReader<Grid,GC,FieldType>;
  public:
    using Base::Base;
  };
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
  // deduction guides
  template <class Grid>
  VtkReader (GridFactory<Grid>&)
    -> VtkReader<Grid, Vtk::ContinuousGridCreator<Grid>>;
 
  template <class GridCreator,
    class = std::void_t<typename GridCreator::Grid>>
  VtkReader (GridCreator&)
    -> VtkReader<typename GridCreator::Grid, GridCreator>;
 
  template <class GridCreator,
    class = std::void_t<typename GridCreator::Grid>>
  VtkReader (std::shared_ptr<GridCreator>)
    -> VtkReader<typename GridCreator::Grid, GridCreator>;
#pragma GCC diagnostic pop
 
} // end namespace Dune
 
#include "vtkreader.impl.hh"