#include "GR_config.h"
#include "GR_events.h"
#include "GR_BFace.h"
#include "GR_Face.h"
#include "GR_Geometry.h"
#include "GR_Mesh.h"
#include "GR_InsertionQueue.h"
#include "GR_WatsonInfo.h"

WatsonInfo::WatsonInfo(const InsertionQueueEntry& IQEIn,
		       const Mesh * const pMIn) :
    spCToRemove(), spBFEncroached(), spFHull(), spFToRemove(), IQE(IQEIn),
    pM(pMIn), vecpCNew(), vecpBFNew() 
{
  SUMAA_LOG_EVENT_BEGIN(WATSON_INFO);
  // Now, for each of the entry types, we need to gather all the Watson
  // insertion information.  The work is done elsewhere, with this
  // constructor basically acting as a wrapper.

  switch (IQE.eType()) {
  case InsertionQueueEntry::eTriCell:
  case InsertionQueueEntry::eTetCell:
    {
      // Identify the cell.
      SimplexCell *pSC = dynamic_cast<SimplexCell*>(IQE.pCCell());

      // First find the new point location
      pSC->vCircumcenter(adNewPoint);
      if (!pSC->qCircumscribes(adNewPoint)) {
	// Are there two faces that are bdry faces?
	assert(InsertionQueueEntry::eTetCell == IQE.eType());
	Face *apF[2];
	int iBF = 0;
	for (int iF = 0; iF < 4; iF++) {
	  Face *pF = pSC->pFFace(iF);
	  if (pF->qIsBdryFace()) 
	    apF[iBF++] = pF;
	}
	// There need to be exactly two bdry faces.
	assert(iBF == 2);
	Vert *pV0 = apF[0]->pVVert(0);
	Vert *pV1 = apF[0]->pVVert(1);
	if (!apF[1]->qHasVert(pV0))
	  pV0 = apF[0]->pVVert(2);
	else if (!apF[1]->qHasVert(pV1))
	  pV1 = apF[0]->pVVert(2);
	assert(apF[0]->qHasVert(pV0) && apF[1]->qHasVert(pV0));
	assert(apF[0]->qHasVert(pV1) && apF[1]->qHasVert(pV1));
	adNewPoint[0] = 0.5 * (pV0->dX() + pV1->dX());
	adNewPoint[1] = 0.5 * (pV0->dY() + pV1->dY());
	adNewPoint[2] = 0.5 * (pV0->dZ() + pV1->dZ());
	assert(pSC->qCircumscribes(adNewPoint));
      }

      // Now find the neighborhood info...
      pM->vGetWatsonData(adNewPoint, pSC, spCToRemove, spFHull,
			 spFToRemove, spBFEncroached);
      break;
    }
  case InsertionQueueEntry::eBdryEdge:
    {
      bool qSplitEdge = false;
      
      // Identify the BFace
      BFace *pBF = IQE.pBFBFace();

      // Find the attached face
      Face *pF = pBF->pFFace(0);
      
      if (pBF->qOffsetInsertionRequested()) {
	const double* adTmpLoc = pBF->adInsertionLocation();
	adNewPoint[0] = adTmpLoc[0];
	adNewPoint[1] = adTmpLoc[1];
      }
      else {
	(dynamic_cast<BdryEdgeBase*>(pBF))->vGetCCSPoint(adNewPoint);
	// Tag this explicitly to avoid circumstances in which this new
	// point falls in the offset insertion zone from causing bad
	// insertion (in the interior, too near the bdry).
	qSplitEdge = true;
      }

      // Find the real cell attached to the face
      Cell *pC = pF->pCCellOpposite(pBF);

      // Now find the neighborhood info...
      spBFEncroached.insert(pBF);
      pM->vGetWatsonData(adNewPoint, pC, spCToRemove, spFHull,
			 spFToRemove, spBFEncroached);

      assert(spBFEncroached.count(pBF) == 1);

      // If there are multiple encroachments, then just split the edge
      // instead.
      if (spBFEncroached.size() > 1) {
	(dynamic_cast<BdryEdgeBase*>(pBF))->vGetCCSPoint(adNewPoint);
	qSplitEdge = true;
	spCToRemove.clear();
	spFHull.clear();
	spFToRemove.clear();
	spBFEncroached.clear();

	spBFEncroached.insert(pBF);
	pM->vGetWatsonData(adNewPoint, pC, spCToRemove, spFHull,
			   spFToRemove, spBFEncroached);
      }
      if (qSplitEdge)
	pBF->vClearOffsetInsertionRequest();
      break;
    }
  case InsertionQueueEntry::eTriBFace:
    {
      bool qSplitFace = false;
      
      // Identify the BFace
      BFace *pBF = IQE.pBFBFace();

      // Find the attached face
      Face *pTF = pBF->pFFace(0);
      
      // Find the new point location
      if (pBF->qOffsetInsertionRequested()) {
	const double* adTmpLoc = pBF->adInsertionLocation();
	adNewPoint[0] = adTmpLoc[0];
	adNewPoint[1] = adTmpLoc[1];
	adNewPoint[2] = adTmpLoc[2];
      }
      else {
	pBF->vCircumcenter(adNewPoint);
	qSplitFace = true;
      }

      // Find the real cell attached to the face
      Cell *pC = pTF->pCCellOpposite(pBF);

      // Now find the neighborhood info...
      spBFEncroached.insert(pBF);
      pM->vGetWatsonData(adNewPoint, pC, spCToRemove, spFHull,
			 spFToRemove, spBFEncroached);

      if (spBFEncroached.size() > 1) {
	pBF->vCircumcenter(adNewPoint);
	qSplitFace = true;
	spCToRemove.clear();
	spFHull.clear();
	spFToRemove.clear();
	spBFEncroached.clear();

	spBFEncroached.insert(pBF);
	pM->vGetWatsonData(adNewPoint, pC, spCToRemove, spFHull,
			   spFToRemove, spBFEncroached);
      }
      if (qSplitFace)
	pBF->vClearOffsetInsertionRequest();
      // This assertion is needed to make sure that VolMesh insertion
      // will work properly.
      assert(spBFEncroached.count(pBF) == 1);

      break;
    }
  case InsertionQueueEntry::eSubSeg:
    {
      Vert *pV0 = IQE.pVVert(0);
      Vert *pV1 = IQE.pVVert(1);
      assert(pV0->iVertType() == Vert::eBdryApex ||
	     pV0->iVertType() == Vert::eBdryCurve);
      assert(pV1->iVertType() == Vert::eBdryApex ||
	     pV1->iVertType() == Vert::eBdryCurve);
      assert(pV0->iSpaceDimen() == 3);

      if (pV0->iVertType() == pV1->iVertType()) {
	// Regardless of whether these are both apexes or both bdry
	// curves, the subsegment should be split exactly in two.
	adNewPoint[0] = (pV0->dX() + pV1->dX()) / 2.;
	adNewPoint[1] = (pV0->dY() + pV1->dY()) / 2.;
	adNewPoint[2] = (pV0->dZ() + pV1->dZ()) / 2.;
      }
      else {
	// In this case, exactly one of the verts is a BdryApex.  To
	// make sure all edges with small incident surface angles get
	// split safely, split all such segments in power of 2 lengths,
	// using the biggest power of 2 less than 2/3 of current edge
	// length.  This will ensure a split between 1/3 and 2/3.
	assert((pV0->iVertType() == Vert::eBdryApex &&
		pV1->iVertType() == Vert::eBdryCurve) ||
	       (pV1->iVertType() == Vert::eBdryApex &&
		pV0->iVertType() == Vert::eBdryCurve));
	double adVec[] = adDIFF3D(pV1->adCoords(), pV0->adCoords());
	double dDistBetween = dMAG3D(adVec);
	vNORMALIZE3D(adVec);
	// Compute the log_e of 2/3 the length, then convert to log_2 by
	// multiplying by log_2(e) (equivalent to dividing by
	// log_e(2)).  Then take the largest integer below this.
	double dSplitPower = floor(log(dDistBetween/1.5) * M_LOG2E);
	double dSplitDist = pow(2., dSplitPower);
	assert(dDistBetween/3 <= dSplitDist &&
	       dSplitDist <= 2*dDistBetween/3);
	if (pV0->iVertType() == Vert::eBdryApex) {
	  // Split at dSplitDist from pV0.
	  adNewPoint[0] = pV0->dX() + dSplitDist * adVec[0];
	  adNewPoint[1] = pV0->dY() + dSplitDist * adVec[1];
	  adNewPoint[2] = pV0->dZ() + dSplitDist * adVec[2];
	}
	else {
	  // Split at dSplitDist from pV1.
	  adNewPoint[0] = pV1->dX() - dSplitDist * adVec[0];
	  adNewPoint[1] = pV1->dY() - dSplitDist * adVec[1];
	  adNewPoint[2] = pV1->dZ() - dSplitDist * adVec[2];
	}
      } // Done with unequal split. 
	

      // Now find an interior cell that contains both pV0 and pV1
      Cell *pC = pCInvalidCell;
      {
	std::set<Cell*> spCInc;
	{
	  std::set<Vert*> spVNeigh;
	  vNeighborhood(pV0, spCInc, spVNeigh);
	  assert(spVNeigh.count(pV1) == 1);
	}
	std::set<Cell*>::iterator iter;
	for (iter = spCInc.begin(); iter != spCInc.end(); iter++) {
	  pC = *iter;
	  assert(pC->qHasVert(pV0));
	  if (pC->qHasVert(pV1)) {
	    break;
	  }
	}
      }
      assert(pC->qHasVert(pV0));
      assert(pC->qHasVert(pV1));

      // Now find the neighborhood info...
      pM->vGetWatsonData(adNewPoint, pC, spCToRemove, spFHull,
			 spFToRemove, spBFEncroached);
      break;
    }
  case InsertionQueueEntry::eInvalid:
    // Should never get here.
    assert2(0, "IQE with no known entry type");
  }
  SUMAA_LOG_EVENT_END(WATSON_INFO);
}

int WatsonInfo::iAddNewEntitiesToQueue(InsertionQueue& IQ) const
  // After insertion, check to see whether any of the new entities are
  // still encroached.  Such encroachment can only (legally) occur when
  // initial bdry encroachment is being cleared up; afterwards, new
  // points that would encroach a bdry entity are never inserted at
  // all, unless the inserted point is also on a bdry.
{
  SUMAA_LOG_EVENT_BEGIN(WATSON_INFO);
  int i, iRetVal = 0;
  for (i = vecpCNew.size() - 1; i >= 0; i--) {
    if (!vecpCNew[i]->qDeleted()) {
      if (IQ.qAddEntry(InsertionQueueEntry(vecpCNew[i]))) {
	iRetVal ++;
      }
    }
  }

  // Add any new bdry faces or bdry segments that happen to be encroached.
  {
    std::set<BFace*> spBF(vecpBFNew.begin(), vecpBFNew.end());
    iRetVal += pM->iQueueEncroachedBdryEntities(IQ, spBF);
  }
  SUMAA_LOG_EVENT_END(WATSON_INFO);
  return iRetVal;
}

