points.h

#ifndef POINTS_H
#define POINTS_H
 
#include "crystal.h"
#include "eigen.h"
#include "exceptions.h"
 
// TODO this should be an enum 
const int crystalCoordinates_ = 0;
const int cartesianCoordinates_ = 1;
 
class Points;
 
class Point {
public:
  Point(Points &points_, const int &index_,
        const Eigen::Vector3d &umklappVector = Eigen::Vector3d::Zero());
 
  Point(const Point &that);
 
  Point &operator=(const Point &that);
 
  Eigen::Vector3d getCoordinates(const int &basis = crystalCoordinates_,
                                 const bool &inWignerSeitz = false) const;
 
  Point operator+(Point &b);
 
  Point operator-(Point &b);
 
  bool hasUmklapp();
 
  int getIndex() const;
 
private:
  Eigen::Vector3d umklappVector;
  int index;
  Points &points;
};
 
class Points {
public:
 
  Points(Crystal &crystalObj_);
 
  Points(Crystal &crystalObj_, const Eigen::Vector3i &mesh_,
         const Eigen::Vector3d &offset_ = Eigen::Vector3d::Zero());
 
  Points(Crystal &crystal_, const Eigen::Tensor<double, 3> &pathExtrema,
         const double &delta);
 
  void setActiveLayer(const Eigen::VectorXi &filter_);
 
  Points(const Points &obj);
 
  Points &operator=(const Points &obj);
 
  std::tuple<Eigen::Vector3i, Eigen::Vector3d> getMesh();
 
  int getNumPoints() const;
 
  Eigen::Vector3d crystalToCartesian(const Eigen::Vector3d &point);
 
  Eigen::Vector3d cartesianToCrystal(const Eigen::Vector3d &point);
 
  Eigen::Vector3d bzToWs(const Eigen::Vector3d &point, const int &basis);
 
  Eigen::Vector3d foldToBz(const Eigen::Vector3d &pointCrystal,
                           const int &basis);
 
  static std::tuple<Eigen::Vector3i, Eigen::Vector3d>
  findMesh(const Eigen::Matrix<double, 3, Eigen::Dynamic> &points);
 
  Crystal &getCrystal();
 
  Point getPoint(const int &index);
 
  Eigen::Vector3d getPointCoordinates(const int &index,
                                      const int &basis = crystalCoordinates);
 
  int getIndex(const Eigen::Vector3d &point);
 
  int isPointStored(const Eigen::Vector3d &crystalCoordinates_);
 
  // note: we could use constexpr to tell the compiler that the class member is
  // available at compilation time. But it wouldn't be compatible on old
  // clusters
  static const int crystalCoordinates;
  static const int cartesianCoordinates;
 
  void
  setIrreduciblePoints(std::vector<Eigen::MatrixXd> *groupVelocities = nullptr);
 
  std::vector<int> irrPointsIterator();
 
  std::vector<int> parallelIrrPointsIterator();
 
  int asIrreducibleIndex(const int &ik);
 
  int asReducibleIndex(const int &ik);
 
  std::tuple<int, Eigen::Matrix3d>
  getRotationToIrreducible(const Eigen::Vector3d &x,
                           const int &basis = cartesianCoordinates);
 
  std::vector<Eigen::Matrix3d> getRotationsStar(const int &ik);
 
  std::vector<int> getReducibleStarFromIrreducible(const int &ik);
 
  //void swapCrystal(Crystal &newCrystal);
 
protected:
  void setMesh(const Eigen::Vector3i &mesh_, const Eigen::Vector3d &offset_);
 
  Crystal &crystalObj;
  Eigen::Vector3i mesh;
  Eigen::Vector3d offset;
  int numPoints = 0;
  // for Wigner Seitz folding
  Eigen::MatrixXd gVectors;
 
  // to store points on a path or active points
  bool explicitlyStored = false;
  bool isPointsListSorted = false;
  Eigen::MatrixXd pointsList;
 
  // for active Points:
  Eigen::VectorXi filteredToFullIndices;
 
  void setupGVectors();
 
  //------------------------------------------
  // Symmetries
 
  std::vector<Eigen::Matrix3d> rotationMatricesCrystal;
  std::vector<Eigen::Matrix3d> rotationMatricesCartesian;
 
  // vector of size 0 to numPoints: given index i, mERI(i) is the index of
  // the rotation R such that R*kIrr = k(i)
  Eigen::VectorXi mapEquivalenceRotationIndex;
 
  // vector of size 0 to numIrrPoints: given an irreducible wavevector i,
  // where i runs from 0 to numIrrPoints,
  // mapIrreducibleToReducibleList(i) is the index of the irreducible point in
  // the reducible list
  Eigen::VectorXi mapIrreducibleToReducibleList;
 
  // vector of size 0 to numPoints. Given a reducible wavevector i,
  // mRTI
  Eigen::VectorXi mapReducibleToIrreducibleList;
  int numIrrPoints = 0;
 
  // vector of size numIrrPoints. Given an irreducible point i
  // (i from 0,numIrrPoints), gets the list of reducible point indices
  // that are symmetry equivalent to point i.
  std::vector<std::vector<int>> irreducibleStars;
 
  // if equiv(i) == i, point is irreducible, otherwise, equiv(i) gives the
  // index of the irreducible equivalent index.
  Eigen::VectorXi equiv;
  //------------------------------------------
};
 
#endif