now properly creates ghost zones & frontiers with new Cactus interp fns

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

2 files changed, 252 insertions, 507 deletions

diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 9fbdfc0..969adaf 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -1,19 +1,14 @@
 // patch_system.cc -- describes the (an) entire system of patches
 // $Id$
 //
-// patch_system::patch_size::grid_array_pars
-// patch_system::patch_size::grid_pars
-//
 // patch_system::patch_system
 // patch_system::~patch_system
+// patch_system::construct_patches
 // patch_system::setup_full_sphere_patch_system
 // patch_system::setup_plus_z_hemisphere_patch_system
 // patch_system::setup_plus_xy_quadrant_patch_system
 // patch_system::setup_plus_xz_quadrant_patch_system
 // patch_system::setup_octant_patch_system
-//
-// patch_system::create_patches
-//
 // patch_system::setup_adjacent_ghost_zones
 // patch_system::setup_adjacent_patch_frontiers
 // patch_system::assert_all_ghost_zones_fully_setup
@@ -35,7 +30,6 @@
 #include "jt/cpm_map.hh"
 #include "jt/linear_map.hh"
 #include "jt/interpolate.hh"
-
 using jtutil::error_exit;
 
 #include "fp.hh"
@@ -54,122 +48,57 @@ using jtutil::error_exit;
 //******************************************************************************
 
 //
-// This function computes, and returns a reference to, a
-// struct grid_arrays::grid_array_pars  from the info in a  struct patch_info .
-//
-// The result refers to an internal static buffer in this function; the
-// usual caveats about lifetimes/overwriting apply.
-//
-// Arguments:
-// N_extend_points = The number of grid points to extend each side of the
-//		     patch beyond its nominal boundary (this is used to
-//		     implement overlapping patches).
-//
-const grid_arrays::grid_array_pars&
-  patch_system::patch_size::grid_array_pars
-	(int N_ghost_points, int N_extend_points, fp delta_drho_dsigma)
-	const
-{
-static struct grid_arrays::grid_array_pars array_pars;
-
-array_pars.min_irho   = round<fp>::to_integer(min_drho  /delta_drho_dsigma);
-array_pars.min_isigma = round<fp>::to_integer(min_dsigma/delta_drho_dsigma);
-array_pars.max_irho
-	= array_pars.min_irho
-	  + round<fp>::to_integer((max_drho  -min_drho  ) / delta_drho_dsigma);
-array_pars.max_isigma
-	= array_pars.min_isigma
-	  + round<fp>::to_integer((max_dsigma-min_dsigma) / delta_drho_dsigma);
-array_pars.min_irho   -= N_extend_points;
-array_pars.min_isigma -= N_extend_points;
-array_pars.max_irho   += N_extend_points;
-array_pars.max_isigma += N_extend_points;
-
-array_pars.min_rho_N_ghost_points = N_ghost_points;
-array_pars.max_rho_N_ghost_points = N_ghost_points;
-array_pars.min_sigma_N_ghost_points = N_ghost_points;
-array_pars.max_sigma_N_ghost_points = N_ghost_points;
-
-return array_pars;
-}
-
-//******************************************************************************
-
-//
-// This function computes, and returns a reference to, a
-// struct grid::grid_pars  from the info in a  struct patch_info .
-//
-// The result refers to an internal static buffer in this function; the
-// usual caveats about lifetimes/overwriting apply.
-//
-const grid::grid_pars&
-  patch_system::patch_size::grid_pars
-	(int N_extend_points, fp delta_drho_dsigma)
-	const
-{
-static struct grid::grid_pars grid_pars_buffer;
-const fp extend_drho_dsigma = fp(N_extend_points) * delta_drho_dsigma;
-
-grid_pars_buffer.  min_drho   = min_drho   - extend_drho_dsigma;
-grid_pars_buffer.delta_drho   = delta_drho_dsigma;
-grid_pars_buffer.  max_drho   = max_drho   + extend_drho_dsigma;
-grid_pars_buffer.  min_dsigma = min_dsigma - extend_drho_dsigma;
-grid_pars_buffer.delta_dsigma = delta_drho_dsigma;
-grid_pars_buffer.  max_dsigma = max_dsigma + extend_drho_dsigma;
-
-return grid_pars_buffer;
-}
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
 // This function constructs a  patch_system  object.
 //
+// Constructor arguments:
+// N_ghost_points = Width in grid points of all ghost zones.
+// N_extend_points = Number of grid points to extend each patch past
+//		     "just touching" so as to overlap neighboring patches.
+//		     Thus patches overlap by
+//			N_overlap_points = 2*N_extend_points + 1
+//		     grid points.  For example, with N_extend_points == 2,
+//		     here are the grid points of two neighboring patches:
+//			x   x   x   x   x   X   X
+//                                      |
+//			        O   O   o   o   o   o   o
+//		     Here | marks the "just touching" boundary,
+//		     x and o the grid points before this extension,
+//		     and X and O the extra grid points added by this
+//		     extension.
+// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
+//
 patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
 			   enum patch_system_type type_in,
-			   int N_ghost_points, int N_overlap_points,
+			   int N_ghost_points, int N_extend_points,
 			   fp delta_drho_dsigma,
 			   int min_gfn_in, int max_gfn_in,
-			   int interpolator_order)
+			   int ghosted_min_gfn_in, int ghosted_max_gfn_in,
+			   int interp_handle_in, int interp_par_table_handle_in)
 
 	: global_coords_(origin_x_in, origin_y_in, origin_z_in),
 	  type_(type_in),
-	  N_patches_(N_patches_of_type(type_in))
+	  N_patches_(N_patches_of_type(type_in)),
+	  all_patches_(N_patches_)
 {
-//
-// At present the implementation restricts N_overlap_points to be even.
-//
-// This is "just" to make it easier to set up the patch sizes in the
-//  serup_*_patch_system()  functions below:  This way we can (and do)
-// extend each side of each patch by N_overlap_points/2 grid points.
-//
-// In contrast, if N_overlap_points were odd, then we'd have to extend
-// different sides by different numbers of points, eg by the floor vs
-// the ceiling of N_overlap_points/2.0, which would require some extra
-// bookkeeping to make sure that for each pair of adjacent patches, one
-// got the floor and the other the ceiling on their adjacent sides.
-//
-if (! jtutil::is_even(N_overlap_points))
-   then error_exit(ERROR_EXIT,
-"***** patch_system::patch_system(): implementation restriction\n"
-"                                    N_overlap_points=%d must be even!\n"
-,
-		   N_overlap_points);				/*NOTREACHED*/
-const int N_extend_points = N_overlap_points/2;
-
-all_patches_ = new array1d<patch *>(0, N_patches()-1);
+const int N_overlap_points = 2*N_extend_points + 1;
 
+// construct/interlink the patches, ghost zones, and frontiers
 switch	(type_in)
 	{
 case full_sphere_patch_system:
-	setup_full_sphere_patch_system(N_ghost_points, N_extend_points,
+	construct_patches(patch_system_info::full_sphere::patch_info_array,
+			  N_ghost_points, N_extend_points,
+			  delta_drho_dsigma);
+	interlink_full_sphere_patch_system(N_ghost_points, N_extend_points,
 				       N_overlap_points, delta_drho_dsigma,
 				       min_gfn_in, max_gfn_in,
 				       interpolator_order);
 	break;
 case plus_z_hemisphere_patch_system:
+	construct_patches(patch_system_info::plus_z_hemisphere
+					   ::patch_info_array,
+			  N_ghost_points, N_extend_points,
+			  delta_drho_dsigma);
 	setup_plus_z_hemisphere_patch_system
 		(N_ghost_points, N_extend_points,
 		 N_overlap_points, delta_drho_dsigma,
@@ -177,19 +106,28 @@ case plus_z_hemisphere_patch_system:
 		 interpolator_order);
 	break;
 case plus_xy_quadrant_patch_system:
+	construct_patches(patch_system_info::plus_xy_quadrant::patch_info_array,
+			  N_ghost_points, N_extend_points,
+			  delta_drho_dsigma);
 	setup_plus_xy_quadrant_patch_system(N_ghost_points, N_extend_points,
 					    N_overlap_points, delta_drho_dsigma,
 					    min_gfn_in, max_gfn_in,
 					    interpolator_order);
 	break;
 case plus_xz_quadrant_patch_system:
+	construct_patches(patch_system_info::plus_xz_quadrant::patch_info_array,
+			  N_ghost_points, N_extend_points,
+			  delta_drho_dsigma);
 	setup_plus_xz_quadrant_patch_system(N_ghost_points, N_extend_points,
 					    N_overlap_points, delta_drho_dsigma,
 					    min_gfn_in, max_gfn_in,
 					    interpolator_order);
 	break;
-case octant_patch_system:
-	setup_octant_patch_system(N_ghost_points, N_extend_points,
+case plus_xyz_octant_patch_system:
+	construct_patches(patch_system_info::plus_xyz_octant::patch_info_array,
+			  N_ghost_points, N_extend_points,
+			  delta_drho_dsigma);
+	setup_plus_xyz_octant_patch_system(N_ghost_points, N_extend_points,
 				  N_overlap_points, delta_drho_dsigma,
 				  min_gfn_in, max_gfn_in,
 				  interpolator_order);
@@ -220,69 +158,71 @@ delete all_patches_;
 //******************************************************************************
 
 //
-// This function constructs and interlinks the patches for a
-// full-sphere patch system.
+// This function constructs a set of patches as described by an array
+// of patch_info structures and associated arguments, and makes these
+// patches members of this patch system.  This function does *NOT*
+// create any of the ghost zones or frontiers, and does *NOT* set up
+// any gridfns.
+//
+void patch_system::construct_patches(const struct patch_info patch_info[],
+				     int N_ghost_points, int N_extend_points,
+				     fp delta_drho_dsigma)
+{
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	const struct patch_info& pi = patch_info[pn];
+
+	struct patch *p;
+	switch	(pi.ctype)
+		{
+	case 'z':
+		p = new z_patch(*this, pn,
+				pi.name, pi.is_plus,
+				pi.grid_array_pars(N_ghost_points,
+						   N_extend_points,
+						   delta_drho_dsigma),
+				pi.grid_pars(N_extend_points,
+					     delta_drho_dsigma));
+	case 'x':
+		p = new x_patch(*this, pn,
+				pi.name, pi.is_plus,
+				pi.grid_array_pars(N_ghost_points,
+						   N_extend_points,
+						   delta_drho_dsigma),
+				pi.grid_pars(N_extend_points,
+					     delta_drho_dsigma));
+	case 'y':
+		p = new y_patch(*this, pn,
+				pi.name, pi.is_plus,
+				pi.grid_array_pars(N_ghost_points,
+						   N_extend_points,
+						   delta_drho_dsigma),
+				pi.grid_pars(N_extend_points,
+					     delta_drho_dsigma));
+		break;
+	default:
+		error_exit(ERROR_EXIT,
+"***** patch_system::construct_patches():\n"
+"        unknown patch_info[pn=%d].ctype=0x%02d='%c'!\n"
+,
+			   pn, pi.ctype, pi.ctype);		/*NOTREACHED*/
+		}
+
+	all_patches_[pn] = p;
+	}
+}
+
+//******************************************************************************
+
+//
+// This function interlinks the patches for a full-sphere patch system.
 //
-void patch_system::setup_full_sphere_patch_system
+void patch_system::interlink_full_sphere_patch_system
 	(int N_ghost_points, int N_extend_points,
 	 int N_overlap_points, fp delta_drho_dsigma,
 	 int min_gfn_in, int max_gfn_in,
-	 int interpolator_order)
+	 int interp_handle, int interp_par_table_handle)
 {
-using namespace patch_system_info::full_sphere_patch_system_info;
-
-// create the patches
-patch& pz = *new z_patch(*this, plus_z.patch_number,
-			 "+z", patch::patch_is_plus,
-			 plus_z.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& px = *new x_patch(*this, plus_x.patch_number,
-			 "+x", patch::patch_is_plus,
-			 plus_x.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& py = *new y_patch(*this, plus_y.patch_number,
-			 "+y", patch::patch_is_plus,
-			 plus_y.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& mx = *new x_patch(*this, minus_x.patch_number,
-			 "-x", patch::patch_is_minus,
-			 minus_x.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& my = *new y_patch(*this, minus_y.patch_number,
-			 "-y", patch::patch_is_minus,
-			 minus_y.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& mz = *new z_patch(*this, minus_z.patch_number,
-			 "-z", patch::patch_is_minus,
-			 minus_z.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-
-// attach the patches to the patch system
-(*all_patches_)(plus_z.patch_number) = & pz;
-(*all_patches_)(plus_x.patch_number) = & px;
-(*all_patches_)(plus_y.patch_number) = & py;
-(*all_patches_)(minus_x.patch_number) = & mx;
-(*all_patches_)(minus_y.patch_number) = & my;
-(*all_patches_)(minus_z.patch_number) = & mz;
-
 // create the ghost zones
 setup_adjacent_ghost_zones(pz, px, N_overlap_points);
 setup_adjacent_ghost_zones(pz, py, N_overlap_points);
@@ -298,18 +238,42 @@ setup_adjacent_ghost_zones(mz, mx, N_overlap_points);
 setup_adjacent_ghost_zones(mz, my, N_overlap_points);
 
 // create the patch frontiers
-setup_adjacent_patch_frontiers(pz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, mx, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, my, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(px, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(py, mx, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mx, my, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(my, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mz, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mz, mx, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mz, my, N_overlap_points, interpolator_order);
+setup_adjacent_patch_frontiers(pz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, mx,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(px, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(py, mx,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mx, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(my, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mz, mx,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mz, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
 
 assert_all_ghost_zones_fully_setup();
 }
@@ -324,53 +288,8 @@ void patch_system::setup_plus_z_hemisphere_patch_system
 	(int N_ghost_points, int N_extend_points,
 	 int N_overlap_points, fp delta_drho_dsigma,
 	 int min_gfn_in, int max_gfn_in,
-	 int interpolator_order)
+	 int interp_handle, int interp_par_table_handle)
 {
-using namespace patch_system_info::plus_z_hemisphere_patch_system_info;
-
-// create the patches
-patch& pz = *new z_patch(*this, plus_z.patch_number,
-			 "+z", patch::patch_is_plus,
-			 plus_z.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& px = *new x_patch(*this, plus_x.patch_number,
-			 "+x", patch::patch_is_plus,
-			 plus_x.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& py = *new y_patch(*this, plus_y.patch_number,
-			 "+y", patch::patch_is_plus,
-			 plus_y.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& mx = *new x_patch(*this, minus_x.patch_number,
-			 "-x", patch::patch_is_minus,
-			 minus_x.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& my = *new y_patch(*this, minus_y.patch_number,
-			 "-y", patch::patch_is_minus,
-			 minus_y.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-
-// attach the patches to the patch system
-(*all_patches_)(plus_z.patch_number) = & pz;
-(*all_patches_)(plus_x.patch_number) = & px;
-(*all_patches_)(plus_y.patch_number) = & py;
-(*all_patches_)(minus_x.patch_number) = & mx;
-(*all_patches_)(minus_y.patch_number) = & my;
 
 // create the ghost zones
 setup_adjacent_ghost_zones(pz, px, N_overlap_points);
@@ -387,14 +306,30 @@ mx.setup_mirror_symmetry_ghost_zone(mx.min_rho_patch_edge());
 my.setup_mirror_symmetry_ghost_zone(my.min_rho_patch_edge());
 
 // create the patch frontiers
-setup_adjacent_patch_frontiers(pz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, mx, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, my, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(px, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(py, mx, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mx, my, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(my, px, N_overlap_points, interpolator_order);
+setup_adjacent_patch_frontiers(pz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, mx,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(px, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(py, mx,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mx, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(my, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
 
 assert_all_ghost_zones_fully_setup();
 }
@@ -409,45 +344,8 @@ void patch_system::setup_plus_xy_quadrant_patch_system
 	(int N_ghost_points, int N_extend_points,
 	 int N_overlap_points, fp delta_drho_dsigma,
 	 int min_gfn_in, int max_gfn_in,
-	 int interpolator_order)
+	 int interp_handle, int interp_par_table_handle)
 {
-using namespace patch_system_info::plus_xy_quadrant_patch_system_info;
-
-// create the patches
-patch& pz = *new z_patch(*this, plus_z.patch_number,
-			 "+z", patch::patch_is_plus,
-			 plus_z.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& px = *new x_patch(*this, plus_x.patch_number,
-			 "+x", patch::patch_is_plus,
-			 plus_x.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& py = *new y_patch(*this, plus_y.patch_number,
-			 "+y", patch::patch_is_plus,
-			 plus_y.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& mz = *new z_patch(*this, minus_z.patch_number,
-			 "-z", patch::patch_is_minus,
-			 minus_z.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-
-// attach the patches to the patch system
-(*all_patches_)(plus_z.patch_number) = & pz;
-(*all_patches_)(plus_x.patch_number) = & px;
-(*all_patches_)(plus_y.patch_number) = & py;
-(*all_patches_)(minus_z.patch_number) = & mz;
 
 // create the ghost zones
 setup_adjacent_ghost_zones(pz, px, N_overlap_points);
@@ -466,11 +364,21 @@ setup_periodic_symmetry_ghost_zones(mz.max_rho_patch_edge(),
 				    true);
 
 // create the patch frontiers
-setup_adjacent_patch_frontiers(pz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(px, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(mz, py, N_overlap_points, interpolator_order);
+setup_adjacent_patch_frontiers(pz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(px, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(mz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
 
 assert_all_ghost_zones_fully_setup();
 }
@@ -485,46 +393,8 @@ void patch_system::setup_plus_xz_quadrant_patch_system
 	(int N_ghost_points, int N_extend_points,
 	 int N_overlap_points, fp delta_drho_dsigma,
 	 int min_gfn_in, int max_gfn_in,
-	 int interpolator_order)
+	 int interp_handle, int interp_par_table_handle)
 {
-using namespace patch_system_info::plus_xz_quadrant_patch_system_info;
-
-// create the patches
-patch& pz = *new z_patch(*this, plus_z.patch_number,
-			 "+z", patch::patch_is_plus,
-			 plus_z.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& px = *new x_patch(*this, plus_x.patch_number,
-			 "+x", patch::patch_is_plus,
-			 plus_x.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& py = *new y_patch(*this, plus_y.patch_number,
-			 "+y", patch::patch_is_plus,
-			 plus_y.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& my = *new y_patch(*this, minus_y.patch_number,
-			 "-y", patch::patch_is_minus,
-			 minus_y.grid_array_pars(N_ghost_points,
-						 N_extend_points,
-						 delta_drho_dsigma),
-			 minus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-
-// attach the patches to the patch system
-(*all_patches_)(plus_z.patch_number) = & pz;
-(*all_patches_)(plus_x.patch_number) = & px;
-(*all_patches_)(plus_y.patch_number) = & py;
-(*all_patches_)(minus_y.patch_number) = & my;
-
 // create the ghost zones
 setup_adjacent_ghost_zones(pz, px, N_overlap_points);
 setup_adjacent_ghost_zones(pz, py, N_overlap_points);
@@ -542,11 +412,21 @@ setup_periodic_symmetry_ghost_zones(py.max_sigma_patch_edge(),
 				    false);
 
 // create the patch frontiers
-setup_adjacent_patch_frontiers(pz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, my, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(px, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(px, my, N_overlap_points, interpolator_order);
+setup_adjacent_patch_frontiers(pz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(px, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(px, my,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
 
 assert_all_ghost_zones_fully_setup();
 }
@@ -555,44 +435,14 @@ assert_all_ghost_zones_fully_setup();
 
 //
 // This function constructs and interlinks the patches for a
-// octant patch system.
+// +xyz octant patch system.
 //
-void patch_system::setup_octant_patch_system
+void patch_system::setup_plus_xyz_octant_patch_system
 	(int N_ghost_points, int N_extend_points,
 	 int N_overlap_points, fp delta_drho_dsigma,
 	 int min_gfn_in, int max_gfn_in,
-	 int interpolator_order)
+	 int interp_handle, int interp_par_table_handle)
 {
-using namespace patch_system_info::octant_patch_system_info;
-
-// create the patches
-patch& pz = *new z_patch(*this, plus_z.patch_number,
-			 "+z", patch::patch_is_plus,
-			 plus_z.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_z.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& px = *new x_patch(*this, plus_x.patch_number,
-			 "+x", patch::patch_is_plus,
-			 plus_x.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_x.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-patch& py = *new y_patch(*this, plus_y.patch_number,
-			 "+y", patch::patch_is_plus,
-			 plus_y.grid_array_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma),
-			 plus_y.grid_pars(N_extend_points, delta_drho_dsigma),
-			 min_gfn_in, max_gfn_in);
-
-// attach the patches to the patch system
-(*all_patches_)(plus_z.patch_number) = & pz;
-(*all_patches_)(plus_x.patch_number) = & px;
-(*all_patches_)(plus_y.patch_number) = & py;
-
 // create the ghost zones
 setup_adjacent_ghost_zones(pz, px, N_overlap_points);
 setup_adjacent_ghost_zones(pz, py, N_overlap_points);
@@ -607,9 +457,15 @@ setup_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(),
 				    true);
 
 // create the patch frontiers
-setup_adjacent_patch_frontiers(pz, px, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(pz, py, N_overlap_points, interpolator_order);
-setup_adjacent_patch_frontiers(px, py, N_overlap_points, interpolator_order);
+setup_adjacent_patch_frontiers(pz, px,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(pz, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
+setup_adjacent_patch_frontiers(px, py,
+			       N_overlap_points,
+			       interp_handle, interp_par_table_handle);
 
 assert_all_ghost_zones_fully_setup();
 }
@@ -619,63 +475,6 @@ assert_all_ghost_zones_fully_setup();
 //******************************************************************************
 
 //
-// This function creates all the patches for a patch_system (with
-// contiguous gridfn arrays), and sets pointers to them in  all_patches_[] .
-//
-// Arguments:
-// all_patch_info[] = An array, subscripted by patch number, which describes
-//		      the information which varies from one patch to another.
-//
-void create_patches(int N_patches,
-		    const struct patch_size all_patch_sizes[],
-		    int N_ghost_points, int N_extend_points,
-		    int N_overlap_points, fp delta_drho_dsigma,
-		    int min_gfn_in, int max_gfn_in,
-		    int interpolator_order)
-{
-//
-// accumulate the total number of grid points in all the patches
-//
-N_grid_points_ = 0;
-ghosted_N_grid_points_ = 0;
-	for (int pn = 0 ; pn < N_patches ; ++pn)
-	{
-	const struct grid_arrays::grid_array_pars gap&
-		= all_patch_sizes[pn].grid_array_pars(N_ghost_points,
-						      N_extend_points,
-						      delta_drho_dsigma);
-	N_grid_points_         += gap.N_grid_points();
-	ghosted_N_grid_points_ += gap.ghosted_N_grid_points();
-	}
-
-//
-// allocate the gridfn storage array
-// ... we use a single array so the patches can be contiguous in it
-//
-fp *gridfn_array = new fp[ghosted_N_grid_points_];
-
-//
-// create the patches, with gridfns stored in contiguous chunks of
-// the storage array
-//
-	for (int pn = 0 ; pn < N_patches ; ++pn)
-	{
-	const struct grid_arrays::grid_array_pars gap&
-		= all_patch_sizes[pn].grid_array_pars(N_ghost_points,
-						      N_extend_points,
-						      delta_drho_dsigma);
-	const struct grid::grid_pars gp&
-		= all_patch_sizes[pn].grid_pars(N_ghost_points,
-						N_extend_points,
-						delta_drho_dsigma);
-	}
-}
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
 // This function creates a pair of periodic-symmetry ghost zones.
 //
 //static
@@ -717,31 +516,21 @@ py.setup_interpatch_ghost_zone(ey, ex, N_overlap_points);
 
 //
 // This function automagically figures out which edges of two adjacent
-// patches are adjacent, then sets up both patches' patch frontiers on
-// those edges, and links these frontiers to the opposite patches'
-// ghost zones.
+// patches are adjacent, then sets up the appropriate pair of patch
+// frontiers and interlinks them with their ghost zones and patches.
 //
 //static
-  void patch_system::setup_adjacent_patch_frontiers(patch &px, patch &py,
-						    int N_overlap_points,
-						    int interpolator_order)
+  void patch_system::setup_adjacent_patch_frontiers
+	(patch &px, patch &py,
+	 int N_overlap_points,
+	 int interp_handle, int interp_par_table_handle)
 {
 const patch_edge& ex = px.edge_adjacent_to_patch(py, N_overlap_points);
 const patch_edge& ey = py.edge_adjacent_to_patch(px, N_overlap_points);
-interpatch_ghost_zone& igx = px.interpatch_ghost_zone_on_edge(ex);
-interpatch_ghost_zone& igy = py.interpatch_ghost_zone_on_edge(ey);
-
-patch_frontier& fx
-	= px.setup_patch_frontier(ex,
-				  igy.other_min_iperp(), igy.other_max_iperp(),
-				  interpolator_order);
-patch_frontier& fy
-	= py.setup_patch_frontier(ey,
-				  igx.other_min_iperp(), igx.other_max_iperp(),
-				  interpolator_order);
-
-igx.set_other_frontier(fy);
-igy.set_other_frontier(fx);
+px.interpatch_ghost_zone_on_edge(ex)
+  .setup_other_frontier(interp_handle, interp_par_table_handle);
+py.interpatch_ghost_zone_on_edge(ey)
+  .setup_other_frontier(interp_handle, interp_par_table_handle);
 }
 
 //******************************************************************************
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 7c2998f..fa3c4d1 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -29,18 +29,17 @@
 
 class	patch_system
 	{
-public:
 	//
 	// ***** static data & functions describing patch systems *****
 	//
-
+public:
 	// what patch-system type are supported?
 	enum patch_system_type {
 			       full_sphere_patch_system,
 			       plus_z_hemisphere_patch_system,
 			       plus_xy_quadrant_patch_system,
 			       plus_xz_quadrant_patch_system,
-			       octant_patch_system
+			       plus_xyz_octant_patch_system
 			       };
 
 	// maximum number of patches in any patch-system type
@@ -58,7 +57,7 @@ public:
 	//
 	// ***** coordinates *****
 	//
-
+public:
 	// global (x,y,z) --> local (x,y,z)
 	fp local_x_of_global_x(fp global_x) const
 		{ return global_coords_.local_x_of_global_x(global_x); }
@@ -84,7 +83,7 @@ public:
 	//
 	// ***** meta-info about the entire patch system *****
 	//
-
+public:
 	// patch-system type
 	enum patch_system_type type() const { return type_; }
 
@@ -106,20 +105,16 @@ public:
 	//
 	// ***** meta-info about gridfns *****
 	//
+public:
 	int min_gfn() const { return ith_patch(0).min_gfn(); }
 	int max_gfn() const { return ith_patch(0).max_gfn(); }
 	int N_gridfns() const { return ith_patch(0).N_gridfns(); }
 
-	//
-	// ***** gridfns *****
-	//
-	// nominal-grid
-
 
 	//
 	// ***** misc stuff *****
 	//
-
+public:
 	// "synchronize" all ghost zones of all patches,
 	// i.e. update the ghost-zone values of the specified gridfns
 	// via the appropriate sequence of symmetry operations
@@ -134,27 +129,24 @@ public:
 			  FILE *output_fp = stdout) const;
 
 
-public:
 	//
 	// ***** constructor, destructor *****
 	//
 	// This constructor doesn't support the full generality of the
 	// patch data structures (which would, eg, allow N_ghost_points
-	// and N_overlap_points and interpolator_order to vary from
+	// and N_extend_points and interpolator_order to vary from
 	// ghost zone to ghost zone, and the grid spacings to vary
 	// from patch to patch.  But in practice we'd probably never
 	// use that generality...
 	//
-	// The present implementation also restricts N_overlap_points
-	// to be even.  N.b. the semantics of N_overlap points are
-	// as documented in patch.hh for patch::edge_adjacent_to_patch().
-	//
+public:
 	patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
 		     enum patch_system_type type_in,
-		     int N_ghost_points, int N_overlap_points,
+		     int N_ghost_points, int N_extend_points,
 		     fp delta_drho_dsigma,
-		     int min_gfn_in, int max_gfn_in,
-		     int interpolator_order);
+		     int         min_gfn_in, int         max_gfn_in,
+		     int ghosted_min_gfn_in, int ghosted_max_gfn_in,
+		     int interp_handle_in, int interp_par_table_handle_in);
 	~patch_system();
 
 
@@ -162,71 +154,40 @@ public:
 	// ***** helper functions for constructor *****
 	//
 private:
-
-	//
-	// this (POD, and hence static-initializable) struct gives a
-	// minimal set of information which varies from one patch to
-	// another
-	//
-	struct	patch_info
-		{
-		const char *name,
-		fp min_drho, max_drho;
-		fp min_dsigma, max_dsigma;
-
-		// compute and return reference to
-		//    struct grid_arrays::grid_array_pars
-		// from the info in this struct
-		// ... result refers to internal static buffer;
-		//     the usual caveats about lifetimes/overwriting apply
-		const grid_arrays::grid_array_pars& grid_array_pars
-			(int N_ghost_points, int N_extend_points,
-			 fp delta_drho_dsigma)
-			const;
-
-		// compute and return reference to  struct grid::grid_pars
-		// from the info in this struct
-		// ... result refers to internal static buffer;
-		//     the usual caveats about lifetimes/overwriting apply
-		const grid::grid_pars& grid_pars
-			(int N_extend_points, fp delta_drho_dsigma)
-			const;
-		};
+	// construct patches as described by patch_info[] array,
+	// and link them into the patch system
+	// does *NOT* create ghost zones or frontiers
+	// does *NOT* set up gridfns
+	void construct_patches(const struct patch_info patch_info[],
+			       int N_ghost_points, int N_extend_points,
+			       fp delta_drho_dsigma);
 
 	// create/interlink patches, ghost zones, and frontiers
-	void setup_full_sphere_patch_system
+	void interlink_full_sphere_patch_system
 		(int N_ghost_points, int N_extend_points,
 		 int N_overlap_points, fp delta_drho_dsigma,
 		 int min_gfn_in, int max_gfn_in,
-		 int interpolator_order);
-	void setup_plus_z_hemisphere_patch_system
+		 int interp_handle, int interp_par_table_handle);
+	void interlink_plus_z_hemisphere_patch_system
 		(int N_ghost_points, int N_extend_points,
 		 int N_overlap_points, fp delta_drho_dsigma,
 		 int min_gfn_in, int max_gfn_in,
-		 int interpolator_order);
-	void setup_plus_xy_quadrant_patch_system
+		 int interp_handle, int interp_par_table_handle);
+	void interlink_plus_xy_quadrant_patch_system
 		(int N_ghost_points, int N_extend_points,
 		 int N_overlap_points, fp delta_drho_dsigma,
 		 int min_gfn_in, int max_gfn_in,
-		 int interpolator_order);
-	void setup_plus_xz_quadrant_patch_system
+		 int interp_handle, int interp_par_table_handle);
+	void interlink_plus_xz_quadrant_patch_system
 		(int N_ghost_points, int N_extend_points,
 		 int N_overlap_points, fp delta_drho_dsigma,
 		 int min_gfn_in, int max_gfn_in,
-		 int interpolator_order);
-	void setup_octant_patch_system
+		 int interp_handle, int interp_par_table_handle);
+	void interlink_plus_xyz_octant_patch_system
 		(int N_ghost_points, int N_extend_points,
 		 int N_overlap_points, fp delta_drho_dsigma,
 		 int min_gfn_in, int max_gfn_in,
-		 int interpolator_order);
-
-	// create patches, set pointers in all_patches_[]
-	void create_patches(int N_patches,
-			    const struct patch_size all_patch_sizes[],
-			    int N_ghost_points, int N_extend_points,
-			    int N_overlap_points, fp delta_drho_dsigma,
-			    int min_gfn_in, int max_gfn_in,
-			    int interpolator_order);
+		 int interp_handle, int interp_par_table_handle);
 
 	// create/interlink a pair of periodic-symmetry ghost zones
 	static
@@ -243,9 +204,10 @@ private:
 	// create/interlink a pair of patch frontiers
 	// ... automagically figures out which edges are adjacent
 	static
-	  void setup_adjacent_patch_frontiers(patch &px, patch &py,
-					      int N_overlap_points,
-					      int interpolator_order);
+	  void setup_adjacent_patch_frontiers
+		(patch &px, patch &py,
+		 int N_overlap_points,
+		 int interp_handle, int interp_par_table_handle);
 
 	// assert() that all ghost zones (and frontiers, where applicable)
 	// of all patches are fully setup
@@ -259,12 +221,6 @@ private:
 	patch_system& operator=(const patch_system &rhs);
 
 private:
-	// actual storage arrays for grid functions
-	// ... array subscripts are [gfn][1d_grid_point_index]
-	// ... FIXME: should be boost::multi_array<> when I get this working
-	array2d<fp>* nominal_gridfns_;
-	array2d<fp>* ghosted_gridfns_;
-
 	// local <--> global coordinate mapping
 	global_coords global_coords_;