partial mods to correspond to recent changes to patch_frontier::

-- still needs more work


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@416 f88db872-0e4f-0410-b76b-b9085cfa78c5

2 files changed, 145 insertions, 101 deletions

diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 6cbdd31..3af6bdc 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -9,12 +9,14 @@
 //
 // interpatch_ghost_zone::interpatch_ghost_zone
 // interpatch_ghost_zone::~interpatch_ghost_zone
-// interpatch_ghost_zone::assert_fully_setup
+// interpatch_ghost_zone::[min,max]_ipar
+// interpatch_ghost_zone::setup_other_patch_frontie
 // interpatch_ghost_zone::synchronize
 //
 
 #include <stdio.h>
 #include <assert.h>
+#include <limits.h>
 #include <math.h>
 
 #include "jt/stdc.h"
@@ -161,9 +163,16 @@ interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in,
 	: 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
+	  other_frontier_(NULL),	// set by setup_other_frontier()
+	  other_iperp_(NULL),		// set in ctor body
+          min_ipar_used_(min_iperp(), max_iperp()),
+          max_ipar_used_(min_iperp(), max_iperp()),
+	  other_par_(min_iperp(), max_iperp(),
+		     min_possible_ipar(), max_possible_ipar()),
+	  interp_result_buffer_
+		(my_patch().ghosted_min_gfn(), my_patch().ghosted_max_gfn(),
+		 min_iperp(), max_iperp(),
+		 min_possible_ipar(), max_possible_ipar())
 {
 //
 // verify that we have the expected relationships between
@@ -247,6 +256,7 @@ if (! local_coords::fuzzy_EQ_ang(my_edge().grid_outer_perp(),
 
 //
 // set up the iperp interpatch coordinate mapping
+// (gives other patch's iperp coordinate for interpolation)
 //
 
 // compute the iperp --> other_iperp mapping for a sample point;
@@ -302,27 +312,31 @@ other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 				       sample_iperp, other_sample_iperp,
 				       iperp_map_is_plus);
 
+
 //
-// set up the par interpatch coordinate mapping
+// set up arrays giving [min,max] ipar that we'll use at each iperp
+//
+	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
+	{
+	min_ipar_used_(iperp) = min_ipar(iperp);
+	max_ipar_used_(iperp) = max_ipar(iperp);
+	}
+
+
+//
+// set up array giving other patch's par coordinate for interpolation
 //
 
 // ... 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
-//	     synchronize()
-//	   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);
+//	   and use the appropriate parts of them in  synchronize()  below
 
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
-	for (int ipar = ghost_zone_min_ipar ;
-	     ipar <= ghost_zone_max_ipar ;
+	for (int ipar = min_possible_ipar() ;
+	     ipar <= max_possible_ipar() ;
 	     ++ipar)
 	{
 	// compute the  other_par corresponding to  (iperp,ipar)
@@ -339,7 +353,7 @@ other_par_ = new array2d<fp>(min_iperp(), max_iperp(),
 	const fp other_par = other_patch()
 			     .modulo_reduce_ang(other_edge().par_is_rho(), tau);
 
-	(*other_par_)(iperp,ipar) = other_par;
+	other_par_(iperp,ipar) = other_par;
 	}
 	}
 }
@@ -351,37 +365,78 @@ other_par_ = new array2d<fp>(min_iperp(), max_iperp(),
 //
 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
+// These functions compute the [min,max] ipar of the ghost zone for
+// a given iperp, taking into account how we treat the corners
+// (cf. the example in the header comments in "ghost_zone.hh"):
 //
-void interpatch_ghost_zone::assert_fully_setup() const
+// 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.  E.g. For the example in the
+//    header comments "ghost_zone.hh", 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 interpatch_ghost_zone::min_ipar(int iperp) const
 {
-assert(other_frontier_ != NULL);
+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);
+}
 
-const ghost_zone& ogz = other_patch().ghost_zone_on_edge(other_edge());
-assert(ogz.is_interpatch());
+int interpatch_ghost_zone::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);
+}
 
-const interpatch_ghost_zone& oigz
-	= static_cast<const interpatch_ghost_zone &>(ogz);
-assert(& oigz.other_patch() == & my_patch());
+//******************************************************************************
+
+//
+// This function constructs and sets up the other patch's patch_frontier::
+// object to handle interpolation from us.
+//
+void interpatch_ghost_zone::setup_other_frontier
+	(int ghosted_min_gfn_in, int ghosted_max_gfn_in,
+	 int interp_handle_in, int interp_par_table_handle_in)
+{
+// set up the interpolator results buffer
+interp_buffer_ = new jtutil::array3d<fp>(ghosted_min_gfn_in, ghosted_max_gfn_in,
+					 min_iperp(), max_iperp(),
+					 min_ipar_used_, max_ipar_used_);
+
+patch_frontier *of
+	= new patch_frontier(other_edge(),
+			     other_min_iperp(), other_max_iperp(),
 }
 
+
 //******************************************************************************
 
 //
 // This function "synchronizes" a ghost zone, i.e. it updates the
 // ghost-zone values of the specified gridfns via the appropriate
-// interpatch interpolations.  The flags specify which part(s) of the
-// ghost zone we want.
+// interpatch interpolations.
+//
+// The flags specify which part(s) of the ghost zone we want, but
+// the present implementation only supports the case where all the
+// flags are  true , i.e. we want the entire ghost zone.
 //
 void interpatch_ghost_zone::synchronize
 	(int ghosted_min_gfn, int ghosted_max_gfn,
@@ -390,25 +445,37 @@ void interpatch_ghost_zone::synchronize
 	 bool want_max_par_corner)
 	const
 {
-// FIXME FIXME
+// make sure the caller wants the entire ghost zone
+if (! (want_min_par_corner && want_non_corner && want_max_par_corner))
+   then error_exit(ERROR_EXIT,
+"***** interpatch_ghost_zone::synchronize(ghosted_[min,max]_gfn=[%d,%d]):\n"
+"        can only handle full ghost zone, not partial,\n"
+"        but not all flags are true!\n"
+"        want_min_par_corner=(int)%d\n"
+"        want_non_corner=(int)%d\n"
+"        want_max_par_corner=(int)%d\n"
+,
+		   ghosted_min_gfn, ghosted_max_gfn,
+		   bool want_min_par_corner,
+		   bool want_non_corner,
+		   bool want_max_par_corner);			/*NOTREACHED*/
+
+// do the interpolation into our result buffer
+other_frontier_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
+			     interp_result_buffer_);
+
+// store the results
+	for (int gfn = ghosted_min_gfn, gfn <= ghosted_max_gfn ; ++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;
+	my_patch().ghosted_gridfn(gfn, irho,isigma)
+		= interp_result_buffer_(gfn, iperp,ipar);
+	}
 	}
 	}
 }
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index f0a270c..342b230 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -282,8 +282,16 @@ public:
 	int outer_iperp() const
 		{ return is_min() ? min_iperp() : max_iperp(); }
 
+	// min/max ipar that might possibly be part of this ghost zone
+	// (derived classes may actually use a subset of this)
+	int min_possible_ipar() const
+		{ return my_edge().min_ipar_with_corners(); }
+	int max_possible_ipar() const
+		{ return my_edge().max_ipar_with_corners(); }
+
 	// assert() that ghost zone is fully setup:
-	// symmetry ghost zone ==> no-op
+	// defined here ==> no-op
+	// symmetry ghost zone ==> unchanged ==> 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
@@ -404,8 +412,8 @@ public:
 		{ 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)
+	// ... we always include 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(); }
 
@@ -444,6 +452,11 @@ private:
 //
 // derived class for interpatch ghost zone of a patch
 //
+// A ghost_zone object maps (my_iperp,my_ipar) coordinates to the other
+// patch's (other_iperp,other_par) coordinates, then calls the other patch's
+// patch_frontier object to interpolate the other patch's data to those
+// coordinates.
+//
 // Note  const  qualifiers refer to the data stored by
 //    synchronize()
 // Since there are no nonconst member functions, once an  interpatch_ghost_zone
@@ -492,35 +505,8 @@ public:
 	// 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);
-		}
+	int min_ipar(int iperp) const;
+	int max_ipar(int iperp) const;
 
 	// 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",
@@ -530,17 +516,10 @@ public:
 	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;
-		}
+	// construct/setup other patch's frontier
+	void setup_other_frontier
+		(int ghosted_min_gfn_in, int ghosted_max_gfn_in,
+		 int interp_handle_in, int interp_par_table_handle_in);
 
 	// constructor, destructor
 	interpatch_ghost_zone(const patch_edge& my_edge_in,
@@ -548,13 +527,6 @@ public:
 			      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
@@ -566,22 +538,27 @@ private:
 	patch& other_patch_;
 	const patch_edge& other_edge_;
 
+	// initialized to NULL in constructor,
+	// set to proper value by  setup_other_frontier()
+	// ... we do *not* own this object (the other patch does)!
+	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
+	// ... we own this object
 	jtutil::cpm_map<fp> *other_iperp_;
 
+	// [min,max]_ipar used at each iperp
+	// = values of [min,max]_ipar() fns
+	// ... index is (my_iperp)
+	jtutil::array1d<int> min_ipar_used_, max_ipar_used_;
+
 	// 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_;
-
-	// initialized to NULL in constructor,
-	// set to proper value by  set_other_frontier()
-	patch_frontier *other_frontier_;
+	jtutil::array2d<fp> other_par_;
 
 	// buffer into which the other patch's patch_frontier object
 	// will store the interpolated gridfn values
 	// ... subscripts are (gfn, my_iperp,my_ipar)
-	jtutil::array3d<fp> 
+	jtutil::array3d<fp> interp_result_buffer_;
 	};