/*-------------------------------------------------------------------------78--
 * Program to read in 3D finite element data from ViTLES simulations
 * and generate output files in Tecplot format.           
 *
 *     *** Fine Resolution Version ***
 *
 * This version creates a linear element associated with each node (thus
 * splits quadratic tetrahedra into 7 linear tetrahedra and interpolates
 * values of pressure at the mid-side nodes).
 *
 * The Tecplot format includes: 
 *               vitles_root.plt     (coordinate, connectivity and solution)
 ----------------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main ( int argc, char *argv[] )
{
  // variables for file names
  FILE  *in_vitles;  // file used to read ViTLES file (to get problem sizes)
  FILE  *in_vitlesx; // file used to read ViTLES geometry
  FILE  *in_vitlesu; // file used to read ViTLES solution
  FILE  *out_plt;    // file used to write Tecplot file

  char  vitles_root[255];

  // variables for filenames
  char  vitles_file[255], geometry_file[255]   , solution_file[255],
        tec_file[255]   ; 

  int   time_iteration;

  // counters
  int    e_counter, e, i, j, n, nz;

  // mesh dimensions
  int    n_elems      , n_nodes        , n_zones      ;

  // "pointers"
  int    local_node   , n_zone_unknowns;
  int    n_dir        , n_elzone       , n_equations  ;

  // fscan temporary variables
  char   temp[128]; 
  double D[11];
  int    I[11];

  // allocatables
//int    *e_conn;
//int    *u_dirichlet_node, *v_dirichlet_node;
//int    *t_dirichlet_node;
//int    *ide_u,  *ide_p;
//double *dir_u;
//int    *node_flag_u, *node_flag_p;
//int    *local_to_global_node;
//int    *element_to_local_u;
//double *dirichlet_valu;
//double *x;

#define T_NODES 40000
#define T_ELEMS 20000

  int    t_nodes = T_NODES;
  int    t_elems = T_ELEMS;
  int    e_conn[T_ELEMS][10];
  int    u_dirichlet_node[T_NODES],     v_dirichlet_node[T_NODES];
  int    w_dirichlet_node[T_NODES];
  int    ide_u[T_NODES][3],             ide_p[T_NODES];
  double dir_u[T_NODES];
  int    node_flag_u[T_NODES],          node_flag_p[T_NODES];
  int    local_to_global_node[T_NODES];
  int    element_to_local_u[T_ELEMS][10], element_to_local_v[T_ELEMS][10];
  int    element_to_local_w[T_ELEMS][10];
  double dirichlet_valu[T_NODES];
  double x[T_NODES][3];
  double u[T_NODES];
  double v[T_NODES];
  double w[T_NODES];
  double p[T_NODES];

#undef T_NODES
#undef T_ELEMS

  int    e0, e1, e2, e3, e4, e5, e6, e7, e8, e9;

  /*-------------------------------------------------------------------------*/
  /*  Open ViTLES geometry and data files                                    */
  /*-------------------------------------------------------------------------*/
  if (argc !=3)
  {
    printf("\n%s%s%s\n\n",
           "Usage:  ", argv[0], "  ViTLES_rootname  time_snapshot_number");
    exit(1);
  }

  strcpy(vitles_root,argv[1]);
  time_iteration = atoi(argv[2]);

  printf("reading solution from time snapshot: %4d\n",time_iteration);

  // get problem size from the vitles_root.vitles file...
  sprintf(vitles_file,"%s.vitles",vitles_root);
  in_vitles = fopen(vitles_file,"r");
  if ( in_vitles==NULL )
  {
    printf("Can't open file %s to get problem size information",vitles_file);
    exit(1);
  }

  fgets(temp,128,in_vitles);
  fscanf(in_vitles,"%d %d %d \n",&I[0],&I[1],&I[2]);
  n_nodes = I[0];
  n_zones = I[2];

  e_counter = 0;         // element counter

  /*-------------------------------------------------------------------------*/
  /*  Eventually Allocate Storage Space (malloc/calloc)                      */
  /*  instead of hard-wiring t_elements and t_nodes                          */
  /*-------------------------------------------------------------------------*/
/*e_conn = calloc(n_elems*6, sizeof(int));
  dirichlet_node = calloc(n_nodes, sizeof(int));
  ide_u =  calloc(n_nodes*3, sizeof(int));
  ide_p =  calloc(n_nodes  , sizeof(int));
  dir_u =  calloc(n_nodes  , sizeof(double));
  node_flag_u = calloc(n_nodes, sizeof(int));
  node_flag_p = calloc(n_nodes, sizeof(int));
  local_to_global_node = calloc(n_nodes, sizeof(int));
  dirichlet_valu = calloc(n_nodes, sizeof(double));
  x                = calloc(n_nodes*2, sizeof(double)); */

  // get total number of elements
  n_elems = 0;
  /*--------------------------------------------------------------------------*/
  /*  Loop over processors and read geometry and solution files               */
  /*--------------------------------------------------------------------------*/
  for (nz=0; nz<n_zones; nz++)
  {
    sprintf(geometry_file,"%s.%03dgeom",vitles_root,nz);
    in_vitlesx = fopen(geometry_file,"rb");
    if (in_vitlesx==NULL)
    {
      printf("Can't open file %s\n",geometry_file);
      exit(1);
    }

    fread(&local_node      , sizeof(int)   , 1         , in_vitlesx);
    fread(&n_elzone        , sizeof(int)   , 1         , in_vitlesx);
    fread(&n_zone_unknowns , sizeof(int)   , 1         , in_vitlesx);
    fread(&n_dir           , sizeof(int)   , 1         , in_vitlesx);

    n_elems = n_elems + n_elzone;

    //--------------------------------------------------------------------------
    // read and write coordinates
    //--------------------------------------------------------------------------
    for (i=0; i<n_nodes; i++)
    {
      fread(&x[i][0]     , sizeof(double), 1 , in_vitlesx);
      fread(&x[i][1]     , sizeof(double), 1 , in_vitlesx);
      fread(&x[i][2]     , sizeof(double), 1 , in_vitlesx);
    }

    //--------------------------------------------------------------------------
    // read elements
    //--------------------------------------------------------------------------
    for (i=0; i<n_elzone; i++)
    {
      fread(&e_conn[i][0]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][1]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][2]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][3]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][4]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][5]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][6]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][7]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][8]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][9]    , sizeof(int), 1, in_vitlesx);
    }
    fclose(in_vitlesx);

    // open ViTLES solution file
    sprintf(solution_file,"%s.%03d.%04d",vitles_root,nz,time_iteration);
    in_vitlesu = fopen(solution_file,"rb");
    if (in_vitlesu==NULL)
    {
      printf("Can't open file %s\n",solution_file);
      exit(1);
    }

    fread(&D, sizeof(double), 1, in_vitlesu);  // read in time

//  fread(&D, sizeof(double), 3, in_vitlesu);  // read over sizes

    for (e=0; e<n_elzone; e++)
    {
      fread(&D, sizeof(double), 10, in_vitlesu);
      for (i=0; i<10; i++)
      {
        u[e_conn[e][i]] = D[i];
      }

      fread(&D, sizeof(double), 10, in_vitlesu);
      for (i=0; i<10; i++)
      {
        v[e_conn[e][i]] = D[i];
      }

      fread(&D, sizeof(double), 10, in_vitlesu);
      for (i=0; i<10; i++)
      {
        w[e_conn[e][i]] = D[i];
      }

      fread(&D, sizeof(double), 4, in_vitlesu);
      for (i=0; i<4; i++)
      {
        p[e_conn[e][i]] = D[i];
      }
      p[e_conn[e][4]] = (D[0]+D[1])/2;
      p[e_conn[e][5]] = (D[1]+D[2])/2;
      p[e_conn[e][6]] = (D[0]+D[2])/2;
      p[e_conn[e][7]] = (D[1]+D[3])/2;
      p[e_conn[e][8]] = (D[2]+D[3])/2;
      p[e_conn[e][9]] = (D[0]+D[3])/2;
    }

    fclose(in_vitlesu);
  }

  /*-------------------------------------------------------------------------*/
  /*  Begin Writing the Tecplot File header  (vitles_root.plt)               */
  /*-------------------------------------------------------------------------*/
  strcpy(tec_file,vitles_root);
  strcat(tec_file,".plt");
  out_plt = fopen(tec_file,"w");
  if (out_plt == NULL) {
    printf("%s: Can't open file %s\n",argv[0],tec_file);
    exit(1);
  }

  fprintf(out_plt,"TITLE = \"Solution: %s, Timestep %d\"\n",
                                                        argv[1],time_iteration);
  fprintf(out_plt,
             "VARIABLES = \"X\", \"Y\", \"Z\", \"U\", \"V\", \"W\", \"P\"\n\n");
 
  fprintf(out_plt,"ZONE N=%5d, E=%5d, DATAPACKING=POINT, ZONETYPE=FETETRAHEDRON\n",
                                                             n_nodes,7*n_elems);

  for (i=0; i<n_nodes; i++)
  {
    fprintf(out_plt,"%8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e\n",
                             x[i][0], x[i][1], x[i][2], u[i], v[i], w[i], p[i]);
  }

  // Send the element connectivity to the Tecplot file 
  for (nz=0; nz<n_zones; nz++)
  {
    sprintf(geometry_file,"%s.%03dgeom",vitles_root,nz);
    in_vitlesx = fopen(geometry_file,"rb");
    if (in_vitlesx==NULL)
    {
      printf("Can't open file %s\n",geometry_file);
      exit(1);
    }

    fread(&local_node      , sizeof(int)   , 1         , in_vitlesx);
    fread(&n_elzone        , sizeof(int)   , 1         , in_vitlesx);
    fread(&n_zone_unknowns , sizeof(int)   , 1         , in_vitlesx);
    fread(&n_dir           , sizeof(int)   , 1         , in_vitlesx);

    for (i=0; i<n_nodes; i++)
    {
      fread(&x[i][0]     , sizeof(double), 1 , in_vitlesx);
      fread(&x[i][1]     , sizeof(double), 1 , in_vitlesx);
      fread(&x[i][2]     , sizeof(double), 1 , in_vitlesx);
    }

    for (i=0; i<n_elzone; i++)
    {
      fread(&e_conn[i][0]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][1]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][2]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][3]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][4]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][5]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][6]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][7]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][8]    , sizeof(int), 1, in_vitlesx);
      fread(&e_conn[i][9]    , sizeof(int), 1, in_vitlesx);
    }
    fclose(in_vitlesx);

    for (e=0; e<n_elzone; e++)
    {
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][4]+1,e_conn[e][1]+1,e_conn[e][5]+1,e_conn[e][7]+1);
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][6]+1,e_conn[e][5]+1,e_conn[e][2]+1,e_conn[e][8]+1);
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][9]+1,e_conn[e][7]+1,e_conn[e][8]+1,e_conn[e][3]+1);
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][4]+1,e_conn[e][5]+1,e_conn[e][0]+1,e_conn[e][7]+1);
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][6]+1,e_conn[e][0]+1,e_conn[e][5]+1,e_conn[e][8]+1);
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][9]+1,e_conn[e][8]+1,e_conn[e][7]+1,e_conn[e][0]+1);
      fprintf(out_plt,"%7d %7d %7d %7d\n",
        e_conn[e][5]+1,e_conn[e][0]+1,e_conn[e][7]+1,e_conn[e][8]+1);
    }
  }

  fclose(out_plt);

  return 0;
}