path: root/src/patch/grid.hh

// 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:
//    <stdio.h>
//    <assert.h>
//    <math.h>		// 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	grid_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_array_pars& grid_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<fp> 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:
	//
	// ***** 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<fp>& rho_map() const { return rho_map_; }
	const linear_map<fp>& sigma_map() const { return sigma_map_; }
	const linear_map<fp>& 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<fp>::noninteger_action
				nia = linear_map<fp>::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<fp>::noninteger_action
				      nia = linear_map<fp>::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<fp>::noninteger_action
				nia = linear_map<fp>::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<fp>::noninteger_action
				  nia = linear_map<fp>::error)
		const
		{ return irho_of_rho(jtutil::radians_of_degrees(drho), nia); }
	int isigma_of_dsigma(fp dsigma, linear_map<fp>::noninteger_action
					nia = linear_map<fp>::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()); }


	//
	// ***** I/O, debugging support, and other misc stuff *****
	//

	// print a gridfn (via C stdio) in ASCII format
	// output is assumed to already be open
	// output format is suitable for gnuplot 'splot'
	void print_gridfn(int gfn, FILE *output_fp = stdout);

	// 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"; }




	//
	// ***** 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_array_pars& grid_array_pars_in,
	     const grid_pars& grid_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(const grid& rhs);
	grid& operator=(const grid& rhs);

private:
	// range of these is the bordered grid
	const linear_map<fp> rho_map_, sigma_map_;

	// angular boundaries of nominal grid
	const fp min_rho_, max_rho_;
	const fp min_sigma_, max_sigma_;
	};