#include "GR_Bdry2D.h"
#include "GR_BFace2D.h"
#include "GR_Cell.h"
#include "GR_Vec.h"

//--------------------------------------------------------
//
// New function for class BdryEdgeBase
// added by Charles Boivin, May 21, 1999
//
//--------------------------------------------------------


void BdryEdgeBase::vGetRatioPoint(double dRatio, double adNewLoc[]) const {
  // This returns the point located at "dRatio" between the first and
  // the second vertex, along the arclength...
  // Get the vertices
  const Vert *pV1 = pVVert(0);
  const Vert *pV2 = pVVert(1);
  double adLoc1[2], adLoc2[2];

	adLoc1[0] = pV1->dX();
	adLoc1[1] = pV1->dY();
	adLoc2[0] = pV2->dX();
	adLoc2[1] = pV2->dY();
	
	pB2D->pBdry(iBdryPatchIndex)->vRatioBetweenTwoPoints(dRatio, adLoc1, adLoc2, adNewLoc);
}

void BdryEdgeBase::vGetBisectionPoint(double adNewLoc[]) const {
  // This returns the point located halfway between pVBegin and
  // pVEnd, along the arclength...
  //    fprintf(stdout, "vGetBisectionPoint: Before the call\n");
  vGetRatioPoint(0.5, adNewLoc);
  //    fprintf(stdout, "vGetBisectionPoint: After the call\n");
}

void BdryEdgeBase::vGetCCSPoint(double adNewLoc[]) const {
  // This returns the point where the segment should be split when using
  // the "Concentric Circular Shells" method, as described in Ruppert's
  // paper "A Delaunay Refinement Algorithm".
  //
  // Basically:
  //   - Return midpoint if both verts are part of the original boundaries
  //   - Return midpoint if both verts are NOT original boundaries
  //   - If one vert is part of the original boundaries:
  //       - Return point at d' = D*2**k'
  //         where - (int) k' = log2(d/D)
  //               - d = 1/2 length of segment
  //               - D = some constant (0.01, for example)

  // Get the vertices
  const Vert *pV1 = pVVert(0);
  const Vert *pV2 = pVVert(1);
  // Const
  double D = 1e-10;
  int    k = 0;
  double dk;
  double dl = dSize();
  double dp;
  bool   qV1, qV2;

//	fprintf(stdout, "Entering CCS\n");

//	fprintf(stdout, "   pV1 type: %d\n", pV1->iVertType());
//	fprintf(stdout, "   pV2 type: %d\n", pV2->iVertType());
//	fprintf(stdout, "   Length: %f\n", dl);

	// If the segment is not connected to a small angle, then there is
	// no need to do anything fancy. Return midpoint
	
	qV1 = pV1->qIsSmallAngleCorner();
	qV2 = pV2->qIsSmallAngleCorner();
	
	if (qV1 || qV2) {
		// There was a small angle. Which one?
		double dAL;		
		BdryPatch2D *pPatch =
		  dynamic_cast<BdryPatch2D*>(pPatchPointer());
		
		// If it's both, cut at midpoint!
		if (qV1 && qV2) {
			vGetRatioPoint(0.5, adNewLoc);
			return;
		}
		
		if (qV1) {
			// pV1 is small		
			dAL = pPatch->dParameterAtPoint(pV1->adCoords());			
		}
		else {
			// It was pV2
			assert(qV2);
			dAL = pPatch->dParameterAtPoint(pV2->adCoords());			
		}

		// Get the distance...
		dk = log(dl / (2.0*D))/log(2.0);
		k = static_cast<int>(dk+0.5);
		dp = D*pow(2.0, k);		
		
		vMessage(4, "Distance is %f\n", dp);
				
		// Distance would be from end or from start?
		if (dAL < (pPatch->dMaxParameter()/2.0)) {
			// From start
			pPatch->vPointAtDistFromStart(dp, adNewLoc);
			return;
		}
		else {
			// From end...
			pPatch->vPointAtDistFromEnd(dp, adNewLoc);
			return;		
		}		
	}
	else {
	 	// No small angles, return midpoint
	 	vGetRatioPoint(0.5, adNewLoc);
	 	return;
	}
	
/*
  if (pV1->iVertType() == Vert::eBdryApex) {
    if (pV2->iVertType() == Vert::eBdryApex) {
      // Both are eBdryApex
      // Return midpoint
      // fprintf(stdout, "   Both are eBdryApex -- bisect\n");
      vGetRatioPoint(0.5, adNewLoc);
    } else {
      // Only one is... pV1
      dk = log(dl / (2.0*D))/log(2.0);
      k = (int)dk;
      if ( (dk - k) > 0.5) k++;
      dp = D*pow(2.0, k);
      //        fprintf(stdout, "   Only pV1 is eBdryApex -- k = %d, dk = %f, dp = %f, dp/dl = %f\n",
      //  	      k, dk, dp, dp/dl);
      vGetRatioPoint(dp/dl, adNewLoc);
    }
  } else {
    if (pV2->iVertType() == Vert::eBdryApex) {
      // Only one is ... pV2
      dk = log(dl / (2.0*D))/log(2.0);
      k = (int)dk;
      if ( (dk - k) > 0.5) k++;
      dp = D*pow(2.0, k);
      // Invert the order, because it's pV2
      // fprintf(stdout, "   Only pV2 is eBdryApex -- k = %d, dk = %f, dp = %f, (dl-dp)/dl = %f\n",
      //	 k, dk, dp, (dl-dp)/dl);
      vGetRatioPoint((dl-dp)/dl, adNewLoc);				
    } else {
      // fprintf(stdout, "   No verts are eBdryApex\n");
      // No verts are BdryApex
      // return midpoint
      vGetRatioPoint(0.5, adNewLoc);
    }
  }
*/

}

enum eEncroachResult BdryEdgeBase::eIsPointEncroachingBall(const double adPoint[],
							   const bool qTieBreak) const {
  // Checks if the point adPoint lies in the circle with diameter defined
  // by pV0 & pV1
  double dX0 = pVVert(0)->dX();
  double dY0 = pVVert(0)->dY();
  double dX1 = pVVert(1)->dX();
  double dY1 = pVVert(1)->dY();

  // Find radius
  double dRadius = dSize() * 0.5;

  // Find the center..
  double adCent[] = {(dX0 + dX1) * 0.5, (dY0 + dY1) * 0.5};

    // Find the distance from center..
  double dDist = dDIST2D(adPoint, adCent);
  switch (iFuzzyComp(dDist, dRadius)) {
  case -1:
    // It encroaches
    return eSplit;
  case 1:
    // It does not
    return eClean;
  case 0:
    // On the border
    if (qTieBreak) return eSplit;
    else           return eClean;
  default:
    assert(0);
    return eClean;
  }					     
}

static enum eEncroachResult
eIsPointEncroachingGeneralLens(const double adPoint[], const bool qTieBreak,
			       const BdryEdgeBase * const pBE,
			       const enum eEncroachType eET)
{
  // Does standard lens checking.  If this is encroached and the
  // encroachment type eET == eNewLens, check to see whether the point
  // happens to hit the jackpot and requires an offset point to be
  // inserted instead of merely splitting the edge.
  const double *adLoc0 = pBE->pVVert(0)->adCoords();
  const double *adLoc1 = pBE->pVVert(1)->adCoords();

    // Find radius
  double dRadius = pBE->dSize() / sqrt(3.0);
  double dDist   = dRadius / 2.0;
  // Find midpoint..
  double adC[2], adC1[2], adC2[2];

  adC[0] = (adLoc0[0] + adLoc1[0]) * 0.5;
  adC[1] = (adLoc0[1] + adLoc1[1]) * 0.5;

  // Get normal..
  double adNormal[2];
  pBE->vUnitNormal(adNormal);

  // Bad idea to use the patch because it can be a CURVE!
//  pPatchPointer()->vNormalAtPoint(adC, adNormal);

    // Find the two centers..
  adC1[0] = adC[0] + adNormal[0]*dDist;
  adC1[1] = adC[1] + adNormal[1]*dDist;
  adC2[0] = adC[0] - adNormal[0]*dDist;
  adC2[1] = adC[1] - adNormal[1]*dDist;

  assert(iFuzzyComp(dDIST2D(adC1, adLoc0), dRadius) == 0);
  assert(iFuzzyComp(dDIST2D(adC1, adLoc1), dRadius) == 0);
  assert(iFuzzyComp(dDIST2D(adC2, adLoc0), dRadius) == 0);
  assert(iFuzzyComp(dDIST2D(adC2, adLoc1), dRadius) == 0);
  assert(iFuzzyComp(dDIST2D(adC1, adC2), dRadius) == 0);

    // Find the distances from centers..
  double dDist1 = dDIST2D(adPoint, adC1);
  double dDist2 = dDIST2D(adPoint, adC2);

  // If both distances are smaller than radius, it encroaches
  int iFC1 = iFuzzyComp(dDist1, dRadius);
  int iFC2 = iFuzzyComp(dDist2, dRadius);

  enum eEncroachResult eRetVal = eClean;
  // If both are -1, then it definitely encroaches
  if (iFC1 == -1 && iFC2 == -1)
    eRetVal = eSplit;
  // If one is negative and one zero, then use the tiebreaker
  if ((iFC1 == -1 && iFC2 ==  0) ||
      (iFC1 ==  0 && iFC2 == -1))
    eRetVal = (qTieBreak ? eSplit : eClean);

  // Note that if an internal bdry edge has had offset insertion done in
  // one region, a request for offset insertion in the other region will
  // instead cause the edge to be split.  This is a theoretical problem,
  // but not expected to be an issue in practice.
  if (eRetVal == eSplit && eET == eNewLens && pBE->qOffsetInsertionAllowed()) {
    // Now check whether to do the offset point thing.  The offset point
    // is inserted iff the proposed new point lies too far (farther than
    // dRadius) from each of the endpoints.  Also, if the offset point
    // lies outside the domain, it isn't inserted.  In principle, the
    // offset point also shouldn't be inserted if it's across an
    // internal bdry from the cell that triggered the insertion.  In
    // practice, once encroached edges are fixed, either the bdry edge
    // with the incorrectly added offset pt will have been split -or-
    // the original triggering cell will still be there.  On the second
    // try, the BFace will be split.

    double dDistA = dDIST2D(adPoint, adLoc0);
    double dDistB = dDIST2D(adPoint, adLoc1);
    // The following gives a provable angle bound of 30 degrees
    // (shortest edge and circumradius equal).
    double dRadius2 = 0.5 * dDIST2D(adLoc0, adLoc1);
    
    double adDiff[] = adDIFF2D(adPoint, adC);
    double dDot = dDOT2D(adNormal, adDiff);

    int iFCA = iFuzzyComp(dDistA, dRadius2);
    int iFCB = iFuzzyComp(dDistB, dRadius2);
    if (iFCA == 1 && iFCB == 1) {
      eRetVal = eOffsetPoint;
      pBE->vRequestOffsetInsertion();
      // Identify which apex is needed; it's the one in the interior, if
      // only one is.
      if (dDot < 0) {
	// Inside domain
	pBE->vSetInsertionLocation(adC2);
      }
      else if (pBE->eType() == Cell::eIntBdryEdge) {
	// Pick the closer apex
	assert(iFuzzyComp(dDist1, dDist2) != 0);
	if (dDist1 < dDist2) {
	  pBE->vSetInsertionLocation(adC1);
	}
	else {
	  pBE->vSetInsertionLocation(adC2);
	}
      }
      else {
	// Don't try to insert this thing outside the domain!
	pBE->vClearOffsetInsertionRequest();
      }
    }
  }
  return eRetVal;
}  

enum eEncroachResult
BdryEdgeBase::eIsPointEncroachingLens(const double adPoint[],
				      const bool qTieBreak) const {
  return eIsPointEncroachingGeneralLens(adPoint, qTieBreak, this, eLens);
}

enum eEncroachResult
BdryEdgeBase::eIsPointEncroachingNewLens(const double adPoint[],
					 const bool qTieBreak) const {
  return eIsPointEncroachingGeneralLens(adPoint, qTieBreak, this, eNewLens);
}

enum eEncroachResult BdryEdgeBase::eIsEncroached(const enum eEncroachType eET) const
{
  // Works for both regular and internal bdry edges.
  for (int iF = 0; iF < iNumFaces(); iF++) {
    const Face *pF = pFFace(iF);
    Cell *pC = pF->pCCellOpposite(this);
    // This cell could be a quad and we can still tell whether the bdry
    // edge is encroached on.
    for (int iV = 0; iV < pC->iNumVerts(); iV++) {
      Vert *pV = pC->pVVert(iV);
      if (!qHasVert(pV)) {
	enum eEncroachResult eER = eIsPointEncroaching(pV->adCoords(), eET,
						       false);
        if (eER != eClean) return eER;
      }
    }
  }
  return eClean;
}
  
void BdryEdgeBase::
vGetOffsetInsertionPoint(const double adTriggerLoc[],
			 double adOffsetLoc[]) const
{
  // Need to return the lens apex on the same side of the bdry edge as
  // the trigger point.

  // Edge endpoint locations
  const double *adLoc0 = pVVert(0)->adCoords();
  const double *adLoc1 = pVVert(1)->adCoords();

    // Find radius of lens arcs
  double dRadius = dSize() / sqrt(3.0);
  double dDist   = dRadius / 2.0;

  // Find midpoint..
  double adC[2];
  adC[0] = (adLoc0[0] + adLoc1[0]) * 0.5;
  adC[1] = (adLoc0[1] + adLoc1[1]) * 0.5;

  // Get normal..
  double adNormal[2];
  vUnitNormal(adNormal);

  // To construct the lens apex that lies on the same side as the
  // trigger point, need a dot product of the normal with the vector
  // from edge center to trigger point.
  double adDiff[] = adDIFF2D(adTriggerLoc, adC);
  double dDot = dDOT2D(adNormal, adDiff);

  if (dDot < 0) {
    // Inside domain
    adOffsetLoc[0] = adC[0] - adNormal[0]*dDist;
    adOffsetLoc[1] = adC[1] - adNormal[1]*dDist;
  }
  else if (eType() == eIntBdryEdge) {
    // Also inside domain, although possibly on the side opposite the
    // triggering cell.  That's not likely to be a problem in practice.
    adOffsetLoc[0] = adC[0] + adNormal[0]*dDist;
    adOffsetLoc[1] = adC[1] + adNormal[1]*dDist;
  }
  else {
    // Don't try to insert this thing outside the domain!
    vClearOffsetInsertionRequest();
    adOffsetLoc[0] = adC[0];
    adOffsetLoc[1] = adC[1];
  }

  assert(iFuzzyComp(dDIST2D(adLoc0, adOffsetLoc), dRadius) == 0);
  assert(iFuzzyComp(dDIST2D(adLoc1, adOffsetLoc), dRadius) == 0);
}
