// patch_interp.cc -- interpolate gridfn data for another patch's ghost zone
// $Id$

//
// patch_interp::patch_interp
// patch_interp::~patch_interp
// patch_interp::[min,max]_ipar_for_gridfn_data
// patch_interp::interpolate
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include "util_Table.h"
#include "cctk.h"

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

#include "fp.hh"
#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"

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

//
// This function constructs an  patch_interp  object.
//
// Note this requires that the adjacent-side  ghost_zone  objects already
// exist, since they're used in determining the ipar range over which the
// interpolation will be done.
//
// It also requires that the patch's gridfns exist, since we size various
// arrays based on the patch's min/max ghosted gfn.
//
patch_interp::patch_interp(const patch_edge& my_edge_in,
			       int min_iperp_in, int max_iperp_in,
			       const jtutil::array1d<int>& min_parindex_used_in,
			       const jtutil::array1d<int>& max_parindex_used_in,
			       const jtutil::array2d<fp>& interp_par_in,
			       int interp_handle_in,
			       int interp_par_table_handle_in)
	: my_patch_(my_edge_in.my_patch()),
	  my_edge_(my_edge_in),
	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
	  min_ipar_(min_ipar_for_gridfn_data()),
	  max_ipar_(max_ipar_for_gridfn_data()),
	  min_parindex_used_(min_parindex_used_in),
	  max_parindex_used_(max_parindex_used_in),
	  interp_par_(interp_par_in),
	  interp_handle_(interp_handle_in),
	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
	  gridfn_coord_origin_(my_edge().par_of_ipar(min_ipar())),
	  gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
	  gridfn_type_codes_ (my_patch().ghosted_min_gfn(),
			      my_patch().ghosted_max_gfn()),
	  gridfn_data_ptrs_  (my_patch().ghosted_min_gfn(),
			      my_patch().ghosted_max_gfn()),
	  result_buffer_ptrs_(my_patch().ghosted_min_gfn(),
			      my_patch().ghosted_max_gfn())
{
int status;

// did interpolator-parameter-table cloning work ok?
status = interp_par_table_handle_;
if (status < 0)
   then error_exit(ERROR_EXIT,
"***** patch_interp::patch_interp():\n"
"        can't clone interpolator parmameter table (handle=%d)!\n"
"        error status=%d\n"
,
		   interp_par_table_handle_in,
		   status);					/*NOTREACHED*/

// set up gridfn storage addressing in interpolator parameter table
const int N_dims = 1;
const CCTK_INT stride = my_patch().ghosted_iang_stride(my_edge().is_rho());
status = Util_TableSetIntArray(interp_par_table_handle_,
			       N_dims, &stride,
			       "input_array_strides");
if (status < 0)
   then error_exit(ERROR_EXIT,
"***** patch_interp::patch_interp():\n"
"        can't set stride in interpolator parmameter table (handle=%d)!\n"
"        error status=%d\n"
,
		   interp_par_table_handle_,
		   status);					/*NOTREACHED*/

// set up gridfn type codes
	for (int gfn = my_patch().ghosted_min_gfn() ;
	     gfn <=    my_patch().ghosted_max_gfn() ;
	     ++gfn)
	{
	gridfn_type_codes_(gfn) = CCTK_VARIABLE_REAL;	// fp == CCTK_REAL
	}
}

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

//
// This function destroys an  patch_interp  object
//
patch_interp::~patch_interp()
{
Util_TableDestroy(interp_par_table_handle_);
}

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

//
// These functions compute the [min,max] range of ipar in this patch
// from which we may use gridfn data in interpolation.  This computation
// depends on our adjacent ghost zones' types as described in the comments
// to our constructor, but needs only my_patch_ and my_edge_ members of
// this object to be valid.
//
int patch_interp::min_ipar_for_gridfn_data() const
{
const ghost_zone& min_par_adjacent_ghost_zone
	= my_patch()
	  .ghost_zone_on_edge(my_edge().min_par_adjacent_edge());
return min_par_adjacent_ghost_zone.is_symmetry()
       ? my_edge().min_ipar_with_corners()
       : my_edge().min_ipar_without_corners();
}

int patch_interp::max_ipar_for_gridfn_data() const
{
const ghost_zone& max_par_adjacent_ghost_zone
	= my_patch()
	  .ghost_zone_on_edge(my_edge().max_par_adjacent_edge());
return max_par_adjacent_ghost_zone.is_symmetry()
       ? my_edge().max_ipar_with_corners()
       : my_edge().max_ipar_without_corners();
}

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

//
// This function interpolates the specified range of ghosted gridfns
// at all the coordinates specified when we were constructed.
// It stores the interpolated data in  result_buffer(gfn, iperp,parindex) .
// It calls  error_exit() if the point is out of range or any other
// interpolator error is detected.
//
void patch_interp::interpolate(int ghosted_min_gfn_to_interp,
				 int ghosted_max_gfn_to_interp,
				 jtutil::array3d<fp>& result_buffer_array)
	const
{
//
// We do a separate interpolation for each iperp.
//
// The interpolator accesses the gridfn data with the generic
// subscripting expression
//	data[offset + i*stride]
// where i is a 0-origin point index for each interpolation, and where
// we get to choose the pointer  data  and the  offset  value for each
// (gridfn,iperp), and the  stride  for each iperp.  Our strategy is to
// leave  offset  unspecified so it defaults to 0, use a common stride
// for all iperp (already set up in our constructor), and specify the
// appropriate data pointers (--> the start of the gridfn data we should
// use) for each (gfn,iperp).  With this strategy we don't need to modify
// the parameter table at all after the constructor.
//
const int N_dims = 1;
const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp,
						ghosted_max_gfn_to_interp);
const CCTK_INT N_gridfn_data_points
	= jtutil::how_many_in_range(min_ipar(), max_ipar());
const int interp_coords_type_code = CCTK_VARIABLE_REAL;	// fp == CCTK_REAL
const CCTK_INT *gridfn_type_codes_ptr
	= & gridfn_type_codes_(ghosted_min_gfn_to_interp);

//
// do the interpolations at each iperp
//
	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
	{
	const int min_parindex = min_parindex_used_(iperp);
	const int max_parindex = max_parindex_used_(iperp);

	// interpolation-point coordinates
	const CCTK_INT N_interp_points
		= jtutil::how_many_in_range(min_parindex, max_parindex);
	const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex);
	const void* const interp_coords[N_dims]
		= { static_cast<const void *>(interp_coords_ptr) };

	// pointers to gridfn data to interpolate, buffer for results
		for (int gfn = ghosted_min_gfn_to_interp ;
		     gfn <= ghosted_max_gfn_to_interp ;
		     ++gfn)
		{
		const int start_irho
			= my_edge().  irho_of_iperp_ipar(iperp, min_ipar());
		const int start_isigma
			= my_edge().isigma_of_iperp_ipar(iperp, min_ipar());
		gridfn_data_ptrs_(gfn)
			= static_cast<const void *>(
				& my_patch()
				  .ghosted_gridfn(gfn, start_irho,start_isigma)
						   );
		result_buffer_ptrs_(gfn)
			= static_cast<void *>(
				& result_buffer_array(gfn, iperp,min_parindex)
					     );
		}
	const void* const* const gridfn_data
		= & gridfn_data_ptrs_(ghosted_min_gfn_to_interp);
	void* const* const result_buffer
		= & result_buffer_ptrs_(ghosted_min_gfn_to_interp);

	// do the actual interpolation
	const int status
		= CCTK_InterpLocalUniform(N_dims,
					  interp_handle_,
					  interp_par_table_handle_,
					  &gridfn_coord_origin_,
					  &gridfn_coord_delta_,
					  N_interp_points,
						  interp_coords_type_code,
						  interp_coords,
					  N_gridfns,
						  &N_gridfn_data_points,
						  gridfn_type_codes_ptr,
						  gridfn_data,
					  N_gridfns,
						  gridfn_type_codes_ptr,
						  result_buffer);
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** patch_interp::interpolate():\n"
"        error return %d from interpolator at iperp=%d of [%d,%d]!\n"
,
			   status, iperp, min_iperp(), max_iperp());
								/*NOTREACHED*/
	}
}