MercuryData.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef TOOLS_MERCURYDATA_H
#define TOOLS_MERCURYDATA_H
 
#include <fstream>
#include <iostream>
#include <string>
#include <sstream>
#include <vector>
 
template<std::size_t NDIMS>
struct MercuryParticle
{
  public:
    double position[NDIMS];
    double velocity[NDIMS];
    double rotation[NDIMS];
    double angularV[NDIMS];
    
    double radius;
    
    std::size_t speciesID;
};
 
template<>
struct MercuryParticle<2>
{
  public:
    double position[3];
    double velocity[3];
    double rotation[3];
    double angularV[3];
    
    double radius;
    
    std::size_t speciesID;
 
};
 
template<std::size_t NDIMS>
std::istream& operator>>(std::istream& in, MercuryParticle<NDIMS>& part)
{
  std::size_t i;
  for (i = 0; i < NDIMS; i++)
    in >> part.position[i];
  
  for (i = 0; i < NDIMS; i++) {
    in >> part.velocity[i];
    // convert velocities to float, because paraview cannot handle doubles
    part.velocity[i] = (float)part.velocity[i];
  }
 
  in >> part.radius;
  
  for (i = 0; i < NDIMS; i++)
    in >> part.rotation[i];
 
  for (i = 0; i < NDIMS; i++) {
    in >> part.angularV[i];
    // convert velocities to float, because paraview cannot handle doubles
    part.angularV[i] = (float)part.angularV[i];
  }
 
  in >> part.speciesID;
  
  return in;
}
 
template<>
std::istream& operator>><2>(std::istream& in, MercuryParticle<2>& part)
{
  in >> part.position[0] >> part.position[2];
  part.position[1] = 0;
  
  in >> part.velocity[0] >> part.velocity[2];
  part.velocity[1] = 0;
  
  in >> part.radius;
 
  in >> part.rotation[1];
  part.rotation[0] = part.rotation[2] = 0;  
 
  in >> part.angularV[1];
  part.angularV[0] = part.angularV[2] = 0;
    
  in >> part.speciesID;
  
  return in;
}
 
 
class MercuryDataFile;
 
template<std::size_t NDIMS>
class MercuryTimeStepIterator;
template<std::size_t NDIMS>
class MercuryTimeStep
{
  public:
    double getTime() const
    {
      return time_;
    }
    
    std::size_t getTimeStepID() const
    {
      return ID_;
    }
    
    std::size_t getNumberOfParticles() const
    {
      return numParticles_;
    }
    
    std::size_t size() const
    {
      return numParticles_;
    }
    
    constexpr std::size_t getNumberOfDimensions() const
    {
      return NDIMS;
    }
    
    typename std::vector< MercuryParticle<NDIMS> >::iterator begin()
    {
      return storage_.begin();
    }
    
    typename std::vector< MercuryParticle<NDIMS> >::const_iterator begin() const
    {
      return storage_.begin();
    }
    
    typename std::vector< MercuryParticle<NDIMS> >::iterator end()
    {
      return storage_.end();
    }
    
    typename std::vector< MercuryParticle<NDIMS> >::const_iterator end() const
    {
      return storage_.end();
    }
    
    MercuryParticle<NDIMS>& operator[](std::size_t idx)
    {
      return storage_[idx];
    }
    
    const MercuryParticle<NDIMS>& operator[](std::size_t idx) const
    {
      return storage_[idx];
    }
    
  private:
    MercuryTimeStep(std::size_t id, MercuryDataFile *pData)
    : time_(0), ID_(id), numParticles_(0), dataFile_(pData)
    {
    }
    
    MercuryTimeStep()
    : time_(0), ID_(0), numParticles_(0), dataFile_(nullptr)
    {
    }
    
    double time_;
    std::size_t ID_;
    std::size_t numParticles_;
    double min_[NDIMS], max_[NDIMS];
    MercuryDataFile *dataFile_;
    
    std::vector< MercuryParticle<NDIMS> > storage_;
    
    template<std::size_t NDIMS2>
    friend std::istream& operator>>(std::istream&, MercuryTimeStep<NDIMS2>&);
    
    friend class MercuryTimeStepIterator<NDIMS>;
    
    friend class MercuryDataFile;
};
 
template<std::size_t NDIMS>
std::istream& operator>>(std::istream& in, MercuryTimeStep<NDIMS>& step)
{
  std::size_t i;
  in >> step.numParticles_ >> step.time_;
  
  for (i = 0; i < NDIMS; i++)
    in >> step.min_[i];
    
  for (i = 0; i < NDIMS; i++)
    in >> step.max_[i];
 
  return in;
}
 
template<std::size_t NDIMS>
class MercuryTimeStepIterator
{
  public:
    bool operator!=(MercuryTimeStepIterator<NDIMS> other) const
    {
      return (isEOFTimeStep_ != other.isEOFTimeStep_);
    }
    
    MercuryTimeStep<NDIMS>& operator*()
    {
      return lastReadTimeStep_;
    }
    
    const MercuryTimeStep<NDIMS>& operator*() const
    {
      return lastReadTimeStep_;
    }
    
    void operator++();
  
  private:
    MercuryTimeStepIterator()
      : isEOFTimeStep_(true), dataFile_(nullptr)
    {
    }
    
    MercuryTimeStepIterator(MercuryDataFile *pData)
      : lastReadTimeStep_(0,pData), isEOFTimeStep_(false), dataFile_(pData)
    {
      ++(*this);
      lastReadTimeStep_.ID_ = 0;
    }
  
    MercuryTimeStep<NDIMS> lastReadTimeStep_;
    bool isEOFTimeStep_;
    MercuryDataFile *dataFile_;
    
    friend class MercuryDataFile;
};
 
 
class MercuryDataFile
{
  public:
    MercuryDataFile(std::string name)
      : file_(name)
    { }
 
    operator bool() const
    {
      return file_.good();
    }
    
    template<std::size_t NDIMS>
    bool isMercuryDataFile()
    {
      //Store the position, so we can jump back at the end of the function..
      std::ios::pos_type currentPosition = file_.tellg();
      //and jump to the start
      file_.seekg(0);
      //get the first line
      std::string line;
      std::getline(file_, line);
      file_.seekg(currentPosition); //and pretend nothing has happened
      
      std::istringstream lineStream(line);
      
      //We'll try to find out if there were exactly enough arguments.
      MercuryTimeStep<NDIMS> step;
      lineStream >> step;
      
      //Did we reach the end yet?
      bool isValid = lineStream.good() || (lineStream.eof() && !lineStream.fail());
      double dummy;
      lineStream >> dummy;
      
      //now we should have reached it.
      isValid = isValid && !lineStream.good();
      return isValid;
    }
    
    template<std::size_t NDIMS>
    class IteratorProxy
    {
    private:
      IteratorProxy(MercuryDataFile* pData)
        : data_(pData)
      { }
      
      MercuryDataFile* data_;
    public:
      MercuryTimeStepIterator<NDIMS> begin()
      { 
        return data_->begin<NDIMS>();
      }
      MercuryTimeStepIterator<NDIMS> end()
      {
        return data_->end<NDIMS>();
      }
      friend class MercuryDataFile;
    };
    
    template<std::size_t NDIMS>
    IteratorProxy<NDIMS> as()
    { 
      return {this};
    };
    
    template<std::size_t NDIMS>
    MercuryTimeStepIterator<NDIMS> begin()
    {
      file_.seekg(0);
      return {this};
    }
    
    template<std::size_t NDIMS>
    MercuryTimeStepIterator<NDIMS> end() const
    {
        return {};
    }
  private:
    std::ifstream file_;
    
    template<std::size_t NDIMS>
    friend class MercuryTimeStep;
    template<std::size_t NDIMS>
    friend class MercuryTimeStepIterator;
};
 
template<std::size_t NDIMS>
void MercuryTimeStepIterator<NDIMS>::operator++()
{
  lastReadTimeStep_.ID_++;
      
  std::string line;
  std::getline(dataFile_->file_, line);
  
  std::istringstream lineStream(line);
  
  lineStream >> lastReadTimeStep_;
  
  //I hope we didn't went beyond end of file...
  if (lineStream.eof())
  {
//   logger(WARN, "The time step header detected an EOF.. Usually this"
//                " means that the format was not what it appeared to be."
//                "\nproceed with caution!");
  }
  //Resize the backing storage container to make sure we can actually
  //fit all the particles in there.
  lastReadTimeStep_.storage_.resize(lastReadTimeStep_.numParticles_);  
  //Well, now that we're set up, read all the particles
  for (MercuryParticle<NDIMS>& part : lastReadTimeStep_)
  {
    //line by line, because no data format can be trusted.
    std::getline(dataFile_->file_, line);
    lineStream.clear();
    lineStream.str(line);
   
    lineStream >> part;
  }
  
  if (dataFile_->file_.eof())
    isEOFTimeStep_ = true;
}
 
#endif