rework again how we set up the gridfn array subscripting for the

interpolator --> now uses   ipar - min_ipar  as the interpolator i


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

2 files changed, 19 insertions, 57 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index ee12e10..421405e 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -61,7 +61,9 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 	  interp_par_(interp_par_in),
 	  interp_handle_(interp_handle_in),
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
-	  gridfn_coord_origin_(my_edge().par_map().origin()),
+	  // note interpolator coordinate origin = par of min_ipar() point
+	  // cf. comments in "patch_interp.hh" on gridfn array subscripting
+	  gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar())),
 	  gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
 	  gridfn_type_codes_ (min_gfn_, max_gfn_),
 	  gridfn_data_ptrs_  (min_gfn_, max_gfn_),
@@ -81,8 +83,8 @@ if (status < 0)
 		   status);					/*NOTREACHED*/
 
 //
-// set up the offset, stride, and min/max subscripts in the
-// parameter table so the interpolator can access the gridfn data arrays
+// set up the gridfn array stride in the parameter table
+// so the interpolator can access the gridfn data arrays
 //
 
 const int N_dims = 1;
@@ -100,49 +102,6 @@ if (status < 0)
 ,
 		   interp_par_table_handle_, status);		/*NOTREACHED*/
 
-// offsets (same for all gfn)
-jtutil::array1d offsets(min_gfn_, max_gfn_);
-	for (int gfn = min_gfn_ ; gfn <= max_gfn_ ; ++gfn)
-	{
-	offsets(gfn) = - min_ipar_ * stride;
-	}
-const int N_gridfns = jtutil::how_many_in_range(min_gfn_, max_gfn_);
-// n.b. Util_TableSetIntArray() *copies* our array,
-//	so it doesn't matter that the array is a local variable
-//	in this function and will disappear when this function returns
-status = Util_TableSetIntArray(interp_par_table_handle_,
-			       N_gridfns, &offsets(min_gfn_),
-			       "input_array_offsets");
-if (status < 0)
-   then error_exit(ERROR_EXIT,
-"***** patch_interp::patch_interp():\n"
-"        can't set gridfn offsets in interpolator parmameter table!\n"
-"        table handle=%d   error status=%d\n"
-,
-		   interp_par_table_handle_, status);		/*NOTREACHED*/
-
-// min/max subscripts
-status = Util_TableSetIntArray(interp_par_table_handle_,
-			       N_dims, &min_ipar_,
-			       "input_array_min_subscripts");
-if (status < 0)
-   then error_exit(ERROR_EXIT,
-"***** patch_interp::patch_interp():\n"
-"        can't set min gridfn subscript in interpolator parmameter table!\n"
-"        table handle=%d   error status=%d\n"
-,
-		   interp_par_table_handle_, status);		/*NOTREACHED*/
-status = Util_TableSetIntArray(interp_par_table_handle_,
-			       N_dims, &max_ipar_,
-			       "input_array_max_subscripts");
-if (status < 0)
-   then error_exit(ERROR_EXIT,
-"***** patch_interp::patch_interp():\n"
-"        can't set max gridfn subscript in interpolator parmameter table!\n"
-"        table handle=%d   error status=%d\n"
-,
-		   interp_par_table_handle_, status);		/*NOTREACHED*/
-
 //
 // set up the gridfn type codes for the interpolator
 //
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index a213b43..c8cf8a0 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -89,17 +89,20 @@
 // the gridfn data is a bit tricky:  The interpolator accesses the
 // gridfn data using the generic (1-D) subscripting expression
 //	data[offset + i*stride]
-// where i is the data array index.  It's convenient to take i=ipar
-// (this way the mapping from i to par is just given by the usual
-// origin and delta values for the par coordinate).  The stride is
-// fixed by the gridfn storage mapping (we get this from the patch).
-// We then calculate the offset such that  offset + min_ipar*stride == 0 ,
-// and then pass the data pointer for each interpolator call as a pointer
-// to the min_ipar grid point at that iperp.  (We also need to set the
-// min and max i in the interpolator parameter table.)  With this strategy
-// we can share the interpolator parameter table across all the iperp
-// values, and we don't need to modify the parameter table after the
-// initial setup.
+// where i is the data array index.  However, we don't want to use
+// offset, because it has to be supplied in the parameter table as
+// an array subscripted by gfn, and so would require changing the
+// parameter table for each call on interpolate() (with potentially
+// different numbers of gridfns being interpolated).  Instead, at each
+// iperp we use i = ipar-min_ipar, so the default offset=0 makes the
+// subscripting expression zero for ipar = min_ipar.  This also makes
+// the interpolator's min_i = 0 and max_i be dims-1 (both the defaults),
+// so those also don't have to be set in the parameter table.  We set
+// the interpolator's data coordinate origin to the par coordinate for
+// min_ipar, so it correctly maps i --> par.  With this strategy we can
+// share the interpolator parameter table across all the iperp values,
+// and we don't need to modify the parameter table at all after the
+// initial setup in our constructor.
 //
 
 class	patch_interp