#include <catalyst.h>
#include <stdio.h>
#include <math.h>
#include "pluto.h"

//-----------------------------------------------------------------------------
/**
 * Initialize Catalyst.
 */
//-----------------------------------------------------------------------------
void do_catalyst_initialization(int scripts, char** pipeline_script)
{
  conduit_node* catalyst_init_params = conduit_node_create();
  // pass scripts pass on command line.
  char buf[256];
  int i;
  for (i=0; i<scripts; i++){
    snprintf(buf, 256, "catalyst/scripts/script%d", i);
    conduit_node_set_path_char8_str(catalyst_init_params, buf, pipeline_script[i]);
  }
  conduit_node_set_path_char8_str(catalyst_init_params, "catalyst_load/implementation", "paraview");
  conduit_node_set_path_char8_str(catalyst_init_params, "catalyst_load/search_paths/paraview", PARAVIEW_IMPL_DIR);
  enum catalyst_status err = catalyst_initialize(catalyst_init_params);
  conduit_node_destroy(catalyst_init_params);
  if (err != catalyst_status_ok)
  {
    printf("Failed to initialize Catalyst: %d\n", err);
  }
}

//-----------------------------------------------------------------------------
/**
 * Execute per cycle
 */
//-----------------------------------------------------------------------------
void do_catalyst_execute(int cycle, double time, Grid* grid, Data* d)
{
  /*   Conversion to Conduit (catalyst) compatible data arrangement   */
  
  unsigned int i, j, k;
  static unsigned int initialized = 0;
  unsigned int counter = 0;
  
  static double *xp, *yp, *zp;
  static uint64_t *CellConn;
  const uint64_t numCells = grid->np_tot[IDIR]*grid->np_tot[JDIR]*grid->np_tot[KDIR];
  const uint64_t numPoints[3] = {grid->np_tot[IDIR]+1, grid->np_tot[JDIR]+1, grid->np_tot[KDIR]+1};
  
  double ***temp  = GetUserVar("temp");
  double ***ndens = GetUserVar("ndens"); 
  double ***mach  = GetUserVar("mach");
  ComputeUserVar (d, grid);  

  if (initialized==0){
    xp = (double *)malloc((grid->np_tot[IDIR]+1) * (grid->np_tot[JDIR]+1) * (grid->np_tot[KDIR]+1) * sizeof(double));
    yp = (double *)malloc((grid->np_tot[IDIR]+1) * (grid->np_tot[JDIR]+1) * (grid->np_tot[KDIR]+1) * sizeof(double));
    zp = (double *)malloc((grid->np_tot[IDIR]+1) * (grid->np_tot[JDIR]+1) * (grid->np_tot[KDIR]+1) * sizeof(double));
    CellConn   = (uint64_t*)malloc(8 * numCells * sizeof(int64_t));
    
    for(k=0; k<grid->np_tot[KDIR]; k++)
    {
      for(j=0; j<grid->np_tot[JDIR]; j++)
      {
        for(i=0; i<grid->np_tot[IDIR]; i++)
        { 
          //cell_count = k * grid->np_tot[IDIR] * grid->np_tot[JDIR] + j * grid->np_tot[IDIR] + i;
          CellConn[counter++] =     k * numPoints[JDIR]*numPoints[IDIR] +     j * numPoints[IDIR] + i; 

          CellConn[counter++] = (k+1) * numPoints[JDIR]*numPoints[IDIR] +     j * numPoints[IDIR] + i;

          CellConn[counter++] = (k+1) * numPoints[JDIR]*numPoints[IDIR] + (j+1) * numPoints[IDIR] + i;

          CellConn[counter++] =     k * numPoints[JDIR]*numPoints[IDIR] + (j+1) * numPoints[IDIR] + i;;

          CellConn[counter++] =     k * numPoints[JDIR]*numPoints[IDIR] +     j * numPoints[IDIR] + i+1;

          CellConn[counter++] = (k+1) * numPoints[JDIR]*numPoints[IDIR] +     j * numPoints[IDIR] + i+1;

          CellConn[counter++] = (k+1) * numPoints[JDIR]*numPoints[IDIR] + (j+1) * numPoints[IDIR] + i+1;

          CellConn[counter++] =     k * numPoints[JDIR]*numPoints[IDIR] + (j+1) * numPoints[IDIR] + i+1;
        }
      }
    }
  }
  
  counter = 0;
  double x1, x2, x3;
  for (k = 0; k<=grid->np_tot[KDIR]; k++){
    for (j = 0; j<=grid->np_tot[JDIR]; j++){
      for (i = 0; i<=grid->np_tot[IDIR]; i++){
        x1 = (i!=grid->np_tot[IDIR])?grid->xl[IDIR][i]:grid->xr[IDIR][i-1];
        x2 = (j!=grid->np_tot[JDIR])?grid->xl[JDIR][j]:grid->xr[JDIR][j-1];
        x3 = (k!=grid->np_tot[KDIR])?grid->xl[KDIR][k]:grid->xr[KDIR][k-1];

        xp[counter]   = x1*sin(x2)*cos(x3);
        yp[counter]   = x1*sin(x2)*sin(x3);
        zp[counter++] = x1*cos(x2);
      }
    }
  }
  counter = 0;
  
  /*   setup conduit for catalyst   */
  conduit_node* catalyst_exec_params = conduit_node_create();
  conduit_node_set_path_int64(catalyst_exec_params, "catalyst/state/timestep", cycle);
  conduit_node_set_path_int64(catalyst_exec_params, "catalyst/state/cycle", cycle);
  conduit_node_set_path_float64(catalyst_exec_params, "catalyst/state/time", time);
  
  conduit_node_set_path_char8_str(catalyst_exec_params, "catalyst/channels/grid/type", "mesh");
  conduit_node* mesh = conduit_node_create();
  
  // add coordsets
  conduit_node_set_path_char8_str(mesh, "coordsets/coords/type", "explicit");
  
  conduit_node_set_path_external_float64_ptr(mesh, "coordsets/coords/values/x", xp, numPoints[IDIR]*numPoints[JDIR]*numPoints[KDIR]);
  conduit_node_set_path_external_float64_ptr(mesh, "coordsets/coords/values/y", yp, numPoints[IDIR]*numPoints[JDIR]*numPoints[KDIR]);
  conduit_node_set_path_external_float64_ptr(mesh, "coordsets/coords/values/z", zp, numPoints[IDIR]*numPoints[JDIR]*numPoints[KDIR]);
  
  // add topologies
  conduit_node_set_path_char8_str(mesh, "topologies/mesh/type", "unstructured");
  conduit_node_set_path_char8_str(mesh, "topologies/mesh/coordset", "coords");
  conduit_node_set_path_char8_str(mesh, "topologies/mesh/elements/shape", "hex");
  conduit_node_set_path_external_int64_ptr(mesh, "topologies/mesh/elements/connectivity", CellConn, numCells * 8);

  // add density (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/density/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/density/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/density/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/density/values", /*Density */(double *)(**d->Vc[RHO]), numCells );
  
  // add pressure (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/pressure/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/pressure/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/pressure/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/pressure/values", /*Pressure */(double *)(**d->Vc[PRS]), numCells );
  
  // add velocityX (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/velocityX/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/velocityX/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/velocityX/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/velocityX/values", /*VX */(double *)(**d->Vc[VX1]), numCells );
  
  // add velocityY (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/velocityY/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/velocityY/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/velocityY/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/velocityY/values", /*VY */(double *)(**d->Vc[VX2]), numCells );
  
  // add velocityZ (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/velocityZ/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/velocityZ/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/velocityZ/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/velocityZ/values", /*VZ */(double *)(**d->Vc[VX3]), numCells );
  
  // add tracer (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/tracer1/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/tracer1/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/tracer1/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/tracer1/values", /*TRC */(double *)(**d->Vc[TRC]), numCells );
  
  // add temperature (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/temperature/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/temperature/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/temperature/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/temperature/values", /*temperature */(double *)(**temp), numCells );
  
  // add ndens (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/ndens/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/ndens/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/ndens/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/ndens/values", /*ndens */(double *)(**ndens), numCells );
  
  // add mach (cell-field)
  conduit_node_set_path_char8_str(mesh, "fields/mach/association", "element");
  conduit_node_set_path_char8_str(mesh, "fields/mach/topology", "mesh");
  conduit_node_set_path_char8_str(mesh, "fields/mach/volume_dependent", "false");
  conduit_node_set_path_external_float64_ptr(mesh, "fields/mach/values", /*mach */(double *)(**mach), numCells );

  // add the mesh info (conduit mesh) to catalyst_exec_params
  conduit_node_set_path_external_node(catalyst_exec_params, "catalyst/channels/grid/data", mesh);
  
  #ifdef CATALYST_DEBUG
  // print for debugging purposes, if needed
  conduit_node_print(catalyst_exec_params);

  // print information with details about memory allocation
  conduit_node* info = conduit_node_create();
  conduit_node_info(catalyst_exec_params, info);
  conduit_node_print(info);
  conduit_node_destroy(info);
  #endif
  
  initialized = 1;
  enum catalyst_status err = catalyst_execute(catalyst_exec_params);
  if (err != catalyst_status_ok)
  {
    printf("Failed to execute Catalyst: %d\n", err);
  }
  conduit_node_destroy(catalyst_exec_params);
  conduit_node_destroy(mesh);
}

//-----------------------------------------------------------------------------
/**
 * Finalize Catalyst.
 */
//-----------------------------------------------------------------------------
void do_catalyst_finalization()
{
  conduit_node* catalyst_fini_params = conduit_node_create();
  enum catalyst_status err = catalyst_finalize(catalyst_fini_params);
  if (err != catalyst_status_ok)
  {
    printf("Failed to execute Catalyst: %d\n", err);
  }
  conduit_node_destroy(catalyst_fini_params);
}