path: root/src/patch/patch_system.hh

// patch_system.hh -- describes the (an) entire system of interlinked patches
// $Id$
//
// patch_system - describes a system of interlinked patches
//

//
// prerequisites:
//	<stdio.h>
//	<assert.h>
//	<math.h>
//	"stdc.h"
//	"config.hh"
//	"../jtutil/util.hh"
//	"../jtutil/array.hh"
//	"../jtutil/linear_map.hh"
//	"coords.hh"
//	"grid.hh"
//	"fd_grid.hh"
//	"patch.hh"
//	"patch_edge.hh"
//	"ghost_zone.hh"
//	"patch_interp.hh"
//	"patch_info.hh"
//

//******************************************************************************

//
// A  patch_system  object describes a system of interlinked patches.
// Its  const  qualifiers refer to the gridfn data.
//

class	patch_system
	{
	//
	// ***** static data & functions describing patch systems *****
	//
public:
	// what patch-system type are supported?
	enum patch_system_type {
			       full_sphere_patch_system,
			       plus_z_hemisphere_patch_system,
			       plus_xy_quadrant_patch_system,
			       plus_xz_quadrant_patch_system,
			       plus_xyz_octant_patch_system
			       };

	// maximum number of patches in any patch-system type
	static const int max_N_patches = 6;

	// decode patch system type into N_patches
	static
	  int N_patches_of_type(enum patch_system_type type_in);

	// patch system type <--> human-readable character-string name
	static
	  const char* name_of_type(enum patch_system_type type_in);
	static
	  enum patch_system_type type_of_name(const char* name_in);


	//
	// ***** coordinates *****
	//
public:
	// global (x,y,z) --> local (x,y,z)
	fp local_x_of_global_x(fp global_x) const
		{ return global_coords_.local_x_of_global_x(global_x); }
	fp local_y_of_global_y(fp global_y) const
		{ return global_coords_.local_y_of_global_y(global_y); }
	fp local_z_of_global_z(fp global_z) const
		{ return global_coords_.local_z_of_global_z(global_z); }

	// local (x,y,z) --> global (x,y,z)
	fp global_x_of_local_x(fp local_x) const
		{ return global_coords_.global_x_of_local_x(local_x); }
	fp global_y_of_local_y(fp local_y) const
		{ return global_coords_.global_y_of_local_y(local_y); }
	fp global_z_of_local_z(fp local_z) const
		{ return global_coords_.global_z_of_local_z(local_z); }

	// ... get global (x,y,z) coordinates of local origin point
	fp origin_x() const { return global_coords_.origin_x(); }
	fp origin_y() const { return global_coords_.origin_y(); }
	fp origin_z() const { return global_coords_.origin_z(); }

 
	//
	// ***** meta-info about the entire patch system *****
	//
public:
	// patch-system type
	enum patch_system_type type() const { return type_; }

	// total number of patches
	int N_patches() const { return N_patches_; }

	// get patches by patch number
	patch &ith_patch(int pn) const
		{ return * all_patches_[pn]; }

	// find a patch by name, return patch number; error_exit() if not found
	int patch_number_of_name(const char* name) const;

	// total number of grid points
	int N_grid_points() const { return N_grid_points_; }
	int ghosted_N_grid_points() const { return ghosted_N_grid_points_; }

	// convert (patch,irho,isigma) <--> 1-D 0-origin grid point number (gpn)
	int gpn_of_patch_irho_isigma(const patch& p, int irho, int isigma)
		const
		{
		return starting_gpn_[p.patch_number()]
		       + p.gpn_of_irho_isigma(irho,isigma);
		}
	int ghosted_gpn_of_patch_irho_isigma(const patch& p,
					     int irho, int isigma)
		const
		{
		return ghosted_starting_gpn_[p.patch_number()]
		       + p.ghosted_gpn_of_irho_isigma(irho,isigma);
		}
	// ... n.b. we return patch as a reference via the function result;
	//     an alternative would be to have a patch*& argument
	patch& patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
		const;
	patch& ghosted_patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
		const;

	//
	// ***** meta-info about gridfns *****
	//
public:
	int min_gfn() const { return ith_patch(0).min_gfn(); }
	int max_gfn() const { return ith_patch(0).max_gfn(); }
	int N_gridfns() const { return ith_patch(0).N_gridfns(); }
	int ghosted_min_gfn() const { return ith_patch(0).ghosted_min_gfn(); }
	int ghosted_max_gfn() const { return ith_patch(0).ghosted_max_gfn(); }
	int ghosted_N_gridfns() const
		{ return ith_patch(0).ghosted_N_gridfns(); }


	//
	// ***** access to gridfns as 1-D arrays[gpn] *****
	// ... n.b. this interface implicitly assumes that gridfn data
	//     arrays are contiguous across patches; this is ensured by
	//     setup_gridfn_storage() (called by our constructor)
	//
public:
	const fp* gridfn_data(int gfn) const
		{ return ith_patch(0).gridfn_data_array(gfn); }
	      fp* gridfn_data(int gfn)
		{ return ith_patch(0).gridfn_data_array(gfn); }
	const fp* ghosted_gridfn_data(int ghosted_gfn) const
		{ return ith_patch(0).ghosted_gridfn_data_array(ghosted_gfn); }
	      fp* ghosted_gridfn_data(int ghosted_gfn)
		{ return ith_patch(0).ghosted_gridfn_data_array(ghosted_gfn); }


	//
	// ***** I/O *****
	//
public:
	// print to a named file (newly (re)created)
	// output format is
	//	dpx	dpy	gridfn
	void print_gridfn(int gfn, const char output_file_name[]) const
		{
		print_unknown_gridfn(false, gfn,
				     false, false, 0,
				     output_file_name, false);
		}
	void print_ghosted_gridfn(int ghosted_gfn,
				  const char output_file_name[],
				  bool want_ghost_zones = true)
		const
		{
		print_unknown_gridfn(true, ghosted_gfn,
				     false, false, 0,
				     output_file_name, want_ghost_zones);
		}

	// print to a named file (newly (re)created)
	// output format is
	//	dpx	dpy	gridfn   global_x   global_y   global_z
	// where global_[xyz} are derived from the angular position
	// and a specified (unknown-grid) radius gridfn
	void print_gridfn_with_xyz
		(int gfn,
		 bool radius_is_ghosted_flag, int unknown_radius_gfn,
		 const char output_file_name[])
		const
		{
		print_unknown_gridfn(false, gfn,
				     true, radius_is_ghosted_flag,
					   unknown_radius_gfn,
				     output_file_name, false);
		}
	void print_ghosted_gridfn_with_xyz
		(int ghosted_gfn,
		 bool radius_is_ghosted_flag, int unknown_radius_gfn,
		 const char output_file_name[],
		 bool want_ghost_zones = true)
		const
		{
		print_unknown_gridfn(true, ghosted_gfn,
				     true, radius_is_ghosted_flag,
					   unknown_radius_gfn,
				     output_file_name, want_ghost_zones);
		}

public:
	// read from a named file
	void read_gridfn(int gfn, const char input_file_name[])
		{ read_unknown_gridfn(false, gfn, input_file_name, false); }
	void read_ghosted_gridfn(int ghosted_gfn,
				 const char input_file_name[],
				 bool want_ghost_zones = true)
		{
		read_unknown_gridfn(true, ghosted_gfn,
				    input_file_name, want_ghost_zones);
		}

private:
	// ... internal worker functions
	void print_unknown_gridfn
		(bool ghosted_flag, int unknown_gfn,
		 bool print_xyz_flag, bool radius_is_ghosted_flag,
				      int unknown_radius_gfn,
		 const char output_file_name[], bool want_ghost_zones)
		const;
	void read_unknown_gridfn(bool ghosted_flag, int unknown_gfn,
				 const char input_file_name[],
				 bool want_ghost_zones);


	//
	// ***** ghost zone operations *****
	//
public:
	// "synchronize" all ghost zones of all patches,
	// i.e. update the ghost-zone values of the specified gridfns
	// via the appropriate sequence of symmetry operations
	// and interpatch interpolations
	void synchronize(int ghosted_min_gfn_to_sync,
			 int ghosted_max_gfn_to_sync);

	// ... do this for all ghosted gridfns
	void synchronize()
		{ synchronize(ghosted_min_gfn(), ghosted_max_gfn()); }


	//
	// compute Jacobian of  synchronize()  into internal buffers,
	// taking into account synchronize()'s full 3-phase algorithm
	//
	void compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync,
					  int ghosted_max_gfn_to_sync)
		const;

	// ... do this for all ghosted gridfns
	void compute_synchronize_Jacobian()
		const
		{
		compute_synchronize_Jacobian(ghosted_min_gfn(),
					     ghosted_max_gfn());
		}

	//
	// The following functions access the Jacobian computed by
	//  compute_synchronize_Jacobian() .  Note this API has the
	// same implicit assumptions on the Jacobian structure documented
	// in the comments in "ghost_zone.hh" immediately following
	//  ghost_zone::compute_Jacobian() .
	//

	// to which patch/edge do the y points in a Jacobian row belong?
	patch& synchronize_Jacobian_y_patch(const ghost_zone& xgz)
		const
		{ return xgz.Jacobian_y_patch(); }
	const patch_edge& synchronize_Jacobian_y_edge (const ghost_zone& xgz)
		const
		{ return xgz.Jacobian_y_edge(); }

	// what is the [min,max] range of m for a Jacobian row?
	int synchronize_Jacobian_min_y_ipar_m(const ghost_zone& xgz)
		const
		{ return xgz.Jacobian_min_y_ipar_m(); }
	int synchronize_Jacobian_max_y_ipar_m(const ghost_zone& xgz)
		const
		{ return xgz.Jacobian_max_y_ipar_m(); }

	// what is the iperp of the Jacobian y points in their (y) patch?
	int synchronize_Jacobian_y_iperp(const ghost_zone& xgz, int x_iperp)
		const
		{ return xgz.Jacobian_y_iperp(x_iperp); }

	// what is the  posn  value of the y points in this Jacobian row?
	int synchronize_Jacobian_y_ipar_posn(const ghost_zone& xgz,
					     int x_iperp, int x_ipar)
		const;

	// what is the Jacobian
	//	partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar)
	//	-------------------------------------------------------------
	//	   partial gridfn(ghosted_gfn, py, y_iperp, y_posn+y_ipar_m)
	// where
	//	y_iperp = Jacobian_y_iperp(x_iperp)
	//	y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar)
	// taking into account synchronize()'s full 3-phase algorithm
	fp synchronize_Jacobian(const ghost_zone& xgz, int x_iperp, int x_ipar,
				int y_ipar_m)
		const;


	//
	// ***** constructor, destructor *****
	//
	// This constructor doesn't support the full generality of the
	// patch data structures (which would, eg, allow N_ghost_points
	// and N_extend_points and the interpolator parameters to vary
	// from ghost zone to ghost zone, and the grid spacings to vary
	// from patch to patch.  But in practice we'd probably never
	// use that generality...
	//
public:
	patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
		     enum patch_system_type type_in,
		     int N_ghost_points, int N_overlap_points,
		     fp delta_drho_dsigma,
		     int min_gfn_in, int max_gfn_in,
		     int ghosted_min_gfn_in, int ghosted_max_gfn_in,
		     int interp_handle_in, int interp_par_table_handle_in);
	~patch_system();


	//
	// ***** helper functions for constructor *****
	//
private:
	// construct patches as described by patch_info[] array,
	// and link them into the patch system
	// does *NOT* create ghost zones
	// does *NOT* set up gridfns
	void create_patches(const struct patch_info patch_info_in[],
			    int N_ghost_points, int N_extend_points,
			    fp delta_drho_dsigma);

	// setup all gridfns with contiguous-across-patches storage
	void setup_gridfn_storage
		(int min_gfn_in, int max_gfn_in,
		 int ghosted_min_gfn_in, int ghosted_max_gfn_in);

	// create/interlink ghost zones
	void interlink_full_sphere_patch_system
		(int N_overlap_points,
		 int interp_handle, int interp_par_table_handle);
	void interlink_plus_z_hemisphere_patch_system
		(int N_overlap_points,
		 int interp_handle, int interp_par_table_handle);
	void interlink_plus_xy_quadrant_patch_system
		(int N_overlap_points,
		 int interp_handle, int interp_par_table_handle);
	void interlink_plus_xz_quadrant_patch_system
		(int N_overlap_points,
		 int interp_handle, int interp_par_table_handle);
	void interlink_plus_xyz_octant_patch_system
		(int N_overlap_points,
		 int interp_handle, int interp_par_table_handle);

	// create/interlink a pair of periodic-symmetry ghost zones
	static
	  void create_periodic_symmetry_ghost_zones
		(const patch_edge& ex, const patch_edge& ey,
		 bool ipar_map_is_plus);

	// construct a pair of interpatch ghost zones
	// ... automagically figures out which edges are adjacent
	static
	  void create_interpatch_ghost_zones(patch &px, patch &py,
					     int N_overlap_points);

	// finish setup of a pair of interpatch ghost zones
	// ... automagically figures out which edges are adjacent
	static
	  void finish_interpatch_setup
		(patch &px, patch &py,
		 int N_overlap_points,
		 int interp_handle, int interp_par_table_handle);

	// assert() that all ghost zones of all patches are fully setup
	void assert_all_ghost_zones_fully_setup() const;


private:
        // we forbid copying and passing by value
        // by declaring the copy constructor and assignment operator
        // private, but never defining them
	patch_system(const patch_system &rhs);
	patch_system& operator=(const patch_system &rhs);

private:
	// local <--> global coordinate mapping
	global_coords global_coords_;

	// meta-info about patch system
	enum patch_system_type type_;
	int N_patches_;
	int N_grid_points_, ghosted_N_grid_points_;

	// [pn] = --> individual patches
	std::vector<patch *> all_patches_;

	// [pn] = starting grid point number of individual patches
	// ... arrays are actually of size N_patches_+1, the [N_patches_]
	//     entries are == N_grid_points_ and ghosted_N_grid_points_
	std::vector<int> starting_gpn_;
	std::vector<int> ghosted_starting_gpn_;

	// pointers to storage blocks for all gridfns
	// ... patches point into these, but we own the storage blocks
	fp* gridfn_storage_;
	fp* ghosted_gridfn_storage_;
	};