// ghost_zone.hh -- fill in gridfn data in patch ghost zones
// $Id$
//
// ***** design notes for ghost zones *****
// ghost_zone - abstract base class to describe ghost zone of patch
// symmetry_ghost_zone - ... derived class for spacetime-symmetry ghost zone
// interpatch_ghost_zone - ... derived class for interpatch ghost zone
//

//
// prerequisites:
//	<stdio.h>
//	<assert.h>
//	<math.h>
//	"jt/util.hh"
//	"jt/array.hh"
//	"jt/cpm_map.hh"
//	"jt/linear_map.hh"
//	fp.hh
//	coords.hh
//	grid.hh
//	fd_grid.hh
//	patch.hh
//	patch_edge.hh
//

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

//
// ***** design notes for ghost zones *****
//

//
// A  ghost_zone  object describes a patch's ghost zone, and knows how
// to compute gridfns there (we usually speak of "extending" a gridfn to
// a ghost zone or zones) based on either the patch system's symmetry
// or interpolation from a neighboring patch.  ghost_zone is an abstract
// base class, from which we derive two classes:
// * A  symmetry_ghost_zone  object describes a ghost zone which
//   is a (discrete) symmetry of spacetime, either mirror-image or
//   periodic.  Such an object knows how to fill in ghost-zone gridfn
//   data from the "other side" of the symmetry.
// * An  interpatch_ghost_zone  object describes a ghost zone which
//   overlaps another patch.  Such an object knows how to get ghost
//   zone gridfn data from the other patch.  More accurately, it gets
//   the data by asking (calling) the appropriate one of the other
//   patch's  patch_frontier  objects.
// Every patch has (points to) 4  ghost_zone  objects, one for each of
// the patch's sides.  See the comments in "patch.hh" for a "big picture"
// discussion of patches, patch edges, ghost zones, and patch frontiers.
//

//
// There are some unobvious complications involved in extending gridfns
// to the ghost zone "corners", i.e. in ghost zone points that are outside
// the nominal grid in *both* coordinates.  There are 3 basic cases here:
// * A corner between two symmetry ghost zones, for example the -x/-y
//   corner in the example below.  In this case it takes *two* sequential
//   symmetry operations to get gridfn data in the corner from the
//   nominal grid.  Symmetry operations commute, so at each point we'll
//   always get the same results independently of in which order we do
//   the symmetry operations.
// * A corner between two interpatch ghost zones, for example the +x/+y
//   corner in the example below.  In this case we could get the gridfn
//   data by either of two distinct interpolation operations (presumably
//   from two distinct patches), which would in general give slightly
//   different results.  In some ideal world we might take the average
//   of these, but at present we split the corner down its diagonal.
//   For the points on the diagonal we make an arbitrary choice; at
//   present this is that they belong to (and get their data via) the
//   rho ghost zone.
//	[For elliptic equations it's important to ensure that
//	each grid point get a value from exactly one computation,
//	and in particular we can't allow two different ghost
//	zones to "own" (assign values) to the same point, since
//	then the first value would be "lost" and wouldn't show
//	up in the Jacobian matrix.]
// * A corner between a symmetry and an interpatch ghost zone, for
//   example the +x/-y or -x/+y corners in the example below.  In this
//   case we need to first do a symmetry operation in the neighboring
//   patch so we don't have to off-center the interpolation in the
//   non-corner part of the interpatch ghost zone.  Then after the
//   interpatch interpolation, we need to do a final symmetry operation
//   in this patch to set up gridfn data in the corner.
//
// To handle all these cases, we use a 3-phase algorithm to extend
// gridfns to the ghost zones:
// Phase 1: Fill in gridfn data at all the non-corner points of symmetry
//	    ghost zones, by using the symmetries to get this data from
//	    its "home patch" nominal grids.
// Phase 2: Fill in gridfn data in all the interpatch ghost zones, by
//	    interpatch interpolating from neighboring patches as described
//	    above.
// Phase 3: Fill in gridfn data at all the corner points of symmetry
//	    ghost zones, by using the symmetries to get this data from
//	    its "home patch" nominal grids or ghost zones.
// Here a given ghost zone corner may be either a full rectangle (so any
// given point is a member of both adjacent corners), or split down its
// diagonal (so any given point is a member of only one corner).  This
// 3-phase algorithm is actually implemented by
//    patch_system::extend_scalar_gridfn_to_all_ghost_zones()
// which in turn calls
//    symmetry_ghost_zone::extend_scalar_gridfn_to_ghost_zone()
//    interpatch_ghost_zone::extend_scalar_gridfn_to_ghost_zone()
//

//
// For example, consider the +z patch in an octant patch system, with
// the ghost zones being 2 points wide.  The following illustration is
// looking down the z axis, and uses (x,y) for the patch coordinates
// for simplicity:
//
//                    #                                                   //
//                   i+y    i+y    i+y    i+y    i+y    i+y    i+y      //
//   (-2,7) (-1,7)  (0,7)  (1,7)  (2,7)  (3,7)  (4,7)  (5,7)  (6,7)  (7,7)
//    <s-x>  <s-x>    #                                              /i+x
//                    #                                            //
//                   i+y    i+y    i+y    i+y    i+y    i+y      //
//   (-2,6) (-1,6)  (0,6)  (1,6)  (2,6)  (3,6)  (4,6)  (5,6)  (6,6)  (7,6)
//    <s-x>  <s-x>    #                                       /i+x    i+x
//                    #                                     //
//                    #                                   //
//   (-2,5) (-1,5)   2,5)--(1,5)--(2,5)--(3,5)--(4,5)--(5,5)  (6,5)  (7,5)
//     s-x    s-x     #                                  |     i+x    i+x
//                    #                                  |
//                    #                                  |
//   (-2,4) (-1,4)  (0,4)  (1,4)  (2,4)  (3,4)  (4,4)  (5,4)  (6,4)  (7,4)
//     s-x    s-x     #                                  |     i+x    i+x
//                    #                                  |
//                    #                                  |
//   (-2,3) (-1,3)  (0,3)  (1,3)  (2,3)  (3,3)  (4,3)  (5,3)  (6,3)  (7,3)
//     s-x    s-x     #                                  |     i+x    i+x
//                    #                                  |
//                    #                                  |
//   (-2,2) (-1,2)  (0,2)  (1,2)  (2,2)  (3,2)  (4,2)  (5,2)  (6,2)  (7,2)
//     s-x    s-x     #                                  |     i+x    i+x
//                    #                                  |
//                    #                                  |
//   (-2,1) (-1,1)  (0,1)  (1,1)  (2,1)  (3,1)  (4,1)  (5,1)  (6,1)  (7,1)
//     s-x    s-x     #                                  |     i+x    i+x
//                    #                                  |
//                    #                                  |
//  #(-2,0)#(-1,0)##(0,0)##(1,0)##(2,0)##(3,0)##(4,0)##(5,0)##(6,0)##(7,0)
//     s-x    s-x     #                                        i+x    i+x
//                    #
//    <s-y>  <s-y>   s-y    s-y    s-y    s-y    s-y    s-y   <s-y>  <s-y>
//   (-2,-1)(-1,-1) (0,-1) (1,-1) (2,-1) (3,-1) (4,-1) (5,-1) (6,-1) (7,-1)
//    <s-x>  <s-x>    #
//                    #
//    <s-y>  <s-y>   s-y    s-y    s-y    s-y    s-y    s-y   <s-y>  <s-y>
//   (-2,-2)(-1,-2) (0,-2) (1,-2) (2,-2) (3,-2) (4,-2) (5,-2) (6,-2) (7,-2)
//    <s-x>  <s-x>    #
//                    #
//
// For this example,
// * The xz plane and yz plane are marked with ### lines
// * The +z patch's nominal grid is ([0,5],[0,5]), i.e. 0 <= x,y <= 5;
//   its boundary lines are shown with single lines --- and | .
// * The diagonal where we've split corners are marked with // lines.
// * The +z patch's ghost zones are
//	-x: (-1,[-1,7]), (-2,[-2,7]) 
//	+x: (6,[-2,6]), (7,[-2,7])
//	-y: ([-2, 7],[-2,-1])
//	+y: ([-2,5],6), ([-2,6],7)
// * The +z patch's frontiers are
//	+x: ([ 3,4],[-2,7])
//	+y: ([-2,7],[ 3,4])
//   Note that in both cases the frontier includes the points computed
//   by symmetry (in phase 1 of our 3-phase algorithm) on the adjacent
//   edges! There are no -x or -y frontiers, since no interpolation is
//   needed across those boundaries of this patch.
//
// Our 3-phase algorithm described above thus becomes:
// Phase 1: Fill in gridfn values at points marked with "s-x" below or
//	    "s-y" above via symmetry mirroring across the -x boundary
//	    (yz plane) or -y boundary (xz plane), as described by the
//	    +z patch's -x or -y  symmetry_ghost_zone  object respectively.
// Phase 2: Fill in gridfn values at points marked with "i+x" below or
//	    "i+y" above via interpatch interpolation from the neighboring
//	    patch across the +z patch's +x or +y boundary, as described
//	    by the +z patch's +x or +y  interpatch_ghost_zone  object
//	    respectively.
// Phase 3: Fill in gridfn values at points marked with "<s-x>" below or
//	    "<s-y>" above via symmetry mirroring across the -x boundary
//	    (yz plane) or -y boundary (xz plane), as described by the
//	    +z patch's -x or -y  symmetry_ghost_zone  object respectively.
//
// The diagonal ** line shows the boundary between the +x and +y ghost
// zones; points there may be interpolated via either of the two possible
// interpatch_ghost_zone  objects.
//	[At present we require all points in a given ghost zone
//	to be interpolated from the same neighboring patch and
//	 patch_frontier  object, so must arbitrarily choose one
//	of the two neighbors for the diagonal points.  In theory
//	it would be better to take the average of the two neighbors,
//	but in practice this doesn't matter for horizon finding
//	or other elliptic stuff (it would matter for stability in
//	time evolutions using patch-angular finite differencing,
//	eg in my "mpe" multipatch black-hole-excision code).]
//

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

//
// ghost_zone - abstract base class to describe ghost zone of patch
//
// This is an abstract base class describing a generic patch ghost zone.
// This might represent either of:
// - a discrete symmetry of spacetime (derived class symmetry_ghost_zone)
// - an overlap with another patch (derived class interpatch_ghost_zone)
//

//
// N.b. const qualifiers on member functions of ghost_zone and its derived
// classes refer to the underlying gridfn data, since this is much more
// useful than applying the qualifiers only to the ghost zone (& derived)
// objects themselves.
// 

// forward declarations
class symmetry_ghost_zone;
class interpatch_ghost_zone;
class patch_system;

class	ghost_zone
	{
public:
	//
	// main client interface: set up gridfn data in a ghost zone
	// (flags specify which part(s) of the ghost zone we want)
	//
	virtual void extend_scalar_gridfn_to_ghost_zone
		(int gfn,
		 bool want_min_par_corner,
		 bool want_noncorner,
		 bool want_max_par_corner)
		const
		= 0;

public:
	// to which patch/edge do we belong?
	patch& my_patch() const { return my_patch_; }
	const patch_edge& my_edge() const { return my_edge_; }

	// what type of ghost zone are we?
	bool is_interpatch() const { return  is_interpatch_; }
	bool is_symmetry()   const { return !is_interpatch_; }

	// convenience forwarding functions down to patch_edge::
	bool is_min() const { return my_edge().is_min(); }
	bool is_rho() const { return my_edge().is_rho(); }

	// adjacent ghost zones to our min/max corners
	const ghost_zone& min_par_adjacent_ghost_zone() const
		{
		return my_patch()
		       .ghost_zone_on_edge( my_edge().min_par_adjacent_edge() );
		}
	const ghost_zone& max_par_adjacent_ghost_zone() const
		{
		return my_patch()
		       .ghost_zone_on_edge( my_edge().max_par_adjacent_edge() );
		}

	// min/max iperp of the ghost zone
	int min_iperp() const
		{
		return my_patch()
		       .minmax_ang_ghost_zone__min_iperp(is_min(), is_rho());
		}
	int max_iperp() const
		{
		return my_patch()
		       .minmax_ang_ghost_zone__max_iperp(is_min(), is_rho());
		}

	// inner/outer iperp of the ghost zone wrt our patch
	int inner_iperp() const
		{ return is_min() ? max_iperp() : min_iperp(); }
	int outer_iperp() const
		{ return is_min() ? min_iperp() : max_iperp(); }

	// assert() that ghost zone is fully setup:
	// symmetry ghost zone ==> no-op
	// interpatch ghost zone ==> assert() that frontier pointer is non-NULL,
	//			     assert() that other patch has interpatch
	//				ghost zone on this edge, and that it
	//				points back to us
	virtual void assert_fully_setup() const { }

protected:
	// ... values for  is_interpatch_in  constructor argument
	static const bool ghost_zone_is_symmetry = false;
	static const bool ghost_zone_is_interpatch = true;

	// constructor
	// ... only used in implementing our derived classes;
	//     the rest of the world constructs our derived classes instead
	ghost_zone(const patch_edge& my_edge_in,
		   bool is_interpatch_in)
		: my_patch_(my_edge_in.my_patch()),
		  my_edge_(my_edge_in),
		  is_interpatch_(is_interpatch_in)
		{ }
public:
	// destructor must be virtual to allow destruction
	// of derived classes via ptr/ref to this class
	virtual ~ghost_zone() { }

private:
	// we forbid copying and passing by value
	// by declaring the copy constructor and assignment operator
	// private, but never defining them (either here or in derived classes)
	ghost_zone(const ghost_zone& rhs);
	ghost_zone& operator=(const ghost_zone& rhs);

private:
	patch& my_patch_;
	const patch_edge& my_edge_;
	const bool is_interpatch_;
	};

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

//
// symmetry_ghost_zone - derived class for spacetime-symmetry ghost zone
//
// In practice, there are two types of spacetime symmetry ghost zone:
// mirror symmetry and periodic symmetry.  However, it turns out that the
// code needed to handle periodic BCs is basically a superset of that
// needed to handle mirror symmetries, so this class represents a generic
// symmetry ghost zone which may be of either type, and once constructed
// doesn't distinguish between the two.
//
// In general, a symmetry ghost zone implies that there's a 1-1 mapping
// between ghost zone points of this patch, and (a subset of the) interior
// points of this or another patch.  If tensors are involved (this isn't
// used at present in the horizon finder), there's also a corresponding
// 1-1 mapping between (angular) tensor components.
//
// A mirror-symmetry ghost zone is specified by (the constructor arguments)
// - a patch edge
// - the (fp) perp coordinate of the mirror plane
// The mapping of ghost zone points is thus "just" the mirror imaging of
// iperp across the symmetry plane within this same patch.  (The mapping
// leaves ipar invariant.)
//
// A periodic-symmetry ghost zone is specified by (the constructor arguments)
// - a patch edge (specifies the ghost zone)
// - the patch edge to which the ghost zone is to be mapped
// - a pair of ipar coordinates, one on this edge and one on the other edge,
//   which map into each other
// - the sign of the ipar mapping (does increasing ipar on this edge map to
//   increasing or decreasing ipar on the other edge?)
// The mapping of ghost zone points is the periodic mapping; this may map
// the ghost zone points to interior points of either this same patch or a
// different one.
//
// In general, the symmetry mapping of ghost zone points is of the form
//	(iperp, ipar) --> (const +/- iperp, const +/- ipar)
// The iperp mapping is always in the direction
//	outside the patch --> inside the patch
// while the ipar mapping might have either sign.
// If there are tensors, the corresponding mapping of tensor components is
//	(index_perp, index_par) --> (+/-) (+/-) (index_perp, index_par)
// (that is, the two +/- signs are multiplied).

//
// Note  const  qualifiers refer to the results of
//	iperp_map_coord()
//	ipar_map_coord()
// Since all the member functions are  const , a  symmetry_ghost_zone
// object is effectively always  const .
//
class	symmetry_ghost_zone
	: public ghost_zone
	{
public:
	// high-level client interface:
	// set up gridfn data in a ghost zone by symmetry mapping
	// (flags specify which part(s) of the ghost zone we want)
	void extend_scalar_gridfn_to_ghost_zone
		(int gfn,
		 bool want_min_par_corner,
		 bool want_noncorner,
		 bool want_max_par_corner)
		const;

	// low-level client interface: symmetry-map coordinates
	const patch& symmetry_patch() const { return symmetry_patch_; }
	const patch_edge& symmetry_edge() const { return symmetry_edge_; }
	int iperp_map_of_iperp(int iperp) const
		{ return iperp_map_->map(iperp); }
	int ipar_map_of_ipar(int ipar) const
		{ return ipar_map_->map(ipar); }
	fp fp_sign_of_iperp_map() const
		{ return iperp_map_->fp_sign(); }
	fp sign_of_ipar_map() const
		{ return ipar_map_->fp_sign(); }

	// min/max/size ipar of the ghost zone
	// (i.e. always including the corners
	// (cf. the example at the start of this file)
	int min_ipar() const { return my_edge().min_ipar_with_corners(); }
	int max_ipar() const { return my_edge().max_ipar_with_corners(); }

	// constructor for mirror-symmetry ghost zone
	symmetry_ghost_zone(const patch_edge& my_edge_in);

	// constructor for periodic-symmetry ghost zone
	// ... ipar mapping specified by giving sample point and mapping sign
	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);

	~symmetry_ghost_zone();

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

private:
	// symmetry mapping --> interior of which patch?  which edge?
	const patch& symmetry_patch_;
	const patch_edge& symmetry_edge_;

	// symmetry mappings for (iperp,ipar)
	// ... we own these objects
	const jtutil::cpm_map<fp> *iperp_map_;
	const jtutil::cpm_map<fp> *ipar_map_;
	};

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

//
// derived class for interpatch ghost zone of a patch
//
// Note  const  qualifiers refer to the data stored by
//    extend_scalar_gridfn_to_ghost_zone()
// Since there are no nonconst member functions, once an  interpatch_ghost_zone
// object is contructed, it's effectively always taken as  const .
//
// Note that as described in the "design notes for ghost zones"
// comments above,  interpatch_ghost_zone  objects are constructed in
// the 2nd and 3rd phase of the overall 3-phase construction process
// described at the comments at the start of this file.
// - first set up the object itslf and its links to/from the patches
//   and their edges
// - then set up and link to the other patch's patch_frontier object
//   from which this object will interpolate
//

class patch_frontier;

class	interpatch_ghost_zone
	: public ghost_zone
	{
public:
	// high-level client interface:
	// set up gridfn data in a ghost zone by symmetry mapping
	// (flags specify which part(s) of the ghost zone we want)
	void extend_scalar_gridfn_to_ghost_zone
		(int gfn,
		 bool want_min_par_corner,
		 bool want_noncorner,
		 bool want_max_par_corner)
		const;

	// basic connectivity info
	const patch& other_patch() const
		{ return other_patch_; }
	const patch_edge& other_edge() const
		{ return other_edge_; }
	const patch_frontier& other_frontier() const
		{
		assert(other_frontier_ != NULL);
		return *other_frontier_;
		}

	// min/max/size ipar of the ghost zone for specified iperp
	// taking into account how we treat the corners
	// (cf. the example at the start of this file)
	// ... that is, if an adjacent ghost zone is mirror-symmetry,
	//		   we do not include that corner;
	//	        if an adjacent ghost zone is interpatch,
	//		   we include up to the diagonal line,
	//		   and if we are a rho ghost zone,
	//		       then also the diagonal line itself
	//		   i.e. in the example at the start of this file,
	//		   the +x ghost zone includes (6,6), (7,6), and (7,7),
	//		   while the +y ghost zone includes (6,7)
	// ... in the following 2 functions,
	//     the  iabs()  term includes the diagonal,
	//     so we must remove the diagonal for !is_rho,
	//     i.e. add 1 to min_ipar and subtract 1 from max_ipar
	int min_ipar(int iperp) const
		{
		return min_par_adjacent_ghost_zone().is_symmetry()
		       ? my_edge().min_ipar_with_corners()
		       : my_edge().min_ipar_without_corners()
			 - iabs(iperp - my_edge().nominal_grid_outer_iperp())
			 + (is_rho() ? 0 : 1);
		}
	int max_ipar(int iperp) const
		{
		return max_par_adjacent_ghost_zone().is_symmetry()
		       ? my_edge().max_ipar_with_corners()
		       : my_edge().max_ipar_without_corners()
			 + iabs(iperp - my_edge().nominal_grid_outer_iperp())
			 - (is_rho() ? 0 : 1);
		}

	// min/max iperp of ghost zone in *other* patch's coordinate system
	// ... if the ghost zone is empty, its min/max won't be "in range",
	//     so we use the *unchecked* map here
	int other_min_iperp() const
		{ return other_iperp_->map_unchecked(min_iperp()); }
	int other_max_iperp() const
		{ return other_iperp_->map_unchecked(max_iperp()); }

	// set "other frontier" pointer -- should be used only by
	//  patch_system::setup_adjacent_patch_frontiers()
	// ... we'd like to declare this private, and make that function
	//     a friend, but at this point we haven't seen the declaration
	//     of class patch_system::, so C++ doesn't grok declaring one
	//     of its member functions as a friend. :(
	void set_other_frontier(patch_frontier &other_frontier_in)
		{
		assert(other_frontier_ == NULL);
		other_frontier_ = & other_frontier_in;
		}

	// constructor, destructor
	interpatch_ghost_zone(const patch_edge& my_edge_in,
			      const patch_edge& other_edge_in,
			      int N_overlap_points);
	~interpatch_ghost_zone();

public:
	// assert() that
	// * 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 assert_fully_setup() const;

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

private:
	patch& other_patch_;
	const patch_edge& other_edge_;

	// initialized to NULL in constructor,
	// set to proper value by  set_other_frontier()
	patch_frontier *other_frontier_;

	// other patch's iperp coordinates of our ghost zone points
	// ... allocated and initialized in constructor (const thereafter)
	// ... maps my_iperp --> other_iperp
	jtutil::cpm_map<fp> *other_iperp_;

	// other patch's (fp) parallel coordinates of our ghost zone points
	// ... allocated and initialized in constructor (const thereafter)
	// ... subscripts are (my_iperp, my_ipar)
	jtutil::array2d<fp> *other_par_;
	};