Grid.hpp

#pragma once
#include <string>
#include <vector>
 
//==============================================================================
enum class GridType { loglinear, logarithmic, linear };
 
struct GridParameters {
  std::size_t num_points;
  double r0, rmax, b;
  GridType type;
 
  GridParameters(std::size_t innum_points = 1, double inr0 = 1.0,
                 double inrmax = 1.0, double inb = 4.0,
                 GridType intype = GridType::loglinear, double indu = 0);
 
  GridParameters(std::size_t innum_points, double inr0, double inrmax,
                 double inb, const std::string &str_type = "loglinear",
                 double indu = 0);
 
  static GridType parseType(const std::string &str_type);
  static std::string parseType(GridType type);
};
 
//==============================================================================
//==============================================================================
class Grid {
 
private:
  // Minimum grid value
  double m_r0;
  // Uniform (u) grid step size
  double m_du;
  // GridType: loglinear, logarithmic, linear [enum]
  GridType gridtype;
  // log-linear grid 'turning point'
  double m_b;
  // Grid values r[i] = r(u)
  std::vector<double> m_r;
  // Convinient: (1/r)*(dr/du)[i]
  std::vector<double> m_drduor;
  // Jacobian (dr/du)[i]
  std::vector<double> m_drdu;
 
public:
  Grid(double in_r0, double in_rmax, std::size_t in_num_points,
       GridType in_gridtype, double in_b = 0);
 
  Grid(const GridParameters &in);
 
  auto r0() const { return m_r.front(); }
  auto front() const { return m_r.front(); }
  auto rmax() const { return m_r.back(); }
  auto back() const { return m_r.back(); }
  auto num_points() const { return m_r.size(); }
  auto size() const { return num_points(); }
  auto du() const { return m_du; }
  auto type() const { return gridtype; }
  auto loglin_b() const { return m_b; }
 
  const std::vector<double> &r() const { return m_r; }
  auto r(std::size_t i) const { return m_r.at(i); };
  auto at(std::size_t i) const { return m_r.at(i); };
  auto operator()(std::size_t i) const { return r(i); }
  const std::vector<double> &drdu() const { return m_drdu; };
  auto drdu(std::size_t i) const { return m_drdu.at(i); };
  const std::vector<double> &drduor() const { return m_drduor; };
  auto drduor(std::size_t i) const { return m_drduor.at(i); };
 
  std::size_t getIndex(double x, bool require_nearest = false) const;
 
  std::string gridParameters() const;
 
  std::vector<double> rpow(double k) const;
 
  auto begin() const { return m_r.cbegin(); }
  auto end() const { return m_r.cend(); }
  auto cbegin() const { return m_r.cbegin(); }
  auto cend() const { return m_r.cend(); }
  auto rbegin() const { return m_r.crbegin(); }
  auto rend() const { return m_r.crend(); }
  auto crbegin() const { return m_r.crbegin(); }
  auto crend() const { return m_r.crend(); }
 
 
  void extend_to(double new_rmax);
 
  GridParameters params() const {
    return GridParameters{num_points(), m_r0, rmax(), m_b, gridtype, m_du};
  }
 
  // Static functions: can be called outside of instantialised object
  static double calc_du_from_num_points(double in_r0, double in_rmax,
                                        std::size_t in_num_points,
                                        GridType in_gridtype, double in_b = 0);
  static std::size_t calc_num_points_from_du(double in_r0, double in_rmax,
                                             double in_du, GridType in_gridtype,
                                             double in_b = 0);
 
private: // helper fns
  static double next_r_loglin(double b, double u, double r_guess);
  static std::vector<double> form_r(const GridType type, const double r0,
                                    const std::size_t num_points,
                                    const double du, const double b);
  static std::vector<double> form_drduor(const GridType type,
                                         const std::vector<double> &in_r,
                                         const double b);
  static std::vector<double> form_drdu(const GridType type,
                                       const std::vector<double> &in_r,
                                       const std::vector<double> &in_drduor);
};