more mods to use Cactus local interpolator

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

2 files changed, 44 insertions, 34 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index e21e09c..c9e170b 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -64,20 +64,27 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in,
 		    ? my_edge_in.max_ipar_with_corners()
 		    : my_edge_in.max_ipar_without_corners())
 	  min_parindex_used_(min_parindex_used_in),
-	  max_parindex_used_(max_parindex_used_in),
+	  N_interp_points_(min_iperp_in, max_iperp_in),
 	  interp_par_(interp_par_in),
 	  ghosted_min_gfn_(ghosted_min_gfn_in),
 	  ghosted_max_gfn_(ghosted_max_gfn_in),
 	  interp_handle_(interp_handle_in),
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
-	  interp_coord_origin_(0.0),	// computed later
-	  interp_coord_delta_(my_patch().delta_ang(my_edge().par_is_rho())),
+	  gridfn_coord_origin_(my_edge().par_of_ipar(min_ipar())),
+	  gridfn_coord_delta_(my_edge().par_map().delta_fp()),
 	  gridfn_data_ptrs_(ghosted_min_gfn_in, ghosted_max_gfn_in,
 			    min_iperp_in, max_iperp_in),
 	  gridfn_offsets_(ghosted_min_gfn_in, ghosted_max_gfn_in,
 			  min_iperp_in, max_iperp_in),
 	  gridfn_par_stride_(0)		// computed later
 {
+// number of points to interpolate at each iperp
+	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
+	{
+	N_interp_points_(iperp)
+		= jtutil::how_many_in_range(min_parindex_used_in(iperp),
+					    max_parindex_used_in(iperp));
+	}
 }
 
 //*****************************************************************************
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 39bdb13..ba62792 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -38,10 +38,13 @@
 //
 // In general we interpolate each gridfn at a number of distinct par
 // for each iperp; the integer "parindex" indexes these points.  We
-// attach no particular semantics to parindex, and it need not be 0-origin.
-// However, we do assume that its range is not too large, and is roughly
-// comparable at all values of iperp, so it's ok to treat (iperp,parindex)
-// as a 2-D rectangular index space for the interpolation points.
+// attach no particular semantics to parindex, and it need not be 0-origin
+// or have the same range for each iperp.  [In practice, parindex will
+// be the other patch's ipar coordinate.]  We use
+//    parsub = parindex - min_parindex_used_(iperp)
+// as a 0-origin array subscript indexing the points for each interpolation.
+// We assume that the range of parsub for each iperp is roughly similar,
+// so it's ok to use (iperp,parsub) as a 2-D rectangular index space.
 //
 // We use the Cactus local interpolator CCTK_InterpLocalUniform()
 // to do the interpolation.  To minimize interpolator overheads, we
@@ -57,7 +60,7 @@ public:
 	const patch& my_patch() const { return my_patch_; }
 	const patch_edge& my_edge() const { return my_edge_; }
 
-	// min/max (iperp,ipar) of the frontier
+	// min/max (iperp,ipar) of the gridfn data we'll use
 	int min_iperp() const { return min_iperp_; }
 	int max_iperp() const { return max_iperp_; }
 	int min_ipar() const { return min_ipar_; }
@@ -88,15 +91,14 @@ public:
 	//			   (for sizing (iperp,parindex) arrays etc).
 	// [min,max]_parindex_used_in(iperp)
 	//	= [min,max] range of parindex actually used at each iperp.
-	//	  See the note below about the required lifetime of this
-	//	  array.
 	// interp_par_in(iperp,parindex)
 	//	= Gives the par coordinates at which we will interpolate;
 	//	  array entries outside the range [min,max]_parindex_in
 	//	  are unused (n.b. this interface implicitly takes the
 	//	  par coordinates to be independent of ghosted_gfn).
-	//	  See the note below about the required lifetime of this
-	//	  array.
+	//	  We keep a reference to the passed-in  interp_par_in
+	//	  array, so this should have a lifetime at last as long
+	//	  as that of this object.
 	// ghosted_[min,max]_gfn = The largest ghosted_gfn range for this
 	//			   frontier; any given interpolate() call
 	//			   may specify a subrange.
@@ -107,11 +109,6 @@ public:
 	//	  This table is not modified by any actions of this
 	//	  class.
 	//
-	// This constructor keeps references to the passed-in
-	// [min,max]_parindex_used_in and interp_par_in arrays, so
-	// these should have lifetimes at least as long as that of
-	// this patch_frontier object.
-	//
 	// Note that this constructor looks at what type of ghost zones
 	// are on the adjacent edges to decide how to handle the corners:
 	// * if an adjacent ghost zone is a symmetry ghost zone
@@ -152,17 +149,23 @@ private:
 	// type as described in comments to constructor (above)
 	int min_ipar_, max_ipar_;
 
-	// [min,max]_parindex_used_(iperp)
-	//	= [min,max] range of parindex actually used at each iperp.
+	// min parindex at each iperp
+	// ... this is a reference to an array passed in to our constructor
+	//     ==> we do *not* own this!
+	const jtutil::array1d<int>& min_parindex_used_;
+
+	// number of points to interpolate at each iperp
+	// == 1 + maximum parsub
+	// ... index is (iperp)
+	const jtutil::array1d<int> N_interp_points_;
+
 	// interp_par_(iperp,parindex)
 	//	= Gives the par coordinates at which we will interpolate;
 	//	  array entries outside the range [min,max]_parindex_in
 	//	  are unused (n.b. this interface implicitly takes the
 	//	  par coordinates to be independent of ghosted_gfn).
-	// ... these are references to arrays passed in constructor
-	//     ==> we do *not* own these
-	const jtutil::array1d<fp>& min_parindex_used_;
-	const jtutil::array1d<fp>& max_parindex_used_;
+	// ... this is a reference to an array passed in to our constructor
+	//     ==> we do *not* own this!
 	const jtutil::array2d<fp>& interp_par_;	// indices (iperp,parindex)
 
 	// range of ghosted gfns for which we may interpolate
@@ -174,31 +177,31 @@ private:
 
 	// Cactus handle to our private Cactus key/value table
 	// giving parameters for the interpolation operator
-	// (we *do* own this -- it's a copy of the passed-in table)
+	// ... this starts out as a copy of the passed-in table,
+	//     then gets extra stuff added to it specific to this
+	//     frontier; it's shared across all iperp
+	// ... we own this table
 	int interp_par_table_handle_;
 
-	// (par) coordinate origin and delta values for the interpolator
-	fp interp_coord_origin_, interp_coord_delta_;
+	// (par) origin and delta values of the gridfn data
+	fp gridfn_coord_origin_, gridfn_coord_delta_;
 
 	//
-	// the interpolator accesses the gridfn data with the generic
+	// The interpolator accesses the gridfn data with the generic
 	// subscripting expression
-	//	data[offset + parindex*stride]
+	//	data[offset + parsub*stride]
 	// where we get to choose the pointer  data  and the  offset
 	// value for each (gridfn,iperp), and the stride for each iperp.
 	//
-	// our strategy will be to choose  offset  so that the []
-	// expression is 0 for parindex = min_parindex_used_(iperp) 
+	// Our strategy will be to leave  offset  unspecified so it
+	// defaults to 0, specify a common stride for all (gfn,iperp),
+	// and specify the appropriate data pointer for each (gfn,iperp).
 	//
 
 	// --> start of gridfn data we can use for each iperp
 	// subscripts are (gfn,iperp)
 	jtutil::array2d<fp *> gridfn_data_ptrs_;
 
-	// offsets for gridfn subscripting by interpolator
-	// subscripts are (gfn,iperp)
-	jtutil::array2d<int> gridfn_offsets_;
-
 	// stride of par coordinate in gridfn data
 	// (same for all gfn and iperp)
 	int gridfn_par_stride_;