// IO_TetraMesh.cpp
// James Liu, ColoState; 2014/07--2016/10

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

#include "silo.h"
// #include "tetgen.h"

#include "PtVec3d.h"
#include "TetraMesh.h"


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

// silo

// Writing scalar data on a tetrahedral mesh into a silo data file
// The tetrahedral mesh is treated as a UCD mesh

int Write_ScalarData_TetraMesh_silo(const TetraMesh &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[4];
  
  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_TET;
  ElemTypeSize[0] = 4;
  ElemTypeCnt[0] = NumEms;
  LenNodeList = NumEms*4;
  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[3];
    NodeList[loc+3] = labelVertex[2];
    loc = loc + 4;
  }
  
  // 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, "TetraMesh", 3, NULL, crd, NumNds, NumEms,
               "ElemList", NULL, DB_DOUBLE, NULL);
  
  // Writing the scalar data into the silo file
  DBPutUcdvar1(fp, silo_varname, "TetraMesh", ScalarData, NumEms,
               NULL, 0, DB_DOUBLE, DB_ZONECENT, NULL);
  
  // Closing the silo file
  DBClose(fp);
  
  delete[] NodeList;

  return(0);
}


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

// TetGen
/*
int Read_MeshData_TetraMesh_TetGen(TetraMesh &mesh, char filenameNode[],
                                   char filenameFace[], char filenameElement[])
{
  FILE *fp0, *fp2, *fp3;
  
  int i, j, tmpn;
  double tmpd;

  int numberNodes, numberFaces, numberElements;
  int (*labelFaceNode)[3], (*labelFaceElement)[2], (*labelElementNode)[4];
  int *boundaryNodeMark, *boundaryFaceMark;
  double (*crd)[3];

  int numberofpoints, numberoftrifaces, numberoftetrahedra;
  int *pointmarkerlist, *trifacelist, *tetrahedronlist;
  double *pointlist;

  // Reading the node file
  if (NULL==(fp0=fopen(filenameNode, "r"))) {
    std::puts("Open node data file failed.  Exit.\n");
    exit(-1);
  }
  //
  fscanf(fp0, "%7d", &numberofpoints);
  fscanf(fp0, "%2d", &tmpn);
  fscanf(fp0, "%2d", &tmpn);
  fscanf(fp0, "%2d", &tmpn);
  //
  pointlist = new double[numberofpoints*3];
  pointmarkerlist = new int[numberofpoints];
  //
  for (i=0; i<numberofpoints; ++i) {
    fscanf(fp0, "7d", &tmpn);
    fscanf(fp0, "20.17f", &(pointlist[3*i+0]));
    fscanf(fp0, "20.17f", &(pointlist[3*i+1]));
    fscanf(fp0, "20.17f", &(pointlist[3*i+2]));
    fscanf(fp0, "3d", &(pointmarkerlist[i]));
  }
  //
  fclose(fp0);
  
  // Reading the face file
  if (NULL==(fp2=fopen(filenameFace, "r"))) {
    std::puts("Open face data file failed.  Exit.\n");
    exit(-1);
  }
  //
  fscanf(fp2, "%7d", &numberoftrifaces);
  fscanf(fp2, "%2d", &tmpn);
  //
  trifacelist = new int[numberoftrifaces];
  for (i=0; i<numberoftrifaces; ++i) {
    fscanf(fp2, "%7d", &tmpn);
    fscanf(fp2, "%7d", &(trifacelist[6*i+0]));
    fscanf(fp2, "%7d", &(trifacelist[6*i+1]));
    fscanf(fp2, "%7d", &(trifacelist[6*i+2]));
    fscanf(fp2, "%2d", &(trifacelist[6*i+3]));
    fscanf(fp2, "%7d", &(trifacelist[6*i+4]));
    fscanf(fp2, "%7d", &(trifacelist[6*i+5]));
  }
  //
  fclose(fp2);

  // Reading the element file
  //
  if (NULL==(fp3=fopen(filenameElement, "r"))) {
    std::puts("Open element data file failed.  Exit.\n");
    exit(-1);
  }
  //
  fscanf(fp3, "%7d", &numberoftetrahedra);
  fscanf(fp3, "%2d", &tmpn);
  fscanf(fp3, "%2d", &tmpn);
  //
  tetrahedronlist = new int[numberoftetrahedra*4];
  //
  for (i=0; i<numberoftetrahedra; ++i) {
    fscanf(fp3, "%7d", &tmpn);
    fscanf(fp3, "%7d", &(tetrahedronlist[4*i+0]));
    fscanf(fp3, "%7d", &(tetrahedronlist[4*i+1]));
    fscanf(fp3, "%7d", &(tetrahedronlist[4*i+2]));
    fscanf(fp3, "%7d", &(tetrahedronlist[4*i+3]));
  }
  //
  fclose(fp3);

  // Counting
  numberNodes = numberofpoints;
  numberFaces = numberoftrifaces;
  numberElements = numberoftetrahedra;

  // Allocatiion for filling-in the tetrahedral mesh
  labelFaceNode = new int[numberFaces][3];
  labelFaceElement = new int[numberFaces][2];
  labelElementNode = new int[numberElements][4];
  boundaryNodeMark = new int[numberNodes];
  boundaryFaceMark = new int[numberFaces];
  crd = new double[numberNodes][3];

  // Converting mesh data
  for (i=0; i<numberofpoints; ++i) {
    crd[i][0] = pointlist[i*3+0];
    crd[i][1] = pointlist[i*3+1];
    crd[i][2] = pointlist[i*3+2];
  }
  //
  for (i=0; i<numberofpoints; ++i)
    boundaryNodeMark[i] = pointmarkerlist[i];
  //
  for (i=0; i<numberoftrifaces; ++i) {
    labelFaceNode[i][0] = trifacelist[6*i+0];
    labelFaceNode[i][1] = trifacelist[6*i+1];
    labelFaceNode[i][2] = trifacelist[6*i+2];
    boundaryFaceMark[i] = trifacelist[6*i+3];
    labelFaceElement[i][0] = trifacelist[6*i+4];
    labelFaceElement[i][1] = trifacelist[6*i+5];
  }
  //
  for (i=0; i<numberoftetrahedra; ++i)
    for (j=0; j<4; ++j)
      labelElementNode[i][j] = tetrahedronlist[4*i+j];

  // Releasing
  delete[] pointmarkerlist, trifacelist, tetrahedronlist;
  delete[] pointlist;

  // Filling-in for the tetrahedral mesh
  mesh.fillNodeInfo(numberNodes, crd, boundaryNodeMark);
  mesh.fillFaceInfo(numberFaces, labelFaceNode, labelFaceElement,
                    boundaryFaceMark);
  mesh.fillElementInfo(numberElements, labelElementNode);
  
  // Releasing
  delete[] boundaryNodeMark, boundaryFaceMark;
  delete[] labelFaceNode, labelFaceElement, labelElementNode;
  delete[] crd;

  return(0);  // if successful
}
*/
// IO_TetraMesh.cpp