finish construction/setup code

--> should all be ready for test compiles now


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

5 files changed, 217 insertions, 126 deletions

diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 3af6bdc..2a5f6c5 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -10,7 +10,7 @@
 // interpatch_ghost_zone::interpatch_ghost_zone
 // interpatch_ghost_zone::~interpatch_ghost_zone
 // interpatch_ghost_zone::[min,max]_ipar
-// interpatch_ghost_zone::setup_other_patch_frontie
+// interpatch_ghost_zone::setup_other_patch_frontier
 // interpatch_ghost_zone::synchronize
 //
 
@@ -80,12 +80,12 @@ const fp fp_symmetry_period_plane_iperp = symmetry_edge().fp_grid_outer_iperp();
 // i.e. if both edges have iperp as the same min/max "direction",
 //	then the mapping is  iperp increasing --> iperp decreasing
 //      (i.e. the map's sign is -1)
-const bool iperp_map_is_plus
+const bool is_iperp_map_plus
 	= ! (my_edge().is_min() == symmetry_edge().is_min());
 iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 				     fp_my_period_plane_iperp,
 				     fp_symmetry_period_plane_iperp,
-				     iperp_map_is_plus);
+				     is_iperp_map_plus);
 
 //
 // parallel map
@@ -164,15 +164,12 @@ interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in,
 	  other_patch_(other_edge_in.my_patch()),
 	  other_edge_(other_edge_in),
 	  other_frontier_(NULL),	// set by setup_other_frontier()
-	  other_iperp_(NULL),		// set in ctor body
-          min_ipar_used_(min_iperp(), max_iperp()),
-          max_ipar_used_(min_iperp(), max_iperp()),
-	  other_par_(min_iperp(), max_iperp(),
-		     min_possible_ipar(), max_possible_ipar()),
-	  interp_result_buffer_
-		(my_patch().ghosted_min_gfn(), my_patch().ghosted_max_gfn(),
-		 min_iperp(), max_iperp(),
-		 min_possible_ipar(), max_possible_ipar())
+	  // remaining subobjects are all set properly in ctor body
+	  other_iperp_(NULL),
+	  other_min_iperp_(0), other_max_iperp_(0),
+          min_ipar_used_(NULL),	max_ipar_used_(NULL),
+	  other_par_(NULL),
+	  interp_result_buffer_(NULL)
 {
 //
 // verify that we have the expected relationships between
@@ -300,26 +297,35 @@ const int other_sample_iperp
 // the overall +/- sign is just the product of the signs of the two individual
 // iperp <--> perp mappings.
 //
-// ... here we encode signs as (floating-point) +/- 1.0
+// ... signs encoded as (floating-point) +/- 1.0
 const double iperp_map_sign_pm1
-	= jtutil::signum(      my_edge().perp_map().delta_fp() )
+	=   jtutil::signum(    my_edge().perp_map().delta_fp() )
 	  * jtutil::signum( other_edge().perp_map().delta_fp() );
-// ... here as standard bool
-const bool iperp_map_is_plus = (iperp_map_sign_pm1 > 0.0);
+// ... signs encoded as is_plus bool flag
+const bool is_iperp_map_plus = (iperp_map_sign_pm1 > 0.0);
 
-// now we finally know enough to set up the iperp interpatch coordinate mapping
+// now we finally know enough to set up the other_iperp(iperp)
+// coordinate mapping
 other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 				       sample_iperp, other_sample_iperp,
-				       iperp_map_is_plus);
+				       is_iperp_map_plus);
+other_min_iperp_ = jtutil::min(other_iperp(min_iperp()),
+			       other_iperp(max_iperp()));
+other_max_iperp_ = jtutil::max(other_iperp(min_iperp()),
+			       other_iperp(max_iperp()));
 
 
 //
-// set up arrays giving [min,max] ipar that we'll use at each iperp
+// set up arrays giving [min,max] ipar that we'll use at each other_iperp
+// ... we will pass these arrays by reference to the other patch's
+//     patch_frontier object, with ipar being parindex there
 //
+min_ipar_used_ = new array1d<int>(other_min_iperp_, other_max_iperp_);
+max_ipar_used_ = new array1d<int>(other_min_iperp_, other_max_iperp_);
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
-	min_ipar_used_(iperp) = min_ipar(iperp);
-	max_ipar_used_(iperp) = max_ipar(iperp);
+	min_ipar_used_(other_iperp(iperp)) = min_ipar(iperp);
+	max_ipar_used_(other_iperp(iperp)) = max_ipar(iperp);
 	}
 
 
@@ -332,11 +338,13 @@ other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 //     (since that depends on the type of our patch's adjacent ghost zones)
 //     ==> we include the corners on the chance we may want them later,
 //	   and use the appropriate parts of them in  synchronize()  below
+other_par_ = new array2d<fp>(other_min_iperp_, other_max_iperp_,
+			     ghost_zone_min_ipar(), ghost_zone_max_ipar());
 
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
-	for (int ipar = min_possible_ipar() ;
-	     ipar <= max_possible_ipar() ;
+	for (int ipar = ghost_zone_min_ipar() ;
+	     ipar <= ghost_zone_max_ipar() ;
 	     ++ipar)
 	{
 	// compute the  other_par corresponding to  (iperp,ipar)
@@ -353,9 +361,19 @@ other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 	const fp other_par = other_patch()
 			     .modulo_reduce_ang(other_edge().par_is_rho(), tau);
 
-	other_par_(iperp,ipar) = other_par;
+	other_par_(other_iperp(iperp),ipar) = other_par;
 	}
 	}
+
+
+//
+// set up interpolation result buffer
+//
+interp_result_buffer_
+	= new array3d<fp>
+		(my_patch().min_ghosted_gfn(), my_patch().max_ghosted_gfn(),
+		 other_min_iperp_, other_max_iperp_,
+		 ghost_zone_min_ipar(), ghost_zone_max_ipar());
 }
 
 //******************************************************************************
@@ -365,6 +383,10 @@ other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 //
 interpatch_ghost_zone::~interpatch_ghost_zone()
 {
+delete interp_result_buffer_;
+deleete other_par_;
+delete max_ipar_used_;
+delete min_ipar_used_;
 delete other_iperp_;
 }
 
@@ -409,21 +431,22 @@ return max_par_adjacent_ghost_zone().is_symmetry()
 //******************************************************************************
 
 //
-// This function constructs and sets up the other patch's patch_frontier::
-// object to handle interpolation from us.
+// This function constructs the *other* patch's patch_frontier:: object,
+// and interlinks it with this ghost zone and the other patch.
+//
+// We use our ipar as the patch_frontier's parindex.
 //
-void interpatch_ghost_zone::setup_other_frontier
-	(int ghosted_min_gfn_in, int ghosted_max_gfn_in,
-	 int interp_handle_in, int interp_par_table_handle_in)
+void interpatch_ghost_zone::setup_other_frontier(int interp_handle_in,
+						 int interp_par_table_handle_in)
 {
-// set up the interpolator results buffer
-interp_buffer_ = new jtutil::array3d<fp>(ghosted_min_gfn_in, ghosted_max_gfn_in,
-					 min_iperp(), max_iperp(),
-					 min_ipar_used_, max_ipar_used_);
-
-patch_frontier *of
+const other_frontier_
 	= new patch_frontier(other_edge(),
-			     other_min_iperp(), other_max_iperp(),
+			     other_min_iperp_, other_max_iperp_,
+			     ghost_zone_min_ipar(), ghost_zone_max_ipar(),
+			     *min_ipar_used_, *max_ipar_used_,
+			     *other_par_,
+			     interp_handle_in, interp_par_table_handle_in);
+other_patch().frontier_ptr_on_edge(other_edge()) = other_frontier_;
 }
 
 
@@ -464,7 +487,7 @@ if (! (want_min_par_corner && want_non_corner && want_max_par_corner))
 other_frontier_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
 			     interp_result_buffer_);
 
-// store the results
+// store the results back into our gridfns
 	for (int gfn = ghosted_min_gfn, gfn <= ghosted_max_gfn ; ++gfn)
 	{
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 342b230..9ce4689 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -284,9 +284,9 @@ public:
 
 	// min/max ipar that might possibly be part of this ghost zone
 	// (derived classes may actually use a subset of this)
-	int min_possible_ipar() const
+	int ghost_zone_min_ipar() const
 		{ return my_edge().min_ipar_with_corners(); }
-	int max_possible_ipar() const
+	int ghost_zone_max_ipar() const
 		{ return my_edge().max_ipar_with_corners(); }
 
 	// assert() that ghost zone is fully setup:
@@ -508,18 +508,17 @@ public:
 	int min_ipar(int iperp) const;
 	int max_ipar(int iperp) const;
 
-	// min/max iperp of ghost zone in *other* patch's coordinate system
-	// ... if the ghost zone is empty, its min/max won't be "in range",
-	//     so we use the *unchecked* map here
-	int other_min_iperp() const
-		{ return other_iperp_->map_unchecked(min_iperp()); }
-	int other_max_iperp() const
-		{ return other_iperp_->map_unchecked(max_iperp()); }
+	// convert our iperp --> other patch's iperp
+	int other_iperp(int iperp) const
+		{
+		assert(other_iperp_ != NULL);
+		return other_iperp_->map(iperp);
+		}
 
-	// construct/setup other patch's frontier
-	void setup_other_frontier
-		(int ghosted_min_gfn_in, int ghosted_max_gfn_in,
-		 int interp_handle_in, int interp_par_table_handle_in);
+	// construct *other* patch's frontier,
+	// interlink it with this ghost zone and the other patch
+	void setup_other_frontier(int interp_handle_in,
+				  int interp_par_table_handle_in);
 
 	// constructor, destructor
 	interpatch_ghost_zone(const patch_edge& my_edge_in,
@@ -541,24 +540,41 @@ private:
 	// initialized to NULL in constructor,
 	// set to proper value by  setup_other_frontier()
 	// ... we do *not* own this object (the other patch does)!
-	patch_frontier *other_frontier_;
+	const patch_frontier *other_frontier_;
+
+	//
+	// our remaining subobjects are all pointed-to because
+	// we won't know the range of other_iperp (which we need
+	// to initialize the subobjects) until partway into the
+	// body of our constructor
+	//
 
 	// other patch's iperp coordinates of our ghost zone points
 	// ... maps my_iperp --> other_iperp
-	// ... we own this object
 	jtutil::cpm_map<fp> *other_iperp_;
 
-	// [min,max]_ipar used at each iperp
-	// = values of [min,max]_ipar() fns
-	// ... index is (my_iperp)
-	jtutil::array1d<int> min_ipar_used_, max_ipar_used_;
+	// other patch's [min,max]_iperp of our ghost zone points
+	int other_min_iperp_, other_max_iperp_;
+
+	// [min,max]_ipar used at each other_iperp
+	// ... we will pass these arrays by reference
+	//     to the other patch's patch_frontier object
+	// ... index is (other_iperp)
+	jtutil::array1d<int>* min_ipar_used_;
+	jtutil::array1d<int>* max_ipar_used_;
 
 	// other patch's (fp) parallel coordinates of our ghost zone points
-	// ... subscripts are (my_iperp, my_ipar)
-	jtutil::array2d<fp> other_par_;
+	// ... we will pass this array by reference
+	//     to the other patch's patch_frontier object
+	//     using my_ipar as the patch_frontier's parindex
+	// ... subscripts are (other_iperp, my_ipar)
+	jtutil::array2d<fp>* other_par_;
 
 	// buffer into which the other patch's patch_frontier object
 	// will store the interpolated gridfn values
-	// ... subscripts are (gfn, my_iperp,my_ipar)
-	jtutil::array3d<fp> interp_result_buffer_;
+	// ... we will pass this array by reference
+	//     to the other patch's patch_frontier object
+	//     using my_ipar as the patch_frontier's parindex
+	// ... subscripts are (gfn, other_iperp,my_ipar)
+	jtutil::array3d<fp>* interp_result_buffer_;
 	};
diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index 6f7b524..d575f70 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -8,12 +8,18 @@
 // x_patch::x_patch
 // y_patch::y_patch
 //
+// patch::minmax_ang_ghost_zone
+// patch::ghost_zone_on_edge
+// patch::ghost_zone_ptr
+// patch::patch_frontier_ptr
 // patch::setup_mirror_symmetry_ghost_zone
 // patch::setup_periodic_symmetry_ghost_zone
 // patch::setup_interpatch_ghost_zone
 // patch::set_ghost_zone
-// patch::edge_adjacent_to_patch
+// patch::minmax_ang_ghost_zone
 // patch::ghost_zone_on_edge
+// patch::frontier_ptr_on_edge
+// patch::edge_adjacent_to_patch
 // patch::interpatch_ghost_zone_on_edge
 // patch::assert_all_ghost_zones_fully_setup
 //
@@ -165,6 +171,80 @@ y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
 //******************************************************************************
 
 //
+// This function returns a reference to the specified ghost zone
+// of this patch.
+//
+ghost_zone& patch::minmax_ang_ghost_zone(bool want_min, bool want_rho)
+	const
+{
+return want_min ? (want_rho ? min_rho_ghost_zone()
+			    : min_sigma_ghost_zone())
+		: (want_rho ? max_rho_ghost_zone()
+			    : max_sigma_ghost_zone());
+}
+
+//******************************************************************************
+
+//
+// This function returns a reference to the specified ghost zone
+// of this patch.
+//
+ghost_zone& patch::ghost_zone_on_edge(const patch_edge &edge)
+const
+{
+assert(& edge.my_patch() == this);
+return minmax_ang_ghost_zone(edge.is_min(), edge.is_rho());
+}
+
+//******************************************************************************
+
+//
+// This function returns a reference to the specified ghost zone
+// pointer data member of this patch, asserting that the current
+// value of this pointer is NULL.  This is used (only) by
+//	patch::setup_*_ghost_zone()
+// to set the ghost zone pointers; the assert() ensures that it can't
+// be used to overwrite a non-NULL pointer.
+//
+ghost_zone*& patch::ghost_zone_ptr(const patch_edge& edge)
+{
+assert(& edge.my_patch() == this);
+
+ghost_zone*& gzp = edge.is_min() ? (edge.is_rho() ? min_rho_ghost_zone_
+						  : min_sigma_ghost_zone_)
+				 : (edge.is_rho() ? max_rho_ghost_zone_
+						  : max_sigma_ghost_zone_);
+assert(gzp == NULL);
+return gzp;
+}
+
+//******************************************************************************
+
+//
+// This function returns a reference to the specified patch frontier
+// pointer data member of this patch, asserting that the current
+// value of this pointer is NULL.  This is used (only) by the other
+// patch's
+//	interpatch_ghost_zone::setup_other_frontier()
+// to set this patch's patch frontier pointers; the assert() ensures that
+// it can't be used to overwrite a non-NULL pointer.
+//
+patch_frontier*& patch::patch_frontier_ptr(const patch_edge& edge)
+{
+assert(& edge.my_patch() == this);
+
+patch_frontier*& pfp
+	= edge.is_min() ? (edge.is_rho() ? min_rho_patch_frontier_
+					 : min_sigma_patch_frontier_)
+			: (edge.is_rho() ? max_rho_patch_frontier_
+					 : max_sigma_patch_frontier_);
+assert(pfp == NULL);
+return pfp;
+}
+
+//******************************************************************************
+
+//
 // This function assert()s that a specified ghost zone of this patch
 // hasn't already been set up, then sets it up as a mirror-symmetry
 // ghost zone.  It returns a reference to the newly-set-up symmetry
@@ -191,22 +271,21 @@ return *temp;
 //
 symmetry_ghost_zone& patch::setup_periodic_symmetry_ghost_zone
 	(const patch_edge& my_edge, const patch_edge& symmetry_edge,
-	 bool ipar_map_is_plus)
+	 bool is_ipar_map_plus)
 {
 // make sure we belong to the right patch
 assert(& my_edge.my_patch() == this);
 
 int my_sample_ipar = my_edge.min_ipar_without_corners();
 int symmetry_sample_ipar
-	= ipar_map_is_plus ? symmetry_edge.min_ipar_without_corners()
+	= is_ipar_map_plus ? symmetry_edge.min_ipar_without_corners()
 			   : symmetry_edge.max_ipar_without_corners();
 
 symmetry_ghost_zone *temp
-    = new symmetry_ghost_zone(my_edge,          symmetry_edge,
-			      my_sample_ipar,	symmetry_sample_ipar,
-			      ipar_map_is_plus);
-set_ghost_zone(temp, my_edge);
-
+	= new symmetry_ghost_zone(my_edge,        symmetry_edge,
+				  my_sample_ipar, symmetry_sample_ipar,
+				  is_ipar_map_plus);
+ghost_zone_ptr(my_edge) = temp;
 return *temp;
 }
 
@@ -228,31 +307,13 @@ assert(& my_edge.my_patch() == this);
 interpatch_ghost_zone *temp
 	= new interpatch_ghost_zone(my_edge, other_edge,
 				    N_overlap_points);
-set_ghost_zone(temp, my_edge);
+ghost_zone_ptr(my_edge) = temp;
 return *temp;
 }
 
 //******************************************************************************
 
 //
-// This is a helper function for  setup_*_ghost_zone().  This function
-// assert()s that one of the ghost zone pointers is NULL, then stores a
-// value in it.
-//
-void patch::set_ghost_zone(ghost_zone* gzp, const patch_edge& my_edge)
-{
-ghost_zone*& my_gzp
-	= my_edge.is_min()
-	  ? (my_edge.is_rho() ? min_rho_ghost_zone_ : min_sigma_ghost_zone_)
-	  : (my_edge.is_rho() ? max_rho_ghost_zone_ : max_sigma_ghost_zone_);
-
-assert(my_gzp == NULL);
-my_gzp = gzp;
-}
-
-//******************************************************************************
-
-//
 // This function finds which patch edge is adjacent to a neighboring
 // patch q, or does an error_exit() if q isn't actually a neighboring patch.
 // The computation is done using only (rho,sigma) coordinate sets and
@@ -359,18 +420,6 @@ return p.minmax_ang_patch_edge(common_is_p_min_q_max, common_is_p_rho);
 //******************************************************************************
 
 //
-// This function returns a reference to the ghost zone on a specified
-// edge of this patch.
-//
-ghost_zone& patch::ghost_zone_on_edge(const patch_edge &e) const
-{
-assert(& e.my_patch() == this);
-return minmax_ang_ghost_zone(e.is_min(), e.is_rho());
-}
-
-//******************************************************************************
-
-//
 // This function verifies that the ghost zone on a specified edge
 // is indeed interpatch, and returns a reference to it as an
 //  interpatch_ghost_zone .
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 536305b..f5d906c 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -281,15 +281,24 @@ public:
 		{ return *min_sigma_ghost_zone_; }
 	ghost_zone& max_sigma_ghost_zone() const
 		{ return *max_sigma_ghost_zone_; }
-	ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho)
-		const
-		{
-		return want_min ? (want_rho ? min_rho_ghost_zone()
-					    : min_sigma_ghost_zone())
-				: (want_rho ? max_rho_ghost_zone()
-					    : max_sigma_ghost_zone());
-		}
-	ghost_zone& ghost_zone_on_edge(const patch_edge &e) const;
+	ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho) const;
+	ghost_zone& ghost_zone_on_edge(const patch_edge &edge) const;
+
+
+	//
+	// ***** access to ghost zone/frontier pointers *****
+	//
+private:
+	// these are used by members/friends only for setup
+	// they return references to our ghost zone/frontier pointers
+	// ... used only from patch::setup_*_ghost_zone()
+	const ghost_zone*&     ghost_zone_ptr    (const patch_edge& edge) const;
+	// ... used only from (friend)
+	//	interpatch_ghost_zone::setup_other_frontier()
+	const patch_frontier*& patch_frontier_ptr(const patch_edge& edge) const;
+	friend interpatch_ghost_zone::setup_other_frontier
+		(int interp_handle_in, int interp_par_table_handle_in);
+
 private:
 	// helper function for setup_*_ghost_zone():
 	// assert() that ghost zone pointer on specified edge is NULL
@@ -304,7 +313,6 @@ public:
 		(const patch_edge &e)
 		const;
 
-
 	//
 	// ***** set up ghost zone and frontier subobjects
 	//
@@ -320,8 +328,8 @@ public:
 	// then set it up as periodic-symmetry
 	// ... return reference to newly-set-up ghost zone object
 	symmetry_ghost_zone& setup_periodic_symmetry_ghost_zone
-	   (const patch_edge& my_edge, const patch_edge& symmetry_edge,
-	    bool ipar_map_is_plus);
+		(const patch_edge& my_edge, const patch_edge& symmetry_edge,
+		 bool ipar_map_is_plus);
 
 	// assert() that this ghost zone hasn't been set up yet,
 	// then set it up as interpatch
@@ -332,13 +340,6 @@ public:
 		(const patch_edge& my_edge, const patch_edge& other_edge,
 		 int N_overlap_points);
 
-	// assert() that this ghost zone is interpatch,
-	// then set up the other patch's frontier
-	// ... return reference to new frontier
-	patch_frontier& setup_other_patch_frontier
-		(const patch_edge& my_edge_in,
-		 int interp_handle_in, int interp_par_table_handle_in);
-
 	// assert() that all ghost zones (and frontiers, where applicable)
 	// are fully setup
 	void assert_all_ghost_zones_fully_setup() const;
@@ -405,11 +406,12 @@ private:
 	ghost_zone* max_sigma_ghost_zone_;
 
 	// frontiers (NULL pointers if no frontier for this ghost zone)
-	// ... pointers are initialized by setup_interpatch_ghost_zone()
-	patch_frontier* min_rho_patch_frontier_,
-	patch_frontier* max_rho_patch_frontier_,
-	patch_frontier* min_sigma_patch_frontier_,
-	patch_frontier* max_sigma_patch_frontier_;
+	// ... pointers are by other patch's
+	//      interpatch_ghost_zone::setup_other_frontier()
+	const patch_frontier* min_rho_patch_frontier_,
+	const patch_frontier* max_rho_patch_frontier_,
+	const patch_frontier* min_sigma_patch_frontier_,
+	const patch_frontier* max_sigma_patch_frontier_;
 	};
 
 //*****************************************************************************
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 61dffe3..c0c6dcb 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -192,18 +192,19 @@ private:
 
 	// gridfn type codes for interpolator
 	// ... must be CCTK_INT so we can pass by reference to interpolator
+	// ... must be mutable so we can set values in ctor
 	// ... index is (gfn)
-	jtutil::array1d<CCTK_INT> gridfn_type_codes_;
+	mutable jtutil::array1d<CCTK_INT> gridfn_type_codes_;
 
 	// --> start of gridfn data to use for interpolation
 	//     (reset for each iperp)
 	// ... we do *not* own the pointed-to data!
 	// ... index is (gfn)
-	jtutil::array1d<fp *> gridfn_data_ptrs_;
+	mutable jtutil::array1d<fp *> gridfn_data_ptrs_;
 
 	// --> start of result buffer for interpolation
 	//     (reset for each iperp)
 	// ... we do *not* own the pointed-to data!
 	// ... index is (gfn)
-	jtutil::array1d<fp *> result_buffer_ptrs_;
+	mutable jtutil::array1d<fp *> result_buffer_ptrs_;
 	}