// gr.hh -- header file for general relativity code
// $Header$

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

//
// number of spatial dimensions in the main Cactus grid
// and in our trial-horizon-surface grid
//
enum	{ N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 };

//
// this enum holds the (a) decoded  geometry_method  parameter,
// i.e. it specifies how we compute the (a) geometry
//
enum	geometry_method
	{
	geometry__global_interp_from_Cactus_grid,
	geometry__local_interp_from_Cactus_grid,
	geometry__Schwarzschild_EF // no comma
	};

//
// this enum holds the (a) decoded  Jacobian_compute_method  parameter,
// i.e. it specifies how we compute the (a) Jacobian matrix
//
enum	Jacobian_compute_method
	{
	Jacobian__numerical_perturbation,
	Jacobian__symbolic_diff_with_FD_dr,
	Jacobian__symbolic_diff // no comma
	};

//
// this enum describes the status of an
//	expansion()
// and/or
//	expansion_Jacobian()
// computation which returns (i.e. which doesn't abort the Cactus run)
//
enum	expansion_status
	{
	// successful computation
	// ... this is also returned for dummy computations
	expansion_success,

	// non-finite() (i.e. +/-infinity or NaN) values found in h
	expansion_failure__surface_nonfinite,

	// geometry interpolator returns CCTK_ERROR_INTERP_POINT_OUTSIDE
	expansion_failure__surface_outside_grid,

	// geometry interpolator returns CCTK_ERROR_INTERP_POINT_EXCISED
	expansion_failure__surface_in_excised_region,	// (not implemented yet)

	// non-finite (i.e. +/-infinity or NaN) values found
	// in interpolated geometry (g_ij, K_ij, partial_k g_ij)
	expansion_failure__geometry_nonfinite,

	// interpolated g_ij isn't positive definite
	expansion_failure__gij_not_positive_definite // no comma
	};

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

//
// This structure holds all the information we need about the Cactus grid
// and gridfns in order to interpolate the geometry information.
//
struct	cactus_grid_info
	{
	cGH *GH;			// --> Cactus grid hierarchy

	// this describes the semantics of the Cactus gij gridfns:
	// true ==> the Cactus g_ij are conformal values as per
	//          CactusEinstein/StaticConformal and the  psi  gridfn
	// false ==> the Cactus g_ij are the physical metric
	bool use_Cactus_conformal_metric;

	// Cactus coordinate system
	fp global_coord_origin[N_GRID_DIMS];	// global (x,y,z)
						//   = global origin + ijk*delta
	fp local_coord_origin[N_GRID_DIMS];	// global (x,y,z)
						// of (i,j,k) = (0,0,0)
						// on this processor
	fp coord_delta[N_GRID_DIMS];		// (x,y,z) grid spacing
	fp mean_coord_delta;		// geometric mean of x,y,z grid spacings

	// dimensions of gridfn data on this processor, viewed as a 3-D array
	// n.b. storage ordering is Fortran,
	//	i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
	CCTK_INT local_gridfn_dims[N_GRID_DIMS];

	//
	// coordinate conversion functions
	// here ijk = gridfn array indices on this processor
	//      xyz = floating-point global coordinates
	//

	// convert integer ijk --> floating-point global xyz coordinates
	fp global_xyz_of_ijk(int axis, int ijk) const
		{ return local_coord_origin[axis] + ijk*coord_delta[axis]; }
	// convert floating-point global xyz --> integer ijk coordinates
	// ... but as a floating-point number
	fp fp_ijk_of_global_xyz(int axis, fp xyz) const
		{ return (xyz - local_coord_origin[axis]) / coord_delta[axis]; }
	// convert floating-point global_xyz --> integer ijk coordinates
	// ... rounding down/up to the next integer (= to the next grid point)
	int ijk_floor_of_global_xyz(int axis, fp xyz) const
		{ return jtutil::ifloor(fp_ijk_of_global_xyz(axis, xyz)); }
	int ijk_ceil_of_global_xyz (int axis, fp xyz) const
		{ return jtutil::iceil(fp_ijk_of_global_xyz(axis, xyz)); }

	//
	// stuff for doing a global interpolation via CCTK_InterpGridArrays()
	// i.e. geometry_info.geometry_method
	//	== geometry__global_interp_from_Cactus_grid
	//

	// Cactus coordinate system
	int coord_system_handle;

	// Cactus variable indices
	int g_dd_11_varindex, g_dd_12_varindex, g_dd_13_varindex,
			      g_dd_22_varindex, g_dd_23_varindex,
						g_dd_33_varindex;
	int K_dd_11_varindex, K_dd_12_varindex, K_dd_13_varindex,
			      K_dd_22_varindex, K_dd_23_varindex,
						K_dd_33_varindex;
	int psi_varindex;	// unused if !use_Cactus_conformal_metric

	//
	// stuff for doing a local interpolation via CCTK_InterpLocalUniform()
	// i.e. geometry_info.geometry_method
	//	== geometry__local_interp_from_Cactus_grid
	//

	// --> Cactus gridfn data for grid posn (0,0,0)
	const fp* g_dd_11_dataptr;
	const fp* g_dd_12_dataptr;
	const fp* g_dd_13_dataptr;
	const fp* g_dd_22_dataptr;
	const fp* g_dd_23_dataptr;
	const fp* g_dd_33_dataptr;
	const fp* K_dd_11_dataptr;
	const fp* K_dd_12_dataptr;
	const fp* K_dd_13_dataptr;
	const fp* K_dd_22_dataptr;
	const fp* K_dd_23_dataptr;
	const fp* K_dd_33_dataptr;
	const fp* psi_dataptr;		// NULL if !use_Cactus_conformal_metric
	};

//
// This structure holds information for computing the spacetime geometry.
// This is normally done by interpolating $g_{ij}$ and $K_{ij}$ from the
// usual Cactus grid, but can optionally instead by done by hard-wiring
// the Schwarzschild geometry in Eddington-Finkelstein coordinates.
//
struct	geometry_info
	{
	enum geometry_method geometry_method;

	// parameters for geometry_method = interpolate_from_Cactus_grid
	int operator_handle;		// Cactus handle to interpolation op
	int param_table_handle;		// Cactus handle to parameter table
					// for the interpolator

	// parameters for geometry_method = Schwarzschild_EF
	fp geometry__Schwarzschild_EF__x_posn;	// x posn of Schwarzschild BH
	fp geometry__Schwarzschild_EF__y_posn;	// y posn of Schwarzschild BH
	fp geometry__Schwarzschild_EF__z_posn;	// z posn of Schwarzschild BH
	fp geometry__Schwarzschild_EF__mass;	// mass of Schwarzschild BH
	fp geometry__Schwarzschild_EF__epsilon;	// use z axis limits if
						// (x^2+y^2)/r^2 <= this
	fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD
						// computation of partial_k g_ij

	// should we check various gridfns for NaNs?
	bool check_that_h_is_finite;
	bool check_that_geometry_is_finite;
	};

//
// This struct holds parameters for computing the Jacobian matrix.
//
struct	Jacobian_info
	{
	enum Jacobian_compute_method     Jacobian_compute_method;
	enum Jacobian_store_solve_method Jacobian_store_solve_method;
	fp perturbation_amplitude;
	};

//
// This struct holds information on what to do if various error/warning
// conditions occur.
//
struct	error_info
	{
	// CCTK_VWarn() level for point outside (or too close to any
	// boundary of) the Cactus grid, i.e. geometry interpolator returns
	// CCTK_ERROR_INTERP_POINT_OUTSIDE
	// ... warning level if error occurs on first Newton iteration,
	//     i.e. when evaluating expansion/Jacobian for initial guess
	int warn_level__point_outside__initial;
	// ... warning level if error occurs on a subsequent Newton iteration
	int warn_level__point_outside__subsequent;

	// CCTK_VWarn() level for the user set
	//  check_that_geometry_is_finite = "true"
	// but the Cactus configure process failed to find the finite()
	// function ==> we have to skip the check
	int warn_level__skipping_finite_check;

	// CCTK_VWarn() level for infinity or NaN in interpolated geometry
	// (g_ij, partial_k g_ij, and/or K_ij)
	int warn_level__nonfinite_geometry;

	// CCTK_VWarn() level for interpolated g_ij isn't positive definite
	// (usually this means we're too close to a singularity)
	// ... warning level if error occurs on first Newton iteration,
	//     i.e. when evaluating expansion/Jacobian for initial guess
	int warn_level__gij_not_positive_definite__initial;
	// ... warning level if error occurs on a subsequent Newton iteration
	int warn_level__gij_not_positive_definite__subsequent;
	};

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

//
// prototypes for functions visible outside their source files
//

// expansion.cc
enum expansion_status
  expansion(patch_system* ps_ptr,
	    const struct cactus_grid_info& cgi,
	    const struct geometry_info& gi,
	    const struct error_info& error_info, bool initial_flag,
	    bool Jacobian_flag = false,
	    bool print_msg_flag = false,
	    jtutil::norm<fp>* H_norms_ptr = NULL);

// expansion_Jacobian.cc
enum expansion_status
  expansion_Jacobian(patch_system* ps_ptr,
		     Jacobian* Jac_ptr,
		     const struct cactus_grid_info& cgi,
		     const struct geometry_info& gi,
		     const struct Jacobian_info& Jacobian_info,
		     const struct error_info& error_info, bool initial_flag,
		     bool print_msg_flag = false);

// Schwarzschild_EF.cc
void Schwarzschild_EF_geometry(patch_system& ps,
			       const struct geometry_info& gi,
			       bool msg_flag);

// misc.cc
enum geometry_method
  decode_geometry_method(const char geometry_method_string[]);
#ifdef NOT_USED
bool geometry_method_is_interp(enum geometry_method geometry_method);
#endif
enum Jacobian_compute_method
  decode_Jacobian_compute_method(const char Jacobian_compute_method_string[]);
const char* expansion_status_string(enum expansion_status status);