// patch.cc -- describes a coordinate/grid patch
// $Id$

//
// patch::patch
// patch::~patch
// z_patch::z_patch
// x_patch::x_patch
// y_patch::y_patch
//
// patch::minmax_ang_ghost_zone
// patch::ghost_zone_on_edge
// patch::interpatch_ghost_zone_on_edge
// patch::create_mirror_symmetry_ghost_zone
// patch::create_periodic_symmetry_ghost_zone
// patch::create_interpatch_ghost_zone
// patch::set_ghost_zone
// patch::edge_adjacent_to_patch
// patch::assert_all_ghost_zones_fully_setup
//
// patch::print_gridfn
// patch::print_ghosted_gridfn
//
// patch_info::grid_array_pars
// patch_info::grid_pars
//

#include <cstdio>
#include <cmath>
#include <assert.h>
using std::fprintf;
using std::printf;

#include "cctk.h"

#include "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
using jtutil::error_exit;

#include "fp.hh"
#include "coords.hh"
#include "grid.hh"
#include "fd_grid.hh"
#include "patch.hh"
#include "patch_edge.hh"
#include "patch_interp.hh"
#include "ghost_zone.hh"

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

//
// This function constructs a  patch  object.
//
patch::patch(patch_system &my_patch_system_in, int patch_number_in,
	     const char name_in[], bool is_plus_in, char ctype_in,
	     local_coords::coords_set coords_set_rho_in,
	     local_coords::coords_set coords_set_sigma_in,
	     local_coords::coords_set coords_set_tau_in,
	     const grid_arrays::grid_array_pars& grid_array_pars_in,
	     const grid::grid_pars& grid_pars_in)

	: fd_grid(grid_array_pars_in, grid_pars_in),

	  my_patch_system_(my_patch_system_in),
	  patch_number_(patch_number_in),
	  name_(name_in),
	  is_plus_(is_plus_in), ctype_(ctype_in),

	  coords_set_rho_  (coords_set_rho_in  ),
	  coords_set_sigma_(coords_set_sigma_in),
	  coords_set_tau_  (coords_set_tau_in  ),

	  min_rho_patch_edge_(*new patch_edge(*this, side_is_min, side_is_rho)),
	  max_rho_patch_edge_(*new patch_edge(*this, side_is_max, side_is_rho)),
	  min_sigma_patch_edge_
	  	(*new patch_edge(*this, side_is_min, side_is_sigma)),
	  max_sigma_patch_edge_
		(*new patch_edge(*this, side_is_max, side_is_sigma)),

	  min_rho_ghost_zone_(NULL),
	  max_rho_ghost_zone_(NULL),
	  min_sigma_ghost_zone_(NULL),
	  max_sigma_ghost_zone_(NULL) // no comma

{ }

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

//
// This function destroys a  patch  object.
//
patch::~patch()
{
// no need to check for null pointers, since  delete NULL  is a silent no-op

delete max_sigma_ghost_zone_;
delete min_sigma_ghost_zone_;
delete max_rho_ghost_zone_;
delete min_rho_ghost_zone_;

delete & max_sigma_patch_edge_;
delete & min_sigma_patch_edge_;
delete & max_rho_patch_edge_;
delete & min_rho_patch_edge_;
}

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

//
// This function constructs a  z_patch  object.
//
z_patch::z_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char* name_in, bool is_plus_in,
		 const grid_arrays::grid_array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'z',
		local_coords::set_mu, local_coords::set_nu,
		local_coords::set_phi,
		grid_array_pars_in, grid_pars_in)
{ }

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

//
// This function constructs an  x_patch  object.
//
x_patch::x_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char* name_in, bool is_plus_in,
		 const grid_arrays::grid_array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'x',
		local_coords::set_nu, local_coords::set_phi,
		local_coords::set_mu,
		grid_array_pars_in, grid_pars_in)
{ }

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

//
// This function constructs a  y_patch  object.
//
y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char* name_in, bool is_plus_in,
		 const grid_arrays::grid_array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'y',
		local_coords::set_mu, local_coords::set_phi,
		local_coords::set_nu,
		grid_array_pars_in, grid_pars_in)
{ }

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

//
// This function returns a reference to the specified ghost zone
// of this patch.
//
ghost_zone& patch::minmax_ang_ghost_zone(bool want_min, bool want_rho)
	const
{
return want_min ? (want_rho ? min_rho_ghost_zone()
			    : min_sigma_ghost_zone())
		: (want_rho ? max_rho_ghost_zone()
			    : max_sigma_ghost_zone());
}

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

//
// This function returns a reference to the specified ghost zone
// of this patch.
//
ghost_zone& patch::ghost_zone_on_edge(const patch_edge &edge)
const
{
assert(& edge.my_patch() == this);
return minmax_ang_ghost_zone(edge.is_min(), edge.is_rho());
}

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

//
// This function verifies that the ghost zone on a specified edge
// is indeed interpatch, and returns a reference to it as an
//  interpatch_ghost_zone .
//
interpatch_ghost_zone& patch::interpatch_ghost_zone_on_edge
	(const patch_edge &e)
	const
{
ghost_zone &g = ghost_zone_on_edge(e);
assert(g.is_interpatch());
return static_cast<interpatch_ghost_zone &>(g);
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then constructs it as a mirror-symmetry
// ghost zone and properly links this to/from the patch.
//
void patch::create_mirror_symmetry_ghost_zone(const patch_edge& my_edge)
{
// make sure we belong to the right patch
assert(& my_edge.my_patch() == this);

symmetry_ghost_zone *temp = new symmetry_ghost_zone(my_edge);
set_ghost_zone(my_edge, temp);
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then creates it as a periodic-symmetry
// ghost zone and properly links this to/from the patch.
//
void patch::create_periodic_symmetry_ghost_zone
	(const patch_edge& my_edge, const patch_edge& symmetry_edge,
	 bool is_ipar_map_plus)
{
// make sure we belong to the right patch
assert(& my_edge.my_patch() == this);

int my_sample_ipar = my_edge.min_ipar_without_corners();
int symmetry_sample_ipar
	= is_ipar_map_plus ? symmetry_edge.min_ipar_without_corners()
			   : symmetry_edge.max_ipar_without_corners();

symmetry_ghost_zone *temp
	= new symmetry_ghost_zone(my_edge,        symmetry_edge,
				  my_sample_ipar, symmetry_sample_ipar,
				  is_ipar_map_plus);
set_ghost_zone(my_edge, temp);
}

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

//
// This function assert()s that a specified ghost zone of this patch
// hasn't already been set up, then creates it as an interpatch ghost
// zone (with lots of NULL pointers for info we can't compute yet)
// and properly links this to/from the patch.
//
void patch::create_interpatch_ghost_zone
	(const patch_edge& my_edge, const patch_edge& other_edge,
	 int N_overlap_points)
{
// make sure we belong to the right patch
assert(& my_edge.my_patch() == this);

interpatch_ghost_zone *temp
	= new interpatch_ghost_zone(my_edge, other_edge,
				    N_overlap_points);
set_ghost_zone(my_edge, temp);
}

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

//
// This is a helper function for  setup_*_ghost_zone().  This function
// assert()s that one of the ghost zone pointers (which one is selected
// by  edge ) is NULL, then stores a value in it.
//
void patch::set_ghost_zone(const patch_edge& edge, ghost_zone* gzp)
{
ghost_zone*& ghost_zone_ptr_to_set
	= edge.is_min()
	  ? (edge.is_rho() ? min_rho_ghost_zone_ : min_sigma_ghost_zone_)
	  : (edge.is_rho() ? max_rho_ghost_zone_ : max_sigma_ghost_zone_);

assert(ghost_zone_ptr_to_set == NULL);
ghost_zone_ptr_to_set = gzp;
}

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

//
// This function finds which patch edge is adjacent to a neighboring
// patch q, or does an error_exit() if q isn't actually a neighboring patch.
// The computation is done using only (rho,sigma) coordinate sets and
// min/max dang bounds ==> it's ok to use this function in setting up
// interpatch ghost zones.
//
// Arguments:
// q = The (supposedly) neighboring patch.
// N_overlap_points = The number of grid points these patches overlap.
//		      If this is nonzero, then these patches must have the
//		      same grid spacing in the perpendicular direction.
//
const patch_edge& patch::edge_adjacent_to_patch(const patch& q,
						int N_overlap_points = 0)
	const
{
const patch& p = *this;

// which (rho,sigma) coordinate do the patches have in common?
// ... this is the perp coordinate for the border
const local_coords::coords_set common_coord_set
	= p.coords_set_rho_sigma() & q.coords_set_rho_sigma();

// is this coordinate rho or sigma in each patch?
const bool common_is_p_rho   = (common_coord_set == p.coords_set_rho  ());
const bool common_is_p_sigma = (common_coord_set == p.coords_set_sigma());
if ((common_is_p_rho ^ common_is_p_sigma) != 0x1)
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        common coordinate isn't exactly one of p.{rho,sigma}!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_p_rho=%d common_is_p_sigma=%d\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_p_rho), int(common_is_p_sigma));
								/*NOTREACHED*/
const bool common_is_q_rho   = (common_coord_set == q.coords_set_rho  ());
const bool common_is_q_sigma = (common_coord_set == q.coords_set_sigma());
if ((common_is_q_rho ^ common_is_q_sigma) != 0x1)
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        common coordinate isn't exactly one of q.{rho,sigma}!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_q_rho=%d common_is_q_sigma=%d\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_q_rho), int(common_is_q_sigma));
								/*NOTREACHED*/

// how much do the patches overlap?
// ... eg N_overlap_points = 3 would be
//	p   p   p   p   p
//		q   q   q   q   q
//     so the overlap would be (N_overlap_points-1) * delta = 2 * delta
if ( (N_overlap_points-1 != 0)
     && jtutil::fuzzy<fp>::NE(p.delta_dang(common_is_p_rho),
			      q.delta_dang(common_is_q_rho)) )
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        N_overlap_points != 0 must have same perp grid spacing in both patches!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_p_rho=%d common_is_q_rho=%d\n"
"        p.delta_dang(common_is_p_rho)=%g\n"
"        q.delta_dang(common_is_q_rho)=%g\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_p_rho), int(common_is_q_rho),
		   double(p.delta_dang(common_is_p_rho)),
		   double(q.delta_dang(common_is_q_rho)));	/*NOTREACHED*/


const fp doverlap = fp(N_overlap_points-1) * p.delta_dang(common_is_p_rho);

// where is the common boundary relative to the min/max sides of each patch?
const bool common_is_p_min_q_max
    = local_coords::fuzzy_EQ_dang(p.min_dang(common_is_p_rho),
				  q.max_dang(common_is_q_rho) - doverlap);
const bool common_is_p_max_q_min
    = local_coords::fuzzy_EQ_dang(p.max_dang(common_is_p_rho),
				  q.min_dang(common_is_q_rho) + doverlap);
if ((common_is_p_min_q_max ^ common_is_p_max_q_min) != 0x1)
   then error_exit(ERROR_EXIT,
"***** patch::edge_adjacent_to_patch():\n"
"        common coordinate isn't exactly one of {pmax/qmin, pmin/qmax}!\n"
"        p.name()=\"%s\" q.name()=\"%s\"\n"
"        common_coord_set=%s\n"
"        common_is_p_rho=%d common_is_q_rho=%d\n"
"        p.delta_dang(common_is_p_rho)=%g\n"
"        q.delta_dang(common_is_q_rho)=%g\n"
"        N_overlap_points=%d doverlap=%g\n"
"        common_is_p_min_q_max=%d common_is_p_max_q_min=%d\n"
,
		   p.name(), q.name(),
		   local_coords::name_of_coords_set(common_coord_set),
		   int(common_is_p_rho), int(common_is_q_rho),
		   double(p.delta_dang(common_is_p_rho)),
		   double(q.delta_dang(common_is_q_rho)),
		   N_overlap_points, double(doverlap),
		   int(common_is_p_min_q_max), int(common_is_p_max_q_min));
								/*NOTREACHED*/

return p.minmax_ang_patch_edge(common_is_p_min_q_max, common_is_p_rho);
}

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

//
// This function verifies (via assert()) that all ghost zones of this
// patch have been fully set up.
//
void patch::assert_all_ghost_zones_fully_setup() const
{
assert(min_rho_ghost_zone_ != NULL);
assert(max_rho_ghost_zone_ != NULL);
assert(min_sigma_ghost_zone_ != NULL);
assert(max_sigma_ghost_zone_ != NULL);

// these calls are no-ops for non-interpatch ghost zones
min_rho_ghost_zone().assert_fully_setup();
max_rho_ghost_zone().assert_fully_setup();
min_sigma_ghost_zone().assert_fully_setup();
max_sigma_ghost_zone().assert_fully_setup();
}

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

//
// This function prints a nominal-grid gridfn in ASCII format to an
// already-open stdio stream.  The output format is suitable for a
// gnuplot 'splot' command.  Individual patches may be selected with
// the  select.patch  program (perl script).
//
void patch::print_gridfn(int gfn, FILE *output_fp = stdout)
	const
{
fprintf(output_fp, "### %s patch\n", name());
fprintf(output_fp, "# gfn=%d\n", gfn);
fprintf(output_fp, "# dpx = %s\n", name_of_dpx());
fprintf(output_fp, "# dpy = %s\n", name_of_dpy());
fprintf(output_fp, "#\n");
fprintf(output_fp, "# dpx\tdpy\tgridfn\tirho\tisigma\n");
	for (int irho = min_irho() ; irho <= max_irho() ; ++irho)
	{
		for (int isigma = min_isigma() ;
		     isigma <= max_isigma() ;
		     ++isigma)
		{
		const fp rho = rho_of_irho(irho);
		const fp sigma = sigma_of_isigma(isigma);
		const fp dpx = dpx_of_rho_sigma(rho, sigma);
		const fp dpy = dpy_of_rho_sigma(rho, sigma);
		fprintf(output_fp,
			"%g\t%g\t%.15g\n",
			dpx, dpy, gridfn(gfn,irho,isigma));
		}
	fprintf(output_fp, "\n");
	}
}

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

//
// This function prints a ghosted-grid gridfn in ASCII format to an
// already-open stdio stream.  The output format is suitable for a
// gnuplot 'splot' command.  Individual patches may be selected with
// the  select.patch  program (perl script).
//
void patch::print_ghosted_gridfn(int ghosted_gfn,
				 FILE *output_fp = stdout,
				 bool want_ghost_zones = true)
	const
{
fprintf(output_fp, "### %s patch\n", name());
fprintf(output_fp, "# ghosted_gfn=%d\n", ghosted_gfn);
fprintf(output_fp, "# dpx = %s\n", name_of_dpx());
fprintf(output_fp, "# dpy = %s\n", name_of_dpy());
fprintf(output_fp, "#\n");
fprintf(output_fp, "# dpx\tdpy\tgridfn\tirho\tisigma\n");
	for (int irho = effective_min_irho(want_ghost_zones) ;
	     irho <= effective_max_irho(want_ghost_zones) ;
	     ++irho)
	{
		for (int isigma = effective_min_isigma(want_ghost_zones) ;
		     isigma <= effective_max_isigma(want_ghost_zones) ;
		     ++isigma)
		{
		const fp rho = rho_of_irho(irho);
		const fp sigma = sigma_of_isigma(isigma);
		const fp dpx = dpx_of_rho_sigma(rho, sigma);
		const fp dpy = dpy_of_rho_sigma(rho, sigma);
		fprintf(output_fp,
			"%g\t%g\t%.15g\n",
			dpx, dpy, ghosted_gridfn(ghosted_gfn,irho,isigma));
		}
	fprintf(output_fp, "\n");
	}
}

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

//
// This function computes, and returns a reference to, a
//  struct grid_arrays::grid_array_pars  from the info in a
//  struct patch_info  and the additional information in the arguments.
//
// The result refers to an internal static buffer in this function; the
// usual caveats about lifetimes/overwriting apply.
//
// Arguments:
// N_ghost_points = Width in grid points of all ghost zones.
// N_extend_points = Number of grid points to extend each patch past
//		     "just touching" so as to overlap neighboring patches.
//		     Thus patches overlap by
//			N_overlap_points = 2*N_extend_points + 1
//		     grid points.  For example, with N_extend_points == 2,
//		     here are the grid points of two neighboring patches:
//			x   x   x   x   x   X   X
//                                      |
//			        O   O   o   o   o   o   o
//		     Here | marks the "just touching" boundary,
//		     x and o the grid points before this extension,
//		     and X and O the extra grid points added by this
//		     extension.
// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
//
const grid_arrays::grid_array_pars&
  patch_info::grid_array_pars(int N_ghost_points, int N_extend_points,
			      fp delta_drho_dsigma)
	const
{
static
  struct grid_arrays::grid_array_pars grid_array_pars_buffer;

grid_array_pars_buffer.min_irho
	= jtutil::round<fp>::to_integer(min_drho  /delta_drho_dsigma);
grid_array_pars_buffer.min_isigma
	= jtutil::round<fp>::to_integer(min_dsigma/delta_drho_dsigma);
grid_array_pars_buffer.max_irho
	= grid_array_pars_buffer.min_irho
	  + jtutil::round<fp>::to_integer(
		(max_drho  -min_drho  ) / delta_drho_dsigma
					 );
grid_array_pars_buffer.max_isigma
	= grid_array_pars_buffer.min_isigma
	  + jtutil::round<fp>::to_integer(
		(max_dsigma-min_dsigma) / delta_drho_dsigma
					 );
grid_array_pars_buffer.min_irho   -= N_extend_points;
grid_array_pars_buffer.min_isigma -= N_extend_points;
grid_array_pars_buffer.max_irho   += N_extend_points;
grid_array_pars_buffer.max_isigma += N_extend_points;

grid_array_pars_buffer.min_rho_N_ghost_points = N_ghost_points;
grid_array_pars_buffer.max_rho_N_ghost_points = N_ghost_points;
grid_array_pars_buffer.min_sigma_N_ghost_points = N_ghost_points;
grid_array_pars_buffer.max_sigma_N_ghost_points = N_ghost_points;

return grid_array_pars_buffer;
}

//******************************************************************************
//
//
// This function computes, and returns a reference to, a
//  struct grid_arrays::grid_pars  from the info in a  struct patch_info
// and the additional information in the arguments.
//
// The result refers to an internal static buffer in this function; the
// usual caveats about lifetimes/overwriting apply.
//
// Arguments:
// N_extend_points = Number of grid points to extend each patch past
//		     "just touching" so as to overlap neighboring patches.
//		     Thus patches overlap by  2*N_extend_points + 1  grid
//		     points.  For example, with N_extend_points == 2, here
//		     are the grid points of two neighboring patches:
//			x   x   x   x   x   X   X
//                                      |
//			        O   O   o   o   o   o   o
//		     Here | marks the "just touching" boundary,
//		     x and o the grid points before this extension,
//		     and X and O the extra grid points added by this
//		     extension.
// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
//
const grid::grid_pars& patch_info::grid_pars(int N_extend_points,
					     fp delta_drho_dsigma)
	const
{
static
  struct grid::grid_pars grid_pars_buffer;

const fp extend_drho_dsigma = fp(N_extend_points) * delta_drho_dsigma;

grid_pars_buffer.  min_drho   = min_drho   - extend_drho_dsigma;
grid_pars_buffer.delta_drho   = delta_drho_dsigma;
grid_pars_buffer.  max_drho   = max_drho   + extend_drho_dsigma;
grid_pars_buffer.  min_dsigma = min_dsigma - extend_drho_dsigma;
grid_pars_buffer.delta_dsigma = delta_drho_dsigma;
grid_pars_buffer.  max_dsigma = max_dsigma + extend_drho_dsigma;

return grid_pars_buffer;
}