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

#include "GR_misc.h"
#include "GR_VolMesh.h"

// The following functions are to convert from little-endian to
// big-endian format.

static bool qFlipEndian = false;

static int convertToInt(const unsigned char raw[4])
{
  // This implementation should be okay even with 64 bit ints, whereas a
  // union might be a little tricky.
//   return ( (reinterpret_cast<int>(raw[0]) << 24)
// 	   + (reinterpret_cast<int>(raw[1]) << 16)
// 	   + (reinterpret_cast<int>(raw[2]) << 8)
// 	   + (raw[3]) );
  static const int is = sizeof(int);
  union {
    unsigned char raw2[is];
    int dummy;
  } a;
  a.raw2[is-1] = raw[0];
  a.raw2[is-2] = raw[1];
  a.raw2[is-3] = raw[2];
  a.raw2[is-4] = raw[3];
  return a.dummy;
}

static double convertToDouble(const unsigned char raw[8])
{
  union {
    unsigned char raw2[8];
    double dummy;
  } a;
  a.raw2[7] = raw[0];
  a.raw2[6] = raw[1];
  a.raw2[5] = raw[2];
  a.raw2[4] = raw[3];
  a.raw2[3] = raw[4];
  a.raw2[2] = raw[5];
  a.raw2[1] = raw[6];
  a.raw2[0] = raw[7];
/*   printf("%02x%02x%02x%02x%02x%02x%02x%02x\n", */
/* 	 raw[0], raw[1], raw[2], raw[3], */
/* 	 raw[4], raw[5], raw[6], raw[7]); */
  return a.dummy;
}

static int readInt(FILE * const fp)
{
  int input = -1;
  if (qFlipEndian) {
    // Just read it
    size_t nRead = fread(&input, sizeof(int), 1, fp);
    assert(nRead == 1);
  }
  else {
    // Read raw data and convert
    unsigned char raw[4];
    size_t nRead = fread(raw, 1, 4, fp);
    input = convertToInt(raw);
    assert(nRead == 4);
  }
  return input;
}

static double readDouble(FILE * const fp)
{
  double input = 0;
  if (qFlipEndian) {
    // Just read it
    size_t nRead = fread(&input, sizeof(double), 1, fp);
    assert(nRead == 1);
  }
  else {
    // Read raw data and convert
    unsigned char raw[8];
    size_t nRead = fread(raw, 1, 8, fp);
    input = convertToDouble(raw);
    assert(nRead == 8);
  }
  return input;
}

// VGrid writes unformatted Fortran files, containing only tets.  The
// info in the .cogsg file is very limited; BC info, for instance, is in
// a separate file, which we don't read at present.
void vReadVGridBinary(const char * const fileName, VolMesh * const pVM)
{
  int i, j;
//   int minVert = 1000000, maxVert = -1;

  int (*tetToVert)[4] = NULL;
//   int (*triToVert)[3] = NULL;

    double (*x)[3];

  FILE *fp = fopen(fileName, "r");

  int lineLength = readInt(fp);
  int junk = readInt(fp);
  int nTet = readInt(fp);
  int nVert = readInt(fp);
  int nBdryPts = readInt(fp);
  int nLayerPts = readInt(fp);
  int nLayerTets = readInt(fp);
  junk = readInt(fp);
  junk = readInt(fp);
  
  int nPredictedLength = (8 + 4*nTet)*4; // That last 4 is sizeof(int)...
  if (lineLength != nPredictedLength) {
    if (qFlipEndian == true) {
      vFatalError("File is neither big nor little endian",
		  "reading VGrid file");
    }
    qFlipEndian = true;
    fclose(fp);
    vReadVGridBinary(fileName, pVM);
    return;
  }    

  vMessage(2, "%8d verts\n", nVert);
  vMessage(2, "%8d bdry verts\n", nBdryPts);
  vMessage(2, "%8d tets\n", nTet);
  vMessage(2, "%8d verts in the bdry layer\n", nLayerPts);
  vMessage(2, "%8d tets in the bdry layer\n", nLayerTets);

  {
    // Estimate GRUMMP data struct sizes; the extra word in cell sizes
    // is a min for the allocated face ptr array overhead.
    int iTetSize = sizeof(TetCell) + 5*sizeof(void*);
    int iVertSize = sizeof(Vert) + 44*sizeof(void*); // Face connectivity
    int iFaceSize = sizeof(Face);

    int iVertMem = nVert * iVertSize/1024;
    int iTetMem = nTet * iTetSize/1024;

    int nTriFace = 2 * nTet;
    int iTriFaceMem = nTriFace * iFaceSize/1024;

    vMessage(2, "Expecting %8d kB of memory for verts\n", iVertMem);
    vMessage(2, "Expecting %8d kB of memory for tets\n", iTetMem);
    vMessage(2, "Expecting %8d kB of memory for tri faces\n", iTriFaceMem);

    vMessage(2, "Total memory approx %d MB + overhead\n",
	     (iVertMem + iTriFaceMem) / 1024);
  }

  /* Now read connectivity. */

  vMessage(2, "Reading tet connectivity\n");
  
  tetToVert = new int[nTet][4];
  
  for (j = 0; j < 4; j++) {
    for (i = 0; i < nTet; i++) {
      tetToVert[i][j] = readInt(fp) - 1;
//       minVert = (tetToVert[i][j] < minVert) ? tetToVert[i][j] : minVert;
//       maxVert = (tetToVert[i][j] > maxVert) ? tetToVert[i][j] : maxVert;
    }
  }

  // minVert should always be 1; this is a Fortran file we're reading.
//   vMessage(2, "Min vert index = %d; max = %d\n", minVert, maxVert);

  int lineCheck = readInt(fp);
  assert(lineCheck == lineLength);

  vMessage(2, "\nReading coordinates\n");

  x = new double[nVert][3];

  lineLength = readInt(fp);
  for (j = 0; j < 3; j++) {
    for (i = 0; i < nVert; i++) {
      x[i][j] = readDouble(fp);
    }
  }
  lineCheck = readInt(fp);
  assert(lineCheck == lineLength);
  
  pVM->vSetupVerts(nVert);

  for (int iV = 0; iV < nVert; iV++) {
    Vert *pV = pVM->pVVert(iV);
    pV->vSetCoords(3, x[iV]);
  }
  vMessage(2, "Set coords of %d vertices\n", nVert);
  delete [] x;

#ifdef BUILD_TETS_SLOWLY
  for (int iT = 0; iT < nTet; iT++) {
    Vert *pV0 = pVM->pVVert(tetToVert[iT][0]);
    Vert *pV1 = pVM->pVVert(tetToVert[iT][1]);
    Vert *pV2 = pVM->pVVert(tetToVert[iT][2]);
    Vert *pV3 = pVM->pVVert(tetToVert[iT][3]);
    
    pVM->createTetCell(pV0, pV1, pV2, pV3, 1);
    if ((iT + 1) % (nTet / 10 + 1) == 0) {
      vMessage(2, "Create %8d tets out of %d\n", iT + 1, nTet);
    }
  }
  vMessage(2, "Created %8d tets\n", nTet);
  delete [] tetToVert;

  vMessage(2, "%8d tri faces.\n",
	   pVM->iNumTriFaces());
  
  // Iterate over all triangles; those with only one tet are bdry faces,
  // so add them.

  int nTri = pVM->iNumTriFaces();
  int nBdryTri = 2*nBdryPts - 4;
  for (int iF = 0; iF < nTri; iF++) {
    Face *pF = pVM->pFFace(iF);
    if (pF->pCCellLeft()->qValid() && pF->pCCellRight()->qValid())
      continue;
    
    BFace *pBF = pVM->createBFace(pF);

    if ((iF + 1) % (nTri / 10 + 1) == 0) {
	vMessage(2, "Creating bdry tris: %8d tris of %d; %6d bdry tris of %d\n",
		 iF + 1, nTri, pVM->iNumBdryFaces(), nBdryTri);
    }
    delete [] triToVert;
  }
  vMessage(2, "Created bdry tris: %6d bdry tris of %d\n",
	   pVM->iNumBdryFaces(), nBdryTri);
#else
  // Use vConvertFromCellVert instead.
  int bdryTriToVert[1][3], bcList1[1], *bcList2;
  pVM->vConvertFromCellVert(nTet, 0, tetToVert, bdryTriToVert, bcList1,
			    bcList2);
#endif

  int iNumOK = 0;
  for (int iF = 0; iF < pVM->iNumFaces(); iF++) {
    Face *pF = pVM->pFFace(iF);
    int iOK = pF->pCCellLeft()->qValid() && pF->pCCellRight()->qValid();
    iNumOK += iOK;
    if (iOK == 0) {
      vMessage(3, "Problem with face: verts %8d %8d %8d (face %d)\n",
	       pVM->iVertIndex(pF->pVVert(0)),
	       pVM->iVertIndex(pF->pVVert(1)),
	       pVM->iVertIndex(pF->pVVert(2)),
	       iF);
    }
  }
  vMessage(2, "%d faces okay out of %d\n", iNumOK, pVM->iNumFaces());
    
  assert(pVM->qValid());
  vMessage(2, "Done reading VGrid file\n");
}
