// horizon_function.cc -- evaluate LHS function H(h)
// $Id$
//
// <<<prototypes for functions local to this file>>>
// horizon_function - top-level driver
/// setup_xyz_posns - setup global xyz posns of grid points
/// interpolate_geometry - interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid
/// compute_H - compute H(h) given earlier setup
//

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

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.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 "../util/coords.hh"
#include "../util/grid.hh"
#include "../util/fd_grid.hh"
#include "../util/patch.hh"
#include "../util/patch_edge.hh"
#include "../util/patch_interp.hh"
#include "../util/ghost_zone.hh"
#include "../util/patch_system.hh"

#include "../elliptic/Jacobian.hh"

#include "gfn.hh"
#include "AHFinderDirect.hh"

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

//
// ***** prototypes for functions local to this file *****
//

namespace {
void setup_xyz_posns(patch_system& ps);
void interpolate_geometry(patch_system& ps,
			  const struct cactus_grid_info& cgi,
			  const struct geometry_interpolator_info& ii);
void compute_H(patch_system& ps, bool Jacobian_flag);
	  }

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

//
// This function computes the LHS function H(h), and optionally also
// its Jacobian coefficients (from which the Jacobian matrix may be
// computed later).
//
// Inputs (angular gridfns, on ghosted grid):
// ... defined on ghosted grid
// ... only values on nominal grid are actually used as input
//	h				# shape of trial surface
//
// Inputs (Cactus 3-D gridfns):
//	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
//	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
//
// Outputs (temporaries computed at each grid point)
//	## computed by hand-written code
//	xx,yy,zz			# xyz positions of grid points
//	X_ud_*, X_udd_*			# xyz derivative coefficients
//	## computed by Maple-generated code
//	g_uu_{11,12,13,22,23,33}	# $g^{ij}$
//	K				# $K$
//	K_dd_{11,12,13,22,23,33}	# $K^{ij}$
//	partial_d_ln_sqrt_g_d		# $\partial_i \ln \sqrt{g}$
//	partial_d_g_uu_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g^{ij}$
//
// Outputs (angular gridfns, all on the nominal grid):
//	## interpolated from 3-D Cactus grid
//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
//	H				# $H = H(h)$
//
// Arguments:
// Jacobian_flag = true to compute the Jacobian coefficients,
//		   false to skip this.
//
void horizon_function(patch_system& ps,
		      const struct cactus_grid_info& cgi,
		      const struct geometry_interpolator_info& ii,
		      bool Jacobian_flag)
{
CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");

// fill in values of all ghosted gridfns in ghost zones
ps.synchronize();

// set up xyz positions of grid points
setup_xyz_posns(ps);

// interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$
interpolate_geometry(ps, cgi, ii);

// compute remaining gridfns --> $H$ and optionally Jacobian coefficients
// by algebraic ops and angular finite differencing
compute_H(ps, Jacobian_flag);
}

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

//
// This function sets up the global xyz positions of the grid points
// in the gridfns global_[xyz].  These will be used by interplate_geometry().
//
namespace {
void setup_xyz_posns(patch_system& ps)
{
CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");

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

		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
		{
		for (int isigma = p.min_isigma() ;
		     isigma <= p.max_isigma() ;
		     ++isigma)
		{
		const fp r = p.ghosted_gridfn(ghosted_gfns::gfn__h,
					      irho,isigma);
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(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 = ps.global_x_of_local_x(local_x);
		const fp global_y = ps.global_y_of_local_y(local_y);
		const fp global_z = ps.global_z_of_local_z(local_z);

		p.gridfn(nominal_gfns::gfn__global_x, irho,isigma) = global_x;
		p.gridfn(nominal_gfns::gfn__global_y, irho,isigma) = global_y;
		p.gridfn(nominal_gfns::gfn__global_z, irho,isigma) = global_z;
		}
		}
	}
}
	  }

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

//
// This function interpolates the Cactus $g_{ij}$ and $K_{ij}$, and
// the spatial derivatives $\partial_k g_{ij}$, to the trial horizon
// surface position given by h.
//
// Inputs (angular gridfns, on ghosted grid):
// ... defined on ghosted grid
// ... only values on nominal grid are actually used as input
//	h				# shape of trial surface
//
// Inputs (angular gridfns, all on the nominal grid):
//	xx,yy,zz			# xyz positions of grid points
//
// Inputs (Cactus 3-D gridfns):
//	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
//	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
//
// Inputs (angular gridfns, all on the nominal grid):
//	## computed directly from h
//	xx,yy,zz			# xyz positions
//
// Outputs (angular gridfns, all on the nominal grid):
//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
//
// On the first call this function also modifies the interpolator
// parameter table.
//
namespace {
void interpolate_geometry(patch_system& ps,
			  const struct cactus_grid_info& cgi,
			  const struct geometry_interpolator_info& gii)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "      interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid");

int status;

const int N_interp_points = ps.N_grid_points();
const int interp_coords_type_code = CCTK_VARIABLE_REAL;
const void* const interp_coords[N_GRID_DIMS]
  = {
    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_x)),
    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_y)),
    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_z)),
    };

const int N_INPUT_ARRAYS = 12;
const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
	= {
	  // $g_{ij}$
	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
						  CCTK_VARIABLE_REAL,
	    // $K_{ij}$
	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
						  CCTK_VARIABLE_REAL,
	  };
const void* const input_arrays[N_INPUT_ARRAYS]
	= {
	  static_cast<const void*>(cgi.g_dd_11_data),
	  static_cast<const void*>(cgi.g_dd_12_data),
	  static_cast<const void*>(cgi.g_dd_13_data),
	  static_cast<const void*>(cgi.g_dd_22_data),
	  static_cast<const void*>(cgi.g_dd_23_data),
	  static_cast<const void*>(cgi.g_dd_33_data),
	  static_cast<const void*>(cgi.K_dd_11_data),
	  static_cast<const void*>(cgi.K_dd_12_data),
	  static_cast<const void*>(cgi.K_dd_13_data),
	  static_cast<const void*>(cgi.K_dd_22_data),
	  static_cast<const void*>(cgi.K_dd_23_data),
	  static_cast<const void*>(cgi.K_dd_33_data),
	  };

const int N_OUTPUT_ARRAYS = 30;

const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
	= {
 // $g_{ij}$         $\partial_x g_{ij}$ $\partial_y g_{ij}$ $\partial_z g_{ij}$
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 // $K_{ij}$
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
		     CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
					 CCTK_VARIABLE_REAL,
	  };
const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
	= {
	  0, 0, 0, 0,		// g_dd_11
	  1, 1, 1, 1,		// g_dd_12
	  2, 2, 2, 2,		// g_dd_13
	  3, 3, 3, 3,		// g_dd_22
	  4, 4, 4, 4,		// g_dd_23
	  5, 5, 5, 5,		// g_dd_33
	   6,  7,  8,		// K_dd_{11,12,13}
	       9, 10,		// K_dd_{22,23}
		  11,		// K_dd_33
	  };
#define DERIV(x)	x
const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
	= {
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_11
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_12
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_13
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_22
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_23
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_33
	  DERIV(0), DERIV(0), DERIV(0),			// K_dd_{11,12,13}
		    DERIV(0), DERIV(0),			// K_dd_{22,23}
			      DERIV(0),			// K_dd_{33}
	  };
void* const output_arrays[N_OUTPUT_ARRAYS]
  = {
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_11)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_111)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_211)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_311)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_12)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_112)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_212)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_312)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_13)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_113)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_213)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_313)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_22)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_122)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_222)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_322)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_23)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_123)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_223)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_323)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__g_dd_33)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_133)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_233)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__partial_d_g_dd_333)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_11)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_12)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_13)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_22)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_23)),
    static_cast<void*>(ps.gridfn_data(nominal_gfns::gfn__K_dd_33)),
    };

bool first_time = true;
if (first_time)
   then {
	first_time = false;

	// store derivative info in interpolator parameter table
	CCTK_VInfo(CCTK_THORNSTRING,
		   "         setting up derivative info for interpolator");

	status = Util_TableSetIntArray(gii.param_table_handle,
				       N_OUTPUT_ARRAYS, operand_indices,
				       "operand_indices");
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        unable to set operand_indices in interpolator parameter table!\n"
"        Util_TableSetIntArray() status=%d\n"
,
			   status);				/*NOTREACHED*/

	status = Util_TableSetIntArray(gii.param_table_handle,
				       N_OUTPUT_ARRAYS, operation_codes,
				       "operation_codes");
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        unable to set operation_codes in interpolator parameter table!\n"
"        Util_TableSetIntArray() status=%d\n"
,
			   status);				/*NOTREACHED*/
	}

CCTK_VInfo(CCTK_THORNSTRING,
	   "         calling interpolator (%d points)",
	   N_interp_points);
status = CCTK_InterpLocalUniform(N_GRID_DIMS,
				 gii.operator_handle, gii.param_table_handle,
				 cgi.coord_origin, cgi.coord_delta,
				 N_interp_points,
				    interp_coords_type_code,
				    interp_coords,
				 N_INPUT_ARRAYS,
				    cgi.gridfn_dims,
				    input_array_type_codes,
				    input_arrays,
				 N_OUTPUT_ARRAYS,
				    output_array_type_codes,
				    output_arrays);
if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
   then {
	// look in interpolator output table entries
	// to see *which* point is out-of-range
	CCTK_INT out_of_range_pt, out_of_range_axis, out_of_range_end;
	if (    (Util_TableGetInt(gii.param_table_handle,
				  &out_of_range_pt,
				  "out_of_range_pt") < 0)
	     || (Util_TableGetInt(gii.param_table_handle,
				  &out_of_range_axis,
				  "out_of_range_axis") < 0)
	     || (Util_TableGetInt(gii.param_table_handle,
				  &out_of_range_end,
				  "out_of_range_end") < 0)    )
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        point out of range when interpolating geometry info from 3-D grid!\n"
"        ==> the trial horizon surface is (at least partially)\n"
"            outside the grid and/or in an excised region!\n"
"        (unable to get info about which point is out of range:\n"
"         maybe an interpolator problem?)\n");			/*NOTREACHED*/

	assert(out_of_range_pt >= 0);
	assert(out_of_range_pt < ps.N_grid_points());
	const double global_x = ps.gridfn_data(nominal_gfns::gfn__global_x)
				   [out_of_range_pt];
	const double global_y = ps.gridfn_data(nominal_gfns::gfn__global_y)
				   [out_of_range_pt];
	const double global_z = ps.gridfn_data(nominal_gfns::gfn__global_z)
				   [out_of_range_pt];

	assert(out_of_range_axis >= 0);
	assert(out_of_range_axis < N_GRID_DIMS);
	const char axis = "xyz"[out_of_range_axis];

	assert((out_of_range_end == -1) || (out_of_range_end == +1));
	const char end = (out_of_range_end == -1) ? '-' : '+';

	error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        the point (%g,%g,%g) on the trial horizon surface\n"
"        is outside the grid in the %c%c direction!\n"
,
		   global_x, global_y, global_z,
		   end, axis);					/*NOTREACHED*/
	}
if (status < 0)
   then error_exit(ERROR_EXIT,
"***** interpolate_geometry(): error return from interpolator!\n"
"        CCTK_InterpLocalUniform() status=%d\n"
,
		   status);					/*NOTREACHED*/
}
	  }

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

//
// This function computes H(h), and optionally its Jacobian coefficients,
// (from which the Jacobian matrix may be computed later).  This function
// uses a mixture of algebraic operations and (rho,sigma) finite differencing.
// The computation is done (entirely) on the nominal angular grid.
//
// Arguments:
// Jacobian_flag = true to compute the Jacobian coefficients,
//		   false to skip this.
//
namespace {
void compute_H(patch_system& ps, bool Jacobian_flag)
{
CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");

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

		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
		{
		for (int isigma = p.min_isigma() ;
		     isigma <= p.max_isigma() ;
		     ++isigma)
		{
		//
		// compute the X_ud and X_udd derivative coefficients
		// ... n.b. this uses the *local* (x,y,z) coordinates
		//
		const fp r = p.ghosted_gridfn(ghosted_gfns::gfn__h,
					      irho,isigma);
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);
		fp xx, yy, zz;
		p.xyz_of_r_rho_sigma(r,rho,sigma, xx, yy, zz);

		// 1st derivative coefficients X_ud
		const fp X_ud_11 = p.partial_rho_wrt_x(xx, yy, zz);
		const fp X_ud_12 = p.partial_rho_wrt_y(xx, yy, zz);
		const fp X_ud_13 = p.partial_rho_wrt_z(xx, yy, zz);
		const fp X_ud_21 = p.partial_sigma_wrt_x(xx, yy, zz);
		const fp X_ud_22 = p.partial_sigma_wrt_y(xx, yy, zz);
		const fp X_ud_23 = p.partial_sigma_wrt_z(xx, yy, zz);

		// 2nd derivative coefficient gridfns X_udd
		const fp X_udd_111 = p.partial2_rho_wrt_xx(xx, yy, zz);
		const fp X_udd_112 = p.partial2_rho_wrt_xy(xx, yy, zz);
		const fp X_udd_113 = p.partial2_rho_wrt_xz(xx, yy, zz);
		const fp X_udd_122 = p.partial2_rho_wrt_yy(xx, yy, zz);
		const fp X_udd_123 = p.partial2_rho_wrt_yz(xx, yy, zz);
		const fp X_udd_133 = p.partial2_rho_wrt_zz(xx, yy, zz);
		const fp X_udd_211 = p.partial2_sigma_wrt_xx(xx, yy, zz);
		const fp X_udd_212 = p.partial2_sigma_wrt_xy(xx, yy, zz);
		const fp X_udd_213 = p.partial2_sigma_wrt_xz(xx, yy, zz);
		const fp X_udd_222 = p.partial2_sigma_wrt_yy(xx, yy, zz);
		const fp X_udd_223 = p.partial2_sigma_wrt_yz(xx, yy, zz);
		const fp X_udd_233 = p.partial2_sigma_wrt_zz(xx, yy, zz);

		//
		// "call" the Maple-generated code
		// ... each cg/*.c file has a separate set of temp variables,
		//     and so must be inside its own set of { } braces
		//

		// gridfn #defins
		#include "cg.hh"

		  {
		// g_uu
		#include "../gr.cg/inverse_metric.c"
		  }

		  {
		// K, K_uu
		#include "../gr.cg/extrinsic_curvature_trace_raise.c"
		  }

		  {
		// partial_d_g_uu
		#include "../gr.cg/inverse_metric_gradient.c"
		  }

		  {
		// partial_d_ln_sqrt_g
		#include "../gr.cg/metric_det_gradient.c"
		  }

		  {
		// HA, HB, HC, HD, H
		#include "../gr.cg/horizon_function.c"
		  }

		if (Jacobian_flag)
		   then {
			// partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
			#include "../gr.cg/horizon_Jacobian.c"
			}
		}
		}
	}
}
	  }