fix subscripting of gridfn arrays by the interpolator

--> now set offset as well as stride and min/max subscript
... also add comments clarifying this


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

2 files changed, 59 insertions, 27 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index be26ac6..ee12e10 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -51,6 +51,8 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 			   int interp_handle_in, int interp_par_table_handle_in)
 	: my_patch_(my_edge_in.my_patch()),
 	  my_edge_(my_edge_in),
+	  min_gfn_(my_patch().ghosted_min_gfn()),
+	  max_gfn_(my_patch().ghosted_max_gfn()),
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
 	  min_ipar_(min_ipar_for_gridfn_data()),
 	  max_ipar_(max_ipar_for_gridfn_data()),
@@ -61,12 +63,9 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
 	  gridfn_coord_origin_(my_edge().par_map().origin()),
 	  gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
-	  gridfn_type_codes_ (my_patch().ghosted_min_gfn(),
-			      my_patch().ghosted_max_gfn()),
-	  gridfn_data_ptrs_  (my_patch().ghosted_min_gfn(),
-			      my_patch().ghosted_max_gfn()),
-	  result_buffer_ptrs_(my_patch().ghosted_min_gfn(),
-			      my_patch().ghosted_max_gfn())
+	  gridfn_type_codes_ (min_gfn_, max_gfn_),
+	  gridfn_data_ptrs_  (min_gfn_, max_gfn_),
+	  result_buffer_ptrs_(min_gfn_, max_gfn_))
 {
 int status;
 
@@ -81,8 +80,14 @@ if (status < 0)
 		   interp_par_table_handle_in,
 		   status);					/*NOTREACHED*/
 
-// set up gridfn storage addressing in interpolator parameter table
+//
+// set up the offset, stride, and min/max subscripts in the
+// parameter table so the interpolator can access the gridfn data arrays
+//
+
 const int N_dims = 1;
+
+// stride
 const CCTK_INT stride = my_edge().ghosted_par_stride();
 status = Util_TableSetIntArray(interp_par_table_handle_,
 			       N_dims, &stride,
@@ -95,7 +100,28 @@ if (status < 0)
 ,
 		   interp_par_table_handle_, status);		/*NOTREACHED*/
 
-// set up patch interpolation region in interpolator parameter table
+// 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");
@@ -117,10 +143,10 @@ if (status < 0)
 ,
 		   interp_par_table_handle_, status);		/*NOTREACHED*/
 
-// set up gridfn type codes
-	for (int gfn = my_patch().ghosted_min_gfn() ;
-	     gfn <=    my_patch().ghosted_max_gfn() ;
-	     ++gfn)
+//
+// set up the gridfn type codes for the interpolator
+//
+	for (int gfn = min_gfn_ ; gfn <= max_gfn_ ; ++gfn)
 	{
 	gridfn_type_codes_(gfn) = CCTK_VARIABLE_REAL;	// fp == CCTK_REAL
 	}
@@ -179,21 +205,6 @@ void patch_interp::interpolate(int ghosted_min_gfn_to_interp,
 				 jtutil::array3d<fp>& result_buffer_array)
 	const
 {
-//
-// We do a separate interpolation for each iperp.
-//
-// The interpolator accesses the gridfn data with the generic
-// subscripting expression
-//	data[offset + i*stride]
-// where i is a 0-origin point index for each interpolation, and where
-// we get to choose the pointer  data  and the  offset  value for each
-// (gridfn,iperp), and the  stride  for each iperp.  Our strategy is to
-// leave  offset  unspecified so it defaults to 0, use a common stride
-// for all iperp (already set up in our constructor), and specify the
-// appropriate data pointers (--> the start of the gridfn data we should
-// use) for each (gfn,iperp).  With this strategy we don't need to modify
-// the parameter table at all after the constructor.
-//
 const int N_dims = 1;
 const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp,
 						ghosted_max_gfn_to_interp);
@@ -223,6 +234,7 @@ const CCTK_INT *gridfn_type_codes_ptr
 		     gfn <= ghosted_max_gfn_to_interp ;
 		     ++gfn)
 		{
+		// set up data pointer to --> (iperp,min_ipar) gridfn
 		const int start_irho
 			= my_edge().  irho_of_iperp_ipar(iperp, min_ipar());
 		const int start_isigma
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 0708b9f..a213b43 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -85,6 +85,22 @@
 // call.  [Different iperp values involve different sets of (1-D)
 // gridfn data, and so inherently require distinct interpolator calls.]
 //
+// Setting up the array subscripting for the interpolator to access
+// 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.
+//
 
 class	patch_interp
 	{
@@ -183,6 +199,10 @@ private:
 	const patch& my_patch_;
 	const patch_edge& my_edge_;
 
+	// range of gfn we can handle
+	// (any given interpolate() call may specify a subrange)
+	const int min_gfn_, max_gfn_;
+
 	// patch interpolation region,
 	// i.e. range of (iperp,ipar) in this patch from which
 	// we will use gridfn data in interpolation