// IO_HexaMesh.cpp
// James Liu, Graham Harper, ColoState; 2014/07--2018/05

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

#include "silo.h"

#include "HexaMesh.h"


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

int Write_ScalarData_HexaMesh_silo(const HexaMesh &mesh,
                                   double *ScalarData,
                                   char filename[], char varname[])
{
  // Declaring
  DBfile *fp = NULL;
  DBfacelist *fl;
  
  char silo_filename[64] = "";
  char silo_varname[64] = "";
  
  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
  
  int NumNds, NumEms;
  
  double *crd[3];
  PtVec3d P;
  
  int i, j, loc;
  int labelVertex[8];
  
  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] = (int)labelVertex[0];
    NodeList[loc+1] = (int)labelVertex[1];
    NodeList[loc+2] = (int)labelVertex[2];
    NodeList[loc+3] = (int)labelVertex[3];
    NodeList[loc+4] = (int)labelVertex[4];
    NodeList[loc+5] = (int)labelVertex[5];
    NodeList[loc+6] = (int)labelVertex[6];
    NodeList[loc+7] = (int)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, "HexaMesh", 3, NULL, crd, NumNds, NumEms,
               "ElemList", NULL, DB_DOUBLE, NULL);
  
  // Writing the scalar data into the silo file
  DBPutUcdvar1(fp, silo_varname, "HexaMesh", ScalarData, NumEms,
               NULL, 0, DB_DOUBLE, DB_ZONECENT, NULL);
  
  // Closing the silo file and cleanup
  DBClose(fp);
  delete[] NodeList;
  
  return(0);
}


// Appending scalar data on a hexahedral mesh to a silo data file

int Append_ScalarData_HexaMesh_silo(const HexaMesh &mesh,
                                    double *ScalarData,
                                    char filename[], char varname[])
{
  // Declaring
  DBfile *fp = NULL;
  
  char silo_filename[64] = "";
  char silo_varname[64] = "";
  
  int NumNds, NumEms;
  
  strcat(silo_filename, filename);
  strcat(silo_filename, ".silo");
  strcat(silo_varname, varname);
  
  // Retrieving mesh info
  NumNds = mesh.numberNodes();
  NumEms = mesh.numberElements();
  
  // Open the Silo file
  fp = DBOpen(silo_filename, DB_PDB, DB_APPEND);
  
  // Writing the scalar data into the silo file
  DBPutUcdvar1(fp, silo_varname, "HexaMesh", ScalarData, NumEms,
               NULL, 0, DB_DOUBLE, DB_ZONECENT, NULL);
  
  // Closing the silo file and cleanup
  DBClose(fp);
  
  return(0);
}


// JL20160728: TO BE REMOVED (after the above version works completely)
// Writing scalar data on a logically rect. hexahedral mesh to a silo data file
/*
int Write_ScalarData_LogiRectHexaMesh_silo_rsv(const LogiRectHexaMesh &mesh,
                                               double *ScalarData,
                                               char filename[], char varname[])
{
  // Declaring
  DBfile *fp = NULL;
  DBfacelist *fl;

  char silo_filename[64] = "";
  char silo_varname[64] = "";
  
  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
  
  int NumNds, NumEms;

  double *crd[3];
  PtVec3d P;

  int i, j, loc;
  int labelVertex[8];

  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] = (int)labelVertex[0];
    NodeList[loc+1] = (int)labelVertex[1];
    NodeList[loc+2] = (int)labelVertex[2];
    NodeList[loc+3] = (int)labelVertex[3];
    NodeList[loc+4] = (int)labelVertex[4];
    NodeList[loc+5] = (int)labelVertex[5];
    NodeList[loc+6] = (int)labelVertex[6];
    NodeList[loc+7] = (int)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, "HexaMesh", 3, NULL, crd, NumNds, NumEms,
               "ElemList", NULL, DB_DOUBLE, NULL);

  // Writing the scalar data into the silo file
  DBPutUcdvar1(fp, silo_varname, "HexaMesh", ScalarData, NumEms,
               NULL, 0, DB_DOUBLE, DB_ZONECENT, NULL);

  // Closing the silo file and cleanup
  DBClose(fp);
  delete[] NodeList;

  return(0);
}
*/

// IO_HexaMesh.cpp
