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

2 files changed, 115 insertions, 96 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index 3335897..e21e09c 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -50,18 +50,33 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in,
 			       int ghosted_min_gfn_in, int ghosted_max_gfn_in,
 			       int interp_handle_in,
 			       int interp_par_table_handle_in)
-	: my_edge_(my_edge_in),
-	  my_patch_(my_edge_in.my_patch()),
+	: my_patch_(my_edge_in.my_patch()),
+	  my_edge_(my_edge_in),
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
+	  min_ipar_(my_patch()
+		    .ghost_zone_on_edge(my_edge_in.min_par_adjacent_edge())
+		    .is_symmetry()
+		    ? my_edge_in.min_ipar_with_corners()
+		    : my_edge_in.min_ipar_without_corners()),
+	  max_ipar_(my_patch()
+		    .ghost_zone_on_edge(my_edge_in.max_par_adjacent_edge())
+		    .is_symmetry()
+		    ? 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),
+	  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_(interp_par_table_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())),
-	  interp_par_coords_(ghosted_min_gfn_, ghosted_max_gfn_,
-			     min_iperp_, max_iperp_,
-			     min_parindex_in, max_parindex_in)
+	  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
 {
 }
 
@@ -72,67 +87,24 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in,
 //
 patch_frontier::~patch_frontier()
 {
-// delete the interpolator objects
-	for (int gfn = my_patch().min_gfn() ;
-	     gfn <= my_patch().max_gfn() ;
-	     ++gfn)
-	{
-	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
-	{
-	delete (*interpolator_)(gfn,iperp);
-	}
-	}
-
-delete interpolator_;
-delete gridfn_buffer_;
-delete ripar_map_;
+Util_TableDestroy(interp_par_table_handle_);
 }
 
 //*****************************************************************************
 
 //
-// This function sets up a specified gridfn's interpolation state, i.e.
-// it copies the gridfn data to the contiguous-in-ipar buffer and sets
-// up any interpolator state
+// This function interpolates the specified range of ghosted gridfns
+// at all the coordinates specified when we were constructed.
+// It stores the interpolated data in  result_buffer(gfn, iperp,parindex) .
+// It calls  error_exit() if the point is out of range or any other
+// interpolator error is detected.
 //
-void patch_frontier::setup_interpolation_for_gridfn(int gfn)
-{
-	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
-	{
-	//
-	// copy this row of gridfn data from the (probably non-contiguous)
-	// gridfn array, into our (contiguous-in-ipar) scratch array
-	//
-		for (int ipar = min_ipar() ; ipar <= max_ipar() ; ++ipar)
-		{
-		int irho   = my_edge().  irho_of_iperp_ipar(iperp, ipar);
-		int isigma = my_edge().isigma_of_iperp_ipar(iperp, ipar);
-		(*gridfn_buffer_)(gfn,iperp,ipar)
-			= my_patch().gridfn(gfn,irho,isigma);
-		}
-
-	//
-	// Now we want to set up the interpolator for the just-copied row
-	// of data.  The interpolator interface requires that we pass a
-	// pointer to the [0] data, where the [] subscripting refers to
-	// the interpolator's integer coordinate.  Since ipar=0 may not
-	// even exist, we introduce  ripar := ipar - min_ipar() and use
-	// this as the interpolator's integer coordinate, so that [0]
-	// refers to the start of this ipar row in the gridfn buffer.
-	//
-	const fp* row_ptr = & (*gridfn_buffer_)(gfn,iperp,min_ipar());
-	interpolator_t_* I = (*interpolator_)(gfn,iperp);
-	I->setup_data_y(row_ptr);
-	}
-}
-
-//******************************************************************************
-
-//
-// This function does the actual interpolation of data for a specified
-// gridfn to a specified (iperp,par) point.  It calls  error_exit()
-// if the point is out of range or any other interpolator error is found.
-//
-fp patch_frontier::interpolate(int gfn, int iperp, fp par) const
+fp patch_frontier::interpolate(int ghosted_min_gfn, ghosted_max_gfn,
+			       jtutil::array3d<fp>& result_buffer) const
 {
 const int N_dims = 1;
+
+const int status
+	= CCTK_InterpLocalUniform(N_dims,
+				  interp_handle_,
+				  interp_par_table_handle_,
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_;
 	};