further work on new design to use Cactus local interp

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

1 files changed, 81 insertions, 34 deletions

diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 19d6f54..39bdb13 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -39,9 +39,9 @@
 // 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 roughly comparable at all
-// values of iperp, so there's no major waste in treating (iperp,parindex)
-// as a 2-D index space for the interpolation points.
+// 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.
 //
 // We use the Cactus local interpolator CCTK_InterpLocalUniform()
 // to do the interpolation.  To minimize interpolator overheads, we
@@ -81,36 +81,47 @@ private:
 public:
 	//
 	// Constructor arguments:
-	// my_edge_in = identifies the patch/edge to which this frontier
-	//		is to belong
-	// [min,max]_iperp_in = the range of iperp for this frontier
-	// [min,max]_parindex_in = extreme [min,max] range of parindex
-	//			   (for sizing (iperp,parindex) arrays etc)
+	// my_edge_in = Identifies the patch/edge to which this frontier
+	//		is to belong.
+	// [min,max]_iperp_in = The range of iperp for this frontier.
+	// [min,max]_parindex_in = Extreme [min,max] range of parindex
+	//			   (for sizing (iperp,parindex) arrays etc).
 	// [min,max]_parindex_used_in(iperp)
-	//	= [min,max] range of parindex actually used at each iperp
-	// interp_par(iperp,parindex)
-	//	= gives the par coordinates at which we will interpolate;
+	//	= [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)
-	// ghosted_[min,max]_gfn = the largest ghosted_gfn range for this
+	//	  par coordinates to be independent of ghosted_gfn).
+	//	  See the note below about the required lifetime of this
+	//	  array.
+	// ghosted_[min,max]_gfn = The largest ghosted_gfn range for this
 	//			   frontier; any given interpolate() call
-	//			   may specify a subrange
-	// interp_handle_in = Cactus handle to the interpolation operator
+	//			   may specify a subrange.
+	// interp_handle_in = Cactus handle to the interpolation operator.
 	// interp_par_table_handle_in
 	//	= Cactus handle to a Cactus key/value table giving
-	//	  parameters (eg order) for the interpolation operator
+	//	  parameters (eg order) for the interpolation operator.
+	//	  This table is not modified by any actions of this
+	//	  class.
 	//
-	// ... 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
-	//	 we assume it's already been mirrored by the time
-	//	 we get set up ==> we can (and do) use the data from it
-	//     * if an adjacent ghost zone is an interpatch ghost zone,
-	//	 we can't assume it's already been interpatch-interpolated
-	//	 by the time we're called ==> we don't use data from it
-	//     ==> constructor requires that the adjacent-side  ghost_zone
-	//         objects already exist!
+	// 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
+	//   we assume it's already been mirrored by the time
+	//   we get set up ==> we can (and do) use the data from it
+	// * if an adjacent ghost zone is an interpatch ghost zone,
+	//   we can't assume it's already been interpatch-interpolated
+	//   by the time we're called ==> we don't use data from it
+	// Thus the constructor requires that the adjacent-side  ghost_zone
+	// objects already exist!
 	//
 	patch_frontier(const patch_edge& my_edge_in,
 		       int min_iperp_in, int max_iperp_in,
@@ -131,11 +142,29 @@ private:
 	patch_frontier& operator=(const patch_frontier& rhs);
 
 private:
-	const patch_edge& my_edge_;
 	const patch& my_patch_;
+	const patch_edge& my_edge_;
 
 	int min_iperp_, max_iperp_;
 
+	// range of ipar in this patch from which we may use gridfn
+	// data in interpolation -- depends on adjacent ghost zones'
+	// 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.
+	// 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_;
+	const jtutil::array2d<fp>& interp_par_;	// indices (iperp,parindex)
+
 	// range of ghosted gfns for which we may interpolate
 	// (an actual interpolation may cover a subrange of this)
 	int ghosted_min_gfn_, ghosted_max_gfn_;
@@ -143,16 +172,34 @@ private:
 	// Cactus handle to the interpolation operator
 	int interp_handle_;
 
-	// Cactus handle to a Cactus key/value table
-	// giving parameters (eg order) for the interpolation operator
+	// 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)
 	int interp_par_table_handle_;
 
 	// (par) coordinate origin and delta values for the interpolator
 	fp interp_coord_origin_, interp_coord_delta_;
 
-	// par coordinates at which we are to interpolate
-	// ... interpolation treats this as a 1-D array;
-	//     we set up values in 3-D
-	// ... subscripts are (gfn, iperp,parindex)
-	jtutil::array3d<fp> interp_par_coords_;
+	//
+	// the interpolator accesses the gridfn data with the generic
+	// subscripting expression
+	//	data[offset + parindex*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) 
+	//
+
+	// --> 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_;
 	};