// patch_system.cc -- describes the (an) entire system of patches
// $Id$
//
// patch_system::patch_system
// patch_system::~patch_system
// patch_system::create_patches
// patch_system::setup_gridfn_storage
// 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_interpatch_ghost_zones
// patch_system::finish_interpatch_setup
// patch_system::assert_all_ghost_zones_fully_setup
//
// patch_system::patch_irho_isigma_of_gpn
// patch_system::ghosted_patch_irho_isigma_of_gpn
//
// patch_system::N_patches_of_type
// patch_system::name_of_type
// patch_system::type_of_name
// patch_system::patch_number_of_name
//
// patch_system::synchronize
// patch_system::compute_synchronize_Jacobian
// patch_system::synchronize_Jacobian_y_ipar_posn
// patch_system::synchronize_Jacobian
//
// patch_system::print_unknown_gridfn
// patch_system::read_unknown_gridfn
//

#include <cstdio>
#include <cmath>
#include <cstring>
#include <vector>
#include <assert.h>
using std::printf;
using std::strcmp;

#include "cctk.h"

#include "stdc.h"
#include "config.hh"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

#include "coords.hh"
#include "grid.hh"
#include "fd_grid.hh"
#include "patch.hh"
#include "patch_edge.hh"
#include "patch_interp.hh"
#include "ghost_zone.hh"
#include "patch_info.hh"
#include "patch_system.hh"
#include "patch_system_info.hh"

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function constructs a  patch_system  object.
//
// Constructor arguments:
// N_ghost_points = Width in grid points of all ghost zones.
// N_overlap_points = Number of grid points that adjacent
//		      nominally-just-touching patches should overlap.
//		      For example, with N_overlap_points == 3, here
//		      are the grid points of two neighboring patches:
//			x   x   x   x   x   X
//                                      |
//			            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.  For this example, the N_extend_points
//		      parameter used by some other functions would
//		      be 1; in general
//			N_overlap_points = 2*N_extend_points + 1
// 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,
			   fp delta_drho_dsigma,
			   int min_gfn_in, int max_gfn_in,
			   int ghosted_min_gfn_in, int ghosted_max_gfn_in,
			   int interp_handle, int interp_par_table_handle)

	: global_coords_(origin_x_in, origin_y_in, origin_z_in),
	  type_(type_in),
	  N_patches_(N_patches_of_type(type_in)),
	  all_patches_(N_patches_),
	  starting_gpn_(N_patches_+1),
	  ghosted_starting_gpn_(N_patches_+1),
	  gridfn_storage_(NULL),	// set in setup_gridfn_storage()
	  ghosted_gridfn_storage_(NULL)	// set in setup_gridfn_storage()
{
if (! jtutil::is_odd(N_overlap_points))
   then error_exit(ERROR_EXIT,
"***** patch_system::patch_system(): implementation restriction:\n"
"        N_overlap_points=%d, but we only support odd values!\n"
,
		   N_overlap_points);				/*NOTREACHED*/
const int N_extend_points = N_overlap_points >> 1;

// construct/interlink the patches and ghost zones
switch	(type_in)
	{
case full_sphere_patch_system:
	create_patches(patch_system_info::full_sphere::patch_info_array,
		       N_ghost_points, N_extend_points,
		       delta_drho_dsigma);
	setup_gridfn_storage(min_gfn_in, max_gfn_in,
			     ghosted_min_gfn_in, ghosted_max_gfn_in);
	interlink_full_sphere_patch_system(N_overlap_points,
					   interp_handle,
					   interp_par_table_handle);
	break;
case plus_z_hemisphere_patch_system:
	create_patches(patch_system_info::plus_z_hemisphere::patch_info_array,
		       N_ghost_points, N_extend_points,
		       delta_drho_dsigma);
	setup_gridfn_storage(min_gfn_in, max_gfn_in,
			     ghosted_min_gfn_in, ghosted_max_gfn_in);
	interlink_plus_z_hemisphere_patch_system(N_overlap_points,
						 interp_handle,
						 interp_par_table_handle);
	break;
case plus_xy_quadrant_patch_system:
	create_patches(patch_system_info::plus_xy_quadrant::patch_info_array,
		       N_ghost_points, N_extend_points,
		       delta_drho_dsigma);
	setup_gridfn_storage(min_gfn_in, max_gfn_in,
			     ghosted_min_gfn_in, ghosted_max_gfn_in);
	interlink_plus_xy_quadrant_patch_system(N_overlap_points,
						interp_handle,
						interp_par_table_handle);
	break;
case plus_xz_quadrant_patch_system:
	create_patches(patch_system_info::plus_xz_quadrant::patch_info_array,
		       N_ghost_points, N_extend_points,
		       delta_drho_dsigma);
	setup_gridfn_storage(min_gfn_in, max_gfn_in,
			     ghosted_min_gfn_in, ghosted_max_gfn_in);
	interlink_plus_xz_quadrant_patch_system(N_overlap_points,
						interp_handle,
						interp_par_table_handle);
	break;
case plus_xyz_octant_patch_system:
	create_patches(patch_system_info::plus_xyz_octant::patch_info_array,
		       N_ghost_points, N_extend_points,
		       delta_drho_dsigma);
	setup_gridfn_storage(min_gfn_in, max_gfn_in,
			     ghosted_min_gfn_in, ghosted_max_gfn_in);
	interlink_plus_xyz_octant_patch_system(N_overlap_points,
					       interp_handle,
					       interp_par_table_handle);
	break;
default:
	error_exit(ERROR_EXIT,
"***** patch_system::patch_system(): bad type_in=(int)%d!\n"
,
		   int(type_in));				/*NOTREACHED*/
	}
}

//******************************************************************************

//
// This function destroys a  patch_system  object.
//
patch_system::~patch_system()
{
delete ghosted_gridfn_storage_;
delete gridfn_storage_;

	for (int pn = N_patches()-1 ; pn >= 0 ; --pn)
	{
	delete &ith_patch(pn);
	}
}

//******************************************************************************

//
// This function is called from the patch_system:: constructor to
// construct a set of patches as described by an array of patch_info
// structures and associated arguments, and make these patches members
// of this patch system.  This function also correctly sets
//	N_grid_points_
//	N_ghosted_grid_points_
//	all_patches_[]
//	starting_gpn_[]
//	ghosted_starting_gpn_[]
// This function does *NOT* create any of the ghost zones, and does
// *NOT* set up any gridfns.
//
void patch_system::create_patches(const struct patch_info patch_info_in[],
				  int N_ghost_points, int N_extend_points,
				  fp delta_drho_dsigma)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "constructing %s patch system",
	   name_of_type(type()));
N_grid_points_ = 0;
ghosted_N_grid_points_ = 0;
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	const struct patch_info& pi = patch_info_in[pn];

	CCTK_VInfo(CCTK_THORNSTRING,
		   "   constructing %s patch",
		   pi.name);

	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));
		break;
	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));
		break;
	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::create_patches():\n"
"        unknown patch_info_in[pn=%d].ctype=0x%02d='%c'!\n"
,
			   pn, pi.ctype, pi.ctype);		/*NOTREACHED*/
		}

	// these record number of grid points in *previous* patches,
	// i.e. they do *not* include the number of grid points in this patch
	starting_gpn_[pn] = N_grid_points_;
	ghosted_starting_gpn_[pn] = ghosted_N_grid_points_;

	N_grid_points_+= p->N_grid_points();
	ghosted_N_grid_points_ += p->ghosted_N_grid_points();

	all_patches_[pn] = p;
	}

starting_gpn_        [N_patches_] = N_grid_points_;
ghosted_starting_gpn_[N_patches_] = ghosted_N_grid_points_;
}

//******************************************************************************

//
// This function is called from the patch_system:: constructor to set
// up the storage for all gridfns in all patches, giving each gridfn a
// contiguous-across-all-patches storage array.  This function also makes
// a number of self-consistency checks to ensure that the gridfn storage
// subscripting is set up properly.
//
// This function assumes that all the patches have already been constructed
// before it is called.
//
// For example, given the patches {x,y,z}, the ghosted gridfns {H,J},
// and the nominal gridfns {a,b,c}, we might picture the storage like
// this:
//
//	xa xa xa ya ya za za za za
//	xb xb xb yb yb zb zb zb zb
//	xc xc xc yc yc zc zc zc zc
//
//	xH xH xH xH yH yH yH zH zH zH zH zH
//	xJ xJ xJ xJ yJ yJ yJ zJ zJ zJ zJ zJ
//
// Here the upper/lower blocks are for nominal/ghosted  gridfns.
// The storage is taken as being contiguous within each row (in fact
// within each block).  Thus the storage for all the nominal gridfns
// (or all the ghosted gridfns) in a single patch is *non*-contiguous.
//
// The creation of patches is done in several phases: first the patches
// are constructed with no gridfn storage, then we are called to set up
// the gridfn storage (taking into account the sizes of the other patches),
// then finally ghost zones are constructed and interlinked.
//
// FIXME: We should pad the gridfn storage as necessary to avoid cache
//        conflicts, but we don't do this at present.
//
void patch_system::setup_gridfn_storage
	(int min_gfn_in, int max_gfn_in,
	 int ghosted_min_gfn_in, int ghosted_max_gfn_in)
{
const int N_gridfns_in = jtutil::how_many_in_range(min_gfn_in, max_gfn_in);
const int ghosted_N_gridfns_in
	= jtutil::how_many_in_range(ghosted_min_gfn_in, ghosted_max_gfn_in);

const int gfn_stride = N_grid_points();
const int ghosted_gfn_stride = ghosted_N_grid_points();

const int N_storage = gfn_stride * N_gridfns_in;
const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in;

CCTK_VInfo(CCTK_THORNSTRING, "setting up gridfn storage");
CCTK_VInfo(CCTK_THORNSTRING,
	   "   gfn=[%d,%d] ghosted_gfn=[%d,%d]",
	   min_gfn_in, max_gfn_in,
	   ghosted_min_gfn_in, ghosted_max_gfn_in);
CCTK_VInfo(CCTK_THORNSTRING,
	   "   N_grid_points()=%d ghosted_N_grid_points()=%d",
	   N_grid_points(), ghosted_N_grid_points());

// storage arrays for all gridfns
gridfn_storage_         = new fp[N_storage];
ghosted_gridfn_storage_ = new fp[ghosted_N_storage];

// divide up the storage array among the patches
// and set up the storage in the individual patches themselves
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	const int         posn =         starting_gpn_[pn];
	const int ghosted_posn = ghosted_starting_gpn_[pn];
	const struct grid_arrays::gridfn_pars gridfn_pars
		= {
		  min_gfn_in, max_gfn_in,
		  & gridfn_storage_[posn],
		  gfn_stride, 0, 0
		  };
	const struct grid_arrays::gridfn_pars ghosted_gridfn_pars
		= {
		  ghosted_min_gfn_in, ghosted_max_gfn_in,
		  & ghosted_gridfn_storage_[ghosted_posn],
		  ghosted_gfn_stride, 0, 0
		  };

	patch& p = ith_patch(pn);
	p.setup_gridfn_storage(gridfn_pars, ghosted_gridfn_pars);
	}

CCTK_VInfo(CCTK_THORNSTRING,
	   "   checking that storage is partitioned properly");

// check to make sure storage for distinct gridfns
// forms a partition of the overall storage array
const patch& pfirst = ith_patch(0);
const patch& plast  = ith_patch(N_patches()-1);
	for (int gfn = min_gfn() ; gfn+1 < max_gfn() ; ++gfn)
	{
	// range of storage occupied by gridfns:
	//	gfn   --> [gfn_first, gfn_last]
	//	gfn+1 --> [gfn1_first, gfn1_last]
	const fp* const gfn_last_ptr
		= & plast.gridfn(gfn, plast.max_irho(),
				      plast.max_isigma());
	const fp* const gfn1_first_ptr
		= & pfirst.gridfn(gfn+1, pfirst.min_irho(),
					 pfirst.min_isigma());
	if (! (gfn1_first_ptr == gfn_last_ptr+1))
	   then error_exit(PANIC_EXIT,
"***** patch_system::setup_gridfn_storage():\n"
"        nominal-grid gridfns don't partition overall storage array!"
"        (this should never happen!)\n"
"        gfn=%d last point at gfn_last_ptr=%p\n"
"	 gfn+1=%d first point at gfn1_first_ptr=%p\n"
"        should have gfn1_first_ptr == gfn_last_ptr+1\n"
,
			   gfn, (const void *) gfn_last_ptr,
			   gfn+1, (const void *) gfn1_first_ptr); /*NOTREACHED*/
	}

	for (int ghosted_gfn = ghosted_min_gfn() ;
	     ghosted_gfn+1 < ghosted_max_gfn() ;
	     ++ghosted_gfn)
	{
	// range of storage occupied by ghosted gridfns:
	//	ghosted_gfn   --> [gfn_first, gfn_last]
	//	ghosted_gfn+1 --> [gfn1_first, gfn1_last]
	const fp* const ghosted_gfn_last_ptr
		= & plast.ghosted_gridfn(ghosted_gfn,
					 plast.ghosted_max_irho(),
					 plast.ghosted_max_isigma());
	const fp* const ghosted_gfn1_first_ptr
		= & pfirst.ghosted_gridfn(ghosted_gfn+1,
					  pfirst.ghosted_min_irho(),
					  pfirst.ghosted_min_isigma());
	if (! (ghosted_gfn1_first_ptr == ghosted_gfn_last_ptr+1))
	   then error_exit(PANIC_EXIT,
"***** patch_system::setup_gridfn_storage():\n"
"        ghosted-grid gridfns don't partition overall storage array!"
"        (this should never happen!)\n"
"        ghosted_gfn=%d last point at ghosted_gfn_last_ptr=%p\n"
"	 ghosted_gfn+1=%d first point at ghosted_gfn1_first_ptr=%p\n"
"        should have ghosted_gfn1_first_ptr == ghosted_gfn_last_ptr+1\n"
,
			   ghosted_gfn, (const void *) ghosted_gfn_last_ptr,
			   ghosted_gfn+1,
			      (const void *) ghosted_gfn1_first_ptr);
								/*NOTREACHED*/
	}

// check to make sure storage for distinct patches
// forms a partition of the storage for each gridfn
	for (int gfn = min_gfn() ; gfn < max_gfn() ; ++gfn)
	{
		for (int pn = 0 ; pn+1 < N_patches() ; ++pn)
		{
		const patch& p = ith_patch(pn);
		const patch& p1 = ith_patch(pn+1);

		// range of storage occupied by gridfn:
		//	p  --> [p_first, p_last]
		//	p1 --> [p1_first, p1_last]
		const fp* const p_last_ptr
			= & p.gridfn(gfn, p.max_irho(), p.max_isigma());
		const fp* const p1_first_ptr
			= & p1.gridfn(gfn, p1.min_irho(), p1.min_isigma());
		if (! (p1_first_ptr == p_last_ptr+1))
		   then error_exit(PANIC_EXIT,
"***** patch_system::setup_gridfn_storage():\n"
"        nominal-grid patches gridfns don't partition storage for gfn=%d!\n"
"        (this should never happen!)\n"
"        gfn=%d %s patch last point at p_last_ptr=%p\n"
"	 gfn=%d %s patch first point at p1_first_ptr=%p\n"
"        should have p1_first_ptr == p_last_ptr+1\n"
,
			   gfn,
			   gfn, p.name(), (const void *) p_last_ptr,
			   gfn+1, p1.name(), (const void *) p1_first_ptr);
								 /*NOTREACHED*/
		}
	}

	for (int ghosted_gfn = ghosted_min_gfn() ;
	     ghosted_gfn < ghosted_max_gfn() ;
	     ++ghosted_gfn)
	{
		for (int pn = 0 ; pn+1 < N_patches() ; ++pn)
		{
		const patch& p = ith_patch(pn);
		const patch& p1 = ith_patch(pn+1);

		// range of storage occupied by ghosted gridfn:
		//	p  --> [p_first, p_last]
		//	p1 --> [p1_first, p1_last]
		const fp* const p_last_ptr
			= & p.ghosted_gridfn(ghosted_gfn,
					     p.ghosted_max_irho(),
					     p.ghosted_max_isigma());
		const fp* const p1_first_ptr
			= & p1.ghosted_gridfn(ghosted_gfn,
					      p1.ghosted_min_irho(),
					      p1.ghosted_min_isigma());
		if (! (p1_first_ptr == p_last_ptr+1))
		   then error_exit(PANIC_EXIT,
"***** patch_system::setup_gridfn_storage():\n"
"        nominal-grid patches gridfns don't partition storage for gfn=%d!\n"
"        (this should never happen!)\n"
"        %s patch (pn=%d) last point at p_last_ptr=%p\n"
"	 %s patch (pn=%d) first point at p1_first_ptr=%p\n"
"        should have p1_first_ptr == p_last_ptr+1\n"
,
			   ghosted_gfn,
			   p.name(), pn, (const void *) p_last_ptr,
			   p1.name(), pn+1, (const void *) p1_first_ptr);
								 /*NOTREACHED*/
		}
	}
}

//******************************************************************************

//
// This function constructs and interlinks the ghost zones for a
// full-sphere patch system.
//
void patch_system::interlink_full_sphere_patch_system
	(int N_overlap_points,
	 int interp_handle, int interp_par_table_handle)
{
CCTK_VInfo(CCTK_THORNSTRING, "interlinking full sphere patch system");

patch& pz = ith_patch(patch_number_of_name("+z"));
patch& px = ith_patch(patch_number_of_name("+x"));
patch& py = ith_patch(patch_number_of_name("+y"));
patch& mx = ith_patch(patch_number_of_name("-x"));
patch& my = ith_patch(patch_number_of_name("-y"));
patch& mz = ith_patch(patch_number_of_name("-z"));

// create the ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
create_interpatch_ghost_zones(pz, px, N_overlap_points);
create_interpatch_ghost_zones(pz, py, N_overlap_points);
create_interpatch_ghost_zones(pz, mx, N_overlap_points);
create_interpatch_ghost_zones(pz, my, N_overlap_points);
create_interpatch_ghost_zones(px, py, N_overlap_points);
create_interpatch_ghost_zones(py, mx, N_overlap_points);
create_interpatch_ghost_zones(mx, my, N_overlap_points);
create_interpatch_ghost_zones(my, px, N_overlap_points);
create_interpatch_ghost_zones(mz, px, N_overlap_points);
create_interpatch_ghost_zones(mz, py, N_overlap_points);
create_interpatch_ghost_zones(mz, mx, N_overlap_points);
create_interpatch_ghost_zones(mz, my, N_overlap_points);

// finish setting up the interpatch ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
finish_interpatch_setup(pz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, mx,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(px, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(py, mx,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mx, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(my, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mz, mx,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mz, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);

assert_all_ghost_zones_fully_setup();
}

//******************************************************************************

//
// This function constructs and interlinks the ghost zones for a
// +z hemisphere patch system.
//
void patch_system::interlink_plus_z_hemisphere_patch_system
	(int N_overlap_points,
	 int interp_handle, int interp_par_table_handle)
{
CCTK_VInfo(CCTK_THORNSTRING, "interlinking +z hemisphere patch system");

patch& pz = ith_patch(patch_number_of_name("+z"));
patch& px = ith_patch(patch_number_of_name("+x"));
patch& py = ith_patch(patch_number_of_name("+y"));
patch& mx = ith_patch(patch_number_of_name("-x"));
patch& my = ith_patch(patch_number_of_name("-y"));

// create the ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
create_interpatch_ghost_zones(pz, px, N_overlap_points);
create_interpatch_ghost_zones(pz, py, N_overlap_points);
create_interpatch_ghost_zones(pz, mx, N_overlap_points);
create_interpatch_ghost_zones(pz, my, N_overlap_points);
create_interpatch_ghost_zones(px, py, N_overlap_points);
create_interpatch_ghost_zones(py, mx, N_overlap_points);
create_interpatch_ghost_zones(mx, my, N_overlap_points);
create_interpatch_ghost_zones(my, px, N_overlap_points);
px.create_mirror_symmetry_ghost_zone(px.max_rho_patch_edge());
py.create_mirror_symmetry_ghost_zone(py.max_rho_patch_edge());
mx.create_mirror_symmetry_ghost_zone(mx.min_rho_patch_edge());
my.create_mirror_symmetry_ghost_zone(my.min_rho_patch_edge());

// finish setting up the interpatch ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
finish_interpatch_setup(pz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, mx,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(px, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(py, mx,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mx, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(my, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);

assert_all_ghost_zones_fully_setup();
}

//******************************************************************************

//
// This function constructs and interlinks the ghost zones for a
// +xy quadrant patch system.
//
void patch_system::interlink_plus_xy_quadrant_patch_system
	(int N_overlap_points,
	 int interp_handle, int interp_par_table_handle)
{
CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xy quadrant patch system");

patch& pz = ith_patch(patch_number_of_name("+z"));
patch& px = ith_patch(patch_number_of_name("+x"));
patch& py = ith_patch(patch_number_of_name("+y"));
patch& mz = ith_patch(patch_number_of_name("-z"));

// create the ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
create_interpatch_ghost_zones(pz, px, N_overlap_points);
create_interpatch_ghost_zones(pz, py, N_overlap_points);
create_interpatch_ghost_zones(px, py, N_overlap_points);
create_interpatch_ghost_zones(mz, px, N_overlap_points);
create_interpatch_ghost_zones(mz, py, N_overlap_points);
create_periodic_symmetry_ghost_zones(pz.min_rho_patch_edge(),
				     pz.min_sigma_patch_edge(),
				     true);
create_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(),
				     py.max_sigma_patch_edge(),
				     true);
create_periodic_symmetry_ghost_zones(mz.max_rho_patch_edge(),
				     mz.max_sigma_patch_edge(),
				     true);

// finish setting up the interpatch ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
finish_interpatch_setup(pz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(px, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(mz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);

assert_all_ghost_zones_fully_setup();
}

//******************************************************************************

//
// This function constructs and interlinks the ghost zones for a
// +xz quadrant patch system.
//
void patch_system::interlink_plus_xz_quadrant_patch_system
	(int N_overlap_points,
	 int interp_handle, int interp_par_table_handle)
{
CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xz quadrant patch system");

patch& pz = ith_patch(patch_number_of_name("+z"));
patch& px = ith_patch(patch_number_of_name("+x"));
patch& py = ith_patch(patch_number_of_name("+y"));
patch& my = ith_patch(patch_number_of_name("-y"));

// create the ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
create_interpatch_ghost_zones(pz, px, N_overlap_points);
create_interpatch_ghost_zones(pz, py, N_overlap_points);
create_interpatch_ghost_zones(pz, my, N_overlap_points);
create_interpatch_ghost_zones(px, py, N_overlap_points);
create_interpatch_ghost_zones(px, my, N_overlap_points);
px.create_mirror_symmetry_ghost_zone(px.max_rho_patch_edge());
py.create_mirror_symmetry_ghost_zone(py.max_rho_patch_edge());
my.create_mirror_symmetry_ghost_zone(my.min_rho_patch_edge());
create_periodic_symmetry_ghost_zones(pz.min_sigma_patch_edge(),
				     pz.min_sigma_patch_edge(),
				     false);
create_periodic_symmetry_ghost_zones(py.max_sigma_patch_edge(),
				     my.min_sigma_patch_edge(),
				     false);

// finish setting up the interpatch ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
finish_interpatch_setup(pz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(px, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(px, my,
			N_overlap_points,
			interp_handle, interp_par_table_handle);

assert_all_ghost_zones_fully_setup();
}

//******************************************************************************

//
// This function constructs and interlinks the ghost zones for a
// +xyz octant patch system.
//
void patch_system::interlink_plus_xyz_octant_patch_system
	(int N_overlap_points,
	 int interp_handle, int interp_par_table_handle)
{
CCTK_VInfo(CCTK_THORNSTRING, "interlinking +xyz octant patch system");

patch& pz = ith_patch(patch_number_of_name("+z"));
patch& px = ith_patch(patch_number_of_name("+x"));
patch& py = ith_patch(patch_number_of_name("+y"));

// create the ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   creating ghost zones");
create_interpatch_ghost_zones(pz, px, N_overlap_points);
create_interpatch_ghost_zones(pz, py, N_overlap_points);
create_interpatch_ghost_zones(px, py, N_overlap_points);
px.create_mirror_symmetry_ghost_zone(px.max_rho_patch_edge());
py.create_mirror_symmetry_ghost_zone(py.max_rho_patch_edge());
create_periodic_symmetry_ghost_zones(pz.min_rho_patch_edge(),
				     pz.min_sigma_patch_edge(),
				     true);
create_periodic_symmetry_ghost_zones(px.min_sigma_patch_edge(),
				     py.max_sigma_patch_edge(),
				     true);

// finish setting up the interpatch ghost zones
CCTK_VInfo(CCTK_THORNSTRING, "   finishing interpatch setup");
finish_interpatch_setup(pz, px,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(pz, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);
finish_interpatch_setup(px, py,
			N_overlap_points,
			interp_handle, interp_par_table_handle);

assert_all_ghost_zones_fully_setup();
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function creates a pair of periodic-symmetry ghost zones.
//
//static
  void patch_system::create_periodic_symmetry_ghost_zones
	(const patch_edge& ex, const patch_edge& ey,
	 bool ipar_map_is_plus)
{
ex.my_patch()
  .create_periodic_symmetry_ghost_zone(ex, ey, ipar_map_is_plus);

if (ex == ey)
   then {
	// ex and ey are the same edge (i.e. the symmetry maps the edge
	// back to itself), so we only want to set up the edge once
	// ==> no-op here
	}
   else ey.my_patch()
	  .create_periodic_symmetry_ghost_zone(ey, ex, ipar_map_is_plus);
}

//******************************************************************************

//
// This function automagically figures out which edges of two adjacent
// patches are adjacent, then creates both patches' ghost zones on those
// edges and interlinks them with their respective patches.
//
//static
  void patch_system::create_interpatch_ghost_zones(patch &px, patch &py,
						   int N_overlap_points)
{
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);

px.create_interpatch_ghost_zone(ex, ey, N_overlap_points);
py.create_interpatch_ghost_zone(ey, ex, N_overlap_points);
}

//******************************************************************************

//
// This function automagically figures out which edges of two adjacent
// patches are adjacent, then finishes setting up both ghost zones.
//
//static
  void patch_system::finish_interpatch_setup
	(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);
px.interpatch_ghost_zone_on_edge(ex)
  .finish_setup(interp_handle, interp_par_table_handle);
py.interpatch_ghost_zone_on_edge(ey)
  .finish_setup(interp_handle, interp_par_table_handle);
}

//******************************************************************************

//
// This function assert()s that all ghost zones of all patches have
// been fully set up.
//
void patch_system::assert_all_ghost_zones_fully_setup() const
{
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	ith_patch(pn).assert_all_ghost_zones_fully_setup();
	}
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function decodes a 0-origin grid point number into a
// (patch,irho,isigma) triple.
//
// Arguments:
// gpn = The grid point number to decode.
// (irho,isigma) = (out) The decoded patch coordinates.
//
// Results:
// This function returns a reference to the decoded patch.  (An alternative
// design would be to return this via a  patch*&  argument, but the design
// used here seems slightly cleaner to use in practice.)
//
patch& patch_system::patch_irho_isigma_of_gpn(int gpn, int& irho, int& isigma)
	const
{
assert( gpn >= 0 );
assert( gpn < N_grid_points() );

	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	// n.b. [pn+1] is ok since starting_gpn_[] has size N_patches()+1
	if ((gpn >= starting_gpn_[pn])  && (gpn < starting_gpn_[pn+1]))
	   then {
		patch& p = ith_patch(pn);
		const int gpn_in_patch = gpn - starting_gpn_[pn];
		assert( gpn_in_patch >= 0 );
		assert( gpn_in_patch < p.N_grid_points() );
		p.irho_isigma_of_gpn(gpn_in_patch, irho,isigma);
		return p;
		}
	}

error_exit(PANIC_EXIT,
"***** patch_system::patch_irho_isigma_of_gpn(gpn=%d):\n"
"        couldn't find any patch! (this should never happen!)\n"
"        N_grid_points()=%d\n"
,
	   gpn,
	   N_grid_points());					/*NOTREACHED*/
}

//******************************************************************************

//
// This function decodes a 0-origin grid point number into a *ghosted*
// (patch,irho,isigma) triple.
//
// Arguments:
// gpn = The grid point number to decode.
// (irho,isigma) = (out) The decoded patch coordinates.
//
// Results:
// This function returns a reference to the decoded patch.  (An alternative
// design would be to return this via a  patch*&  argument, but the design
// used here seems slightly cleaner to use in practice.)
//
patch& patch_system::ghosted_patch_irho_isigma_of_gpn
	(int gpn, int& irho, int& isigma)
	const
{
assert( gpn >= 0 );
assert( gpn < ghosted_N_grid_points() );

	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	// n.b. [pn+1] is ok since ghosted_starting_gpn_[]
	//      has size N_patches()+1
	if (    (gpn >= ghosted_starting_gpn_[pn])
	     && (gpn <  ghosted_starting_gpn_[pn+1]))
	   then {
		patch& p = ith_patch(pn);
		const int gpn_in_patch = gpn - ghosted_starting_gpn_[pn];
		assert( gpn_in_patch >= 0 );
		assert( gpn_in_patch < p.ghosted_N_grid_points() );
		p.ghosted_irho_isigma_of_gpn(gpn_in_patch, irho,isigma);
		return p;
		}
	}

error_exit(PANIC_EXIT,
"***** patch_system::ghosted_patch_irho_isigma_of_gpn(gpn=%d):\n"
"        couldn't find any patch! (this should never happen!)\n"
"        ghosted_N_grid_points()=%d\n"
,
	   gpn,
	   ghosted_N_grid_points());				/*NOTREACHED*/
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function decodes a patch system's type into N_patches.
//
//static
  int patch_system::N_patches_of_type(enum patch_system_type type_in)
{
switch	(type_in)
	{
case full_sphere_patch_system:
	return patch_system_info::full_sphere::N_patches;
case plus_z_hemisphere_patch_system:
	return patch_system_info::plus_z_hemisphere::N_patches;
case plus_xy_quadrant_patch_system:
	return patch_system_info::plus_xy_quadrant::N_patches;
case plus_xz_quadrant_patch_system:
	return patch_system_info::plus_xz_quadrant::N_patches;
case plus_xyz_octant_patch_system:
	return patch_system_info::plus_xyz_octant::N_patches;
default:
	error_exit(PANIC_EXIT,
"***** patch_system::N_patches_of_type(): bad type=(int)%d!\n"
"                                         (this should never happen!)\n",
		   int(type_in));				/*NOTREACHED*/
	}
}

//******************************************************************************

//
// This function decodes a patch system's type into a human-readable
// type name (a C string).
//
//static
  const char* patch_system::name_of_type(enum patch_system_type type_in)
{
switch	(type_in)
	{
case full_sphere_patch_system:		return "full sphere";
case plus_z_hemisphere_patch_system:	return "+z hemisphere";
case plus_xy_quadrant_patch_system:	return "+xy quadrant";
case plus_xz_quadrant_patch_system:	return "+xz quadrant";
case plus_xyz_octant_patch_system:	return "+xyz octant";
default:
	error_exit(PANIC_EXIT,
"***** patch_system::name_of_type(): bad type=(int)%d!\n"
"                                    (this should never happen!)\n",
		   int(type_in));				/*NOTREACHED*/
	}
}

//******************************************************************************

//
// This function encodes a human-readable type name (a C string) into
// a patch system's type into.
//
//static
  enum patch_system::patch_system_type
    patch_system::type_of_name(const char* name_in)
{
if	(STRING_EQUAL(name_in, "full sphere"))
   then return full_sphere_patch_system;
else if (STRING_EQUAL(name_in, "+z hemisphere"))
   then return plus_z_hemisphere_patch_system;
else if (STRING_EQUAL(name_in, "+xy quadrant"))
   then return plus_xy_quadrant_patch_system;
else if (STRING_EQUAL(name_in, "+xz quadrant"))
   then return plus_xz_quadrant_patch_system;
else if (STRING_EQUAL(name_in, "+xyz octant"))
   then return plus_xyz_octant_patch_system;
else	error_exit(PANIC_EXIT,
"***** patch_system::type_of_name(): unknown name=\"%s\"!",
		   name_in);					/*NOTREACHED*/
}

//******************************************************************************

//
// This function finds a patch from its human-readable name, and returns
// the patch number, or does an error_exit() if no patch is found with
// the specified name.
//
int patch_system::patch_number_of_name(const char* name) const
{
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	if (STRING_EQUAL(ith_patch(pn).name(), name))
	   then return pn;
	}

error_exit(ERROR_EXIT,
"***** patch_system::patch_number_of_name():\n"
"        no patch with name \"%s\"!\n",
	   name);						/*NOTREACHED*/
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function "synchronizes" all ghost zones of all patches, i.e. it
// update the ghost-zone values of the specified gridfns via the appropriate
// sequence of symmetry operations and interpatch interpolations.  This
// process is described in detail in the header comments in "ghost_zone.hh".
//
void patch_system::synchronize(int ghosted_min_gfn_to_sync,
			       int ghosted_max_gfn_to_sync)
{
//
// Phase 1:
// Fill in gridfn data at all the non-corner points of symmetry ghost
// zones, using the symmetries to get this data from its "home patch"
// nominal grids.
//
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	patch& p = ith_patch(pn);
		// n.b. these loops must use _int_ variables for the loop
		//      to terminate!
		for (int want_min = false ; want_min <= true ; ++want_min)
		{
		for (int want_rho = false ; want_rho <= true ; ++want_rho)
		{
		ghost_zone& gz = p.minmax_ang_ghost_zone(want_min, want_rho);
		if (gz.is_symmetry())
		   then gz.synchronize(ghosted_min_gfn_to_sync,
				       ghosted_max_gfn_to_sync,
				       false,	// want min-par corner?
				       true,	// want non-corner?
				       false);	// want max-par corner?
		}
		}
	}

//
// Phase 2:
// Fill in gridfn data in all the interpatch ghost zones, using interpatch
// interpolation from neighboring patches as described above.
//
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	patch& p = ith_patch(pn);
		// n.b. these loops must use _int_ variables for the loop
		//      to terminate!
		for (int want_min = false ; want_min <= true ; ++want_min)
		{
		for (int want_rho = false ; want_rho <= true ; ++want_rho)
		{
		ghost_zone& gz = p.minmax_ang_ghost_zone(want_min, want_rho);
		if (gz.is_interpatch())
		   then gz.synchronize(ghosted_min_gfn_to_sync,
				       ghosted_max_gfn_to_sync);
		}
		}
	}

//
// Phase 3:
// Fill in gridfn data at all the corner points of symmetry ghost zones,
// using the symmetries to get this data from its "home patch" nominal
// grids or ghost zones.
//
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	patch& p = ith_patch(pn);
		// n.b. these loops must use _int_ variables for the loop
		//      to terminate!
		for (int want_min = false ; want_min <= true ; ++want_min)
		{
		for (int want_rho = false ; want_rho <= true ; ++want_rho)
		{
		ghost_zone& gz = p.minmax_ang_ghost_zone(want_min, want_rho);
		if (gz.is_symmetry())
		   then gz.synchronize(ghosted_min_gfn_to_sync,
				       ghosted_max_gfn_to_sync,
				       true,	// want min-par corner?
				       false,	// want non-corner?
				       true);	// want max-par corner?
		}
		}
	}
}

//******************************************************************************

//
// This function computes the Jacobian of  synchronize()  into internal
// buffers, taking into account synchronize()'s full 3-phase algorithm.
//
// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just
//		pass all the Jacobian operations down to the ghost zones
//
void patch_system::compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync,
						int ghosted_max_gfn_to_sync)
	const
{
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	patch& p = ith_patch(pn);
		// n.b. these loops must use _int_ variables for the loop
		//      to terminate!
		for (int want_min = false ; want_min <= true ; ++want_min)
		{
		for (int want_rho = false ; want_rho <= true ; ++want_rho)
		{
		ghost_zone& gz = p.minmax_ang_ghost_zone(want_min, want_rho);
		gz.compute_Jacobian(ghosted_min_gfn_to_sync,
				    ghosted_max_gfn_to_sync);
		}
		}
	}
}

//******************************************************************************

//
// Given that  compute_synchronize_Jacobian()  has been called, this
// function computes the  posn  value of the y points in a Jacobian row.
//
// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just
//		pass all the operation down to the ghost zones
//
int patch_system::synchronize_Jacobian_y_ipar_posn(const ghost_zone& xgz,
						   int x_iperp, int x_ipar)
	const
{
return xgz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
}

//******************************************************************************

//
// Given that  compute_synchronize_Jacobian()  has been called, this
// function computes the Jacobian
//	partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar)
//	-------------------------------------------------------------
//	   partial gridfn(ghosted_gfn, py, y_iperp, y_posn+y_ipar_m)
// where
//	y_iperp = synchronize_Jacobian_y_iperp(xgz, x_iperp)
//	y_posn = synchronize_Jacobian_y_ipar_posn(xgz, x_iperp, x_ipar)
// taking into account synchronize()'s full 3-phase algorithm
//
// FIXME FIXME: at the moment we ignore the 3-phase algorithm and just
//		pass all the operation down to the ghost zones
//
fp patch_system::synchronize_Jacobian
	(const ghost_zone& xgz, int x_iperp, int x_ipar, 
	 int y_ipar_m)
	const
{
return xgz.Jacobian(x_iperp, x_ipar, y_ipar_m);
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function prints an unknown-grid gridfn in ASCII format to a
// named output file.  The output format is suitable for a gnuplot
// 'splot' command.  (Individual patches may be selected with the
//  select.patch  program (perl script).)  The output format is either
//	# print_xyz_flag == false
//	dpx	dpy	gridfn
// or
//	# print_xyz_flag == true
//	dpx	dpy	gridfn   global_x   global_y   global_z
// where global_[xyz} are derived from the angular position and a
// specified (unknown-grid) radius gridfn.
//
void patch_system::print_unknown_gridfn
	(bool ghosted_flag, int unknown_gfn,
	 bool print_xyz_flag, bool radius_is_ghosted_flag,
			      int unknown_radius_gfn,
	 const char output_file_name[], bool want_ghost_zones)
	const
{
if (want_ghost_zones && !ghosted_flag)
   then error_exit(PANIC_EXIT,
"***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n"
"        can't have want_ghost_zones && !ghosted_flag !\n"
,
		   unknown_gfn);				/*NOTREACHED*/
if (want_ghost_zones && print_xyz_flag && !radius_is_ghosted_flag)
   then error_exit(PANIC_EXIT,
"***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n"
"        can't have want_ghost_zones && print_xyz_flag\n"
"        && !radius_is_ghosted_flag!\n"
"        unknown_radius_gfn=%d\n"
,
		   unknown_gfn,
		   unknown_radius_gfn);				/*NOTREACHED*/

FILE *output_fp = std::fopen(output_file_name, "w");
if (output_fp == NULL)
   then error_exit(ERROR_EXIT,
"***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n"
"        can't open output file \"%s\"\n!"
,
		   unknown_gfn,
		   output_file_name);				/*NOTREACHED*/

	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	const patch& p = ith_patch(pn);

	fprintf(output_fp, "### %s patch\n", p.name());
	fprintf(output_fp, "# %s_gfn=%d\n",
		(ghosted_flag ? "ghosted" : "nominal"), unknown_gfn);
	fprintf(output_fp, "# dpx = %s\n", p.name_of_dpx());
	fprintf(output_fp, "# dpy = %s\n", p.name_of_dpy());
	fprintf(output_fp, "#\n");
	fprintf(output_fp,
		print_xyz_flag
		? "# dpx\tdpy\tgridfn\tglobal_x\tglobal_y\tglobal_z\n"
		: "# dpx\tdpy\tgridfn\n");

		for (int irho = p.effective_min_irho(want_ghost_zones) ;
		     irho <= p.effective_max_irho(want_ghost_zones) ;
		     ++irho)
		{
		for (int isigma = p.effective_min_isigma(want_ghost_zones) ;
		     isigma <= p.effective_max_isigma(want_ghost_zones) ;
		     ++isigma)
		{
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);
		const fp dpx = p.dpx_of_rho_sigma(rho, sigma);
		const fp dpy = p.dpy_of_rho_sigma(rho, sigma);
		fprintf(output_fp,
			"%g\t%g\t%#.15g",
			dpx, dpy, p.unknown_gridfn(ghosted_flag,
						   unknown_gfn, irho,isigma));
		if (print_xyz_flag)
		   then {
			const fp r = p.unknown_gridfn(radius_is_ghosted_flag,
						      unknown_radius_gfn,
						      irho,isigma);
			fp local_x, local_y, local_z;
			p.xyz_of_r_rho_sigma(r,rho,sigma,
					     local_x,local_y,local_z);
			const fp global_x = global_x_of_local_x(local_x);
			const fp global_y = global_y_of_local_y(local_y);
			const fp global_z = global_z_of_local_z(local_z);
			fprintf(output_fp,
				"\t%#g\t%#g\t%#g",
				global_x, global_y, global_z);
			}
		fprintf(output_fp, "\n");
		}
		fprintf(output_fp, "\n");
		}
	fprintf(output_fp, "\n");
	}

std::fclose(output_fp);
}

//******************************************************************************

//
// This function reads an unknown-grid gridfn in ASCII format from
// a named input file.  Comments ('#' in column 1) and blank lines
// are ignored, otherwise the input format matches that written by
// print_unknown_gridfn(): the first 3 numbers on each line are taken
// to be dpx, dpy, and the gridfn value; anything else on the line is
// ignored.
//
void patch_system::read_unknown_gridfn(bool ghosted_flag, int unknown_gfn,
			       const char input_file_name[],
			       bool want_ghost_zones)
{
if (want_ghost_zones && !ghosted_flag)
   then error_exit(PANIC_EXIT,
"***** patch_system::read_unknown_gridfn(unknown_gfn=%d):\n"
"        can't have want_ghost_zones && !ghosted_flag !\n"
,
		   unknown_gfn);				/*NOTREACHED*/

FILE *input_fp = std::fopen(input_file_name, "r");
if (input_fp == NULL)
   then error_exit(ERROR_EXIT,
"***** patch_system::read_unknown_gridfn(unknown_gfn=%d):\n"
"        can't open input file \"%s\"\n!"
,
		   unknown_gfn,
		   input_file_name);				/*NOTREACHED*/

int line_number = 1;
	for (int pn = 0 ; pn < N_patches() ; ++pn)
	{
	patch& p = ith_patch(pn);

	for (int irho = p.effective_min_irho(want_ghost_zones) ;
	     irho <= p.effective_max_irho(want_ghost_zones) ;
	     ++irho)
	{
	for (int isigma = p.effective_min_isigma(want_ghost_zones) ;
	     isigma <= p.effective_max_isigma(want_ghost_zones) ;
	     ++isigma)
	{
	const fp rho   = p.rho_of_irho(irho);
	const fp sigma = p.sigma_of_isigma(isigma);
	const fp dpx = p.dpx_of_rho_sigma(rho, sigma);
	const fp dpy = p.dpy_of_rho_sigma(rho, sigma);

	const int N_buffer = 100;
	char buffer[N_buffer];
	// read/discard comments and blank lines
		do
		{
		if (std::fgets(buffer, N_buffer, input_fp) == NULL)
		   then error_exit(ERROR_EXIT,
"***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n"
"        I/O error or unexpected end-of-file on input!\n"
"        at irho=%d of [%d,%d], isigma=%d of [%d,%d]\n"
"        dpx=%g dpy=%g\n"
,
				   p.name(), unknown_gfn,
				   irho, p.effective_min_irho(want_ghost_zones),
					 p.effective_max_irho(want_ghost_zones),
				   isigma,
				      p.effective_min_isigma(want_ghost_zones),
				      p.effective_max_isigma(want_ghost_zones),
				   dpx, dpy);			/*NOTREACHED*/
		++line_number;
		} while ((buffer[0] == '#') || (buffer[0] == '\n'));

	double read_dpx, read_dpy, read_gridfn_value;
	if (sscanf(buffer, "%lf %lf %lf",
		   &read_dpx, &read_dpy, &read_gridfn_value) != 3)
	   then error_exit(ERROR_EXIT,
"***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n"
"        bad input data at input line %d!\n"
,
			   p.name(), unknown_gfn,
			   line_number);			/*NOTREACHED*/
	if (! (    jtutil::fuzzy<fp>::EQ(read_dpx,dpx)
		&& jtutil::fuzzy<fp>::EQ(read_dpy,dpy)    ) )
	   then error_exit(ERROR_EXIT,
"***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n"
"        wrong (dpx,dpy) at input line %d!\n"
"        expected (%g,%g)\n"
"        read     (%g,%g)\n"
,
			   p.name(), unknown_gfn,
			   line_number,
			   dpx, dpy,
			   read_dpx, read_dpy);			/*NOTREACHED*/

	p.unknown_gridfn(ghosted_flag,
			 unknown_gfn, irho,isigma) = read_gridfn_value;
	}
	}

	}

std::fclose(input_fp);
}