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

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

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

void vReadUGridBinary(const char * const fileName, VolMesh * const pVM)
{
  int nTet, nPyr, nPrism, nHex;
  int nVert, nBdryTri, nBdryQuad;
  int i, j;
  int minVert = 1000000, maxVert = -1;

  int (*tetToVert)[4] = NULL;
  int (*pyrToVert)[5] = NULL;
  int (*prismToVert)[6] = NULL;
  int (*hexToVert)[8] = NULL;

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

  unsigned char raw[36];

  double (*x)[3];

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

  size_t nRead = fread(raw, 1, 28, fp);

  nVert = convertToInt(raw);
  nBdryTri = convertToInt(raw + 4);
  nBdryQuad = convertToInt(raw + 8);
  nTet = convertToInt(raw + 12);
  nPyr = convertToInt(raw + 16);
  nPrism = convertToInt(raw + 20);
  nHex = convertToInt(raw + 24);

  vMessage(2, "%8d verts\n", nVert);
  vMessage(2, "%8d bdry tris\n", nBdryTri);
  vMessage(2, "%8d bdry quads\n", nBdryQuad);
  vMessage(2, "%8d tets\n", nTet);
  vMessage(2, "%8d pyrs\n", nPyr);
  vMessage(2, "%8d prisms\n", nPrism);
  vMessage(2, "%8d hexes\n\n", nHex);

  {
    // 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 iPyrSize = sizeof(Cell) + 6*sizeof(void*);
    int iPrismSize = sizeof(Cell) + 7*sizeof(void*);
    int iHexSize = sizeof(Cell) + 9*sizeof(void*);
    int iVertSize = sizeof(Vert) + 35*sizeof(void*); // Face connectivity
    int iFaceSize = sizeof(Face);

    int iVertMem = nVert * iVertSize/1024;
    int iTetMem = nTet * iTetSize/1024;
    int iPyrMem = nPyr * iPyrSize/1024;
    int iPrismMem = nPrism * iPrismSize/1024;
    int iHexMem = nHex * iHexSize/1024;

    int nTriFace = (4 * (nTet + nPyr) + 2 * nPrism) / 2 - nBdryTri;
    int nQuadFace = (3 * nPrism + nPyr + 6 * nHex) / 2 - nBdryQuad;

    int iTriFaceMem = nTriFace * iFaceSize/1024;
    int iQuadFaceMem = nQuadFace * 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 pyrs\n", iPyrMem);
    vMessage(2, "Expecting %8d kB of memory for prisms\n", iPrismMem);
    vMessage(2, "Expecting %8d kB of memory for hexes\n", iHexMem);
    vMessage(2, "Expecting %8d kB of memory for tri faces\n", iTriFaceMem);
    vMessage(2, "Expecting %8d kB of memory for quad faces\n",
	     iQuadFaceMem);

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

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

  x = new double[nVert][3];

  for (i = 0; i < nVert; i++) {
    nRead = fread(raw, 1, 24, fp);
    x[i][0] = convertToDouble(raw);
    x[i][1] = convertToDouble(raw+8);
    x[i][2] = convertToDouble(raw+16);
  }

  /* Now read connectivity. */

  if (nBdryTri > 0) {
    vMessage(2, "Reading bdry tri connectivity\n");
    triToVert = new int[nBdryTri][3];
    for (i = 0; i < nBdryTri; i++) {
      for (j = 0; j < 3; j++) {
	nRead = fread(raw, 1, 4, fp);
	triToVert[i][j] = convertToInt(raw);
      }
    }
  }    

  if (nBdryQuad > 0) {
    vMessage(2, "Reading bdry quad connectivity\n");
    quadToVert = new int[nBdryQuad][4];
    for (i = 0; i < nBdryQuad; i++) {
      for (j = 0; j < 4; j++) {
	nRead = fread(raw, 1, 4, fp);
	quadToVert[i][j] = convertToInt(raw);
      }
    }
  }    

  if (nBdryTri + nBdryQuad > 0) {
    // Read bdry condition info and throw it away.
    for (i = 0; i < nBdryTri + nBdryQuad; i++) {
      nRead = fread(raw, 1, 4, fp);
    }
  }
  
  if (nTet > 0) {
    vMessage(2, "Reading tet connectivity\n");
    
    tetToVert = new int[nTet][4];
    
    for (i = 0; i < nTet; i++) {
      for (j = 0; j < 4; j++) {
	nRead = fread(raw, 1, 4, fp);
	tetToVert[i][j] = convertToInt(raw);
	minVert = (tetToVert[i][j] < minVert) ? tetToVert[i][j] : minVert;
	maxVert = (tetToVert[i][j] > maxVert) ? tetToVert[i][j] : maxVert;
      }
    }
  }

  if (nPyr > 0) {
    vMessage(2, "Reading pyramid connectivity\n");
    
    pyrToVert = new int[nPyr][5];
    
    for (i = 0; i < nPyr; i++) {
      for (j = 0; j < 5; j++) {
	nRead = fread(raw, 1, 4, fp);
	pyrToVert[i][j] = convertToInt(raw);
	minVert = (pyrToVert[i][j] < minVert) ? pyrToVert[i][j] : minVert;
	maxVert = (pyrToVert[i][j] > maxVert) ? pyrToVert[i][j] : maxVert;
      }
    }
  }
  
  if (nPrism > 0) {
    vMessage(2, "Reading prism connectivity\n");
    
    prismToVert = new int[nPrism][6];
    
    for (i = 0; i < nPrism; i++) {
      for (j = 0; j < 6; j++) {
	nRead = fread(raw, 1, 4, fp);
	prismToVert[i][j] = convertToInt(raw);
	minVert = (prismToVert[i][j] < minVert) ? prismToVert[i][j] : minVert;
	maxVert = (prismToVert[i][j] > maxVert) ? prismToVert[i][j] : maxVert;
      }
    }
  }

  if (nHex > 0) {
    vMessage(2, "Reading hex connectivity\n");
    
    hexToVert = new int[nHex][8];
    
    for (i = 0; i < nHex; i++) {
      for (j = 0; j < 8; j++) {
	nRead = fread(raw, 1, 4, fp);
	hexToVert[i][j] = convertToInt(raw);
	minVert = (hexToVert[i][j] < minVert) ? hexToVert[i][j] : minVert;
	maxVert = (hexToVert[i][j] > maxVert) ? hexToVert[i][j] : maxVert;
      }
    }
  }

  vMessage(2, "Min vert index = %d; max = %d\n", minVert, maxVert);
  
  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;

  // Now, if minVert = 1 (Fortran numbering), decrement the indices in
  // all the connectivity tables.
  assert(minVert == 0 || minVert == 1);
  if (minVert == 1) {
    for (int iT = 0; iT < nBdryTri; iT++) {
      triToVert[iT][0] --;
      triToVert[iT][1] --;
      triToVert[iT][2] --;
    }

    for (int iQ = 0; iQ < nBdryQuad; iQ++) {
      quadToVert[iQ][0] --;
      quadToVert[iQ][1] --;
      quadToVert[iQ][2] --;
      quadToVert[iQ][3] --;
    }

    for (int iT = 0; iT < nTet; iT++) {
      tetToVert[iT][0] --;
      tetToVert[iT][1] --;
      tetToVert[iT][2] --;
      tetToVert[iT][3] --;
    }

    for (int iP = 0; iP < nPyr; iP++) {
      pyrToVert[iP][0] --;
      pyrToVert[iP][1] --;
      pyrToVert[iP][2] --;
      pyrToVert[iP][3] --;
      pyrToVert[iP][4] --;
    }

    for (int iP = 0; iP < nPrism; iP++) {
      prismToVert[iP][0] --;
      prismToVert[iP][1] --;
      prismToVert[iP][2] --;
      prismToVert[iP][3] --;
      prismToVert[iP][4] --;
      prismToVert[iP][5] --;
    }

    for (int iH = 0; iH < nHex; iH++) {
      hexToVert[iH][0] --;
      hexToVert[iH][1] --;
      hexToVert[iH][2] --;
      hexToVert[iH][3] --;
      hexToVert[iH][4] --;
      hexToVert[iH][5] --;
      hexToVert[iH][6] --;
      hexToVert[iH][7] --;
    }
  }
  // Done decrementing all the connectivity tables.

  // Now do the conversion
  pVM->vConvertFromCellVert(nTet, nPyr, nPrism, nHex, nBdryTri, nBdryQuad,
			    tetToVert, pyrToVert, prismToVert, hexToVert,
			    triToVert, quadToVert);
  if (tetToVert)   free(tetToVert);
  if (pyrToVert)   free(pyrToVert);
  if (prismToVert) free(prismToVert);
  if (hexToVert)   free(hexToVert);
  if (triToVert)   free(triToVert);
  if (quadToVert)  free(quadToVert);
}
