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

#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "GR_config.h"
#include "GR_AdaptPred.h"
#include "GR_Bdry2D.h"
#include "GR_InsertionQueue.h"
#include "GR_Mesh2D.h"
#include "GR_events.h"
#include "GR_misc.h"
#include "GR_ADT.h"
#include <vector>
#include <set>

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

typedef struct _refinementInfo_ {
  double dX, dY, dLen;
} RefinementInfo;

static char GR_tri_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]",
  "    -c    Use diametral lenses instead of diametral circles [1]",
  "          This will push the minimum angle to 30 degrees (instead of 20)",
  "    -j    Coarsen quasi-structured mesh anisotropically [0]",
  "  ",
  "  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]",
  " ",
  "    -q #  Mesh quality measure used [2]",
  "         0  Maximum angle",
  "         1  Minimum angle",
  "         2  All angles",
  "         6  Aspect ratio I (inscribed to circum radius)",
  "         9  Aspect ratio II (based on area versus perimeter)",
  "        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_tri_Usage[i]) > 0; i++) {
    vMessage(0, "%s\n", GR_tri_Usage[i]);
  }
  exit(1);
}

int main(int iNArg, char *apcArgs[]) {
  char strBaseFileName[FILE_NAME_LEN], strQualFileName[FILE_NAME_LEN];
  char strBdryFileName[FILE_NAME_LEN]; // Used when reading both a mesh
				       // file and a bdry file for
				       // refinement. 
  char strBaseOutFileName[FILE_NAME_LEN];
  bool qReadOutFileName = false;
  bool qReadAuxBdryFileName = false;
  bool qLengthScaleGiven = false;
  char strLengthFileName[FILE_NAME_LEN];
  double dScale = 1.;
  double dGrading = 1.;
  enum eEncroachType eET = eLens;

  char strExecName[FILE_NAME_LEN];
  char strCanonicalName[20];
  enum {UNKNOWN, TRI, MESHOPT, COARSEN, HEAL} iExec = UNKNOWN;

  qAdaptPred = false;

  sprintf(strExecName, "%s", apcArgs[0]);
  if (strstr(strExecName, "coarsen2d")) {
    iExec = COARSEN;
    sprintf(strCanonicalName, "coarsen2d");
  }
  else if (strstr(strExecName, "meshopt2d")) {
    iExec = MESHOPT;
    sprintf(strCanonicalName, "meshopt2d");
  }
  else if (strstr(strExecName, "tri")) {
    iExec = TRI;
    sprintf(strCanonicalName, "tri");
  }
  else if (strstr(strExecName, "heal2d")) {
    iExec = HEAL;
    sprintf(strCanonicalName, "heal2d");
  }
  else {
    vFatalError("Invalid executable name.  Be less creative renaming files.\n",
		"tri.cxx:main()");
  }

  vGRUMMPInit(strCanonicalName);

  int iQualMeasure = 2, iPreviousPasses = 0;
  bool qAllowBdryChanges = true;
  bool qAnisoCoarsening = false;

  int iPChar, iError = 1;
  while ((iPChar = getopt(iNArg, apcArgs, "A:b:B:c:g:i:j:l:o:q:r:")) != EOF) {
    switch (iPChar) {
    case 'A':
      qAdaptPred = (optarg[0] == '1');
      break;
    case 'b':
      qAllowBdryChanges = (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 'c':
      switch (optarg[0]) {
      case '0':
	eET = eBall;
	break;
      case '1':
      default:
	eET = eLens;
	break;
      case '2':
	eET = eNewLens;
	break;
      }
      break;
    case 'g':
      sscanf(optarg, "%lf", &dGrading);
      if (dGrading < 1) dGrading = 1;
      break;
    case 'i':
      strncpy(strBaseFileName, 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("Input file name too long", "command line");
      }
      iError--;
      break;
    case 'j':
      qAnisoCoarsening = (optarg[0] == '1');
      break;
    case 'l':
      if (iExec == MESHOPT || iExec == TRI) {
	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 {
	iError += 10; // Be sure to crash out.
      }
      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 'q':
      {
	int iTmp;
	sscanf(optarg, "%d", &iTmp);
	if (iTmp == 0 || iTmp == 1 || iTmp == 2 || iTmp == 6 || iTmp == 9
	    || iTmp == 10)
	  iQualMeasure = iTmp;
	break;
      }
    case 'r':
      sscanf(optarg, "%lf", &dScale);
      break;
    default:
      iError += 10;
      break;
    }
  }
  if (iError)
    vUsage(strCanonicalName);

  // File name manipulation
  vMakeFileName(strQualFileName, "%s.qual",
		strBaseFileName, "main() [tri.C]");
  if (!qReadOutFileName) {
    strncpy(strBaseOutFileName, strBaseFileName, FILE_NAME_LEN-1);
  }

  // Open input and log files
  vOpenMessageFile(strBaseFileName);
  vMessage(1, "Mesh quality output file: %s\n", strQualFileName);

  Mesh2D *pM2D;
  // Read and validate 2D mesh
  if (iExec == TRI) {
    // Reading a bdry file and creating a mesh from that
    Bdry2D B2D(strBaseFileName, false);
    pM2D = new Mesh2D(B2D, iQualMeasure);

    vMessage(2, "Number of points in bdry file: %d\n",
	     B2D.iNumPoints());
    vMessage(2, "Number of verts in initial triangulation:  %d\n",
	     pM2D->iNumVerts());
    vMessage(1, "Number of triangles in initial triangulation:  %d\n",
	     pM2D->iNumCells());
  }
  else { // coarsen2d or meshopt2d
    Bdry2D *pB2D = NULL;
    if (qReadAuxBdryFileName) {
      pB2D = new Bdry2D(strBdryFileName);
    }
    // Reading an existing mesh file
    pM2D = new Mesh2D(strBaseFileName, iQualMeasure, pB2D);
    vMessage(1, "Number of cells read: %d\n", pM2D->iNumCells());
    if (pB2D != NULL) delete pB2D;
  }
  assert(pM2D->qValid());

  pM2D->vEvaluateQuality();

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

  int iNumPasses = 0;
  switch (iExec) {
  case MESHOPT:
    {
      pM2D->vInitLengthScale(1./dScale, 1./dGrading);
      pM2D->vSetSwapType(eDelaunay);
      int iNSwaps;
//       clock_t clBegin = clock();
      do {
	iNSwaps = pM2D->iSwap();
	vMessage(1, "Swapped %d times.\n", iNSwaps);
      } while (iNSwaps != 0);
//       clock_t clEnd = clock();
//       vMessage(1, "Elapsed time %d clocks at %d per second.\n",
// 	       clEnd - clBegin, CLOCKS_PER_SEC);
      pM2D->vEvaluateQuality();
      
      pM2D->vCheckForSmallAngles();
      pM2D->vSetLSChecking(true);
      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 = pM2D->iNumVerts();
	while (iDataRead == 2) {
	  if (iV < 0 || iV > iNV) {
	    vFatalError("Vertex index out of range",
			"length scale input");
	  }
	  Vert *pV = pM2D->pVVert(iV);
	  pV->vSetLS(dLen);
	  pM2D->vMarkForLengthScaleCheck(pV);
	  vMessage(4, "Set length scale for vert %d.\n", iV);
	  iDataRead = fscanf(pFLenFile, "%d%lf", &iV, &dLen);
	};
	pM2D->vUpdateLengthScale();
	pM2D->vCheckEdgeLengths();
	pM2D->vAdaptToLengthScale(eET);
      }
      break;
    }
  case TRI:
    {
      pM2D->vCheckForSmallAngles();

      pM2D->vSetEncroachmentType(eET);
      pM2D->vCheckForLargeAngles();
      pM2D->vInitLengthScale(1.0/dScale, 1.0/dGrading);
      pM2D->vSetLSChecking(true);

      InsertionQueue IQ(pM2D);
      IQ.vQualityRefine();

      vMessage(1,"Number of triangles after initial refinement:  %d\n",
	       pM2D->iNumCells());

      // Do it again, because our Shewchuk implementation doesn't queue
      // tri's that are created when removing verts, with the result
      // that some get missed in the final mesh.
      InsertionQueue IQ2(pM2D);
      IQ2.vQualityRefine();

      vMessage(1,"Number of triangles after initial refinement:  %d\n",
	       pM2D->iNumCells());

      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");

	
	RefinementInfo RefInfo;
	std::vector<RefinementInfo> vecRefinementInfo;
	int iDataRead;

	while(1) {
	  iDataRead = fscanf(pFLenFile, "%lf%lf%lf",
			     &(RefInfo.dX), &(RefInfo.dY), &(RefInfo.dLen));
	  if (iDataRead != 3) break;
	  vecRefinementInfo.push_back(RefInfo);
	}
	fclose(pFLenFile);
	vMessage(1,"Read Length Scale data for %zd locations\n", vecRefinementInfo.size());

	int iNumNeedMoreRefinement;
	do {
	  iNumNeedMoreRefinement = 0;

	  vMessage(1, "Smoothing vertices...\n");
	  pM2D->vSetSmoothingThreshold(40.);
	  pM2D->iSmooth(1);
 	  pM2D->iSwap();
	  // 	  pM2D->iSmooth(1);

	  std::set<int> setAlreadyRefined;
	  ADT VTree(pM2D->pVAGetVerts());
	  std::vector<int> veciQueryResult;
	  for (std::vector<RefinementInfo>::iterator Iter = vecRefinementInfo.begin(); 
	       Iter != vecRefinementInfo.end();
	       Iter++) {
	    double dEps = (Iter->dLen)/2.;
	    do {
	      double adBB[] = {Iter->dX - dEps, Iter->dX + dEps,
			       Iter->dY - dEps, Iter->dY + dEps};
	      veciQueryResult = VTree.veciRangeQuery(adBB);
	      dEps *= 1.5;
	    } while (veciQueryResult.empty());
	    for (int i = veciQueryResult.size() - 1; i >= 0; i--) {
	      int iV = veciQueryResult[i];
	      if (setAlreadyRefined.find(iV) != setAlreadyRefined.end()) {
		continue;
	      }
	      setAlreadyRefined.insert(iV);
	      Vert *pV = pM2D->pVVert(iV);
	      if (pV->dLS() > Iter->dLen) {
		double dNewLS = Iter->dLen;
 		if (dNewLS < (pV->dLS()/1.5)) { // reduce it and flag it as needing more refinement
 		  dNewLS = pV->dLS()/1.5;
 		  iNumNeedMoreRefinement++;
 		}
		vMessage(2,"Updating Length Scale of Vert %d was %f now %f\n",
			 veciQueryResult[i], pV->dLS(), dNewLS);
		pV->vSetLS(dNewLS);
		pM2D->vMarkForLengthScaleCheck(pV);
	      }
	    }
	  }
	  pM2D->vUpdateLengthScale();
	  InsertionQueue IQ3(pM2D);
	  IQ3.vQualityRefine();
	  vMessage(2,"Number of triangles after refinement step:  %d\n",
		   pM2D->iNumCells());
	  vMessage(2,"Number of Verts not meeting final length scale:  %d\n", iNumNeedMoreRefinement);

	} while (iNumNeedMoreRefinement != 0);
      }

      pM2D->vEvaluateQuality();

      // Smoothing parameters for later
      pM2D->vSetSmoothingThreshold(40.);
//       iNumPasses = 1;
    }
    vMessage(1,"Number of triangles after refinement:  %d\n",
	     pM2D->iNumCells());
    break;
  case COARSEN:
    pM2D->vInitLengthScale(1.0/dScale, 1.0/dGrading);
    pM2D->vSetSmoothingThreshold(25.);
    if (qAnisoCoarsening) {
      pM2D->vSetInteriorSkip(3);
      pM2D->vSetBoundaryFill(true);
    }
    else {
      pM2D->vSetInteriorSkip(1);
      pM2D->vSetBoundaryFill(false);
    }
    pM2D->vCoarsen();
    vMessage(1, "Number of cells after coarsening: %d\n", pM2D->iNumCells());

    pM2D->vSetSwapType(eMaxDihed);
    pM2D->iSwap();
    pM2D->vEvaluateQuality();

    iNumPasses = 5;
    iPreviousPasses = 2;
    break;
  case HEAL:
    pM2D->vSetLengthScale(7.5);

//     {
//       // Total hack job for simulating mesh adaptation.
//       for (int ii = 0; ii < pM2D->iNumVerts(); ii++) {
//	Vert *pV = pM2D->pVVert(ii);
//	double dX = pV->dX();
//	double dY = pV->dY();
//	double dX_star = 250 - (dY-541)*(dY-541) / 130;
//	double dLen;
//	if (dX < dX_star) {
//	  dLen = 1. / (0.125 + 10*exp(-0.5*(dX_star - dX)));
//	}
//	else {
//	  dLen = 1. / (0.05 + 10*exp(-(dX-dX_star)/3.));
//	}
//	pV->vSetLS(dLen);
//       }
//     }

    vMessage(1, "Read mesh, preparing to heal it\n");
    pM2D->vCheckEdgeLengths();

    pM2D->vAdaptToLengthScale(eET);

    pM2D->vEvaluateQuality();
    pM2D->vCheckEdgeLengths();
    // vWriteFile_Mesh2D(*pM2D, strBaseOutFileName, ".insert");

    vMessage(1, "Smoothing\n");
    // Smoothing parameters for later
    pM2D->vSetSmoothingThreshold(40.);
    iNumPasses = 4;
    break;
  default:
    assert(0);
    break;
  }

  if (iNumPasses > 0) {
    // Do some smoothing if requested
    vMessage(1, "Smoothing vertices...\n");
    for (int ii = iPreviousPasses; ii < iNumPasses; ii++) {
      vMessage(2, "Pass %d of %d...\n", ii+1, iNumPasses);
      pM2D->iSmooth(1);
      if ((ii+1)%3 == 0) {
	// After a few smoothing passes, try swapping to see if that
	// helps any.
	pM2D->iSwap();
      }
    }
    pM2D->vEvaluateQuality();
  }

//   // Hack for ubercoarse meshes
//   pM2D->vSetSmoothingThreshold(25.);
//   pM2D->vSetInteriorSkip(1);
//   pM2D->vSetBoundaryFill(false);
//   pM2D->vCoarsen();
//   vMessage(1, "Number of cells after coarsening: %d\n", pM2D->iNumCells());

//   pM2D->vSetSwapType(eMaxDihed);
//   pM2D->iSwap();
//   pM2D->iSmooth(1);

    // Quality output
  pM2D->vWriteQualityToFile(strQualFileName);
  // Mesh output
  vMessage(1,"Writing mesh to file(s)...\n");
  if (iExec == TRI) {
    pM2D->vReorder();
    vWriteNative(*pM2D, strBaseOutFileName, true);
//     vWriteFEA(*pM2D, strBaseOutFileName, 5);
  }
  else {
    //    pM2D->vReorder();
    vWriteNative(*pM2D, strBaseOutFileName, true, ".out");
//     vWriteFEA(*pM2D, strBaseOutFileName, 5, ".out");
  }

  vMessage(1,"Done\n");

//   pM2D->vCheckEdgeLengths();

  assert(pM2D != NULL);
  delete pM2D;

  return (0);
}
