// grid.hh -- classes for a 2D uniform tensor-product grid // $Id$ // // grid_arrays - data arrays for a 2D tensor-product grid // grid - uniform 2D tensor-product grid // // // prerequisites: // // // // for M_PI (used by degree/radian conversions) // "jt/util++.hh" // jtutil:: stuff: // // how_many_in_range(), // // degrees_of_radians(), radians_of_degrees(), // "jt/array.hh" // "jt/linear_map.hh" // fp.hh // coords.hh // //***************************************************************************** // // grid_arrays - data arrays for a 2D tensor-product grid // // This is a helper class for class grid (below). This class stores // most of the actual grid function (gridfn) data arrays for a uniform // tensor-product 2D grid. // // The integer grid coordinates are (irho,isigma). This class deals // with the grid solely at the level of arrays with integer subscripts; // the derived class grid deals with the floating-point coordinates // related to those subscripts. // // The grid has a nominal extent, surrounded by a "ghost zone" for // finite differencing purposes. For simplicity, all grid functions // (gridfns) are stored on the full grid (nominal + ghost zone), even // though we actually only need this for finite differencing targets. // // We identify a gridfn by a small-integer "grid function number", // a.k.a. "gfn". Each gridfn is stored contiguously. // class grid_arrays { public: // // ***** {min,max}_{rho,sigma} "sides" of grid ***** // // // A grid has 4 (angular) "sides", which we identify as // {min,max}_{rho,sigma}. Given a side, we define coordinates // (perpendicular,parallel) to it, normally abbreviated to // (perp,par). // // As well as functions directly referring to a specific side, // we also support referring to one of these chosen at run-time, // via Boolean flags: // // // generic (irho,isigma) coordinate // iang = want_rho ? irho : isigma // // // opposite (irho,isigma) coordinate // ixang = want_rho ? isigma : irho // // // generic (min,max) direction // minmax = want_min ? min : max // // FIXME: This system of Boolean flags works ok, but it requires // a lot of repetitive code conditional-expression functions // in this class. Is there a cleaner solution? // there are precisely this many possible sides static const int N_possible_sides = 4; // we specify {min,max} with a Boolean want_min static const bool side_min = true; // for want_min static const bool side_max = false; // for want_min // we specify {rho,sigma} with a Boolean want_rho static const bool side_rho = true; // for want_rho static const bool side_sigma = false; // for want_rho // human-readable names for the sides (for debugging) static const char *minmax_name(bool minmax) { return minmax ? "min" : "max"; } static const char *iang_name(bool want_rho) { return want_rho ? "irho" : "isigma"; } // // ***** gfn-checking and gridfn access functions ***** // int N_gridfns() const { return N_gridfns_; } bool is_valid_gfn(int gfn) const { return (gfn >= 0) && (gfn < N_gridfns()); } // access to gridfn data // ... rvalue (may be slightly faster since it's a const function // and it returns by value rather than reference fp gridfn(int gfn, int irho, int isigma) const { return gridfn_data_(gfn, irho, isigma); } // ... lvalue (must return a reference) fp& gridfn(int gfn, int irho, int isigma) { return gridfn_data_(gfn, irho, isigma); } // // ***** array info ***** // // nominal-grid min/max/sizes int min_irho() const { return min_irho_; } int max_irho() const { return max_irho_; } int min_isigma() const { return min_isigma_; } int max_isigma() const { return max_isigma_; } int minmax_iang(bool want_min, bool want_rho) const { return want_min ? (want_rho ? min_irho() : min_isigma()) : (want_rho ? max_irho() : max_isigma()); } int N_irho() const { return jtutil::how_many_in_range(min_irho(), max_irho()); } int N_isigma() const { return jtutil::how_many_in_range(min_isigma(), max_isigma()); } int N_grid_points() const { return N_irho() * N_isigma(); } // full-grid min/max (don't need sizes) int full_grid__min_irho() const { return full_grid__min_irho_; } int full_grid__max_irho() const { return full_grid__max_irho_; } int full_grid__min_isigma() const { return full_grid__min_isigma_; } int full_grid__max_isigma() const { return full_grid__max_isigma_; } int full_grid__minmax_iang(bool want_min, bool want_rho) const { return want_min ? (want_rho ? full_grid__min_irho() : full_grid__min_isigma()) : (want_rho ? full_grid__max_irho() : full_grid__max_isigma()); } // // ***** ghost zones ***** // // ghost zone min/max perpendicular coordinates int min_rho_ghost_zone__min_iperp() const { return full_grid__min_irho(); } int min_rho_ghost_zone__max_iperp() const { return min_irho() - 1; } int max_rho_ghost_zone__min_iperp() const { return max_irho() + 1; } int max_rho_ghost_zone__max_iperp() const { return full_grid__max_irho(); } int min_sigma_ghost_zone__min_iperp() const { return full_grid__min_isigma(); } int min_sigma_ghost_zone__max_iperp() const { return min_isigma() - 1; } int max_sigma_ghost_zone__min_iperp() const { return max_isigma() + 1; } int max_sigma_ghost_zone__max_iperp() const { return full_grid__max_isigma(); } int minmax_ang_ghost_zone__min_iperp(bool want_min, bool want_rho) const { return want_min ? (want_rho ? min_rho_ghost_zone__min_iperp() : min_sigma_ghost_zone__min_iperp()) : (want_rho ? max_rho_ghost_zone__min_iperp() : max_sigma_ghost_zone__min_iperp()); } int minmax_ang_ghost_zone__max_iperp(bool want_min, bool want_rho) const { return want_min ? (want_rho ? min_rho_ghost_zone__max_iperp() : min_sigma_ghost_zone__max_iperp()) : (want_rho ? max_rho_ghost_zone__max_iperp() : max_sigma_ghost_zone__max_iperp()); } // ghost zone min/max parallel coordinates // ... not including corners int rho_ghost_zone_without_corners__min_ipar() const { return min_isigma(); } int rho_ghost_zone_without_corners__max_ipar() const { return max_isigma(); } int sigma_ghost_zone_without_corners__min_ipar() const { return min_irho(); } int sigma_ghost_zone_without_corners__max_ipar() const { return max_irho(); } int ang_ghost_zone_without_corners__min_ipar(bool want_rho) const { return want_rho ? rho_ghost_zone_without_corners__min_ipar() : sigma_ghost_zone_without_corners__min_ipar(); } int ang_ghost_zone_without_corners__max_ipar(bool want_rho) const { return want_rho ? rho_ghost_zone_without_corners__max_ipar() : sigma_ghost_zone_without_corners__max_ipar(); } // ... including corners int rho_ghost_zone_with_corners__min_ipar() const { return full_grid__min_isigma(); } int rho_ghost_zone_with_corners__max_ipar() const { return full_grid__max_isigma(); } int sigma_ghost_zone_with_corners__min_ipar() const { return full_grid__min_irho(); } int sigma_ghost_zone_with_corners__max_ipar() const { return full_grid__max_irho(); } int ang_ghost_zone_with_corners__min_ipar(bool want_rho) const { return want_rho ? rho_ghost_zone_with_corners__min_ipar() : sigma_ghost_zone_with_corners__min_ipar(); } int ang_ghost_zone_with_corners__max_ipar(bool want_rho) const { return want_rho ? rho_ghost_zone_with_corners__max_ipar() : sigma_ghost_zone_with_corners__max_ipar(); } // // ***** membership predicates ***** // bool is_in_nominal_grid(int irho, int isigma) const { return (irho >= min_irho()) && (irho <= max_irho()) && (isigma >= min_isigma()) && (isigma <= max_isigma()); } bool is_in_full_grid(int irho, int isigma) const { return (irho >= full_grid__min_irho()) && (irho <= full_grid__max_irho()) && (isigma >= full_grid__min_isigma()) && (isigma <= full_grid__max_isigma()); } bool is_in_ghost_zone(int irho, int isigma) const { return is_in_full_grid(irho, isigma) && !is_in_nominal_grid(irho, isigma); } // // ***** argument structure for constructor ***** // // this structure bundles related arguments together so we don't // have 20+ (!) separate arguments to our top-level constructors struct array_pars { int min_irho, max_irho; int min_isigma, max_isigma; int min_rho_N_ghost_zones, max_rho_N_ghost_zones; int min_sigma_N_ghost_zones, max_sigma_N_ghost_zones; }; // // ***** constructor, destructor ***** // grid_arrays(const grid_arrays::array_pars &array_pars_in, int N_gridfns_in); // compiler-generated default destructor is ok private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator // private, but never defining them grid_arrays(const grid_arrays &rhs); grid_arrays &operator=(const grid_arrays &rhs); private: // // ***** the actual gridfn storage arrays ***** // // n.b. this array is *first* data member in this class // ==> possibly slightly faster access (0 offset from pointer) array3d gridfn_data_; // indices are (gfn, irho, isigma) // number of gridfns const int N_gridfns_; // nominal grid min/max bounds const int min_irho_, max_irho_; const int min_isigma_, max_isigma_; // full grid min/max bounds const int full_grid__min_irho_, full_grid__max_irho_; const int full_grid__min_isigma_, full_grid__max_isigma_; }; //****************************************************************************** // // grid - uniform 2D tensor-product grid // // The grid is uniform in the floating point grid coordinates (rho,sigma). // There is also some (limited) support for expressing these coordinates // in degrees (drho,dsigma); this is useful for humans trying to specify // things in parameter files. // // The nominal (non-bordered) angular grid boundaries may coincide with // grid points, or they may be at "half-integer" grid coordinates. That // is, suppose we have a unit grid spacing, and a boundary at an angular // coordinate of 0; then the grid may be either 0, 1, 2, ..., or // 0.5, 1.5, 2.5, ... . // class grid : public grid_arrays { public: // human-readable names for the sides (for debugging) static const char *ang_name(bool want_rho) { return want_rho ? "rho" : "sigma"; } static const char *dang_name(bool want_rho) { return want_rho ? "drho" : "dsigma"; } // // ***** low-level access to coordinate maps ***** // // direct (read-only) access to the underlying linear_map objects // ... useful for (eg) passing to interpolators const linear_map &rho_map() const { return rho_map_; } const linear_map &sigma_map() const { return sigma_map_; } const linear_map &ang_map(bool want_rho) const { return want_rho ? rho_map() : sigma_map(); } // // ***** single-axis coordinate conversions ***** // // ... angules in radians fp rho_of_irho(int irho) const { return rho_map().fp_of_int(irho); } fp sigma_of_isigma(int isigma) const { return sigma_map().fp_of_int(isigma); } fp ang_of_iang(bool want_rho, int iang) const { return want_rho ? rho_of_irho(iang) : sigma_of_isigma(iang); } fp fp_irho_of_rho(fp rho) const { return rho_map().fp_int_of_fp(rho); } int irho_of_rho(fp rho, linear_map::noninteger_action nia = linear_map::error) const { return rho_map().int_of_fp(rho, nia); } fp fp_isigma_of_sigma(fp sigma) const { return sigma_map().fp_int_of_fp(sigma); } int isigma_of_sigma(fp sigma, linear_map::noninteger_action nia = linear_map::error) const { return sigma_map().int_of_fp(sigma, nia); } fp fp_iang_of_ang(bool want_rho, fp ang) const { return want_rho ? fp_irho_of_rho(ang) : fp_isigma_of_sigma(ang); } int iang_of_ang(bool want_rho, fp ang, linear_map::noninteger_action nia = linear_map::error) const { return want_rho ? irho_of_rho(ang, nia) : isigma_of_sigma(ang, nia); } // ... angles in degrees fp drho_of_irho(int irho) const { return jtutil::degrees_of_radians(rho_of_irho(irho)); } fp dsigma_of_isigma(int isigma) const { return jtutil::degrees_of_radians(sigma_of_isigma(isigma)); } int irho_of_drho(fp drho, linear_map::noninteger_action nia = linear_map::error) const { return irho_of_rho(jtutil::radians_of_degrees(drho), nia); } int isigma_of_dsigma(fp dsigma, linear_map::noninteger_action nia = linear_map::error) const { return isigma_of_sigma(jtutil::radians_of_degrees(dsigma), nia); } // // ***** grid info ***** // // grid spacings fp delta_rho() const { return rho_map().delta_fp(); } fp delta_sigma() const { return sigma_map().delta_fp(); } fp delta_drho() const { return jtutil::degrees_of_radians(delta_rho()); } fp delta_dsigma() const { return jtutil::degrees_of_radians(delta_sigma()); } // inverse grid spacings fp inverse_delta_rho() const { return rho_map().inverse_delta_fp(); } fp inverse_delta_sigma() const { return sigma_map().inverse_delta_fp(); } // nominal grid min/max fp min_rho() const { return min_rho_; } fp max_rho() const { return max_rho_; } fp min_sigma() const { return min_sigma_; } fp max_sigma() const { return max_sigma_; } fp min_drho() const { return jtutil::degrees_of_radians(min_rho()); } fp max_drho() const { return jtutil::degrees_of_radians(max_rho()); } fp min_dsigma() const { return jtutil::degrees_of_radians(min_sigma()); } fp max_dsigma() const { return jtutil::degrees_of_radians(max_sigma()); } // full-grid min/max fp full_grid__min_rho() const { return rho_of_irho(full_grid__min_irho()); } fp full_grid__max_rho() const { return rho_of_irho(full_grid__max_irho()); } fp full_grid__min_sigma() const { return sigma_of_isigma(full_grid__min_isigma()); } fp full_grid__max_sigma() const { return sigma_of_isigma(full_grid__max_isigma()); } // // ***** argument structure for constructor ***** // // this structure bundles related arguments together so we don't // have 20+ (!) separate arguments to our top-level constructors struct grid_pars // *** note angles in degrees *** { fp min_drho, delta_drho, max_drho; fp min_dsigma, delta_dsigma, max_dsigma; }; // // ***** constructor, destructor ***** // grid(const grid_arrays::array_pars &grid_array_pars, const grid::grid_pars &grid_pars_in, int N_gridfns_in, bool verbose_flag = false); // compiler-generated default destructor is ok private: // we forbid copying and passing by value // by declaring the copy constructor and assignment operator // private, but never defining them grid(const grid &rhs); grid &operator=(const grid &rhs); private: // range of these is the bordered grid const linear_map rho_map_, sigma_map_; // angular boundaries of nominal grid const fp min_rho_, max_rho_; const fp min_sigma_, max_sigma_; };