#include "GR_config.h"
#include "GR_misc.h"
#include "GR_Geometry.h"
#include "GR_Vertex.h"
#include "GR_Face.h"
#include "GR_Vec.h"
#include "GR_BFace.h"
#include "GR_VolMesh.h"

void QuadFace::vNormal(double * const adNorm) const
     // Checks to see whether verts are in 2D or 3D.  For 2D, returns
     // with the normal (directed area) in adNorm[0] (this allows
     // vNormal to be called with the address of a double as its arg).
     // For 3D, returns a true vector normal.
{
  assert(qValid());
  assert(apVVerts[0]->iSpaceDimen() == 3);

  double adV1[3], adV2[3];
  ::vNormal(apVVerts[0]->adCoords(), apVVerts[1]->adCoords(),
	    apVVerts[2]->adCoords(), adV1);
  ::vNormal(apVVerts[0]->adCoords(), apVVerts[2]->adCoords(),
	    apVVerts[3]->adCoords(), adV2);
  adNorm[0] = (adV1[0] + adV2[0]) * 0.5;
  adNorm[1] = (adV1[1] + adV2[1]) * 0.5;
  adNorm[2] = (adV1[2] + adV2[2]) * 0.5;
}

void QuadFace::vUnitNormal(double * const adUnitNorm) const
     // Checks to see whether verts are in 2D or 3D.  For 2D, returns
     // with the unitNormal (directed area) in adUnitNorm[0] (this allows
     // vUnitNormal to be called with the address of a double as its arg).
     // For 3D, returns a true vector unitNormal.
{
  assert(qValid());
  assert(apVVerts[0]->iSpaceDimen() == 3);

  double adV1[3], adV2[3];
  ::vNormal(apVVerts[0]->adCoords(), apVVerts[1]->adCoords(),
	    apVVerts[2]->adCoords(), adV1);
  ::vNormal(apVVerts[0]->adCoords(), apVVerts[2]->adCoords(),
	    apVVerts[3]->adCoords(), adV2);
  adUnitNorm[0] = (adV1[0] + adV2[0]);
  adUnitNorm[1] = (adV1[1] + adV2[1]);
  adUnitNorm[2] = (adV1[2] + adV2[2]);
  vNORMALIZE3D(adUnitNorm);
}

double QuadFace::dSize() const
{
  double adTmp1[] = adDIFF3D(apVVerts[0]->adCoords(),
			     apVVerts[2]->adCoords());
  double adTmp2[] = adDIFF3D(apVVerts[1]->adCoords(),
			     apVVerts[3]->adCoords());
  double adTmp[3];
  vCROSS3D(adTmp1, adTmp2, adTmp);
  return (0.5*dMAG3D(adTmp));
}

int QuadFace::iFullCheck() const
{
  if (!Face::iFullCheck()) return 0;

  // Check whether face location is correct.
  assert(pCR->qValid() && pCL->qValid());
  switch (uiLoc) {
  case eBdryFace:
    if (pCR->eType() != Cell::eQuadBFace &&
	pCL->eType() != Cell::eQuadBFace) {
      return 0;
    }
    break;
  case eBdryTwoSide:
    if (pCR->eType() != Cell::eIntQuadBFace &&
	pCL->eType() != Cell::eIntQuadBFace) {
      return 0;
    }
    break;
  case eExterior:
  case eInterior:
    if ( (pCR->eType() != Cell::eHex &&
	  pCR->eType() != Cell::ePyr &&
	  pCR->eType() != Cell::ePrism)
	 ||
	 (pCL->eType() != Cell::eHex &&
	  pCL->eType() != Cell::ePyr &&
	  pCL->eType() != Cell::ePrism)
	 ||
	 pCL->iRegion() != pCR->iRegion()
	 ) {
      return 0;
    }
    break;
  default:
    // For all other face types, something is wrong.
    // Catches eUnknown, ePseudoSurface, and eInvalidFaceLoc.
    return 0;
  }
  // Exit is here rather than case-by-case so that other checks can be
  // added easily at a later date.
  return 1;
}

void TriFace::vNormal(double * const adNorm) const
     // Checks to see whether verts are in 2D or 3D.  For 2D, returns
     // with the normal (directed area) in adNorm[0] (this allows
     // vNormal to be called with the address of a double as its arg).
     // For 3D, returns a true vector normal.
{
  assert(qValid());
  assert(apVVerts[0]->iSpaceDimen() == 3);

  ::vNormal(apVVerts[0]->adCoords(), apVVerts[1]->adCoords(),
	    apVVerts[2]->adCoords(), adNorm);
  vSCALE3D(adNorm, 0.5);
}

void TriFace::vUnitNormal(double * const adUnitNorm) const
{
  assert(qValid());
  assert(apVVerts[0]->iSpaceDimen() == 3);

  ::vNormal(apVVerts[0]->adCoords(), apVVerts[1]->adCoords(),
	    apVVerts[2]->adCoords(), adUnitNorm);
  vNORMALIZE3D(adUnitNorm);
}

void TriFace::vCircumcenter(double adCircCent[]) const
{
  assert(qValid());
#ifndef NDEBUG
  double adA[3] = {0,0,0};
#endif
  double adB[3] = {pVVert(1)->dX() - pVVert(0)->dX(),
		   pVVert(1)->dY() - pVVert(0)->dY(),
		   pVVert(1)->dZ() - pVVert(0)->dZ()};
  double adC[3] = {pVVert(2)->dX() - pVVert(0)->dX(),
		   pVVert(2)->dY() - pVVert(0)->dY(),
		   pVVert(2)->dZ() - pVVert(0)->dZ()};


  double adRow1[3], adRow2[3], adRow3[3];

  adRow1[0] = adB[0];
  adRow1[1] = adB[1];
  adRow1[2] = adB[2];

  adRow2[0] = adC[0];
  adRow2[1] = adC[1];
  adRow2[2] = adC[2];

  vCROSS3D(adRow1, adRow2, adRow3);

  double dRHS1 = 0.5 * dDOT3D(adB, adB);
  double dRHS2 = 0.5 * dDOT3D(adC, adC);
  double dRHS3 = 0;

  vSolve3By3(adRow1, adRow2, adRow3, dRHS1, dRHS2, dRHS3, adCircCent);

#ifndef NDEBUG
  double dMagA = dDIST3D(adCircCent, adA);
  double dMagB = dDIST3D(adCircCent, adB);
  double dMagC = dDIST3D(adCircCent, adC);
  assert( 0 == iFuzzyComp(dMagA, dMagB));
  assert( 0 == iFuzzyComp(dMagA, dMagC));
#endif

  adCircCent[0] += pVVert(0)->dX();
  adCircCent[1] += pVVert(0)->dY();
  adCircCent[2] += pVVert(0)->dZ();
}

double TriFace::dCircumradius() const
{
  assert(qValid());
  double adCent[3];
  vCircumcenter(adCent);

  double dDistance = dDIST3D(pVVert(0)->adCoords(), adCent);
#ifndef NDEBUG
  double dDist1 = dDIST3D(pVVert(1)->adCoords(), adCent);
  double dDist2 = dDIST3D(pVVert(2)->adCoords(), adCent);
  assert( 0 == iFuzzyComp(dDistance, dDist1));
  assert( 0 == iFuzzyComp(dDistance, dDist2));
#endif
  return dDistance;
}

double TriFace::dSize() const
{
  double adTmp1[] = adDIFF3D(apVVerts[0]->adCoords(),
			     apVVerts[1]->adCoords());
  double adTmp2[] = adDIFF3D(apVVerts[0]->adCoords(),
			     apVVerts[2]->adCoords());
  double adTmp[3];
  vCROSS3D(adTmp1, adTmp2, adTmp);
  return (0.5*dMAG3D(adTmp));
}

enum eEncroachResult
TriBFaceBase::eIsEncroached(const enum eEncroachType eET) const
  // In a constrained Delaunay tetrahedralization, the point opposite
  // from a bdry face must encroach on it if any point is going to, so
  // check this point only.
{
  const Face *pF = pFFace();
  Cell *pC = pF->pCCellOpposite(this);
  assert(pC->eType() == Cell::eTet);
  Vert *pVOpp = dynamic_cast<TetCell*>(pC)->pVVertOpposite(pF);
  return eIsPointEncroaching(pVOpp->adCoords(), eET, false);
}

enum eEncroachResult
TriBFaceBase::eIsPointEncroachingBall(const double adPoint[3],
				      const bool qTieBreaker) const
  // This is the easy version of encroachment: a bdry face is encroached
  // iff a point lies within its circumsphere.
{
  const Face *pF = pFFace();
  double adCircCent[3];
  dynamic_cast<const TriFace*>(pF)->vCircumcenter(adCircCent);
  double dRadius = dDIST3D(adCircCent, pF->pVVert(0)->adCoords());
  double dDistance = dDIST3D(adCircCent, adPoint);
  switch (iFuzzyComp(dDistance, dRadius)) {
  case -1:
    // Distance smaller than radius -> inside
    return eSplit;
  case 1:
    // Distance larger than radius -> outside
    return eClean;
  case 0:
    return qTieBreaker ? eSplit : eClean;
  default:
    assert(0);
    return eClean;
  }
}

static enum eEncroachResult
eIsPointEncroachingGeneralLens(const double adPoint[3],
			       const bool qTieBreaker,
			       const TriBFaceBase * const pBTF,
			       const enum eEncroachType eET)
  // This version of encroachment is more complicated, because we need
  // to find the centers and radii of the spheres that make up the lens.
  // Then the point encroaches iff it's inside both spheres.  Also, if
  // 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 Face *pF = pBTF->pFFace();
  double adCircCent[3];
  dynamic_cast<const TriFace*>(pF)->vCircumcenter(adCircCent);
  double dFaceRadius = dDIST3D(adCircCent, pF->pVVert(0)->adCoords());
  double dSphereRadius = dFaceRadius*2 / sqrt(3.);
  double dEncDist = dSphereRadius / 2.;

  double adNorm[3], adCent1[3], adCent2[3];
  vUnitNormal(pF->pVVert(0)->adCoords(), pF->pVVert(1)->adCoords(),
	      pF->pVVert(2)->adCoords(), adNorm);
  adCent1[0] = adCircCent[0] + adNorm[0]*dEncDist;
  adCent1[1] = adCircCent[1] + adNorm[1]*dEncDist;
  adCent1[2] = adCircCent[2] + adNorm[2]*dEncDist;

  adCent2[0] = adCircCent[0] - adNorm[0]*dEncDist;
  adCent2[1] = adCircCent[1] - adNorm[1]*dEncDist;
  adCent2[2] = adCircCent[2] - adNorm[2]*dEncDist;

  double dDist1 = dDIST3D(adCent1, adPoint);
  double dDist2 = dDIST3D(adCent2, adPoint);

  int iFC1 = iFuzzyComp(dDist1, dSphereRadius);
  int iFC2 = iFuzzyComp(dDist2, dSphereRadius);

  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 = (qTieBreaker ? eSplit : eClean);

  // Note that if an internal bdry face has had offset insertion done in
  // one region, a request for offset insertion in the other region will
  // instead cause the face to be split.  This is a theoretical problem,
  // but not expected to be an issue in practice.
  if (eRetVal == eSplit && eET == eNewLens && pBTF->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 corners.  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 = dDIST3D(adPoint, pF->pVVert(0)->adCoords());
    double dDistB = dDIST3D(adPoint, pF->pVVert(1)->adCoords());
    double dDistC = dDIST3D(adPoint, pF->pVVert(2)->adCoords());
    // The following gives a provable angle bound of 30 degrees
    // (shortest edge and circumradius equal).

    int iFCA = iFuzzyComp(dDistA, dFaceRadius);
    int iFCB = iFuzzyComp(dDistB, dFaceRadius);
    int iFCC = iFuzzyComp(dDistC, dFaceRadius);
    if (iFCA == 1 && iFCB == 1 && iFCC == 1) {
      eRetVal = eOffsetPoint;
      pBTF->vRequestOffsetInsertion();

      double adDiff[] = adDIFF3D(adPoint, adCircCent);
      double dDot = dDOT3D(adNorm, adDiff);

      // Identify which apex is needed; it's the one in the interior, if
      // only one is.
      if (dDot < 0) {
	// Inside domain
	pBTF->vSetInsertionLocation(adCent2);
      }
      else if (pBTF->eType() == Cell::eIntTriBFace) {
	// Pick the closer apex
	assert(iFuzzyComp(dDist1, dDist2) != 0);
	if (dDist1 < dDist2) {
	  pBTF->vSetInsertionLocation(adCent1);
	}
	else {
	  pBTF->vSetInsertionLocation(adCent2);
	}
      }
      else {
	// Don't try to insert this thing outside the domain!
	pBTF->vClearOffsetInsertionRequest();
      }
    }
  }
  return eRetVal;
}

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

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

bool TriBFaceBase::qPointProjectsInside(const double adPointIn[3]) const
{
  // Find out where the given point projects onto the plane of the bdry
  // face.
  const Face *pF = pFFace();
  double adPoint[] = {adPointIn[0], adPointIn[1], adPointIn[2]};
  pF->vProjOntoFace(adPoint);

  // Now find the barycentric coordinates of the projection in the tet
  // opposite this bdry face
  assert(pF->pCCellOpposite(this)->eType() == Cell::eTet);
  TetCell *pTC = dynamic_cast<TetCell*>(pF->pCCellOpposite(this));
  double adBary[4];
  pTC->vBarycentrics(adPoint, adBary);

  // If all the barycentrics are zero or positive, great.  If one (or
  // more) is (are) negative, then the projection falls outside the
  // face.
  int aiBarySign[] = {iFuzzyComp(adBary[0], 0),
		      iFuzzyComp(adBary[1], 0),
		      iFuzzyComp(adBary[2], 0),
		      iFuzzyComp(adBary[3], 0)};
  int iNumNeg = ((aiBarySign[0] == -1 ? 1 : 0) +
		 (aiBarySign[1] == -1 ? 1 : 0) +
		 (aiBarySign[2] == -1 ? 1 : 0) +
		 (aiBarySign[3] == -1 ? 1 : 0));
  if (iNumNeg > 0) return false;

  // Also, projection onto a corner is worthless...
  int iNumZero = ((aiBarySign[0] == 0 ? 1 : 0) +
		  (aiBarySign[1] == 0 ? 1 : 0) +
		  (aiBarySign[2] == 0 ? 1 : 0) +
		  (aiBarySign[3] == 0 ? 1 : 0));
  if (iNumZero == 3) return false;
  return true;
}

bool TriFace::qIsLocallyDelaunay() const
  // Return true if you would -not- swap this face when doing Delaunay
  // face swapping.
{
  Vert *pVA, *pVB, *pVC, *pVD, *pVE, *pVPivot0, *pVPivot1, *pVOther;
  TetCell* apTCTets[4];
  int iNTets;
  eFaceCat eFC = eCategorizeFace(pVA, pVB, pVC, pVD, pVE, apTCTets, iNTets,
				 pVPivot0, pVPivot1, pVOther);
  switch (eFC) {
  case eT23:
  case eT32:
  case eT22:
  case eT44:
    {
      int iRes = iInsphere(pVA, pVB, pVC, pVD, pVE);
      if (eFC == eT32) iRes *= -1;
      return (iRes != 1);
    }
  default:
    return true;
  }
}

int TriFace::iFullCheck() const
{
  if (!Face::iFullCheck()) return 0;

  // Check whether face location is correct.
  assert(pCR->qValid() && pCL->qValid());
  switch (uiLoc) {
  case eBdryFace:
    if (pCR->eType() != Cell::eTriBFace &&
	pCL->eType() != Cell::eTriBFace) {
      return 0;
    }
    break;
  case eBdryTwoSide:
    if (pCR->eType() != Cell::eIntTriBFace &&
	pCL->eType() != Cell::eIntTriBFace) {
      return 0;
    }
    break;
  case eExterior:
  case eInterior:
    if ( (pCR->eType() != Cell::eTet &&
	  pCR->eType() != Cell::ePyr &&
	  pCR->eType() != Cell::ePrism)
	 ||
	 (pCL->eType() != Cell::eTet &&
	  pCL->eType() != Cell::ePyr &&
	  pCL->eType() != Cell::ePrism)
	 ||
	 pCL->iRegion() != pCR->iRegion()
	 ) {
      return 0;
    }
    break;
  default:
    // For all other face types, something is wrong.
    // Catches eUnknown, ePseudoSurface, and eInvalidFaceLoc.
    return 0;
  }
  // Exit is here rather than case-by-case so that other checks can be
  // added easily at a later date.
  return 1;
}
