// IO_BrickMesh.cpp
// James Liu, ColoState; 2014/07--2016/05

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

#include "silo.h"

#include "BrickMesh.h"
#include "PtVec3d.h"


// Writing scalar data on a brick mesh to BOV data files
// (coupled: bov-header & value)

int Write_ScalarData_BrickMesh_BOV(const BrickMesh &mesh, double *ScalarData,
  char filename[], char varname[])
{
  FILE *fph, *fpv;
  char BOV_header_filename[64] = "";
  char BOV_value_filename[64] = "";

  double xa = mesh.xa();
  double xb = mesh.xb();
  double yc = mesh.yc();
  double yd = mesh.yd();
  double ze = mesh.ze();
  double zf = mesh.zf();

  int nx = mesh.nx();
  int ny = mesh.ny();
  int nz = mesh.nz();

  strcat(BOV_header_filename, filename);
  strcat(BOV_header_filename, ".bov");

  strcat(BOV_value_filename, filename);
  strcat(BOV_value_filename, ".value");

  // Writing the BOV header file
  fph = fopen(BOV_header_filename, "wt");
  fprintf(fph, "TIME: 0.0\n");
  fprintf(fph, "DATA_FILE: %s\n", BOV_value_filename);
  fprintf(fph, "DATA_SIZE: %d  %d  %d\n", nx, ny, nz);
  fprintf(fph, "DATA_FORMAT: DOUBLE\n");
  fprintf(fph, "VARIABLE: %s\n", varname);
#ifdef WORDS_BIGENDIAN
  fprintf(fph, "DATA_ENDIAN: BIG\n");
#else
  fprintf(fph, "DATA_ENDIAN: LITTLE\n");
#endif
  fprintf(fph, "CENTERING: zonal\n");
  fprintf(fph, "BRICK_ORIGIN: %lf  %lf  %lf\n", xa, yc, ze);
  fprintf(fph, "BRICK_SIZE: %lf  %lf  %lf\n", xb-xa, yd-yc, zf-ze);
  fclose(fph);

  // Writing the BOV value file
  fpv = fopen(BOV_value_filename, "wb");
  fwrite(ScalarData, sizeof(double), nx*ny*nz, fpv);
  fclose(fpv);

  return(0);
}


// Writing scalar data on a brick mesh to a silo data file

int Write_ScalarData_BrickMesh_silo(const BrickMesh &mesh, double *ScalarData,
  char filename[], char varname[])
{
  // Declaring
  DBfile *fp = NULL;
  char silo_filename[64] = "";
  // char silo_varname[64] = "";
  char *DimName[3];
  int DimSz[3];
  float *crd[3];

  strcat(silo_filename, filename);
  strcat(silo_filename, ".silo");

  // Setting the dimension names, e.g., "X", "Y", "Z"
  DimName[0] = strdup("X");
  DimName[1] = strdup("Y");
  DimName[2] = strdup("Z");

  // Setting the dimension sizes (number of nodes in each dimension)
  DimSz[0] = (int)mesh.nx() + 1;
  DimSz[1] = (int)mesh.ny() + 1;
  DimSz[2] = (int)mesh.nz() + 1;

  // Allocating ...
  float *nodex = new float[DimSz[0]];
  float *nodey = new float[DimSz[1]];
  float *nodez = new float[DimSz[2]];

  // Retrieving x,y,z-partition coordinates
  for (int i=0; i<DimSz[0]; ++i)  nodex[i] = (float)mesh.xCrd(i);
  for (int j=0; j<DimSz[1]; ++j)  nodey[j] = (float)mesh.yCrd(j);
  for (int k=0; k<DimSz[2]; ++k)  nodez[k] = (float)mesh.zCrd(k);

  // Assembling the array of coordinate arrays
  crd[0] = nodex;
  crd[1] = nodey;
  crd[2] = nodez;

  // Creating a Silo file
  fp = DBCreate(silo_filename, DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);

  // Writing the mesh data into the silo file
  DBPutQuadmesh(fp, "BrickMesh", DimName, crd, DimSz, 3,
                DB_FLOAT, DB_COLLINEAR, NULL);

  // Resizing
  DimSz[0] = (int)mesh.nx();
  DimSz[1] = (int)mesh.ny();
  DimSz[2] = (int)mesh.nz();

  // Writing the scalar data into the silo file
  DBPutQuadvar1(fp, "NumerPres", "BrickMesh", ScalarData,
                DimSz, 3,
                NULL, 0,
                DB_DOUBLE, DB_ZONECENT, NULL);

  // Closing the silo file
  DBClose(fp);

  return(0);
}


int Write_VectorData_BrickMesh_silo(const BrickMesh &mesh, double *VectorData[3],
  char filename[], char varname[])
{
  // Declaring
  DBfile *fp = NULL;
  char silo_filename[64] = "";
  // char silo_varname[64] = "";
  char *DimName[3], *VarName[3];
  int DimSz[3];
  float *crd[3];
  
  strcat(silo_filename, filename);
  strcat(silo_filename, ".silo");
  
  // Setting the dimension names, e.g., "X", "Y", "Z"
  DimName[0] = strdup("X");
  DimName[1] = strdup("Y");
  DimName[2] = strdup("Z");
  
  // Setting the dimension names, e.g., "u", "v", "w"
  VarName[0] = strdup("u");
  VarName[1] = strdup("v");
  VarName[2] = strdup("w");

  // Setting the dimension sizes (number of nodes in each dimension)
  DimSz[0] = (int)mesh.nx() + 1;
  DimSz[1] = (int)mesh.ny() + 1;
  DimSz[2] = (int)mesh.nz() + 1;
  
  // Allocating ...
  float *nodex = new float[DimSz[0]];
  float *nodey = new float[DimSz[1]];
  float *nodez = new float[DimSz[2]];
  
  // Retrieving x,y,z-partition coordinates
  for (int i=0; i<DimSz[0]; ++i)  nodex[i] = (float)mesh.xCrd(i);
  for (int j=0; j<DimSz[1]; ++j)  nodey[j] = (float)mesh.yCrd(j);
  for (int k=0; k<DimSz[2]; ++k)  nodez[k] = (float)mesh.zCrd(k);
  
  // Assembling the array of coordinate arrays
  crd[0] = nodex;
  crd[1] = nodey;
  crd[2] = nodez;
  
  // Creating a Silo file
  fp = DBCreate(silo_filename, DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
  
  // Writing the mesh data into the silo file
  DBPutQuadmesh(fp, "BrickMesh", DimName, crd, DimSz, 3,
                DB_FLOAT, DB_COLLINEAR, NULL);
  
  // Resizing
  DimSz[0] = (int)mesh.nx();
  DimSz[1] = (int)mesh.ny();
  DimSz[2] = (int)mesh.nz();
  
  // Writing the vector data into the silo file
  DBPutQuadvar(fp, "NumerVel", "BrickMesh", 3, VarName, VectorData,
               DimSz, 3,
               NULL, 0,
               DB_DOUBLE, DB_ZONECENT, NULL);

  // Closing the silo file
  DBClose(fp);
  
  return(0);
}


// Writing scalar data on a brick mesh (TREATED AS A UCD MESH!!)
// to a silo data file

int Write_ScalarData_BrickMeshUcd_silo(const BrickMesh &mesh,
  double *ScalarData, char filename[], char varname[])
{
  // Declaring
  DBfile *fp = NULL;
  // DBfacelist *fl;
  char silo_filename[64] = "";
  char silo_varname[64] = "";
  int NumNds, NumEms;
  int i, loc;
  int labelVertex[8];
  double *crd[3];
  PtVec3d P;

  int NumElemTypes = 1;
  int ElemType[1], ElemTypeCnt[1], ElemTypeSize[1];
  int LenNodeList, *NodeList;
  
/*
  int nfaces, nzones, nnodes;
  int lfacelist, lzonelist;
  int fshapesize, fshapecnt, zshapesize, zshapecnt;
  int facelist[11200];  // JL20150321: TO BE REVISED
  int zoneno[2800];  // JL20150321: TO BE REVISED
*/
  
  strcat(silo_filename, filename);
  strcat(silo_filename, ".silo");
  strcat(silo_varname, varname);
  
  // Retrieving mesh info
  NumNds = mesh.numberNodes();
  NumEms = mesh.numberElements();
  crd[0] = new double[NumNds];
  crd[1] = new double[NumNds];
  crd[2] = new double[NumNds];
  for (i=0; i<NumNds; ++i) {
    P = mesh.node(i+1);
    crd[0][i] = P.xCrd();
    crd[1][i] = P.yCrd();
    crd[2][i] = P.zCrd();
  }
  
  // Creating a "zone"/element/cell list
  ElemType[0] = DB_ZONETYPE_HEX;
  ElemTypeSize[0] = 8;
  ElemTypeCnt[0] = NumEms;
  LenNodeList = NumEms*8;
  NodeList = new int[LenNodeList];
  loc = 0;
  for (i=1; i<=NumEms; ++i) {
    mesh.getElementNode(i,labelVertex);
    NodeList[loc+0] = labelVertex[0];
    NodeList[loc+1] = labelVertex[1];
    NodeList[loc+2] = labelVertex[2];
    NodeList[loc+3] = labelVertex[3];
    NodeList[loc+4] = labelVertex[4];
    NodeList[loc+5] = labelVertex[5];
    NodeList[loc+6] = labelVertex[6];
    NodeList[loc+7] = labelVertex[7];
    loc = loc + 8;
  }

/*
  // Calculating the external face list
  fl = DBCalcExternalFacelist(NodeList, NumNds, 1,
                              ElemTypeSize, ElemTypeCnt, NumElemTypes,
                              NULL, 0);
  
  nfaces = fl->nfaces;
  fshapecnt = fl->nfaces;
  fshapesize = 4;
  lfacelist = fl->lnodelist;
  for (i=0; i<lfacelist; i++)  facelist[i] = fl->nodelist[i];
  for (i=0; i<nfaces; i++)  zoneno[i] = fl->zoneno[i];
  DBFreeFacelist(fl);
*/
  
  // Creating a Silo file
  fp = DBCreate(silo_filename, DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
  
  // Writing the zone/element/cell list into the silo file
  DBPutZonelist2(fp, "ElemList", NumEms, 3, NodeList, LenNodeList, 1, 0, 0,
                 ElemType, ElemTypeSize, ElemTypeCnt, NumElemTypes, NULL);

/*
  // Writing the face list into the silo file
  DBPutFacelist(fp, "ExtFaceList", fl->nfaces, 3,
                facelist, lfacelist, 0, zoneno,
                &fshapesize, &fshapecnt, 1, NULL, NULL, 0);
*/
  
  // Writing the mesh data into the silo file
  DBPutUcdmesh(fp, "BrickMeshUcd", 3, NULL, crd, NumNds, NumEms,
               "ElemList", NULL, DB_DOUBLE, NULL);
  
  // Writing the scalar data into the silo file
  DBPutUcdvar1(fp, silo_varname, "BrickMeshUcd", ScalarData, NumEms,
               NULL, 0, DB_DOUBLE, DB_ZONECENT, NULL);
  
  // Closing the silo file
  DBClose(fp);
  
  delete[] NodeList;

  return(0);
}

// IO_BrickMesh.cpp