// ghost_zone.cc -- fill in gridfn data in patch ghost zones
// $Id$

//
// symmetry_ghost_zone::symmetry_ghost_zone (mirror symmetry)
// symmetry_ghost_zone::symmetry_ghost_zone (periodic BC)
// symmetry_ghost_zone::~symmetry_ghost_zone
// symmetry_ghost_zone::extend_scalar_gridfn_to_ghost_zone
//
// interpatch_ghost_zone::interpatch_ghost_zone
// interpatch_ghost_zone::~interpatch_ghost_zone
// interpatch_ghost_zone::assert_fully_setup
// interpatch_ghost_zone::extend_scalar_gridfn_to_ghost_zone
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
#include "jt/interpolate.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 "ghost_zone.hh"
#include "patch_frontier.hh"

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

//
// This function constructs a mirror-symmetry ghost zone object
//
symmetry_ghost_zone::symmetry_ghost_zone(const patch_edge& my_edge_in)
	: ghost_zone(my_edge_in, ghost_zone_is_symmetry),
	  symmetry_patch_(my_edge_in.my_patch()),
	  symmetry_edge_(my_edge_in)
{
// iperp_map: i --> (i of ghost zone) - i
iperp_map_ = new cpm_map<fp>(min_iperp(), max_iperp(),
			     my_edge_in.fp_grid_outer_iperp());

// ipar_map_: identity map
ipar_map_ = new cpm_map<fp>(min_ipar(), max_ipar());
}

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

//
// This function constructs a periodic-symmetry ghost zone object.
//
symmetry_ghost_zone::symmetry_ghost_zone
	(const patch_edge& my_edge_in, const patch_edge& symmetry_edge_in,
	 int my_edge_sample_ipar,      int symmetry_edge_sample_ipar,
	 bool ipar_map_is_plus)
	: ghost_zone(my_edge_in, ghost_zone_is_symmetry),
	  symmetry_patch_(symmetry_edge_in.my_patch()),
	  symmetry_edge_(symmetry_edge_in)
{
//
// perpendicular map
//
const fp fp_my_period_plane_iperp = my_edge().fp_grid_outer_iperp();
const fp fp_symmetry_period_plane_iperp = symmetry_edge().fp_grid_outer_iperp();

// iperp mapping must be outside --> inside
// i.e. if both edges have iperp as the same min/max "direction",
//	then the mapping is  iperp increasing --> iperp decreasing
//      (i.e. the map's sign is -1)
const bool iperp_map_is_plus
	= ! (my_edge().is_min() == symmetry_edge().is_min());
iperp_map_ = new cpm_map<fp>(min_iperp(), max_iperp(),
			     fp_my_period_plane_iperp,
				     fp_symmetry_period_plane_iperp,
			     iperp_map_is_plus);

//
// parallel map
//
ipar_map_ = new cpm_map<fp>(min_ipar(), max_ipar(),
			    my_edge_sample_ipar, symmetry_edge_sample_ipar,
			    ipar_map_is_plus);
}

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

//
// This function destroys a  symmetry_ghost_zone  object.
//
symmetry_ghost_zone::~symmetry_ghost_zone()
{
delete ipar_map_;
delete iperp_map_;
}

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

//
// This function symmetry-maps a scalar gridfn in a specified part
// of the ghost zone.
//
// N.b. "scalar" here refers solely to the symmetry properties of the
// gridfn under the ghost zone's "const +/-" coordinate transformations;
// it need not be an actual scalar under more general coordinate
// transformations.
//
void symmetry_ghost_zone::extend_scalar_gridfn_to_ghost_zone
	(int gfn,
	 bool want_min_par_corner,
	 bool want_non_corner,
	 bool want_max_par_corner)
	const
{
	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	for (int ipar = min_ipar() ; ipar <= max_ipar() ; ++ipar)
	{
	// do we want to do this point?
	if (! my_edge().ipar_is_in_selected_corner(want_min_par_corner,
						   want_non_corner,
						   want_max_par_corner,
						   ipar) )
	   then continue;				// *** LOOP CONTROL ***

	const int sym_iperp = iperp_map_of_iperp(iperp);
	const int sym_ipar  = ipar_map_of_ipar  (ipar );
	const int sym_irho
		= symmetry_edge().irho_of_iperp_ipar  (sym_iperp,sym_ipar);
	const int sym_isigma
		= symmetry_edge().isigma_of_iperp_ipar(sym_iperp,sym_ipar);
	const fp sym_gridfn = symmetry_patch().gridfn(gfn, sym_irho,sym_isigma);

	const int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
	const int isigma = my_edge().isigma_of_iperp_ipar(iperp,ipar);
	my_patch().gridfn(gfn, irho,isigma) = sym_gridfn;
	}
	}
}

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

//
// This function constructs an  interpatch_ghost_zone  object.
//
interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in,
					     const patch_edge& other_edge_in,
					     int N_overlap_points)
	: ghost_zone(my_edge_in, ghost_zone_is_interpatch),
	  other_patch_(other_edge_in.my_patch()),
	  other_edge_(other_edge_in),
	  other_frontier_(NULL),	// other frontier object
					// may not exist yet
	  other_iperp_(NULL), other_par_(NULL)	// initialized in ctor body
{
//
// verify that we have the expected relationships between
// this and the other patch's (mu,nu,phi) coordinates:
//

// perp coordinate is common to us and the other patch, so
// ghost zone/frontier must be min in one patch, max in the other
if (my_edge().is_min() == other_edge().is_min())
   then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        ghost zone/frontier must be min in one patch, max in the other!\n"
,
		   my_patch().name(), my_edge().name(),
		   other_patch().name(), other_edge().name());	/*NOTREACHED*/

// coord in common between the two patches must be perp coord in both patches
// and this patch's tau coordinate must be other edge's parallel coordinate
const local_coords::coords_set common_coords_set
	= local_coords::coords_set_not(my_patch().coords_set_rho_sigma()
				       ^
				       other_patch().coords_set_rho_sigma());
if (! (    (common_coords_set == my_edge().coords_set_perp())
	&& (common_coords_set == other_edge().coords_set_perp())
	&& (my_patch().coords_set_tau() == other_edge().coords_set_par())    ) )
   then error_exit(PANIC_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        (rho,sigma,tau) coordinates don't match up properly\n"
"        between this patch/edge and the other patch/edge!\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        my_patch().coords_set_{rho,sigma,tau}={%s,%s,%s}\n"
"        my_edge().coords_set_{perp,par}={%s,%s}\n"
"        other_patch().coords_set_{rho,sigma,tau}={%s,%s,%s}\n"
"        other_edge().coords_set_{perp,par}={%s,%s}\n"
,
	my_patch().name(), my_edge().name(),
	other_patch().name(), other_edge().name(),
	local_coords::name_of_coords_set(my_patch().coords_set_rho()),
	local_coords::name_of_coords_set(my_patch().coords_set_sigma()),
	local_coords::name_of_coords_set(my_patch().coords_set_tau()),
	local_coords::name_of_coords_set(my_edge().coords_set_perp()),
	local_coords::name_of_coords_set(my_edge().coords_set_par()),
	local_coords::name_of_coords_set(other_patch().coords_set_rho()),
	local_coords::name_of_coords_set(other_patch().coords_set_sigma()),
	local_coords::name_of_coords_set(other_patch().coords_set_tau()),
	local_coords::name_of_coords_set(other_edge().coords_set_perp()),
	local_coords::name_of_coords_set(other_edge().coords_set_par()));
								/*NOTREACHED*/

// perp coordinate must match (mod 2*pi) across the two patches
// after taking into account any overlap
const fp other_overlap = N_overlap_points * other_edge().perp_map().delta_fp();
const fp other_outer_perp_minus_overlap	// move back inwards into other patch
					// by overlap distance, to get a value
					// that should match our own
					// grid_outer_perp() value
	= other_edge().grid_outer_perp()
	  + (other_edge().is_min() ? + other_overlap : - other_overlap);
if (! local_coords::fuzzy_EQ_ang(my_edge().grid_outer_perp(),
				 other_outer_perp_minus_overlap))
   then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        perp coordinate doesn't match (mod 2*pi) across the two patches!\n"
"        my_edge().grid_outer_perp()=%g   <--(compare this)\n"
"        N_overlap_points=%d other_overlap=%g\n"
"        other_edge.grid_outer_perp()=%g\n"
"        other_outer_perp_minus_overlap=%g   <--(against this)\n"
,
		   my_patch().name(), my_edge().name(),
		   other_patch().name(), other_edge().name(),
		   double(my_edge().grid_outer_perp()),
		   N_overlap_points, double(other_overlap),
		   double(other_edge().grid_outer_perp()),
		   double(other_outer_perp_minus_overlap));	/*NOTREACHED*/


//
// set up the iperp interpatch coordinate mapping
//

// compute the iperp --> other_iperp mapping for a sample point;
// ... if the ghost zone is empty, then the sample point will necessarily
//     be out-of-range in the ghost zone, so we use the *unchecked*
//     conversions to avoid errors in this case
// ... we do the computation using the fact that  perp  is the same
//     coordinate in both patches (modulo 2*pi radians = 360 degrees)
const int sample_iperp = outer_iperp();
const fp sample_perp = my_edge().perp_map()
				.fp_of_int_unchecked(sample_iperp);
						// unchecked conversion here!
const fp other_sample_perp = other_patch()
			     .modulo_reduce_ang(other_edge().perp_is_rho(),
						sample_perp);
const fp fp_other_sample_iperp = other_edge()
				 .fp_iperp_of_perp(other_sample_perp);

// verify that this is fuzzily a grid point
if (! fuzzy<fp>::is_integer(fp_other_sample_iperp))
   then error_exit(ERROR_EXIT,
"***** interpatch_ghost_zone::interpatch_ghost_zone:\n"
"        my_patch().name()=\"%s\" my_edge().name()=%s\n"
"        other_patch().name()=\"%s\" other_edge().name()=%s\n"
"        sample_iperp=%d sample_perp=%g\n"
"        other_sample_perp=%g fp_other_sample_iperp=%g\n"
"        ==> fp_other_sample_iperp isn't fuzzily an integer!\n"
"        ==> patches aren't commensurate in the perpendicular coordinate!\n"
,
		   my_patch().name(), my_edge().name(),
		   other_patch().name(), other_edge().name(),
		   sample_iperp, double(sample_perp),
		   double(other_sample_perp),
		   double(fp_other_sample_iperp));		/*NOTREACHED*/
const int other_sample_iperp = round<fp>::to_integer(fp_other_sample_iperp);

// compute the +/- sign (direction) of the iperp --> other_iperp mapping
//
// Since perp is the same in both patches (mod 2*pi radians = 360 degrees),
// the overall +/- sign is just the product of the signs of the two individual
// iperp <--> perp mappings.
//
// ... here we encode signs as (floating-point) +/- 1.0
const double iperp_map_sign_pm1
	= jtutil::signum(my_edge().perp_map().delta_fp())
	  * jtutil::signum(other_edge().perp_map().delta_fp());
// ... here as standard bool
const bool iperp_map_is_plus = (iperp_map_sign_pm1 > 0.0);

// now we finally know enough to set up the iperp interpatch coordinate mapping
other_iperp_ = new cpm_map<fp>(min_iperp(), max_iperp(),
			       sample_iperp, other_sample_iperp,
			       iperp_map_is_plus);

//
// set up the par interpatch coordinate mapping
//

// ... note that we can't yet count on any other ghost zone existing,
//     so we don't yet know whether or not we'll want our corners
//     (since that depends on the type of our patch's adjacent ghost zones)
//     ==> we include the corners on the chance we may want them later,
//	   and use the appropriate parts of them in
//	     extend_scalar_gridfn_to_ghost_zone()
//	   below
const int ghost_zone_min_ipar = my_edge().min_ipar_with_corners();
const int ghost_zone_max_ipar = my_edge().max_ipar_with_corners();

other_par_ = new array2d<fp>(min_iperp(), max_iperp(),
			     ghost_zone_min_ipar, ghost_zone_max_ipar);

	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	for (int ipar = ghost_zone_min_ipar ;
	     ipar <= ghost_zone_max_ipar ;
	     ++ipar)
	{
	// compute the  other_par corresponding to  (iperp,ipar)
	// ... here we use the fact (which we verified above) that
	//     other edge's parallel coordinate == our tau coordinate
	//     (at least modulo 2*pi radians = 360 degrees)
	const fp perp  = my_edge().perp_of_iperp(iperp);
	const fp par   = my_edge().par_of_ipar(ipar);

	const fp rho   = my_edge().  rho_of_perp_par(perp, par);
	const fp sigma = my_edge().sigma_of_perp_par(perp, par);

	const fp tau   = my_patch().tau_of_rho_sigma(rho, sigma);
	const fp other_par = other_patch()
			     .modulo_reduce_ang(other_edge().par_is_rho(), tau);

	(*other_par_)(iperp,ipar) = other_par;
	}
	}
}

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

//
// this function destroys an  interpatch_ghost_zone  object.
//
interpatch_ghost_zone::~interpatch_ghost_zone()
{
delete other_par_;
delete other_iperp_;
}

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

//
// This function assert()s that
// * our frontier pointer is non-NULL,
// * the other patch has an interpatch ghost zone on this edge
// * that interpatch ghost zone points back to us
//
void interpatch_ghost_zone::assert_fully_setup() const
{
assert(other_frontier_ != NULL);

const ghost_zone& ogz = other_patch().ghost_zone_on_edge(other_edge());
assert(ogz.is_interpatch());

const interpatch_ghost_zone& oigz
	= static_cast<const interpatch_ghost_zone &>(ogz);
assert(& oigz.other_patch() == & my_patch());
}

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

//
// This function extends a scalar gridfn defined on the other patch's
// nominal grid, to a specified part of this ghost zone, by interpatch
// interpolating it from the neighboring patch's frontier.
//
// N.b. "scalar" here refers solely to the symmetry properties of the
// gridfn under the interpatch coordinate transformations; it need not
// be an actual scalar under more general coordinate transformations.
//
void interpatch_ghost_zone::extend_scalar_gridfn_to_ghost_zone
	(int gfn,
	 bool want_min_par_corner,
	 bool want_non_corner,
	 bool want_max_par_corner)
	const
{
other_frontier_->setup_interpolation_for_gridfn(gfn);

	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	for (int ipar = min_ipar(iperp) ; ipar <= max_ipar(iperp) ; ++ipar)
	{
	// do we want to do this point?
	if (! my_edge().ipar_is_in_selected_corner(want_min_par_corner,
						   want_non_corner,
						   want_max_par_corner,
						   ipar) )
	   then continue;				// *** LOOP CONTROL ***

	int oiperp = other_iperp_->map(iperp);
	fp opar = (*other_par_)(iperp,ipar);
	fp ogridfn = other_frontier_->interpolate(gfn, oiperp, opar);

	int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
	int isigma = my_edge().isigma_of_iperp_ipar(iperp,ipar);
	my_patch().gridfn(gfn, irho,isigma) = ogridfn;
	}
	}
}