/***************************************************************************
			  InsertWatson2D.cxx
			     -------------------
    begin                : Thu Jul 15 1999
    copyright            : (C) 1999 by Charles Boivin
    email                : cboivin@mech.ubc.ca

    Heavily modified for better compatibility with InsertionQueue by
    Carl Ollivier-Gooch, October, 2003.

    Also, modified to allow offset edge insertion (and therefore a
    better angle bound) by Carl Ollivier-Gooch, February, 2004.
 ***************************************************************************/

#include <algorithm>
#include <queue>

#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"
#include "GR_WatsonInfo.h"

void Mesh2D::vGetWatsonData(const double adNewPoint[3],
			    Cell* const pCSeed,
			    std::set<Cell*>& spCSearch,
			    std::set<Face*>& spFHull,
			    std::set<Face*>& spFToRemove,
			    std::set<BFace*>& spBFEncroached) const
  // The cell, face and bdry face listings from this routine are only
  // guaranteed to be correct for interior vertices.  However, boundary
  // vertices will be identifiable as such, because spBFEncroached will
  // not be empty on return.  In this case, normal encroachment handling
  // (i.e., calling iInsertPointOnBdryWatson appropriately) will
  // straighten out the situation.  Ordinarily, the InsertionQueue code
  // will call this routine only when the vertex is in the interior, so
  // no problems will arise.
{
  static int iHullsChecked = 0, iTotalHullSize = 0;
  // This code used to be part of Mesh2D::iInsertPointWatson(Cell *
  // pCInitialCell, const double adNewpoint[]).

  assert(pCSeed->iFullCheck());
  assert(pCSeed->eType() == Cell::eTriCell);

  spCSearch.clear();
  spFHull.clear();
  spFToRemove.clear();
  spBFEncroached.clear();

  // Define the new vertex
  Vert VDummy;
  VDummy.vSetCoords (2, adNewPoint);

  // Make sure point is really in the circumcircle of the given cell
  assert (iIncircle (pCSeed->pVVert (0),
		     pCSeed->pVVert (1),
		     pCSeed->pVVert (2),
		     &VDummy) >= 0
	  );

  // First item on the list is the initial cell..
  spCSearch.insert (pCSeed);
  std::queue<Cell*> qpC;
  qpC.push(pCSeed);

  // This list keeps track of the BFaces encroached on by the new vertex
  std::set<BdryEdge*> spBEEncroachedBFaces;

  // Must now determine a list of cells in the mesh with the new point
  // in their circumcircle

  SUMAA_LOG_EVENT_BEGIN (WATSON_NEIGHBORS);

  int iF;
  for (iF = 0; iF < 3; iF++) {
    Face *pF = pCSeed->pFFace(iF);
    spCSearch.insert(pF->pCCellOpposite(pCSeed));
    qpC.push(pF->pCCellOpposite(pCSeed));
  }

  int iInd = 1;
  while (!qpC.empty()) {
    // Go through the list and find the cells with the point in their
    // circumcircle 
    Cell *pCCand = qpC.front();
    qpC.pop();
    if (pCCand->eType() == Cell::eTriCell) {
      // For triangles:
      if (iIncircle(pCCand->pVVert(0), pCCand->pVVert(1),
		    pCCand->pVVert(2), &VDummy) != -1) {
	// If new point is in or on circumcircle, add neighbors to search list
	for (iF = 0; iF < 3; iF++) {
	  Face *pF = pCCand->pFFace(iF);
	  Cell * pCTmp = pF->pCCellOpposite(pCCand);
	  if (spCSearch.count(pCTmp) == 0) {
	    spCSearch.insert(pCTmp);
	    qpC.push(pCTmp);
	  }
	}
	iInd++; // Go on to the next candidate.
      }
      else {
	// Point not in circumcircle; delete from the list
	spCSearch.erase(pCCand);
      }
    } // Done with triangles
    else if (pCCand->eType() == Cell::eQuadCell ||
	     pCCand->qIsQuasiStructured()) {
      // .clear and return is also an alternative here, resulting in no
      // insertion, but this is quite restrictive.
      spCSearch.erase(pCCand);
    }
    else {
      assert(pCCand->eType() == Cell::eBdryEdge ||
	     pCCand->eType() == Cell::eIntBdryEdge);
      BdryEdgeBase* pBECand = dynamic_cast<BdryEdgeBase*> (pCCand);
      // For bdry faces:
      // If new point encroaches, add bface to encroached list.  In this
      // case, we need to tag the BFace as encroached in case of ties.
      enum eEncroachResult eER =
	pBECand->eIsPointEncroaching(adNewPoint, eEncroachmentType(), true);
      if (eER != eClean) {
	spBFEncroached.insert(pBECand);
      }
      // Delete bface from search list
      spCSearch.erase(pCCand);
      // If the vertex falls on an internal bdry edge, add the cells
      // adjacent to the BFace on each side; this is needed to check for
      // other encroached faces.
      if (pBECand->eType() == Cell::eIntBdryEdge &&
	  !pBECand->pPatchPointer()->qDisjointFrom(&VDummy)) {
	Cell *pCTmp = pBECand->pFFace(0)->pCCellOpposite(pCCand);
	if (spCSearch.count(pCTmp) == 0) {
	  spCSearch.insert(pCTmp);
	  qpC.push(pCTmp);
	}
	pCTmp = pBECand->pFFace(1)->pCCellOpposite(pCCand);
	if (spCSearch.count(pCTmp) == 0) {
	  spCSearch.insert(pCTmp);
	  qpC.push(pCTmp);
	}
      }
    }
  } // Queue is now empty!

  // Ensure that there will be no offset insertion when multiple BFaces
  // are encroached, as this can lead to odd exterior insertion cases.
  switch (spBFEncroached.size()) {
  case 0:
    break;
  case 1:
    {
      BFace *pBECand = *(spBFEncroached.begin());
      if (pBECand->eType() == Cell::eIntBdryEdge &&
	  pBECand->qOffsetInsertionRequested()) {
	// For this case, make certain that there won't be insertion in
	// the wrong region, as this will create inverted cells and other
	// disasters.  The "wrongness" test makes sure the new point and
	// the centroid of the cell that started it all are on the same
	// side of an internal bdry edge.
	double adCent[2];
	pCSeed->vCentroid(adCent);
	if (iOrient2D(pBECand->pVVert(0)->adCoords(),
		      pBECand->pVVert(1)->adCoords(), adCent) !=
	    iOrient2D(pBECand->pVVert(0)->adCoords(),
		      pBECand->pVVert(1)->adCoords(), adNewPoint)) {
	  // Rig things so that this internal bdry edge will be split.
	  double adTmp[2];
	  (dynamic_cast<BdryEdgeBase*>(pBECand))->vGetCCSPoint(adTmp);
	  pBECand->vSetInsertionLocation(adTmp);
	  pBECand->vForbidOffsetInsertion();
	}
      }
    }
    break;
  default:
    std::set<BFace*>::iterator iterBF;
    for (iterBF = spBFEncroached.begin(); iterBF != spBFEncroached.end();
	 iterBF++) {
      (*iterBF)->vForbidOffsetInsertion();
    }
    break;
  }

  std::multiset<Face*> mspFTemp;
  std::set<Cell*>::iterator iterC;
  for (iterC = spCSearch.begin(); iterC != spCSearch.end(); iterC++) {
    Cell *pC = *iterC;
    for (iF = 0; iF < 3; iF ++) {
      mspFTemp.insert(pC->pFFace(iF));
    }
  }

  // The following loop contains some extra insert calls, but so what?
  std::multiset<Face*>::iterator iterF;
  for (iterF = mspFTemp.begin(); iterF != mspFTemp.end(); iterF++) {
    Face *pF = *iterF;
    int iCount = mspFTemp.count(pF);
    assert(iCount == 1 || iCount == 2);
    if (iCount == 1) spFHull.insert(pF);
    else             spFToRemove.insert(pF);
  }      

  vMessage(4, "Watson: Done sorting faces into hull and interior.\n");
  vMessage(4, "Watson: # cells to delete: %zd\n", spCSearch.size());
  vMessage(4, "Watson: # faces to delete: %zd\n", spFToRemove.size());
  vMessage(4, "Watson: # faces in hull: %zd\n", spFHull.size());
  iTotalHullSize += spCSearch.size();
  iHullsChecked ++;
  vMessage(3, "Watson: Total cells deleted: %5d in %5d tries.\n",
	   iTotalHullSize, iHullsChecked);
  assert(3*spCSearch.size() - 2*spFToRemove.size() ==
	 spFHull.size());
  SUMAA_LOG_EVENT_END (WATSON_NEIGHBORS);
}

int Mesh2D::iInsertPointWatson(WatsonInfo& WI)
{
  int iTriCellsBefore = iNumTriCells();
  int iBFacesBefore = iNumBdryFaces();
  int iVertsBefore = iNumVerts();
  int iVertsInserted;

  List<Cell*> LpCToRemove;
  List<Face*> LpFToRemove;
  List<Face*> LpFHull;
  List<BFace*> LpBFEncroached;

  // Copy data into the archaic lists.  This horrible hack will
  // eventually vanish anyway when Serge cleans up the Watson insertion
  // code.
  {
    std::set<Cell*>::iterator iterC;
    for (iterC = WI.spCToRemove.begin(); iterC != WI.spCToRemove.end(); iterC++)
      LpCToRemove.vAppendItem(*iterC);

    std::set<Face*>::iterator iterF;
    for (iterF = WI.spFToRemove.begin(); iterF != WI.spFToRemove.end(); iterF++)
      LpFToRemove.vAppendItem(*iterF);

    for (iterF = WI.spFHull.begin(); iterF != WI.spFHull.end(); iterF++)
      LpFHull.vAppendItem(*iterF);
  }

  switch (WI.eType()) {
  case InsertionQueueEntry::eTriCell:
    {
      Vert *pV = pVNewVert();
      pV->vSetCoords(2, WI.adCoords());
      vInsertWatsonInterior(pV, LpCToRemove, LpFToRemove, LpFHull);
      iVertsInserted = 1;
      break;
    }
  case InsertionQueueEntry::eBdryEdge:
    {
      BFace *pBF = WI.IQE.pBFBFace();
      LpBFEncroached.vAppendItem(pBF);

      std::set<BFace*>::iterator iterBF;
      for (iterBF = WI.spBFEncroached.begin(); iterBF != WI.spBFEncroached.end();
	   iterBF++) {
	BFace *pBFTmp = *iterBF;
	if (pBFTmp == pBF) continue;
	LpBFEncroached.vAppendItem(pBFTmp);
      }

      if (pBF->qOffsetInsertionRequested()  && pBF->qOffsetInsertionAllowed()) {
	Vert *pV = pVNewVert();
	pV->vSetCoords(2, WI.adCoords());
	vMessage(3, "Attempting offset point insertion at (%f, %f).\n",
		 pV->dX(), pV->dY());
	// First, delete all verts inside the circle with center at pV
	// passing through the endpoints of the edge.
	int iNDeleted = 0;
	bool qSuccess = qDeleteInteriorPointsInBall
	  (pBF, pV->adCoords(),
	   dDIST2D(pV->adCoords(), pBF->pVVert(0)->adCoords()),
	   iNDeleted);
	assert(qSuccess);
	vMessage(3, "Deleted %d points to clear ball for offset bdry insertion.\n",
		 iNDeleted);
	if (0 == iNDeleted) {
	  // If nothing was deleted, just go ahead and insert the point.
	  vInsertWatsonInterior(pV, LpCToRemove, LpFToRemove,
				LpFHull);
	}
	else {
	  // Need to re-compute the neighborhood of this point, because
	  // the old one just got trashed.  Fortunately, the interior
	  // cell opposite the BFace that we're trying to insert for is
	  // guaranteed to circumscribe the new vertex.
	  Face *pFBdry = pBF->pFFace();
	  Cell *pCSeed = pFBdry->pCCellOpposite(pBF);
	  assert(iIncircle(pCSeed->pVVert(0), pCSeed->pVVert(1),
			   pCSeed->pVVert(2), pV) == 1);
	  std::set<BFace*> spBFDummy;
	  vGetWatsonData(pV->adCoords(), pCSeed, WI.spCToRemove,
			 WI.spFHull, WI.spFToRemove, spBFDummy);
	  {
	    std::set<Cell*>::iterator iterC;
	    for (iterC = WI.spCToRemove.begin(); iterC != WI.spCToRemove.end();
		 iterC++)
	      LpCToRemove.vAppendItem(*iterC);
	    
	    std::set<Face*>::iterator iterF;
	    for (iterF = WI.spFToRemove.begin(); iterF != WI.spFToRemove.end();
		 iterF++)
	      LpFToRemove.vAppendItem(*iterF);
	    
	    for (iterF = WI.spFHull.begin(); iterF != WI.spFHull.end(); iterF++)
	      LpFHull.vAppendItem(*iterF);
	  }
	  assert(int(spBFDummy.size()) == LpBFEncroached.iLength());
	  // The encroached BFaces will the same as they were the last
	  // time we checked for this point, and they've already been
	  // queued, so ignore them.
	  vInsertWatsonInterior(pV, LpCToRemove, LpFToRemove, LpFHull);
	}
	iVertsInserted = 1;
	// The following line looks like a bug at this point.  CFO-G 9/1/04
//	pBF->vForbidOffsetInsertion();
      }
      else {
	// Split the edge instead.
	Face *pF = pBF->pFFace();
	Cell *pCInit = pF->pCCellOpposite(WI.IQE.pBFBFace());

	int iFaceIdx;
	for (iFaceIdx = 0; (iFaceIdx < pCInit->iNumFaces() &&
			    pCInit->pFFace(iFaceIdx) != pF); iFaceIdx++)
	  {}
	// Now copy data into a List<BdryEdge*> ...
	List<BdryEdgeBase*> LpBE;
	for (int ii = LpBFEncroached.iLength()-1; ii >= 0; ii--) {
	  BdryEdgeBase *pBE =
	    dynamic_cast<BdryEdgeBase*>(LpBFEncroached[ii]);
	  assert(pBE->qValid());
	  LpBE.vAppendItem(pBE);
	}
	// Must be sure to get the right insertion point here.
	double adNewLoc[2];
	(dynamic_cast<BdryEdgeBase*>(pBF))->vGetCCSPoint(adNewLoc);
	iVertsInserted =
	  iInsertPointOnBdryWatson(pCInit, iFaceIdx, adNewLoc,
				   LpBE, false);
      }
      break;
    }
  default:
    // Should never have any other cases for Mesh2D insertion.
    vFatalError("impossible insertion attempt",
		"Mesh2D::iInsertPointWatson(WatsonInfo&)");
  }
  int iVertsAfter = iNumVerts();
  int iTriCellsAfter = iNumTriCells();
  int iBFacesAfter = iNumBdryFaces();

  vMessage(4, "Started with %d verts,  now have %d.\n", iVertsBefore,
	   iVertsAfter);
  vMessage(4, "Started with %d tri cells,  now have %d.\n", iTriCellsBefore,
	   iTriCellsAfter);
  vMessage(4, "Started with %d bfaces, now have %d.\n", iBFacesBefore,
	   iBFacesAfter);

  vUpdateLengthScale();

  assert(WI.vecpCNew.size() == 0);
  assert(WI.vecpBFNew.size() == 0);
  for (int iC = iTriCellsBefore; iC < iTriCellsAfter; iC++) {
    WI.vecpCNew.push_back(pCCell(iC));
  }

  for (int iBF = iBFacesBefore; iBF < iBFacesAfter; iBF++) {
    WI.vecpBFNew.push_back(pBFBFace(iBF));
  }

  // The following assertion can't always be an equality because under
  // some circumstances, a vert is created, then never inserted into the
  // mesh because (for example) it lies outside the domain.
  assert(iVertsAfter - iVertsBefore >= iVertsInserted);

  // We might have deleted some verts to clear out space before
  // splitting an edge, but this return value ignores those, counting
  // only verts added to the mesh.  By implication, keeping track of
  // total vertex count by adding up the number of verts inserted is not
  // accurate.
  return iVertsInserted;
}

// The following class is used primarily for building face-cell
// connectivity from cell-vertex data.  Can be used for both 2D and 3D
// applications.
class VertSort {
private:
  Vert *pV;
  int iCLeft;
public:
  VertSort() : pV(pVInvalidVert), iCLeft(-1) {}
  VertSort(Vert * const pVIn, const int iCIn) : pV(pVIn), iCLeft(iCIn) {}
  VertSort operator=(const VertSort& VS)
    {
      if (this != &VS) {
	pV = VS.pV;
	iCLeft = VS.iCLeft;
      }
      return (*this);
    }
  Vert *pVVert() const
    {
      return pV;
    }
  int iLeftCell() const {return iCLeft;}
  bool operator<(const VertSort& VS) const
    {
      return (pV < VS.pV ||
	      (pV == VS.pV && iCLeft > VS.iCLeft));
    }
  bool operator==(const VertSort& VS) const
    {
      return (pV == VS.pV && iCLeft == VS.iCLeft);
    }
}; // End of class VertSort

static int iMaxNumEdges;
static VertSort* aVSVerts;

static void vSetUpVertList(const List<Face*> LpFHull)
{
  int iNumHullFaces = LpFHull.iLength();
  // Figure out the connectivity of the hull so that new faces can be
  // built.  This requires making and sorting a list of edges.

  // First make sure there's enough space to store them all; if not,
  // make it so there is.  The last 'new' here leaks memory, but we
  // aren't talking about very much, and it only happens once per
  // invocation of the -program- (not once per subroutine call).  So who
  // cares?
  if (iNumHullFaces*2 > iMaxNumEdges) {
    if (aVSVerts != NULL)
      delete [] aVSVerts;
    aVSVerts = new VertSort[iNumHullFaces*3];
    iMaxNumEdges = iNumHullFaces*3;
  }
  int i;
  for (i = 0; i < iNumHullFaces; i++) {
    Face *pF = LpFHull[i];
    aVSVerts[2*i  ] = VertSort(pF->pVVert(0), i);
    aVSVerts[2*i+1] = VertSort(pF->pVVert(1), i);
  }
  std::sort(aVSVerts, aVSVerts + 2*iNumHullFaces);
}


//@ Add a vertex using Watson's method when only interior cells are involved
void Mesh2D::vInsertWatsonInterior(Vert* const pVNew,
				   const List<Cell*>& LpCToRemove,
				   const List<Face*>& LpFToRemove,
				   const List<Face*>& LpFHull)
{
  SUMAA_LOG_EVENT_BEGIN(WATSON_INSERT);
  int i, iRegion = LpCToRemove[0]->iRegion();

  // Disconnect deleted cells from all their faces and mark for later
  // cleanup.
  for (i = 0; i < LpCToRemove.iLength(); i++) {
    Cell *pC = LpCToRemove[i];
    for (int iF = 0; iF < 3; iF++) {
      Face *pF = pC->pFFace(iF);
      pF->vRemoveCell(pC);
    }
    pC->vMarkDeleted();
  }

  // Mark deleted faces for later cleanup.
  for (i = 0; i < LpFToRemove.iLength(); i++) {
    LpFToRemove[i]->vMarkDeleted();
  }

  // Create a new cell for every hull face and connect to hull face.
  // Reset vertex hints for hull face verts.
  int iNumHullFaces = LpFHull.iLength();
  Cell **apCNew = new Cell*[iNumHullFaces];
  for (i = 0; i < iNumHullFaces; i++) {
    Face *pF = LpFHull[i];
    Cell *pC = apCNew[i] = pCNewCell();
    pC->vClearFaces();
    pC->vAddFace(pF);
    pC->vSetRegion(iRegion);
    pF->vAddCell(pC);
    pF->pVVert(0)->vSetHintFace(pF);
    pF->pVVert(1)->vSetHintFace(pF);
  }

  // Now create a new face for every hull vertex, and connect to new
  // cells.
  vSetUpVertList(LpFHull);
  int iNumNewFaces = iNumHullFaces;
  for (i = 0; i < iNumNewFaces; i++) {
    Face *pF = pFNewFace();
    Vert *pV = aVSVerts[2*i].pVVert();
    int iFLeft, iFRight;

    // Now determine which face is on which side of the -new- face
    // pVNew, pV0, pV1.
    {
      int iFA = aVSVerts[2*i + 0].iLeftCell();
      int iFB = aVSVerts[2*i + 1].iLeftCell();
      Face *pFA = LpFHull[iFA];
      Vert *pVOther;
      if (pFA->pVVert(0) != pV)
	pVOther = pFA->pVVert(0);
      else {
	assert(pFA->pVVert(1) != pV);
	pVOther = pFA->pVVert(1);
      }
      int iOrient = iOrient2D(pVNew, pV, pVOther);
      switch (iOrient) {
      case 1:
	// iFA is on the -right- side of the new face.
	iFRight = iFA;
	iFLeft  = iFB;
	break;
      case -1:
	// iFA is on the -left- side of the new face.
	iFRight = iFB;
	iFLeft  = iFA;
	break;
      default:
	vFoundBug("colinear points");
      }
    }

    pF->vAssign(apCNew[iFLeft], apCNew[iFRight], pV, pVNew);
    apCNew[iFLeft ]->vAddFace(pF);
    apCNew[iFRight]->vAddFace(pF);
    assert(pF->iFullCheck());
  }

  // Set vertex hint for new vertex.
  pVNew->vSetHintFace(pFFace(iNumFaces() - 1));
  spVUpdateLS.insert(pVNew);

  // Validity checks on new cells.
#ifndef NDEBUG
  SUMAA_LOG_EVENT_BEGIN(WATSON_CHECK);
  for (i = 0; i < iNumHullFaces; i++) {
    Cell *pC = apCNew[i];
    assert(pC->iFullCheck());
    assert(iOrient2D(pC->pVVert(0), pC->pVVert(1), pC->pVVert(2)) == 1);
  }
  // assert(qIsLocallyDelaunay());
  SUMAA_LOG_EVENT_END(WATSON_CHECK);
#endif
  delete [] apCNew;
  SUMAA_LOG_EVENT_END(WATSON_INSERT);
}

static int
iOtherFaceIndexOnVertex (List < structFaceList > LsF,
			 List < int >LiIndex,
			 int iActualF,
			 Vert * pV)
{
  // Use to return the index (in the list LiIndex) of the
  // face (other than iActualF, in LsF) connected to the vertex
  // pV. It will only search using the indices in LiIndex.

  int i;
  bool qFaceFound = false;

  i = -1;

  while (!qFaceFound) {
    i++;
    // Make sure we're not beyond the range of the list
    assert (i < LiIndex.iLength ());
    // check if any of the two vertices of the face are the same as pV
    if ((LsF[LiIndex[i]].pV0 == pV) || (LsF[LiIndex[i]].pV1 == pV)) {
      // found one.. make sure it's not the same as the actual one..
      if (i != iActualF) {
	// it's not the same...
	qFaceFound = true;
      }
    }
  }

  return i;
}

bool Mesh2D::
qDeleteInteriorPointsInBall(const Cell *pCBdry, const double adLoc[],
			    const double dRadius, int& iNDeleted)
{
  // Delete all the points that are located inside the given
  // ball. pBFBdry is a BFace that contains or is very close to adLoc,
  // so it can be used as a starting point.

  // Here is the idea:
  // 1) Find the cell on the other side of the face, put it in the list of
  //    cells to check.
  // 2) Get all the neighbor cells too, put them in the list.
  // 3) Go through the list. Check the verts. Mark the ones to be
  //    deleted as we go. Put the ones outside the diametral circle in
  //    another list.
  // 4) If a vertex is inside the diametral circle (to be deleted), the
  //    cells next to the cell who owned the vertex have to be put in
  //    the list too.
  // 5) Stops when there are no more cells to check.
  // 6) Delete all the vertex that were marked for deletion.
  // 7) Done.

  iNDeleted = 0;

  assert(pCBdry->eType() == Cell::eBdryEdge ||
	 pCBdry->eType() == Cell::eIntBdryEdge);
  const Face *pF = pCBdry->pFFace(0);

  List<Cell*> LpCCheck;
  List<Vert*> LpVDelete;
  List<const Vert*> LpVOk;

  int i, j, k; // counters

  Cell *pCOpp, *pCOther, *pCCheck; // Cell pointers
  Vert *pVCheck;

  const Vert *pV0 = pF->pVVert(0); // The face endpoints -- never put them
  const Vert *pV1 = pF->pVVert(1); // in the list

  double adTestPt[2]; // Used to calculate

  // These are not considered part of the diametral circle
  LpVOk.vAddItem(pV0);
  LpVOk.vAddItem(pV1);

  pCOpp = pF->pCCellOpposite(pCBdry);

  // Add the first cells to check to the list
  LpCCheck.vAddItem(pCOpp);

  // Add its neighbors too...
  for (i = 0; i < 3; i++) {
    Face *pFCell = pCOpp->pFFace(i);
    if (pFCell != pF) {
      pCOther = pFCell->pCCellOpposite(pCOpp);
      // Add to list if it's not a boundary face
      if (pCOther->eType() == Cell::eTriCell) {
	LpCCheck.vAddItem(pCOther);
      }
    }
  }

  // LpVOk contains two items: pV0 and pV1
  // LpVDelete is still empty
  assert(LpVOk.iLength() == 2);
  assert(LpVDelete.iLength() == 0);

  // Now go through the list of cells, and check the verts..
  for (i = 0; i < LpCCheck.iLength(); i++) {
    // Get the cell
    pCCheck = LpCCheck[i];
    // Check the three verts
    for (j = 0; j < 3; j++) {
      // Get the vert
      pVCheck = pCCheck->pVVert(j);
      if ( (pVCheck->iVertType() == Vert::eInterior) &&
	   (!LpVDelete.qInList(pVCheck)) &&
	   (!LpVOk.qInList(pVCheck)) ) {
				// If the vert has not been checked yet..
	adTestPt[0] = pVCheck->dX();
	adTestPt[1] = pVCheck->dY();
	if ( iFuzzyComp(dDIST2D(adTestPt, adLoc), dRadius) < 0 ) {
	  // The point is inside the diametral circle.
	  // 1) add the vert to the delete list
	  // 2) add the cell's neighbors to the list of cells..
	  LpVDelete.vAddItem(pVCheck);
	  for (k = 0; k < 3; k++) {
	    Face *pFCell = pCCheck->pFFace(k);
	    pCOpp = pFCell->pCCellOpposite(pCCheck);
	    if (pCOpp->eType() == Cell::eTriCell) {
	      LpCCheck.vAddItem(pCOpp);
	    }
	  } // Loop over all the faces of the cell
	} // inside diametral circle
	else {
	  // Vert is outside. Add it to the "OK" list
	  LpVOk.vAddItem(pVCheck);
	} // if it's inside the diametral circle
      } // if the vert has not been checked yet
    } // loop over the three verts of a cell
  } // loop over the list of cells

  // Now delete all the verts that are in the list of verts to be deleted
  if (LpVDelete.iLength() > 0) {
    //		vMessage(2, "Have %d verts to delete\n", LpVDelete.iLength());
    for (i = 0; i < LpVDelete.iLength(); i++) {
      // Get the vert
      pVCheck = LpVDelete[i];
      // Make sure we can delete it
      assert(pVCheck->iVertType() == Vert::eInterior);
      // Don't delete any quasi-structured verts
      if (pVCheck->qIsStructured()) continue;
      // Delete it
      int iNSwaps;
#ifndef NDEBUG
      bool qSuccess =
#endif
	qRemoveVertByContraction(pVCheck, iNSwaps);
#ifndef NDEBUG
      assert(qSuccess);
#endif
      iNDeleted ++;
    }

  } // there are verts to delete
  // Always return true, because in 2D, the scheme can't fail.
  return true;
}

static void
vFindBdryCircumcircleCellsAndOrderThem
(
 Cell * pInitialCell,	        // Cell next to boundary
 int iFaceIdx,		        // BFace index..
 const double adNewPoint[2],	// New point on the bdry
 Vert * pVNew,		        // The new vert..
 List<BdryEdgeBase*>& LpBEEncroachedBFaces,	// List of encroached BFaces
 List<Cell*>& LpCDelCells,	// List of Deletable Cells
 List<structFaceList>& LsAllFaces,	// Keeps track of all faces (hull & interior)
 List<int>& LiHullFaces,	// Index to hull faces in LsAllFaces
 List<int>& LiIntFaces,	        // Index to int faces in LsAllFaces
 List<structHullOrder>& LsHullOrder,	// Has the hull faces in CCW order
 eEncroachType eET	// Are we using circles or diametral lenses?
 )
{
  int iDC;
  int i, j;
  Face *pInitialF;

  // Make sure cell is good..
  assert (pInitialCell->iFullCheck ());
  assert (pInitialCell->eType () == Cell::eTriCell);

  // Make sure face is good...
  assert (pInitialCell->pFFace (iFaceIdx)->iFullCheck ());
  pInitialF = pInitialCell->pFFace (iFaceIdx);

  /* Not necessarily always true...
     // Make sure point is really in the circumcircle of the given cell
     // It should, if it's on it..
     assert (iIncircle (pInitialCell->pVVert (0),
     pInitialCell->pVVert (1),
     pInitialCell->pVVert (2),
     pVNew) >= 0
     );
  */
  // Make sure the face is a boundary face..
  assert (pInitialF->qIsBdryFace ());


#ifndef NDEBUG
  // Check that it's really on it..
  // Get the BdryEdge
  // BdryEdge *pBE = ((BdryEdge *) pInitialF->pCCellOpposite (pInitialCell));
  // Get the patch
  // BdryPatch2D *pPatch = (BdryPatch2D*)pBE->pPatchPointer ();
  // Check if the distance is small...
  // double dDist = pPatch->dDistFromBdry (adNewPoint);
  // assert(iFuzzyComp (dDist/1e6, 0.0) == 0);
#endif

  // First item on the list is the initial cell..
  LpCDelCells.vAddItem (pInitialCell);
  // First item on the Encroached list is the BdryEdge..
  LpBEEncroachedBFaces.vAddItem
    (dynamic_cast<BdryEdgeBase*>(pInitialF->pCCellOpposite(pInitialCell)));
  // Where is it?
  // int iEdgeIndex = LpBEEncroachedBFaces.iSeek(((BdryEdge*)pInitialF->pCCellOpposite(pInitialCell)));

  // Must now determine a list of cells in the mesh with the new point
  // in their circumcircle

  SUMAA_LOG_EVENT_BEGIN (WATSON_BDRY_NEIGHBORS);

  for (iDC = 0; iDC < LpCDelCells.iLength (); iDC++) {
    // Go through the list and find the cells with the point in their circumcircle
    Cell *pC;
    Face *pF;
    // Find neighbors
    // 3 neighbors..
    for (i = 0; i < 3; i++) {
      structFaceList sF;
      bool qAlreadyChecked;
      // Since we're already going through the faces, might as
      // well keep track of them here...

      // Get the face
      pF = LpCDelCells[iDC]->pFFace (i);

      // Make sure it's good
      assert (pF->iFullCheck ());
      // Get the opposite cell
      pC = pF->pCCellOpposite (LpCDelCells[iDC]);
      // Make sure it's good
      assert (pC->iFullCheck ());
      // Set the sF stuff to keep track of faces
      sF.pF = pF;
      sF.pC = LpCDelCells[iDC];
      sF.pCOpp = pC;
      sF.pV0 = pF->pVVert (0);
      sF.pV1 = pF->pVVert (1);

      // Make sure it hasn't been visited already
      qAlreadyChecked = false;
      if (iDC > 0) {
	for (j = 0; j < iDC; j++) {
	  if (LpCDelCells[j] == pC) {
	    // Checked already, skipping
	    qAlreadyChecked = true;
	  }
	}
      }

      if (!qAlreadyChecked) {
	// only check if it hasn't already been added
	// Make sure it's really a cell, not a BFace
	if (!pF->qIsBdryFace ()) {
	  if ((pC->eType () == Cell::eTriCell) &&
	      iIncircle (pC->pVVert (0), pC->pVVert (1), pC->pVVert (2),
			 pVNew) == 1) {
	    // It's in it, add cell to list and increment counter
	    // No need to check for duplicates again, it's been done.
	    LpCDelCells.vAppendItem (pC);
	    // the face visited is definitely an interior face..
	    int iIns = LsAllFaces.iLength();
	    LiIntFaces.vAppendItem (iIns);
	    LsAllFaces.vAppendItem (sF);
	  }
	  else {
	    // not in circle, so the face is a hull face..
	    // keep track
	    int iIns = LsAllFaces.iLength();
	    LiHullFaces.vAppendItem (iIns);
	    // add it to faces list
	    LsAllFaces.vAppendItem (sF);
	  }
	}
	else {
	  // do something for a BFace.. should check encroachment
	  // The opposite CELL is the BFACE ...
	  if (dynamic_cast<BdryEdgeBase *>(sF.pCOpp)
	      ->eIsPointEncroaching (adNewPoint, eET) != eClean) {
	    LpBEEncroachedBFaces.vAddItem (dynamic_cast<BdryEdgeBase *>
					   (sF.pCOpp));
	  }
	  // However, a BFace is a Hull face for sure: add it
	  // This will be the bface..
	  int iIns = LsAllFaces.iLength();
	  LiHullFaces.vAppendItem (iIns);
	  LsAllFaces.vAppendItem (sF);
	}			// end if it's a BFace or not

      }				// if not checked already..

    }				// end for loop over neighbors

  }				// end for loop over list of deletable cells

  SUMAA_LOG_EVENT_END (WATSON_BDRY_NEIGHBORS);

  // This should match up
  assert ((LiHullFaces.iLength () + 2 * LiIntFaces.iLength ()) == 3 * LpCDelCells.iLength ());

  // Alright. We have 1) a list of deletable cells, 2) a list of hull faces
  // and 3) a list of deletable faces.
  // We now need to ORDER the hull faces so that they're right-handed.

  SUMAA_LOG_EVENT_BEGIN (WATSON_BDRY_ORDER);

  Vert *pVSearch, *pVPrev, *pV0, *pV1;

  structHullOrder *asHullOrder = new structHullOrder[LiHullFaces.iLength ()];

  // We want to make sure that the first face in the ordered list is
  // the boundary face that we are inserting on. This will make it easier
  // later on..

  bool qFaceFound;

  qFaceFound = false;

  i = -1;
  while (!qFaceFound) {
    i++;
    if (LsAllFaces[LiHullFaces[i]].pF == pInitialF)
      qFaceFound = true;
  }

  bool qIsRight;

  // i now contains the index of the face that needs to be split...
  pV0 = LsAllFaces[LiHullFaces[i]].pV0;
  pV1 = LsAllFaces[LiHullFaces[i]].pV1;

  // Will try to find the face connected to pV1
  pVSearch = pV1;
  pVPrev = LsAllFaces[LiHullFaces[i]].pV0;

  qIsRight = true;

  // Add this to the list, as the first item...
  asHullOrder[0].iIndex = i;
  asHullOrder[0].qIsRight = qIsRight;

  // Make an ordered list of all hull faces. Then, once the whole
  // list is done, check to see if it's right-handed. If it's not,
  // reverse it

  for (j = 1; j < LiHullFaces.iLength (); j++) {
    // Must find other face connected to pVSearch
    int iNewIdx;

    iNewIdx = iOtherFaceIndexOnVertex (LsAllFaces,
				       LiHullFaces,
				       asHullOrder[j - 1].iIndex,
				       pVSearch);

    if (pVSearch == LsAllFaces[LiHullFaces[iNewIdx]].pV0) {
      pV0 = LsAllFaces[LiHullFaces[iNewIdx]].pV0;
      pV1 = LsAllFaces[LiHullFaces[iNewIdx]].pV1;
      qIsRight = true;
    }
    else {
      pV1 = LsAllFaces[LiHullFaces[iNewIdx]].pV0;
      pV0 = LsAllFaces[LiHullFaces[iNewIdx]].pV1;
      qIsRight = false;
    }

    // If the above is true, then this will work
    pVSearch = pV1;
    // Add to list..
    asHullOrder[j].iIndex = iNewIdx;
    asHullOrder[j].qIsRight = qIsRight;
  }

  // Checking if right-handed

  if (asHullOrder[LiHullFaces.iLength () - 1].qIsRight) {
    pVPrev = LsAllFaces[LiHullFaces[asHullOrder[LiHullFaces.iLength () - 1].iIndex]].pV0;
  }
  else {
    pVPrev = LsAllFaces[LiHullFaces[asHullOrder[LiHullFaces.iLength () - 1].iIndex]].pV1;
  }

  pV0 = LsAllFaces[LiHullFaces[asHullOrder[0].iIndex]].pV0;
  pV1 = LsAllFaces[LiHullFaces[asHullOrder[0].iIndex]].pV1;

  if (iOrient2D (pVPrev, pV0, pV1) < 0) {
    // if not, Reverse the connectivity...
    for (j = 0; j < LiHullFaces.iLength (); j++) {
      asHullOrder[j].qIsRight = !asHullOrder[j].qIsRight;
    }

    for (j = 1; j < (LiHullFaces.iLength () + 1) / 2; j++) {
      // Switch entries around
      structHullOrder sH = asHullOrder[j];
      asHullOrder[j] = asHullOrder[LiHullFaces.iLength () - j];
      asHullOrder[LiHullFaces.iLength () - j] = sH;
    }
  }

  // Transform the array in a list.. so we can pass it down..
  for (i = 0; i < LiHullFaces.iLength (); i++) {
    LsHullOrder.vAppendItem (asHullOrder[i]);
  }

  delete[]asHullOrder;

  SUMAA_LOG_EVENT_END (WATSON_BDRY_ORDER);

  // We now have a list hull faces that is right-handed, and we have the
  // list of deletable cells, and the list of deletable faces.

}

Face *Mesh2D::
pFBdryCreateNewCellsAndAssignThem (
				   List < Cell * >LpCDelCells,
				   List < structFaceList > LsAllFaces,
				   List < int >LiHullFaces,
				   List < int >LiIntFaces,
				   List < structHullOrder > LsHullOrder,
				   Vert * pVNew,
				   BFace * pNewBF,
				   int,
				   int,
				   int iSide
				   )
{
  SUMAA_LOG_EVENT_BEGIN (WATSON_BDRY_INSERT);
  // iSide is 1 for the first side in an interior bdry edge insertion, 2
  // for the second side, and 0 for a regular bdry edge insertion.
  assert(iSide == 0 || iSide == 1 || iSide == 2);

  // Will have (DeletableCells + 1) cells in total..
  // will have (DeletableCells) int faces in total..
  // Will have 1 new hull faces in total
  // (including the two new faces on the boundary)
  // Also one NEW Boundary face.. (actually a cell)..
  Cell **pCNew = new Cell *[LpCDelCells.iLength () + 1];
  Face **pFNew = new Face *[LpCDelCells.iLength ()];
  Face **pHFNew = new Face *[LpCDelCells.iLength () + 3];

  int i;
  bool qNewItem;

  Face *pFReturn = pFInvalidFace;

  for (i = 0; i < LpCDelCells.iLength () + 1; i++) {
    pCNew[i] = pCNewCell ();
    // All cells will be in the same region
    pCNew[i]->vSetRegion (LpCDelCells[0]->iRegion ());
  }

  for (i = 0; i < LpCDelCells.iLength (); i++) {
    pFNew[i] = pFNewFace ();
    // These are interior faces..
    pFNew[i]->vSetFaceLoc (Face::eInterior);
  }

  for (i = 0; i < LpCDelCells.iLength () + 3; i++) {
    pHFNew[i] = pFNewFace ();
  }

  // This should be true for this to work..
  assert (LiHullFaces.iLength () == LpCDelCells.iLength () + 2);
  // This too
  assert (LiHullFaces.iLength () == LiIntFaces.iLength () + 3);

  // Assign the hull faces...
  for (i = 0; i < LiHullFaces.iLength () + 1; i++) {
    // There is 1 more hull face than at beginning, so (+1)..
    // The face data
    structFaceList sF;
    // Vertices
    Vert *pV0;
    Vert *pV1;
    // Cells
    Cell *pCLeft;
    Cell *pCRight;
    // Boolean
    bool qRight;

    qNewItem = false;		// by default

    if (i < LiHullFaces.iLength () - 1) {
      sF = LsAllFaces[LiHullFaces[LsHullOrder[i + 1].iIndex]];
      qRight = LsHullOrder[i + 1].qIsRight;
      pCLeft = pCNew[i];
      pCRight = sF.pF->pCCellOpposite (sF.pC);
      if (qRight) {
	pV0 = sF.pV0;
	pV1 = sF.pV1;
      }
      else {
	pV0 = sF.pV1;
	pV1 = sF.pV0;
      }
    }
    else {
      // These are the last 2 hull faces.
      // These are the hull faces that come from the boundary face that
      // is split.

      // They both come from the same old face...
      sF = LsAllFaces[LiHullFaces[LsHullOrder[0].iIndex]];
      qRight = LsHullOrder[0].qIsRight;


      if (i < LiHullFaces.iLength ()) {
	// it's the second to last one

	if (iSide == 1 || iSide == 0) {
	  pCLeft = pCNew[i - 1];
	  pCRight = sF.pCOpp;
	}
	else {
	  pCLeft = pCNew[i - 1];
	  pCRight = pNewBF;
	  // This is the one we want..
	  pFReturn = pHFNew[i];
	}
	if (iSide == 0 || iSide == 2) qNewItem = true;
	if (qRight) {
	  pV0 = sF.pV0;
	  pV1 = pVNew;
	}
	else {
	  pV0 = sF.pV1;
	  pV1 = pVNew;
	}
      }
      else {
	// it's the last one
	if (iSide != 2) {
	  pCLeft = pCNew[0];
	  // This is going to be the new Bface.
	  pCRight = pNewBF;
	  // This is the one we want..
	  pFReturn = pHFNew[i];
	  qNewItem = true;
	}
	else {
	  pCLeft = pCNew[0];
	  pCRight = sF.pCOpp;
	}
	if (qRight) {
	  pV0 = pVNew;
	  pV1 = sF.pV1;
	}
	else {
	  pV0 = pVNew;
	  pV1 = sF.pV0;
	}
      }
    }

    // Assign them...
    pHFNew[i]->vAssign (pCLeft, pCRight, pV0, pV1);

    if (qNewItem) {
      // pCRight is a new BFace... can not use replace face..
      if ((pCRight->eType () == Cell:: eBdryEdge) ||
	  (pCRight->eType () == Cell::eIntBdryEdge)) {
	// Set the face type
	if (pCRight->eType () == Cell::eBdryEdge) {
	  pHFNew[i]->vSetFaceLoc (Face::eBdryFace);
	  // Can assign the face now too..
	  dynamic_cast<BdryEdgeBase*>(pCRight)->vAssign (pHFNew[i]);
	  // And the bdry patch...
	  dynamic_cast<BdryEdgeBase*>(pCRight)->
	    vSetPatch(dynamic_cast<BdryEdgeBase*>(sF.pCOpp)->iPatchIndex(),
		      &B2D);
	}
	else {
	  pHFNew[i]->vSetFaceLoc (Face::eBdryTwoSide);
	  // Need to know both face and assign at the same time..
	  // This is done later..
	}
      }
      else {
	// Should never be different than the two types above
	assert (0);
      }

    }
    else {
      // pCRight is not new.. can just replace stuff...
      // set face type
      switch (pCRight->eType ()) {
      case Cell::eBdryEdge:
	pHFNew[i]->vSetFaceLoc (Face::eBdryFace);
	// Out with the old face, in with the new, if it's periodic.
	{
	  BFace *pBF = dynamic_cast<BFace*>(pCRight);
	  assert(pBF); // Should never fail.
	  BdryPatch *pBP = pBF->pPatchPointer();
	  if (pBP->qIsPeriodic()) {
	    Periodic *pP = pBP->pPGetPeriodicData();
	    pP->vRemoveFace(pBP, pBF->pFFace());
	    pP->vAddFace(pBP, pHFNew[i]);
	  }
	}
	break;
      case Cell::eIntBdryEdge:
	pHFNew[i]->vSetFaceLoc (Face::eBdryTwoSide);
	break;
      case Cell::eTriCell:
	pHFNew[i]->vSetFaceLoc (Face::eInterior);
	break;
      }
      sF.pCOpp->vReplaceFace (sF.pF, pHFNew[i]);
    }

    // set a hint face too
    pV0->vSetHintFace (pHFNew[i]);
  }

  // Now do the interior faces...

  // This must be true
  assert (LpCDelCells.iLength () <= LiHullFaces.iLength ());

  for (i = 0; i < LpCDelCells.iLength (); i++) {
    Cell *pCLeft, *pCRight;
    Vert *pV0, *pV1;
    structFaceList sF;

    sF = LsAllFaces[LiHullFaces[LsHullOrder[i + 1].iIndex]];

    pCLeft = pCNew[i];
    pCRight = pCNew[i + 1];
    pV1 = pVNew;
    if (LsHullOrder[i + 1].qIsRight) {
      pV0 = sF.pV1;
    }
    else {
      pV0 = sF.pV0;
    }

    pFNew[i]->vAssign (pCLeft, pCRight, pV0, pV1);
  }

  // All that's left to do is to assign the cells..

  for (i = 0; i < LpCDelCells.iLength () + 1; i++) {
    Face *pF0, *pF1, *pF2;

    if (i == 0) {
      // First cell
      pF0 = pHFNew[LpCDelCells.iLength () + 2];
      pF1 = pHFNew[i];
      pF2 = pFNew[i];
    }
    else {
      if (i == LpCDelCells.iLength ()) {
	// Last cell...
	pF0 = pFNew[i - 1];
	pF1 = pHFNew[i];
	pF2 = pHFNew[i + 1];
      }
      else {
	// Intermediate cells..
	pF0 = pFNew[i - 1];
	pF1 = pHFNew[i];
	pF2 = pFNew[i];
      }
    }

    pCNew[i]->vAssign (pF0, pF1, pF2);
  }

  SUMAA_LOG_EVENT_END (WATSON_BDRY_INSERT);
  SUMAA_LOG_EVENT_BEGIN (WATSON_BDRY_CHECK);

  // Mark the old cells for deletion...
  for (i = 0; i < LpCDelCells.iLength (); i++) {
    LpCDelCells[i]->vMarkDeleted ();
  }

  // Mark the old INTERIOR faces for deletion...
  for (i = 0; i < LiIntFaces.iLength (); i++) {
    LsAllFaces[LiIntFaces[i]].pF->vMarkDeleted ();
  }

  // Mark the old HULL faces for deletion..
  for (i = 0; i < LiHullFaces.iLength (); i++) {
    LsAllFaces[LiHullFaces[i]].pF->vMarkDeleted ();
  }

  // Vertex hint
  pVNew->vSetHintFace (pHFNew[LiHullFaces.iLength ()]);


  // make sure the new cells are valid...
  for (i = 0; i < LpCDelCells.iLength () + 1; i++) {
    assert (pCNew[i]->iFullCheck ());
  }

  // Make sure the int faces are valid
  for (i = 0; i < LpCDelCells.iLength (); i++) {
    // interior face
    assert (pFNew[i]->iFullCheck ());
  }

  // Make sure the hull faces are valid
//   for (i = 0; i < LiHullFaces.iLength () + 1; i++) {
//     assert (pHFNew[i]->iFullCheck ());
//   }

  // Should not check mesh now, in case the insertion is on both sides of the
  // boundary face...
  // make sure mesh is valid
  //    assert(qValid());

  SUMAA_LOG_EVENT_END (WATSON_BDRY_CHECK);

  delete [] pCNew;
  delete [] pFNew;
  delete [] pHFNew;

  return pFReturn;
}

bool Mesh2D::qEnsurePointIsInterior(List <BdryEdgeBase*>& LpBEEncroachedBFaces,
				    int iWhichEdge,
				    Vert *pVNew,
				    bool qFix  // to fix or not
				    )
{
  // Returns whether or not a point was inserted outside..
  bool qExterior = false;

//	vMessage(2, "Content of the encroached BFaces list -- Edge: %d -- point to check: (%f, %f)\n", iWhichEdge, pVNew->dX(), pVNew->dY());
//	for (int iBF = 0; iBF < LpBEEncroachedBFaces.iLength(); iBF++) {
//		BdryEdge *pBE = LpBEEncroachedBFaces[iBF];
//		vMessage(2, "%d - Patch: %d\n", pBE, pBE->pPatchPointer());
//		vMessage(2, "   pV0: (%f, %f) pV1 (%f, %f)\n",
//			pBE->pVVert(0)->dX(), pBE->pVVert(0)->dY(), pBE->pVVert(1)->dX(), pBE->pVVert(1)->dY());
//	}

	// Get the original patch...
  BdryPatch2D *pInitialPatch = dynamic_cast<BdryPatch2D*>
    (LpBEEncroachedBFaces[iWhichEdge]->pPatchPointer());

  // First item is the edge that is about to be split so skip it
  for (int iBF = 0; iBF < LpBEEncroachedBFaces.iLength(); iBF++) {
    double adMid[2];
    BdryEdgeBase *pBE = LpBEEncroachedBFaces[iBF];
    BdryPatch2D *pPatch = dynamic_cast<BdryPatch2D*>(pBE->pPatchPointer());
    pBE->vGetRatioPoint(0.5, adMid);

    if ( (pPatch != pInitialPatch) && (pPatch->eGeom2DType != BdryPatch2D::ePolyLine) ) {
      // No need to look at an edge that is located on the same patch...
      // This problem does not occur with lines either, so no need to look at them
      Vert *pV0, *pV1;
      double adNew[2];
      double adEdgeMid[2];
      int iOrientMid, iOrientNew;

      // Get the edge verts
      pV0 = pBE->pVVert(0);
      pV1 = pBE->pVVert(1);

      // Get the midpoint of the segment
      adEdgeMid[0] = (pV0->dX() + pV1->dX())/2.0;
      adEdgeMid[1] = (pV0->dY() + pV1->dY())/2.0;

      // Set the new point too
      adNew[0] = pVNew->dX();
      adNew[1] = pVNew->dY();

      // Get the orientation result...
      iOrientMid = iOrient2D(pV0->adCoords(), pV1->adCoords(), adMid);

      // Do the same for the new point...
      iOrientNew = iOrient2D(pV0->adCoords(), pV1->adCoords(), adNew);

      if ( iOrientMid == iOrientNew) {
				// Both points are located on the same side of the edge...
				// Check distances...
	double dNewDist = dDIST2D(adEdgeMid, adNew);
	double dMidDist = dDIST2D(adEdgeMid, adMid);

	//				vMessage(2, "dNewDist: %f  dMidDist: %f\n", dNewDist, dMidDist);
	if (iFuzzyComp(dNewDist, dMidDist) < 0) {
	  vMessage(2, "NEW METHOD: Point is located at the exterior\n");
	  //					qExterior = true;
	}
				// Do another kind of check...
	int iOEdgeMid = iOrient2D(pV0->adCoords(), pVNew->adCoords(), adEdgeMid);
	int iOMid     = iOrient2D(pV0->adCoords(), pVNew->adCoords(), adMid);
	if (iOEdgeMid != iOMid) {
	  vMessage(2, "EVEN NEWER METHOD: Point is located outside\n");
	  qExterior = true;
	}
      }

    }

    if (qExterior && qFix) {

      pVNew->vMarkDeleted();
      // Point is located on the exterior of the mesh.  Let's fix this
      // then...
      Face *pF = pBE->pFFace(0);
      Cell *pC = pF->pCCellOpposite(pBE);
      List <BdryEdgeBase*> LpBEDummy;
      int kk = 0;
      double adPt[2];

      pBE->vGetCCSPoint(adPt);

      while (pC->pFFace(kk) != pF) {
	kk++;
      }

      vMessage(2, "Attempting insertion to fix exterior bug at (%f, %f)\n", adPt[0], adPt[1]);
				// insert point
      iInsertPointOnBdryWatson(pC, kk, adPt, LpBEDummy, false);
      LpBEDummy.vClear();

    } // if it was exterior

  } // for loop over BFaces other than the first one..

  return qExterior;
}

// CHANGE ALERT
// My inclination is to leave this one intact, called through the
// wrapper iInsertPointWatson(WatsonInfo&).  I figure changing it is too
// likely to break some of the curved bdry stuff, and I'm not interested
// in re-testing that thoroughly.  Besides, the only interrupt here is
// for insertion outside the domain, and handling that interrupt inline
// is a pretty minor issue from the elegance viewpoint.

int Mesh2D::
iInsertPointOnBdryWatson (Cell *pInitialCell, int iFaceIdx,
			  const double adNewPoint[2],
			  List <BdryEdgeBase*>& LpBEEncroachedBFaces,
			  bool qApex)
{
  // Insert a point on the boundary using Watson's insertion method.
  // The initial cell given should have the point in its circumcircle
  // as the point is supposed to be its circumcenter...!
  // The face index given by iFaceIdx should be the index of the face
  // (of pInitialCell) on which the point to be inserted lies.


  //
  //            3                         2
  //       __________               ___________
  //      |          /             |       __--|
  //  B   |        / |          3  |   __--    |
  //  d  0|      /   | 2  ==>      |_--        |  1
  //  r   |    /     |             | --__      |
  //  y   |  /       |          4  |     --__  |
  //      |/_________|             |___________|
  //            1                         0
  //
  //
  //
  //
  //

  SUMAA_LOG_EVENT_BEGIN (WATSON_BDRY);

  int iRetVal = 0;
  int iPatchIndex;
  bool qExtAtAll;
  double adNew[2];

//  vMessage(2, "BdryInsert: About to insert (%f, %f)\n", adNewPoint[0], adNewPoint[1]);


  Face *pF = pInitialCell->pFFace (iFaceIdx);
  Cell *pCBdry = pF->pCCellOpposite (pInitialCell);
  BFace *pNewBF;

  iPatchIndex = dynamic_cast<BdryEdgeBase *>(pCBdry)->iPatchIndex ();
  qExtAtAll = false;

  // Define the new vertex
  Vert *pVNew = pVNewVert ();
  pVNew->vSetCoords (2, adNewPoint);

  adNew[0] = pVNew->dX();
  adNew[1] = pVNew->dY();


  // Add a vertex type for the new vert
  if (qApex) {
    pVNew->vSetType (Vert::eBdryApex);
  }
  else {
    pVNew->vSetType (Vert::eBdry);
  }


  List <Cell*> LpCDelCells, LpCDelCells2;	// List of deletable cells

  List <structFaceList> LsAllFaces, LsAllFaces2;	// All faces (hull & int)

  List <int> LiHullFaces, LiHullFaces2;	// Index of hull faces

  List <int> LiIntFaces, LiIntFaces2;	// Index of int faces

  List <structHullOrder> LsHullOrder, LsHullOrder2;	// Order of the hull faces (CCW)

  if (pCBdry->eType () == Cell::eIntBdryEdge) {

    // This is an interior bEdge... Must do things differently..
    // Get the other face...
    const Face *pF2 = dynamic_cast<IntBdryEdge*>(pCBdry)->pFOtherFace (pF);
    // Get the other cell too
    Cell *pCInitial2 = pF2->pCCellOpposite (pCBdry);
    // Find the index..
    int iFaceIdx2 = 0;
    while (pCInitial2->pFFace (iFaceIdx2) != pF2)
      iFaceIdx2++;
    /*
      vMessage(2, "Verts for initial cell: (%f, %f) (%f, %f) (%f, %f)\n",
      pInitialCell->pVVert(0)->dX(), pInitialCell->pVVert(0)->dY(),
      pInitialCell->pVVert(1)->dX(), pInitialCell->pVVert(1)->dY(),
      pInitialCell->pVVert(2)->dX(), pInitialCell->pVVert(2)->dY());
      vMessage(2, "Verts for other cell: (%f, %f) (%f, %f) (%f, %f)\n",
      pCInitial2->pVVert(0)->dX(), pCInitial2->pVVert(0)->dY(),
      pCInitial2->pVVert(1)->dX(), pCInitial2->pVVert(1)->dY(),
      pCInitial2->pVVert(2)->dX(), pCInitial2->pVVert(2)->dY());
    */

    Face *pOtherF1, *pOtherF2;

    // Get the cells first

    vFindBdryCircumcircleCellsAndOrderThem
      (pInitialCell, iFaceIdx, adNewPoint, pVNew, LpBEEncroachedBFaces,
       LpCDelCells, LsAllFaces, LiHullFaces, LiIntFaces, LsHullOrder,
       this->eEncroachmentType());

    vFindBdryCircumcircleCellsAndOrderThem
      (pCInitial2, iFaceIdx2, adNewPoint, pVNew, LpBEEncroachedBFaces,
       LpCDelCells2, LsAllFaces2, LiHullFaces2, LiIntFaces2, LsHullOrder2,
       this->eEncroachmentType());

    // And create the new stuff..

    if (LpBEEncroachedBFaces.iLength() > 1) {
      BdryEdgeBase *pBE = dynamic_cast<BdryEdgeBase*>
	(pInitialCell->pFFace(iFaceIdx)->pCCellOpposite(pInitialCell));
      int iEdgeIdx = LpBEEncroachedBFaces.iSeek(pBE);
      qExtAtAll = qEnsurePointIsInterior(LpBEEncroachedBFaces, iEdgeIdx, pVNew, true);
    }

    if (!qExtAtAll) {
      // Must delete points inside the diametral circle on both sides now..
      //			vDeleteDiametralCirclePoints(pF, pCBdry);
      //			vDeleteDiametralCirclePoints(pF2, pCBdry);

      // Create new face...
      pNewBF = pIBENewIntBdryEdge ();

      pOtherF1 = pFBdryCreateNewCellsAndAssignThem
	(LpCDelCells, LsAllFaces, LiHullFaces, LiIntFaces, LsHullOrder,
	 pVNew, pNewBF, Cell::eIntBdryEdge, Face::eBdryTwoSide, 1);

      pOtherF2 = pFBdryCreateNewCellsAndAssignThem
	(LpCDelCells2, LsAllFaces2, LiHullFaces2, LiIntFaces2, LsHullOrder2,
	 pVNew, pNewBF, Cell::eIntBdryEdge, Face::eBdryTwoSide, 2);

      // Do the face assignment ...
      // Must invert the second face...
      Cell *pCLeft, *pCRight;
      Vert *pV0, *pV1;

      pCLeft = pOtherF2->pCCellLeft ();
      pCRight = pOtherF2->pCCellRight ();
      pV0 = pOtherF2->pVVert (0);
      pV1 = pOtherF2->pVVert (1);


      pOtherF2->vAssign (pCRight, pCLeft, pV1, pV0);

      // Do the assignments..
      dynamic_cast<IntBdryEdge*>(pNewBF)->vAssign   (pOtherF1, pOtherF2);
      dynamic_cast<IntBdryEdge*>(pNewBF)->vSetPatch (iPatchIndex, &B2D);
      // Set face type
      pOtherF1->vSetFaceLoc (Face::eBdryTwoSide);
      pOtherF2->vSetFaceLoc (Face::eBdryTwoSide);

      pVNew->vSetType(Vert::eBdryTwoSide);

      // Make sure to set a length scale for this vertex.
      spVUpdateLS.insert(pVNew);

      assert(!pCBdry->qDeleted());
      BdryEdgeBase* pBE = dynamic_cast<BdryEdgeBase*>(pCBdry);
      pBE->vClearOffsetInsertionRequest();
      pBE->vAllowOffsetInsertion();
      iRetVal ++;

      assert(pNewBF->iFullCheck());
    }
    else {
      // Some points were about to be inserted to the exterior...
      // Try to insert the "old" point again

      // Shouldn't we delete the "old" pVNew here? It will never be used...


      vMessage(2, "And trying to insert at (%f, %f) again...\n", adNew[0], adNew[1]);
      // And try to insert the old one too after
      List <BdryEdgeBase*> LpBEDummy;
      BdryEdgeBase *pBE = LpBEEncroachedBFaces[0];
      pF = pBE->pFFace(0);
      Cell *pC = pF->pCCellOpposite(pBE);
      int kk = 0;

      while (pC->pFFace(kk) != pF) {
	kk++;
      }
      SUMAA_LOG_EVENT_END(WATSON_BDRY);
      iRetVal += iInsertPointOnBdryWatson(pC, kk, adNew, LpBEDummy, false);
      SUMAA_LOG_EVENT_BEGIN(WATSON_BDRY);
    }
  } // done with case of a interior boundary edge
  else {

    // Get the list of cells with the new point in their circumcircle
    // LsHullOrder will have the list of hull faces in CCW order

    vFindBdryCircumcircleCellsAndOrderThem
      (pInitialCell, iFaceIdx, adNewPoint, pVNew, LpBEEncroachedBFaces,
       LpCDelCells, LsAllFaces, LiHullFaces, LiIntFaces, LsHullOrder,
       this->eEncroachmentType());

    if (LpBEEncroachedBFaces.iLength() > 1) {
      BdryEdgeBase *pBE = dynamic_cast<BdryEdgeBase*>
	(pInitialCell->pFFace(iFaceIdx)->pCCellOpposite(pInitialCell));
      int iEdgeIdx = LpBEEncroachedBFaces.iSeek(pBE);
      qExtAtAll = qEnsurePointIsInterior(LpBEEncroachedBFaces, iEdgeIdx, pVNew, true);
    } // if number of encroached bfaces is not 1

    if (!qExtAtAll) {

      // Check for points inside the diametral circle of the
      // boundary face and delete them

      // BFace is pF ...
      int iNDeleted;
      {
	double adLoc[2], dRadius;
	// For now, only diametral circle deletion.
	BdryEdgeBase* pBE = dynamic_cast<BdryEdgeBase*>(pCBdry);
	pBE->vCircumcenter(adLoc);
	dRadius = dDIST2D(adLoc, pBE->pVVert(0)->adCoords());
	bool qSuccess = qDeleteInteriorPointsInBall(pCBdry, adLoc,
						    dRadius, iNDeleted);
	assert(qSuccess);
      }

      if (0 != iNDeleted) {
	// Does pF still exist?
	assert(!pF->qDeleted());
	assert(!pCBdry->qDeleted());

	pInitialCell = pF->pCCellOpposite(pCBdry);

				// Delete the vertex -- not used
	pVNew->vMarkDeleted();

	iFaceIdx = 0;
	while (pInitialCell->pFFace(iFaceIdx) != pF) iFaceIdx++;

	List <BdryEdgeBase*> LpBEDummy;

	SUMAA_LOG_EVENT_END(WATSON_BDRY);
	iRetVal += iInsertPointOnBdryWatson(pInitialCell, iFaceIdx, adNew,
					    LpBEDummy, false);
	SUMAA_LOG_EVENT_BEGIN(WATSON_BDRY);
      }
      else {
	// Now insert...
	pNewBF = pBFNewBFace ();
	Face *pOtherF = pFBdryCreateNewCellsAndAssignThem
	  (LpCDelCells, LsAllFaces, LiHullFaces, LiIntFaces, LsHullOrder,
	   pVNew, pNewBF, Cell::eBdryEdge, Face::eBdryFace, 0);

	// Do the assignments..
	dynamic_cast<BdryEdgeBase *>(pNewBF)->vAssign (pOtherF);
// 	dynamic_cast<BdryEdgeBase *>(pNewBF)->vSetPatch (iPatchIndex, &B2D);
	// Set face type
	pOtherF->vSetFaceLoc (Face::eBdryFace);

	// Make sure to set a length scale for this vertex.
	spVUpdateLS.insert(pVNew);

	assert(!pCBdry->qDeleted());
	BdryEdgeBase* pBE = dynamic_cast<BdryEdgeBase*>(pCBdry);
	pBE->vClearOffsetInsertionRequest();
	pBE->vAllowOffsetInsertion();
	iRetVal ++;
	assert (pNewBF->iFullCheck ());
      }

    }
    else {
      // Some points were inserted outside... Try to insert the
      // point again...

      vMessage(2, "And trying to insert at (%f, %f) again...\n", adNew[0], adNew[1]);
      // And try to insert the old one too after
      List <BdryEdgeBase*> LpBEDummy;
      BdryEdgeBase *pBE = LpBEEncroachedBFaces[0];
      pF = pBE->pFFace(0);
      Cell *pC = pF->pCCellOpposite(pBE);
      int kk = 0;

      while (pC->pFFace(kk) != pF) {
	kk++;
      }

      SUMAA_LOG_EVENT_END(WATSON_BDRY);
      iRetVal += iInsertPointOnBdryWatson(pC, kk, adNew, LpBEDummy, false);
      SUMAA_LOG_EVENT_BEGIN(WATSON_BDRY);

    } // If any exterior fix was done or not..

  } // if it's an interior bdry or not

  SUMAA_LOG_EVENT_END (WATSON_BDRY);

  return iRetVal;
}

// This call is never used, so nuke it.

// int Mesh2D::
// iInsertPointOnBdryWatson (const double adNewPt[2],
//			  List <BdryEdge*>& LpBEEncroachedBFaces)
// {

//   // This is used when the cells / faces kept in memory to determine which ones
//   // the new point lies on do not exist anymore...
//   // Try to find them...

//   double dSmallestDist =1e20;

//   int    iClosestFace = -1;
//   BdryEdgeBase *pBE;

//   SUMAA_LOG_EVENT_BEGIN (WATSON_BDRY_UNKNOWN);

//   for (int i = 0; i < iNumTotalBdryFaces(); i++) {
//     if (i < iNumBdryFaces()) {
//       pBE = dynamic_cast<BdryEdgeBase*>(pBFBFace (i));
//     } else {
//       pBE = dynamic_cast<BdryEdgeBase*>(pIBEIntBdryEdge(i-iNumBdryFaces()));
//     }

//     vMessage(2, "Looking at edge %d : (%f,%f)->(%f,%f)\n",
//	     i, pBE->pVVert(0)->dX(), pBE->pVVert(0)->dY(),
//	     pBE->pVVert(1)->dX(), pBE->pVVert(1)->dY());

//     if (pBE->qIsPointEncroaching (adNewPt, qCircleEncroaching)) {
//       // Find distance from boundary.
//       double dDistance;

//       dDistance = fabs(pBE->pPatchPointer()->dDistFromBdry(adNewPt));

//       vMessage(2, "   Distance: %f\n", dDistance);

//       if (dDistance < dSmallestDist) {
//	dSmallestDist = dDistance;
//	iClosestFace = i;
//       }
//     }
//   }

//   if (iClosestFace > -1) {
//     if (iClosestFace < iNumBdryFaces()) {
//       pBE = ((BdryEdge*) pBFBFace (iClosestFace));
//     } else {
//       pBE = ((BdryEdge*) pIBEIntBdryEdge(iClosestFace-iNumBdryFaces()));
//     }

//     vMessage(2, "Using this edge: %d\n", iClosestFace);

//     Face *pF = pBE->pFFace (0);
//     Cell *pC = pF->pCCellOpposite (pBE);
//     int j = 0;
//     // Find the face index...
//     while (pC->pFFace (j) != pF)
//       j++;

//     SUMAA_LOG_EVENT_END (WATSON_BDRY_UNKNOWN);

//     return iInsertPointOnBdryWatson (pC, j, adNewPt, LpBEEncroachedBFaces, false);
//   }
//   else {
//     fprintf (stdout, "ERROR: iInsertPointOnBdryWatson: Could not find a BFace with the point on it\n");
//     assert (false);
//     return -1;
//   }
// }
