#ifndef WIN32
#include <unistd.h>
#else /* WIN32 defined */
int getopt(int argc, char *argv[], char *opstring);
#endif

#include <stdlib.h>

#include <math.h>
#include <string.h>
#include "GR_AdaptPred.h"
#include "GR_Bdry3D.h"
#include "GR_InsertionQueue.h"
#include "GR_SurfMesh.h"
#include "GR_VolMesh.h"
#include "GR_events.h"
#include "GR_misc.h"

#ifdef IRIX
#include <getopt.h>
#endif

static char GR_tetra_Usage[][80] = {
  "    -i    Base file name for input file",
  " ",
  "  Boolean options (argument of 1 for true, 0 for false):",
  "    -A    Use adaptive precision geometric predicates [0]",
  "    -b    Allow boundary changes [1]",
  "    -e    Allow edge swapping [1]",
  "    -m    Merge co-planar boundary patches as appropriate [1]",
  "  ",
  "  Numeric options",
  "    -a #  Max allowable angle for surface reconfiguration [5]",
  "    -c    Use diametral lenses instead of diametral circles [0]",
  "          0 -> circumradius < shortest edge length * 2 [good]",
  "          1 -> circumradius < shortest edge length * sqrt(2) [better]",
  "          2 -> circumradius < shortest edge length [best known result]",
  " ",
  "  Length scale options (values < 1 don't make sense):",
  "    -g #  Grading (rate of cell size change); larger is slower [1]",
  "    -r #  Resolution (feature size / cell size); larger is finer [1]",
  " ",
  "  Adaptation option: (meshopt2d only)",
  "    -l filename  Refine to length scale specified (per vertex) in file",
  " ",
  "    -o    Base file name for output file [input file base name]",
  " ",
  "  Mesh optimization string, in shorthand; see manual:",
  "    -O    [wsmf5]",
  " ",
  "    -q #  Mesh quality measure used [2]",
  "         0  Maximum dihedral angle",
  "         1  Minimum dihedral angle",
  "         2  All dihedral angles",
  "         3  Maximum solid angle",
  "         4  Minimum solid angle",
  "         5  All solid angles",
  "         6  Aspect ratio I (inscribed to circum radius)",
  "         7  Sliver quality measure",
  "         8  Skewness quality measure",
  "         9  Aspect ratio II (based on volume versus surface area)",
  "        10  Aspect ratio III (shortest edge to circumradius)",
  ""  // Must end with a null string.
};
  
static void vUsage(const char strCanonName[])
{
  vMessage(0, "Usage:  %s -i basefilename [options] \n", strCanonName);
  for (int i = 0; strlen(GR_tetra_Usage[i]) > 0; i++) {
    vMessage(0, "%s\n", GR_tetra_Usage[i]);
  }
  exit(1);
}

int main(int iNArg, char *apcArgs[]) {
  char strBaseFileName[FILE_NAME_LEN], strQualFileName[FILE_NAME_LEN];
  char strBaseOutFileName[FILE_NAME_LEN];
  char strBdryFileName[FILE_NAME_LEN];
  bool qReadOutFileName = false;
  bool qReadAuxBdryFileName = false;
  bool qLengthScaleGiven = false;
  char strLengthFileName[FILE_NAME_LEN];
  int iBdryOptLevel = 1;

  char strOptParams[256];
  char strExecName[FILE_NAME_LEN], strCanonicalName[20];
  enum {UNKNOWN, TETRA, MESHOPT, COARSEN, TRANSLATE} iExec = UNKNOWN;
  
  sprintf(strExecName, "%s", apcArgs[0]);
  if (strstr(strExecName, "coarsen3d")) {
    iExec = COARSEN;
    sprintf(strCanonicalName, "coarsen3d");
  }
  else if (strstr(strExecName, "meshopt3d")) {
    iExec = MESHOPT;
    sprintf(strCanonicalName, "meshopt3d");
  }
  else if (strstr(strExecName, "trans3d")) {
    iExec = TRANSLATE;
    sprintf(strCanonicalName, "trans3d");
  }
  else if (strstr(strExecName, "tetra")) {
    iExec = TETRA;
    sprintf(strCanonicalName, "tetra");
  }
  else {
    vFatalError("Invalid executable name.  Be less creative renaming files.\n",
		"tetra.cxx:main()");
  }
  
  // Set default values for arguments
  bool qDoEdgeSwapping = true;
  bool qAllowBdryChanges = true;
  double dGrading = 1;
  double dResolution = 1;
  double dMaxAngle = 5; // In degrees
  int iQualMeasure = 2;
  eEncroachType eET = eBall;
  // The following options are currently inactive in 3D.
//    bool qAnisoCoarsening = false;
  // This one is global, declared in GR_AdaptPred.h
  qAdaptPred = false;
  
  bool qFile = false, qError = false, qParams = false;
  int iPChar;
  while ((iPChar =
	  getopt(iNArg, apcArgs, "A:a:B:b:e:g:i:j:l:m:o:O:q:r:")) != EOF) {
    switch (iPChar) {
    case 'a':
      sscanf(optarg, "%lf", &dMaxAngle);
      if (dMaxAngle < 0) dMaxAngle = 0;
      if (dMaxAngle > 180) dMaxAngle = 180;
      break;
    case 'A':
      qAdaptPred = (optarg[0] == '1');
      break;
    case 'B':
      strncpy(strBdryFileName, optarg, FILE_NAME_LEN-1);
      qReadAuxBdryFileName = true;
      {
	int iDum;
	for (iDum = 0; iDum < FILE_NAME_LEN && strBdryFileName[iDum];
	     iDum++) {}
	if (iDum == FILE_NAME_LEN) // No null; file name too long
	  vFatalError("Boundary file name too long", "command line");
      }
      break;
    case 'b':
      qAllowBdryChanges = (optarg[0] == '1');
      break;
//     case 'c':
//       switch (optarg[0]) {
//       case '0':
// 	eET = eBall;
// 	break;
//       case '1':
//       default:
// 	eET = eLens;
// 	break;
//       case '2':
// 	eET = eNewLens;
// 	break;
//       }
//       break;
    case 'e':
      qDoEdgeSwapping = (optarg[0] == '1');
      break;
    case 'g':
      sscanf(optarg, "%lf", &dGrading);
      if (dGrading < 1) dGrading = 1;
      break;
    case 'i':
      if (qFile) qError = true;
      strncpy(strBaseFileName, optarg, FILE_NAME_LEN-1);
      {
	size_t iDum;
	for (iDum = 0; iDum < FILE_NAME_LEN && strBaseFileName[iDum];
	     iDum++) {}
	if (iDum == FILE_NAME_LEN) // No null; file name too long
	  vFatalError("File name too long", "command line");
      }
      qFile = true;
      break;
    case 'j':
//        qAnisoCoarsening = (optarg[0] == '1');
      break;
    case 'l':
      if (iExec == MESHOPT) {
	strncpy(strLengthFileName, optarg, FILE_NAME_LEN-1);
	{
	  int iDum;
	  for (iDum = 0; iDum < FILE_NAME_LEN && strLengthFileName[iDum];
	       iDum++) {}
	  if (iDum == FILE_NAME_LEN) // No null; file name too long
	    vFatalError("Lengthscale file name too long", "command line");
	}
	qLengthScaleGiven = true;
      }
      else {
	qError = true; // Be sure to crash out.
      }
      break;      
    case 'm':
      sscanf(optarg, "%d", &iBdryOptLevel);
      if (iBdryOptLevel != 0 && iBdryOptLevel != 1)
	qError = true;
      break;
    case 'o':
      strncpy(strBaseOutFileName, optarg, FILE_NAME_LEN-1);
      {
	int iDum;
	for (iDum = 0; iDum < FILE_NAME_LEN && strBaseFileName[iDum];
	     iDum++) {}
	if (iDum == FILE_NAME_LEN) // No null; file name too long
	  vFatalError("Output file name too long", "command line");
      }
      qReadOutFileName = true;
      break;
    case 'O':
      sscanf(optarg, "%s", strOptParams);
      qParams = true;
      break;
    case 'q':
      sscanf(optarg, "%d", &iQualMeasure);
      break;
    case 'r':
      sscanf(optarg, "%lf", &dResolution);
      break;
    default:
      qError = true;
      break;
    }
  }

  if (qError || !qFile) {
    vUsage(strCanonicalName);
  }
  if (!qReadOutFileName) {
    strncpy(strBaseOutFileName, strBaseFileName, FILE_NAME_LEN-1);
  }    

  // Do this after option parsing so that the usage message doesn't have
  // other stuff around it.
  vGRUMMPInit(strCanonicalName);

  // Open log file and name quality file
  vOpenMessageFile(strBaseFileName);
  vMakeFileName(strQualFileName, "%s.qual", strBaseFileName,
		"main() [tetra.C]");  
  vMessage(1, "Mesh quality output file: %s\n", strQualFileName);

  VolMesh *pVMesh;
  if (TETRA == iExec) {
    // tetra initialization code

    // Read, validate, and optimize boundary rep
    Bdry3D *pB3D = new Bdry3D(strBaseFileName, iBdryOptLevel);
    
    // Generate a surface mesh from the boundary rep
    SurfMesh SM(pB3D, iQualMeasure);

    // SM.vMarkSmallAngles();
    vWriteFile_Surface(SM, strBaseOutFileName, ".surf");

    // Create an initial tetrahedralization of that surface mesh
    pVMesh = new VolMesh(SM, iQualMeasure, 1./dResolution, 1./dGrading);

    assert(pVMesh->qValid());
    vMessage(1, "Number of tets generated initially:  %d\n",
	     pVMesh->iNumCells());
  }
  else {
    // meshopt3d, trans3d and coarsen3d initialization code
    Bdry3D *pB3D = NULL;
    if (qReadAuxBdryFileName) {
      pB3D = new Bdry3D(strBdryFileName);
    }
    
    // Read volume mesh
    pVMesh = new VolMesh(strBaseFileName, iQualMeasure, dMaxAngle, pB3D);
    vMessage(1, "Number of cells read: %d\n", pVMesh->iNumCells());
  }

  if (qDoEdgeSwapping)
    pVMesh->vAllowEdgeSwapping();
  else
    pVMesh->vDisallowEdgeSwapping();

  if (qAllowBdryChanges)
    pVMesh->vAllowBdryChanges();
  else
    pVMesh->vDisallowBdryChanges();    

  pVMesh->vSetSmoothingThreshold(25.);

  switch (iExec) {
  case MESHOPT:
    dResolution /= sqrt(2.);
    if (!qParams) {
      sscanf("wsmf5mf5rmf5mf5", "%s", strOptParams);
      qParams = true;
    }
    pVMesh->vInitLengthScale(1./dResolution, 1./dGrading);
    if (qLengthScaleGiven) {
      // Read it and store it.
      FILE *pFLenFile = fopen(strLengthFileName, "r");
      if (NULL == pFLenFile)
	vFatalError("Couldn't open file for reading",
		    "length scale input");
      double dLen;
      int iV, iDataRead = fscanf(pFLenFile, "%d%lf", &iV, &dLen);
      int iNV = pVMesh->iNumVerts();
      while (iDataRead == 2) {
	if (iV < 0 || iV > iNV) {
	  vFatalError("Vertex index out of range",
		      "length scale input");
	}
	Vert *pV = pVMesh->pVVert(iV);
	pV->vSetLS(dLen);
	pVMesh->vMarkForLengthScaleCheck(pV);
	vMessage(4, "Set length scale for vert %d.\n", iV);
	iDataRead = fscanf(pFLenFile, "%d%lf", &iV, &dLen);
      };
      // Now use it.
      pVMesh->vUpdateLengthScale();
      pVMesh->vAdaptToLengthScale(eET);
    }
    pVMesh->vClassifyTets();
    pVMesh->vEvaluateQuality();
    break;
  case TETRA:
    {
      for (int iV = 0; iV < pVMesh->iNumVerts(); iV++)
	vMessage(4, "Vert: %3d @ (%.5G %.5G %.5G) Length scale: %.5G\n",
		 iV, pVMesh->pVVert(iV)->dX(), pVMesh->pVVert(iV)->dY(),
		 pVMesh->pVVert(iV)->dZ(), pVMesh->pVVert(iV)->dLS());
      // pVMesh->vProtectSmallAngles();
      pVMesh->vSetSwapType(eDelaunay);
      pVMesh->vSetEncroachmentType(eET);
      InsertionQueue IQ(pVMesh);
      IQ.vQualityRefine();
    }
    pVMesh->vClassifyTets();
    pVMesh->vEvaluateQuality();
    if (!qParams)
      sscanf("wsmf5", "%s", strOptParams);
    break;
  case COARSEN:
    // The constant in the length scale initialization is an empirical
    // constant, chosen to approximately match the edge-length
    // statistics and to get about a factor of 8 reduction in mesh size.
    vMessage(2, "Improving mesh quality before coarsening...\n");
    pVMesh->vSetSmoothingThreshold(20.);
    pVMesh->iSmooth(2);
    pVMesh->vSetSwapType(eMinSine);
    pVMesh->iSwap();
    pVMesh->vInitLengthScale(1.2/dResolution, 1./dGrading);
    pVMesh->vEvaluateQuality();
    // The following line leads to minor disasters in terms of (for
    //    instance) applying T22 swaps on the bdry.
    //    pVMesh->vSetStrictPatchChecking();

    // Set the initial optimization-based smoothing threshold fairly
    // low, since the mesh is likely to be fairly poor at that point.
    pVMesh->vSetSmoothingThreshold(10.);

    pVMesh->vCoarsen();
    vMessage(1, "Number of tets after coarsening: %d\n", pVMesh->iNumCells());

    pVMesh->vClassifyTets();
    pVMesh->vEvaluateQuality();

    if (!qParams) {
      sscanf("wsmf5mf5wsrmf5mf5", "%s", strOptParams);
      qParams = true;
    }
    break;
  case TRANSLATE:
    // Do nothing
    break;
  default:
    break;
  }
  
  if (iExec != TRANSLATE) {
    size_t iLen = strlen(strOptParams), i = 0;

    if (qDoEdgeSwapping)
      pVMesh->vAllowEdgeSwapping();
    else
      pVMesh->vDisallowEdgeSwapping();
    
    while (i < iLen) {
      char cType = strOptParams[i++];
      assert(cType);
      switch (cType) {
      case 'W':
      case 'w':  // Do some swapping
	{
	  char cMeasure = strOptParams[i++];
	  pVMesh->vAllowSwapRecursion();
	  switch (cMeasure) {
	  case 'I':
	  case 'i': // Insphere criterion
	    vMessage(1, "Doing Delaunay face swapping...\n");
	    pVMesh->vSetSwapType(eDelaunay);
	    break;
	  case 'M':
	  case 'm': // Minmax dihedral criterion
	    vMessage(1, "Doing face swapping to reduce maximum dihedral angle...\n"); 
	    pVMesh->vSetSwapType(eMaxDihed);
	    break;
	  case 'S':
	  case 's': // Maxmin sine criterion
	    vMessage(1, "Doing face swapping to increase minimum sine of dihedral angle...\n");
	    pVMesh->vSetSwapType(eMinSine);
	    break;
	  case 'q':
	  case 'Q':
	    vMessage(1, "Doing face swapping based on user-specified quality criterion...\n");
	    pVMesh->vSetSwapType(eFromQual);
	    break;
	  default:
	    vFatalError("Bad option for swapping criterion", "main()");
	  }
	  //	clock_t clBegin = clock();
	  int iNSwaps = pVMesh->iSwap();
	  //	clock_t clEnd = clock();
	  vMessage(1, "Number of swaps: %d\n", iNSwaps);
	  vMessage(1, "Number of tets after optimization:  %d\n",
		   pVMesh->iNumCells());
	  //	vMessage(1, "Number of clock ticks: %d\n", clEnd - clBegin);
	  break;
	}
      case 'M':
      case 'm': // Do some smoothing
	{
	  char cTech = strOptParams[i++];
	  int iFunc = atoi(&strOptParams[i++]);
	  pVMesh->vSetSmoothingTechnique(cTech);
	  if (iFunc == 0) iFunc = 10;
	  pVMesh->vSetSmoothingGoal(iFunc);
	  vMessage(1, "Smoothing vertices...\n");
	  pVMesh->iSmooth(1);
	}
	break;
      case 'N':
      case 'n':
	// No-op
	break;
      case 'R':
	pVMesh->iRepairBadCells(true);
	break;
      case 'r':
	pVMesh->iRepairBadCells(false);
	break;
      case 'B':
	pVMesh->iBreakBadBdryCells(true);
	break;
      case 'b':
	pVMesh->iBreakBadBdryCells(false);
	break;
      default:
	vFatalError("Bad option for mesh optimization operation", "main()");
      }
      pVMesh->vClassifyTets();
      pVMesh->vEvaluateQuality();
    }
  }

  pVMesh->vPurge();
  pVMesh->vWriteQualityToFile(strQualFileName);

  if ((MESHOPT != iExec) || qLengthScaleGiven) {
    pVMesh->vCheckEdgeLengths();
  }

  vMessage(1,"Writing volume mesh file...\n");
  pVMesh->vReorder();
  if (TETRA == iExec) 
    vWriteFile_VolMesh(*pVMesh, strBaseOutFileName);
  else
    vWriteFile_VolMesh(*pVMesh, strBaseOutFileName, ".out");

//   // Extra diagnostic stuff, to test the theory that bad tets cluster at
//   // the boundary.
//   for (int ii = 0; ii < pVMesh->iNumCells(); ii++) {
//     Cell *pC = pVMesh->pCCell(ii);
//     double dQual = pVMesh->dEvaluateQuality(pC);
//     vMessage(2, "Cell %5d.  Qual = %10.6f. Vert types: %2d %2d %2d %2d.\n",
// 	     ii, dQual, pC->pVVert(0)->iVertType(),
// 	     pC->pVVert(1)->iVertType(),
// 	     pC->pVVert(2)->iVertType(),
// 	     pC->pVVert(3)->iVertType());
//   }

  delete pVMesh;
  
  vMessage(1,"Done\n");

  return (0);
}
