#include <math.h>
#include "GR_assert.h"
#include "GR_events.h"
#include "GR_BFace.h"
#include "GR_Geometry.h"
#include "GR_List.h"
#include "GR_Mesh2D.h"
#include "GR_Vec.h"

#ifdef IRIX6
#include <values.h>
#endif

/*
   ------------------------------------------------

   This file includes the functions necessary for
   Ruppert's algorithm in 2D.

   ------------------------------------------------
 */

void Mesh2D::
vSetLSChecking (bool qCheck)
{

  if (qCheck == true) {
    // Should mark all cells as potentially bad again..

    for (int i = 0; i < iNumCells (); i++) {
      pCCell (i)->vMarkIllSized ();
    }
  }
}


static double
dReturnAngleAtMiddle (const Vert * pV1, const Vert * pV0, const Vert * pV2)
{
  // calculates and returns the angle formed by pV1->pV0->pV2
  // This code creates
  //   a) two vectors normal to the sides
  //   b) a vector joining the two end vectors
  // It uses those vectors to determine if the angle is smaller or greater than 180 degrees..

  double a2dV1[2], a2dV2[2], dLen1, dLen2;
  // Create vectors..
  // Normal to pV0->pV1, pointing to the "right" of pV0->pV1
  a2dV1[0] = pV1->dX () - pV0->dX ();
  a2dV1[1] = pV1->dY () - pV0->dY ();
  dLen1 = sqrt (a2dV1[0] * a2dV1[0] + a2dV1[1] * a2dV1[1]);
  assert (dLen1 > 0.0);
  a2dV2[0] = pV2->dX () - pV0->dX ();
  a2dV2[1] = pV2->dY () - pV0->dY ();
  dLen2 = sqrt (a2dV2[0] * a2dV2[0] + a2dV2[1] * a2dV2[1]);
  assert (dLen2 > 0.0);

  // Calculate dot product
  double dDotResult = dDOT2D (a2dV1, a2dV2);
  // Divide by lengths
  dDotResult /= (dLen1 * dLen2);
  // dDotResult now has cos(angle)
  return acos (dDotResult);
}

void Mesh2D::
vCheckForLargeAngles (const double dMaxAngle)
{


  SUMAA_LOG_EVENT_BEGIN (RUPPERT_INIT_LARGE);

  int i = -1;
  bool qFixedSomething = true;

  vMessage(1, "Checking for large angles..\n");

  // This is because, at 90 degrees, the circumcenter is
  // ON the boundary. Anything greater will make the circumcenter
  // fall beyond the boundary.. And you don't want that.

  // Now check for edges encroaching on verts already..

  while (qFixedSomething) {
    qFixedSomething = false;
    i = -1;

    while (i < iNumBdryFaces () - 1) {
      i++;

//      vMessage(3, "Looking at face %d -- Max is %d\n", i, iNumBdryFaces());

      Vert *pV0, *pV1, *pVOpp;
      BFace *pBF = pBFBFace (i);

      if (!pBF->qDeleted ()) {

	Face *pF = pBF->pFFace (0);
	Cell *pC = pF->pCCellOpposite (pBF);
	if (pC->eType() == Cell::eQuadCell) continue;

	pVOpp = dynamic_cast<TriCell*>( pC)->pVVertOpposite (pF);
	pV0 = pF->pVVert (0);
	pV1 = pF->pVVert (1);

	assert (pV0 != pVOpp);
	assert (pV1 != pVOpp);

	vMessage(4, "pV0: (%f, %f), pV1 (%f, %f), pVOpp (%f, %f)\n",
		 pV0->dX(), pV0->dY(), pV1->dX(), pV1->dY(), pVOpp->dX(), pVOpp->dY());

	double dAngle = dReturnAngleAtMiddle (pV0, pVOpp, pV1);

	vMessage(3, "Angle found here: %f\n", dAngle);

	if (dAngle > dMaxAngle) {
	  // The angle is greater than the maximum allowed
	  // Must split the edge.

	  qFixedSomething = true;

	  vMessage(4, "   Splitting - %d = %d...\n", i, iNumBdryFaces ());
	  vMessage(4, "   Angle was %f\n", dAngle);
	  vMessage(4, "   pV0 (%f, %f), pV1 (%f, %f), pVOpp (%f, %f)\n",
		   pV0->dX (), pV0->dY (), pV1->dX (), pV1->dY (), pVOpp->dX (), pVOpp->dY ());

	  double adMid[2];
	  dynamic_cast<BdryEdgeBase*>( pBF)->vGetBisectionPoint (adMid);
//				  vMessage(3, "   Midpoint is (%f, %f)\n", adMid[0], adMid[1]);
	  // Find the face index...
	  int j = -1;
	  bool qFaceFound = false;

	  while (!qFaceFound) {
	    j++;
	    assert (j < 3);
	    if (pC->pFFace (j) == pF)
	      qFaceFound = true;
	  }

	  // Insert the point on the boundary...
	  List < BdryEdgeBase * >LpBEList;
	  iInsertPointOnBdryWatson (pC, j, adMid, LpBEList, true);
	  if (qLengthScaleInitialized) {
	    vSetLengthScale(pVVert(iNumVerts()-1));
	  }

	  // Reset the counter.. start from beginning..
	  i = -1;
	}				// end if angle > 90.0

      }				// end if the BFace was not deleted

    }				// end while i < iNumBdryFaces()

    // Now do the same thing for internal boundary faces...

    i = -1;
//    vMessage(3, "Now looking at internal boundary faces\n");

    while (i < iNumIntBdryFaces () - 1) {
      i++;

//      vMessage(3, "Looking at intbface %d -- Max is %d\n", i, iNumIntBdryFaces () - 1);

      Vert *pV0, *pV1, *pVOpp;
      IntBdryEdge *pBF = pIBEIntBdryEdge (i);

      if (!pBF->qDeleted ()) {

	// must check both sides..
	for (int k = 0; k < 2; k++) {
	  Face *pF = pFInvalidFace;
	  if (k == 0) {
	    pF = pBF->pFFace (0);
	  }
	  else {
	    pF = pBF->pFFace (1);
	  }
	  Cell *pC = pF->pCCellOpposite (pBF);
	  if (pC->eType() == Cell::eQuadCell) continue;

	  pVOpp = dynamic_cast<TriCell*>( pC)->pVVertOpposite (pF);
	  pV0 = pF->pVVert (0);
	  pV1 = pF->pVVert (1);

	  assert (pV0 != pVOpp);
	  assert (pV1 != pVOpp);

	  double dAngle = dReturnAngleAtMiddle (pV0, pVOpp, pV1);

//	  vMessage(3, "Angle found here: %f\n", dAngle * 180.0 / M_PI);

	  if (dAngle > dMaxAngle) {
	    // The angle is greater than the maximum allowed
	    // Must split the edge.

	    qFixedSomething = true;

	    vMessage(4, "   Splitting - %d = %d...\n", i, iNumBdryFaces ());
	    //        fprintf(stdout, "   Angle was %f\n", dAngle);
	    //        fprintf(stdout, "   pV0 (%f, %f), pV1 (%f, %f), pVOpp (%f, %f)\n",
	    //                pV0->dX(), pV0->dY(), pV1->dX(), pV1->dY(), pVOpp->dX(), pVOpp->dY());

	    double adMid[2];
	    dynamic_cast<BdryEdgeBase*>( pBF)->vGetBisectionPoint (adMid);
//	    vMessage(3, "   Midpoint is (%f, %f)\n", adMid[0], adMid[1]);
	    // Find the face index...
	    int j = -1;
	    bool qFaceFound = false;

	    while (!qFaceFound) {
	      j++;
	      assert (j < 3);
	      if (pC->pFFace (j) == pF)
		qFaceFound = true;
	    }

	    // Insert the point on the boundary...
	    List < BdryEdgeBase * >LpBEList;
	    iInsertPointOnBdryWatson (pC, j, adMid, LpBEList, true);
	    if (qLengthScaleInitialized) {
	      vSetLengthScale(pVVert(iNumVerts()-1));
	    }

	    // Reset the counter.. start from beginning..
	    i = -1;
	  }			// end if angle > 90.0

	}				// checking both sides..

      }				// end if the BFace was not deleted

    }				// end while i < iNumBdryFaces()

  }// while there are no more large angles to inspect ANYWHERE
  // (both internal and external bfaces..


  if (qLengthScaleInitialized) {
    vUpdateLengthScale();
  }

  vPurge();
  // If you don't clear this, you run the risk of having out-of-range
  // pointers, or at least pointers that aren't to the verts you
  // thought.  And InitLengthScale gets called after anyway, so it makes
  // to difference if there's anything in here...
  spVUpdateLS.clear();
  vMessage(2, "Done checking for large angles\n");

  assert (qValid ());

  SUMAA_LOG_EVENT_END (RUPPERT_INIT_LARGE);
}

void Mesh2D::
vCheckForSmallAngles ()
{
  //
  // This functions checks for small angles in corners of the original
  // geometry and marks the vertices in those corners when it finds some.
  // It also isolates that corner by inserting at the shortest edge's midpoint
  // distance, on both segments
  //

  SUMAA_LOG_EVENT_BEGIN (RUPPERT_INIT_SMALL);

  int i;

  vMessage(1, "Looking for small angles in the geometry\n");

  // We already know where there are small angles, from the B2D object...
  for (i = 0; i < B2D.iNumPoints(); i++) {
    if (B2D.qIsBdryPoint(i)) {
      if (B2D.qVertSmallAngle(i)) {
	vMessage(2, "Point %d has a small angle\n", i);
	// Now find the vert that is located there
	double adPt[2];
	bool   qFoundPt = false;
	B2D.vGetPoint(i, adPt);
	int j = -1;
	Vert *pV;

	do {
	  j++;
	  pV = pVVert(j);
	  if ((!pV->qDeleted()) && (pV->qIsSmallAngleCorner())) {
	    if (iFuzzyComp(dDIST2D(pV->adCoords(), adPt), 0) == 0) {
	      // We have the right vert!
	      qFoundPt = true;
	    }
	  }
	} while (!qFoundPt);

	vMessage(2, "Point is located at (%f, %f)\n", pV->dX(), pV->dY());

	// We know it will have B2D.iNumPatchesConnectedToPoint(i)
	// BFaces around it
	List < Face * >LpFNext;
	List < Vert * >LpVNext;

	vGetNeighborVertsAndFaces(pV, LpVNext, LpFNext);

	int iF;
	for (iF = 0; iF < LpFNext.iLength(); iF++) {
	  Face *pF = LpFNext[iF];
	  Cell *pCL = pF->pCCellLeft();
	  Cell *pCR = pF->pCCellRight();

	  bool qBdryL, qBdryR;
	  qBdryL = ((pCL->eType() == Cell::eBdryEdge) || (pCL->eType() == Cell::eIntBdryEdge));
	  qBdryR = ((pCR->eType() == Cell::eBdryEdge) || (pCR->eType() == Cell::eIntBdryEdge));

	  if ( (!qBdryL) && (!qBdryR) ) {
	    // delete it from list
	    LpFNext.vDeleteByIndex(iF);
	    LpVNext.vDeleteByIndex(iF);
	    iF--;
	  }
	} // for loop over neighbor faces

	assert (LpFNext.iLength() == B2D.iNumPatchesConnectedToPoint(i));

	// now find the smallest and next smallest distance of verts to this
	// point...

	double dSmallestDist = 1e10;
	double dNextSmallest = 1e10;
	for (int k = 0; k < LpFNext.iLength(); k++) {
	  // Find smallest distance
	  double dDistance = dDIST2D(LpVNext[k]->adCoords(), pV->adCoords());
	  if (dDistance < dSmallestDist) {
	    dNextSmallest = dSmallestDist;
	    dSmallestDist = dDistance;
	  }
	  else {
	    if (dDistance < dNextSmallest) {
	      dNextSmallest = dDistance;
	    }
	  }
	} // for loop over all neighbor faces

	vMessage(2, "Smallest distance found to that corner is %f\n", dSmallestDist);
	vMessage(2, "Will split at half that distance...\n");

				// We want only half of this distance...
	dSmallestDist /= 2.0;


	// Ok. Now need to insert points on each BFace.
	// Good news is we already have access to the patch pointers!

	List <BdryPatch2D*> LpBP2DDone;
	LpBP2DDone.vDisallowSorting();
	List <double> LdLocs;
	double adMid[2];
	LpBP2DDone.vClear();
	List <int> LiWhichLeft;

	for (j = 0; j < B2D.iNumPatchesConnectedToPoint(i); j++) {
	  BdryPatch2D *pPatch = B2D.pB2DPatchConnectedToPoint(i, j);
	  // Must figure out from which end to find point...
	  if (pPatch->iBeginPointIndex() == pPatch->iEndPointIndex()) {
	    // Special treatment for closed loops
	    if (LpBP2DDone.qInList(pPatch)) {
	      pPatch->vPointAtDistFromEnd(dSmallestDist, adMid);
	      // add the coords to the list
	      LdLocs.vAppendItem(adMid[0]);
	      LdLocs.vAppendItem(adMid[1]);
	    }
	    else {
	      pPatch->vPointAtDistFromStart(dSmallestDist, adMid);
	      // add the coords to the list
	      LdLocs.vAppendItem(adMid[0]);
	      LdLocs.vAppendItem(adMid[1]);
	      // add patch to list
	      LpBP2DDone.vAddItem(pPatch);
	    }
	  } // for special loops
	  else {
	    if (pPatch->iBeginPointIndex() == i) {
	      pPatch->vPointAtDistFromStart(dSmallestDist, adMid);
	      // add the coords to the list
	      LdLocs.vAppendItem(adMid[0]);
	      LdLocs.vAppendItem(adMid[1]);
	      // add patch to list
	      LpBP2DDone.vAddItem(pPatch);
	    }
	    else {
	      assert(pPatch->iEndPointIndex() == i);

	      pPatch->vPointAtDistFromEnd(dSmallestDist, adMid);
	      // add the coords to the list
	      LdLocs.vAppendItem(adMid[0]);
	      LdLocs.vAppendItem(adMid[1]);
	      // add patch to list
	      LpBP2DDone.vAddItem(pPatch);
	    }
	  } // if the patch is not a closed loop
	  // Now figure out which edge was connected to that patch
	  int iSmallestTotal = -1;
	  for (iF = 0; iF < LpFNext.iLength(); iF++) {
	    Face *pF = LpFNext[iF];
	    Cell *pCL = pF->pCCellLeft();
	    Cell *pCR = pF->pCCellRight();
	    BdryEdgeBase *pBE;
	    if ( (pCL->eType() == Cell::eBdryEdge) ||
		 (pCL->eType() == Cell::eIntBdryEdge) ) {
	      // it's pCL
	      pBE = dynamic_cast<BdryEdgeBase*>(pCL);
	    }
	    else {
	      // it's pCR
	      pBE = dynamic_cast<BdryEdgeBase*>(pCR);
	    }
	    BdryPatch2D *pEdgePatch =
	      dynamic_cast<BdryPatch2D*>(pBE->pPatchPointer());
	    if (pEdgePatch == pPatch) {
	      iSmallestTotal = iF;
	    }
	  }
	  LiWhichLeft.vAppendItem(iSmallestTotal);
	} // for loop over patches connected to the point

	vMessage(2, "Here are the points that will be inserted to fix the small angle\n");
	int iLoc;
	for (iLoc = 0; iLoc < LdLocs.iLength()/2; iLoc++) {
	  vMessage(2, "Edge %d: (%f, %f)\n",
		   LiWhichLeft[iLoc],
		   LdLocs[iLoc*2], LdLocs[(iLoc*2) + 1]);
	}

	// Must figure out in which order to insert them...
	List <BdryEdgeBase *> LpBEdge;
	for (iF = 0; iF < LpFNext.iLength(); iF++) {
	  Face *pF = LpFNext[iF];
	  Cell *pCL = pF->pCCellLeft();
	  Cell *pCR = pF->pCCellRight();

	  bool qBdryL, qBdryR;
	  qBdryL = ((pCL->eType() == Cell::eBdryEdge) || (pCL->eType() == Cell::eIntBdryEdge));
	  qBdryR = ((pCR->eType() == Cell::eBdryEdge) || (pCR->eType() == Cell::eIntBdryEdge));

	  if (qBdryL) {
	    LpBEdge.vAppendItem(dynamic_cast<BdryEdgeBase*>(pCL));
	  }

	  if (qBdryR) {
	    LpBEdge.vAppendItem(dynamic_cast<BdryEdgeBase*>(pCR));
	  }
	}

	// Now have a list of points to insert...
	int iWhich = 0;
	while (LiWhichLeft.iLength() > 0) {
	  // Dummy list...
	  iLoc = LiWhichLeft[iWhich];
	  List < BdryEdgeBase * >LpEBFDummy;
	  adMid[0] = LdLocs[2*iLoc];
	  adMid[1] = LdLocs[(2*iLoc)+1];
	  Vert VDum;
	  VDum.vSetCoords(2, adMid);
	  // Find the cell, and the face index..
	  BdryEdgeBase *pBE = LpBEdge[iWhich];
	  Face *pF = pBE->pFFace(0);
	  Cell *pC = pF->pCCellOpposite(pBE);
	  // Find the face index now..
	  int iFaceIdx = 0;
	  while (pC->pFFace(iFaceIdx) != pF) iFaceIdx++;
	  // Check if the point is in the exterior.. don't fix it.
	  bool qExterior = qEnsurePointIsInterior(LpBEdge, iWhich, &VDum, false);
	  if (!qExterior) {
	    vMessage(2, "Inserting at (%f, %f)...\n", adMid[0], adMid[1]);
	    iInsertPointOnBdryWatson(pC, iFaceIdx, adMid, LpEBFDummy, false);
	    if (qLengthScaleInitialized) {
	      vSetLengthScale(pVVert(iNumVerts()-1));
	    }
	    LiWhichLeft.vDelete(iLoc);
	    LpBEdge.vDelete(pBE);
	  }
	  else {
	    vMessage(2, "Skipping Edge %d because it's outside..\n", iWhich);
	  }
	  iWhich++;
	  if (iWhich >= LiWhichLeft.iLength()) {
	    iWhich = 0;
	  }
	} // loop over list of points to insert

      } // small angle

    } // boundary point

  } // all points

  //  vUpdateLengthScale();

  vMessage(2, "Done looking for small angles\n");

  assert (qValid ());

  SUMAA_LOG_EVENT_END (RUPPERT_INIT_SMALL);
}
