# Interpolation
#TODO: This will be needed for interface with SymmetryInterpolate
#CCTK_INT FUNCTION SymmetryInterpolate(
#  CCTK_POINTER_TO_CONST IN cctkGH,
#  CCTK_INT IN N_dims,
#  CCTK_INT IN local_interp_handle,
#  CCTK_INT IN param_table_handle,
#  CCTK_INT IN coord_system_handle,
#  CCTK_INT IN N_interp_points,
#  CCTK_INT IN interp_coords_type,
#  CCTK_POINTER_TO_CONST ARRAY IN interp_coords,
#  CCTK_INT IN N_input_arrays,
#  CCTK_INT ARRAY IN input_array_indices,
#  CCTK_INT IN N_output_arrays,
#  CCTK_INT ARRAY IN output_array_types,
#  CCTK_POINTER ARRAY IN output_arrays)
#USES FUNCTION SymmetryInterpolate

CCTK_INT FUNCTION DriverInterpolate(
  CCTK_POINTER_TO_CONST IN cctkGH,
  CCTK_INT IN N_dims,
  CCTK_INT IN local_interp_handle,
  CCTK_INT IN param_table_handle,
  CCTK_INT IN coord_system_handle,
  CCTK_INT IN N_interp_points,
  CCTK_INT IN interp_coords_type,
  CCTK_POINTER_TO_CONST ARRAY IN interp_coords,
  CCTK_INT IN N_input_arrays,
  CCTK_INT ARRAY IN input_array_indices,
  CCTK_INT IN N_output_arrays,
  CCTK_INT ARRAY IN output_array_types,
  CCTK_POINTER ARRAY IN output_arrays)
PROVIDES FUNCTION DriverInterpolate WITH CarpetX_DriverInterpolate LANGUAGE C

#ifndef INTERP_H
#define INTERP_H
#include "cctk.h"
/* Scheduled functions */
int CarpetX_InterpStartup(void);
int CarpetX_InterpGridArrays(
    cGH const *const cGH, int const N_dims, int const local_interp_handle,
    int const param_table_handle, int const coord_system_handle,
    int const N_interp_points, int const interp_coords_type_code,
    void const *const interp_coords[], int const N_input_arrays,
    CCTK_INT const input_array_variable_indices[], int const N_output_arrays,
    CCTK_INT const output_array_type_codes[], void *const output_arrays[]);
extern "C" CCTK_INT CarpetX_DriverInterpolate(
    CCTK_POINTER_TO_CONST const cctkGH, CCTK_INT const N_dims,
    CCTK_INT const local_interp_handle, CCTK_INT const param_table_handle,
    CCTK_INT const coord_system_handle, CCTK_INT const N_interp_points,
    CCTK_INT const interp_coords_type_code,
    CCTK_POINTER_TO_CONST const coords_list[], CCTK_INT const N_input_arrays,
    CCTK_INT const input_array_variable_indices[],
    CCTK_INT const N_output_arrays, CCTK_INT const output_array_type_codes[],
    CCTK_POINTER const output_arrays[]);
#endif /* !defined(INTERP_H) */

#include "interp.hxx"
#include "util_Table.h"
#include "util_ErrorCodes.h"

extern "C" CCTK_INT CarpetX_InterpGridArrays(
    cGH const *const cctkGH, int const N_dims, int const local_interp_handle,
    int const param_table_handle, int const coord_system_handle,
    int const N_interp_points, int const interp_coords_type_code,
    void const *const coords[], int const N_input_arrays,
    CCTK_INT const input_array_variable_indices[], int const N_output_arrays,
    CCTK_INT const output_array_type_codes[], void *const output_arrays[]) {
/* TODO: verify that the interface with SymmetryInterpolate can be simply copied from
   Carpet like below */
//  if (CCTK_IsFunctionAliased("SymmetryInterpolate")) {
//    return SymmetryInterpolate(
//        cctkGH, N_dims, local_interp_handle, param_table_handle,
//        coord_system_handle, N_interp_points, interp_coords_type_code, coords,
//        N_input_arrays, input_array_variable_indices, N_output_arrays,
//        output_array_type_codes, output_arrays);
//  } else {
    return CarpetX_DriverInterpolate(
        cctkGH, N_dims, local_interp_handle, param_table_handle,
        coord_system_handle, N_interp_points, interp_coords_type_code, coords,
        N_input_arrays, input_array_variable_indices, N_output_arrays,
        output_array_type_codes, output_arrays);
//  }
}
extern "C" CCTK_INT CarpetX_InterpStartup() {
  CCTK_OverloadInterpGridArrays(CarpetX_InterpGridArrays);
  return 0;
}
extern "C" CCTK_INT CarpetX_DriverInterpolate(
    CCTK_POINTER_TO_CONST const cctkGH, CCTK_INT const N_dims,
    CCTK_INT const local_interp_handle, CCTK_INT const param_table_handle,
    CCTK_INT const coord_system_handle, CCTK_INT const N_interp_points,
    CCTK_INT const interp_coords_type_code,
    CCTK_POINTER_TO_CONST const coords[], CCTK_INT const N_input_arrays,
    CCTK_INT const input_array_variable_indices[],
    CCTK_INT const N_output_arrays, CCTK_INT const output_array_type_codes[],
    CCTK_POINTER const output_arrays[]) {
  DECLARE_CCTK_PARAMETERS;
// This verifies that the order in param_table_handle matches the order of the
// runtime parameter from CarpetX
  CCTK_INT order;
  int n_elems = Util_TableGetInt(param_table_handle, &order, "order");
  assert(n_elems == 1);
  assert(order == interpolation_order);

  vector<CCTK_INT> varinds;
  varinds.resize(N_output_arrays);
  n_elems = Util_TableGetIntArray
    (param_table_handle, N_output_arrays, varinds.data(), "operand_indices");
  assert(n_elems == N_output_arrays);
  vector<CCTK_INT> operations;
  operations.resize(N_output_arrays);
  n_elems = Util_TableGetIntArray
    (param_table_handle, N_output_arrays, operations.data(), "operation_codes");
  assert(n_elems == N_output_arrays);
  const CCTK_POINTER resultptrs = (const CCTK_POINTER)output_arrays;
  CarpetX_Interpolate(cctkGH, N_interp_points, static_cast<const CCTK_REAL*>(coords[0]),
              static_cast<const CCTK_REAL*>(coords[1]), static_cast<const CCTK_REAL*>(coords[2]),
              N_output_arrays, varinds.data(), operations.data(), resultptrs);
  return 0;
}