/*-------------------------------------------------------------------------78--
 * Program to read in 3D finite element data from ViTLES simulations
 * and generate output files in VTK format.
 *
 *     *** Coarse Linear Element Version ***
 *
 * This version creates 4 noded linear elements (vtk type=10) from each
 * quadratic element.
 *
 * The VTK format includes:
 *               vitles_root.vtk     (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_vtk;    // file used to write VTK file

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

  // input arguments
  char  vitles_root[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;
  double *x;
  double *u;
  double *v;
  double *w;
  double *p;

  /*-------------------------------------------------------------------------*/
  /*  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

  /*--------------------------------------------------------------------------*/
  // get total number of elements
  /*--------------------------------------------------------------------------*/
  n_elems   = 0;

  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);

    n_elems = n_elems + n_elzone;

    fclose(in_vitlesx);
  }

  //----------------------------------------------------------------------------
  //  Allocate Storage Space (malloc/calloc)
  //----------------------------------------------------------------------------
  e_conn         = malloc(n_elems*10* sizeof(int));
  x              = malloc(n_nodes*3 * sizeof(double));
  u              = malloc(n_nodes   * sizeof(double));
  v              = malloc(n_nodes   * sizeof(double));
  w              = malloc(n_nodes   * sizeof(double));
  p              = malloc(n_nodes   * sizeof(double));

  //----------------------------------------------------------------------------
  //  Begin Writing the VTK File header  (vitles_root.vtk)
  //----------------------------------------------------------------------------
  sprintf(vtk_file,"%s_%04d.vtk",vitles_root,time_iteration);

  out_vtk = fopen(vtk_file,"w");
  if (out_vtk == NULL) {
    printf("%s: Can't open file %s\n",argv[0],vtk_file);
    exit(1);
  }

  fprintf(out_vtk,"# vtk DataFile Version 2.0\n");
  fprintf(out_vtk,"Solution: %s, Timestep %d\n",argv[1],time_iteration);
  fprintf(out_vtk,"ASCII\n\n");

  /*--------------------------------------------------------------------------*/
  /*  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);

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

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

    if (nz==0)
    {
      fprintf(out_vtk,"DATASET UNSTRUCTURED_GRID\n");
      fprintf(out_vtk,"POINTS %8d double\n",n_nodes);

      for (i=0; i<n_nodes; i++)
      {
        fprintf(out_vtk,"%12.4e %12.4e %12.4e\n",x[0+i*3],x[1+i*3],x[2+i*3]);
      }
      fprintf(out_vtk,"\n");
      fprintf(out_vtk,"CELLS %7d %7d\n",n_elems,5*n_elems);
//    fprintf(out_vtk,"CELLS %7d %7d\n",n_elems,11*n_elems);
    }

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

    //--------------------------------------------------------------------------
    // read 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[i+e*10]] = D[i];
      }

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

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

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

    fclose(in_vitlesu);
  }

  //----------------------------------------------------------------------------
  // Since all of the solution variables are known, continue writing
  // the vtk file.
  //----------------------------------------------------------------------------

  fprintf(out_vtk,"\nCELL_TYPES %d\n",n_elems);
  for (e=0; e<n_elems; e++)
  {
    fprintf(out_vtk,"10 \n");
    // fprintf(out_vtk,"24\n");
  }

  fprintf(out_vtk,"\nPOINT_DATA %d\n",n_nodes);
  fprintf(out_vtk,"\nSCALARS pressure double 1\n");
  fprintf(out_vtk,"LOOKUP_TABLE default\n");
  for (i=0; i<n_nodes; i++)
  {
    fprintf(out_vtk," %12.4e\n", p[i]);
  }

  fprintf(out_vtk,"VECTORS velocity double\n");
  for (i=0; i<n_nodes; i++)
  {
    fprintf(out_vtk," %14.6e %14.6e %14.6e\n", u[i],v[i],w[i]);
  }


  fclose(out_vtk);

  return 0;
}