// 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::decode_integration_method
// patch::integrate_gridfn
// patch::integration_coeff
// patch::ghost_zone_on_edge
// patch::corner_ghost_zone_containing_point
// patch::ghost_zone_containing_point
// 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
//

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

#include "cctk.h"

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

#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 decodes the character-string name of an integration method
// into an  enum integration_method .  See the comments in "patch.hh" on the
// declaration of  enum integration_method  for details on the methods and
// their character-string names.
//
//static
  enum patch::integration_method
    patch::decode_integration_method(const char method_string[])
{
if	(    STRING_EQUAL(method_string, "trapezoid")
	  || STRING_EQUAL(method_string, "trapezoid rule")    )
   then return integration_method__trapezoid;
else if (    STRING_EQUAL(method_string, "Simpson")
	  || STRING_EQUAL(method_string, "Simpson's rule")    )
   then return integration_method__Simpson;
else if (    STRING_EQUAL(method_string, "Simpson (variant)")
	  || STRING_EQUAL(method_string, "Simpson's rule (variant)")    )
   then return integration_method__Simpson_variant;
else	error_exit(ERROR_EXIT,
"***** patch::decode_integration_method():\n"
"        unknown method_string=\"%s\"!\n"
,
		   method_string);				/*NOTREACHED*/
}

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

//
// This function computes an approximation to the (surface) integral of
// a gridfn over the patch's nominal area,
//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
// where $f(\rho,\sigma)$ is the gridfn and $dA$ is the area element in
// $(\rho,sigma)$ coordinates.
//
// The integration scheme is selected based on the  method  argument.
//
// Bugs:
// The way the integration coefficients are computed is somewhat inefficient.
//
fp patch::integrate_gridfn(int src_gfn,
			   bool area_weighting_flag,
			   enum integration_method method)
	const
{
fp sum = 0.0;
	for (int irho = min_irho() ; irho <= max_irho() ; ++irho)
	{
	for (int isigma = min_isigma() ; isigma <= max_isigma() ; ++isigma)
	{
	sum += gridfn(src_gfn, irho,isigma)
	       * integration_coeff(method,
				   max_irho()-min_irho(),
				   irho      -min_irho())
	       * integration_coeff(method,
				   max_isigma()-min_isigma(),
				   isigma      -min_isigma());
	}
	}
return delta_rho() * delta_sigma() * sum;
}

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

//
// This function computes the integration coefficients for
//  integrate_gridfn() .  That is, if we write
//	$\int_{x_0}^{x_N} f(x) \, dx
//		\approx \Delta x \, \sum_{i=0}^N c_i f(x_i)$
// then this function computes $c_i$.
//
// Arguments:
// method = Specifies the integration method.
// N = The number of integration intervals.
// i = Specifies the point at which the coefficient is desired.
//
//static
  fp patch::integration_coeff(enum integration_method method, int N, int i)
{
assert(i >= 0);
assert(i <= N);

switch	(method)
	{
case integration_method__trapezoid:
	if ((i == 0) || (i == N))
	   then return 0.5;
	   else return 1.0;

case integration_method__Simpson:
	if ((N % 2) != 0)
	   then error_exit(ERROR_EXIT,
"***** patch::integration_coeff():\n"
"        Simpson's rule requires N to be even, but N=%d!\n",
			   N);					/*NOTREACHED*/
	if	((i == 0) || (i == N))
	   then return 1.0/3.0;
	else if ((i % 2) == 0)
	   then return 2.0/3.0;
	else	return 4.0/3.0;

case integration_method__Simpson_variant:
	if (N < 7)
	   then error_exit(ERROR_EXIT,
"***** patch::integration_coeff():\n"
"        Simpson's rule (variant) requires N >= 7, but N=%d!\n",
			   N);					/*NOTREACHED*/
	if	((i == 0) || (i == N))
	   then return 17.0/48.0;
	else if ((i == 1) || (i == N-1))
	   then return 59.0/48.0;
	else if ((i == 2) || (i == N-2))
	   then return 43.0/48.0;
	else if ((i == 3) || (i == N-3))
	   then return 49.0/48.0;
	else	return 1.0;

default:
	error_exit(ERROR_EXIT,
"***** patch::integration_coeff(): unknown method=(int)%d!\n"
"                                  (this should never happen!)\n"
,
		   int(method));				/*NOTREACHED*/
	}
}

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

//
// This function returns a reference to the ghost zone on a specified
// edge, after first assert()ing that the edge belongs to this patch.
//
// N.b. This function can't be inline in "patch.hh" because it needs
//	member functions of class patch_edge, which comes after class patch
//	in our #include order.
//
ghost_zone& patch::ghost_zone_on_edge(const patch_edge& e)
	const
{
assert(& e.my_patch() == this);
return minmax_ang_ghost_zone(e.is_min(), e.is_rho());
}

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

//
// This function determines which of the two adjacent ghost zones meeting
// at a specified corner, contains a specified point.  If the point isn't
// in either ghost zone, an abort(0) is done.
//
// Arguments:
// {rho,sigma}_is_min = Specify the corner (and implicitly the ghost zones).
// irho,isigma = Specify the point.
//
// Results:
// This function returns (a reference to) the desired ghost zone.
ghost_zone& patch::corner_ghost_zone_containing_point
	(bool rho_is_min, bool sigma_is_min,
	 int irho, int isigma)
	const
{
ghost_zone&   rho_gz =   minmax_rho_ghost_zone(  rho_is_min);
ghost_zone& sigma_gz = minmax_sigma_ghost_zone(sigma_is_min);

const patch_edge&   rho_edge =   rho_gz.my_edge();
const patch_edge& sigma_edge = sigma_gz.my_edge();

const int   rho_iperp =   rho_edge.iperp_of_irho_isigma(irho, isigma);
const int   rho_ipar  =   rho_edge. ipar_of_irho_isigma(irho, isigma);
const int sigma_iperp = sigma_edge.iperp_of_irho_isigma(irho, isigma);
const int sigma_ipar  = sigma_edge. ipar_of_irho_isigma(irho, isigma);

const bool is_in_rho_ghost_zone
	=   rho_gz.is_in_ghost_zone(  rho_iperp,   rho_ipar);
const bool is_in_sigma_ghost_zone
	= sigma_gz.is_in_ghost_zone(sigma_iperp, sigma_ipar);

// check that point is in exactly one ghost zone
assert(is_in_rho_ghost_zone ^ is_in_sigma_ghost_zone);

return is_in_rho_ghost_zone ? rho_gz : sigma_gz;
}

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

//
// This function determines which ghost zone contains a specified point.
// For a corner point between two symmetry ghost zones, it's unspecified
// which ghost zone will be chosen.
//
// If the point isn't in any ghost zone of this patch, an error_exit()
// is done.
//
// Arguments:
// irho,isigma = Specify the point.
//
// Results:
// This function returns (a reference to) the desired ghost zone.
ghost_zone& patch::ghost_zone_containing_point(int irho, int isigma)
	const
{
	// n.b. these loops must use _int_ variables for the loop
	//      to terminate!
	for (int want_min = false ; want_min <= true ; ++want_min)
	{
	for (int want_rho = false ; want_rho <= true ; ++want_rho)
	{
	const patch_edge& e = minmax_ang_patch_edge(want_min, want_rho);
	const int iperp = e.iperp_of_irho_isigma(irho, isigma);
	const int ipar  = e.ipar_of_irho_isigma (irho, isigma);

	ghost_zone& gz = minmax_ang_ghost_zone(want_min, want_rho);
	if (gz.is_in_ghost_zone(iperp, ipar))
	   then return gz;
	}
	}

error_exit(ERROR_EXIT,
"***** patch::ghost_zone_containing_point():\n"
"        no ghost zone contains point (this should never happen)!\n"
"        patch=%s irho=%d isigma=%d\n"
,
	   name(), irho, isigma);				/*NOTREACHED*/
}

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

//
// 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& other_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 other_sample_ipar = is_ipar_map_plus
			? other_edge.min_ipar_without_corners()
			: other_edge.max_ipar_without_corners();

symmetry_ghost_zone *temp
	= new symmetry_ghost_zone(my_edge,        other_edge,
				  my_sample_ipar, other_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();
}