#include <queue>

#include <math.h>
#include "GR_assert.h"
#include "GR_events.h"
#include "GR_misc.h"
#include "GR_BFace.h"
#include "GR_Geometry.h"
#include "GR_Quality.h"
#include "GR_Sort.h"
#include "GR_Vec.h"
#include "GR_VolMesh.h"
#include "GR_Watson.h"
#include "GR_WatsonInfo.h"

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

//  bool operator==(const SEdgeSort& SA, const SEdgeSort& SB) {
//    return (SA.pV0 == SB.pV0 &&
//	  SA.pV1 == SB.pV1 &&
//	  SA.iF  == SB.iF );
//  }

static FaceSort* aFSEdges = static_cast<FaceSort*> (NULL);
static int iMaxNumEdges = 0;

//  //@@ A function to compare one edge to another, with no ties allowed
//  extern "C" {
//    static int iEdgeComp(const void * pSEdge1, const void * pSEdge2)
//    {
//      struct SEdgeSort {Vert *pV0, *pV1; int iF;};
//      const SEdgeSort *pSE1 = (const SEdgeSort*) pSEdge1;
//      const SEdgeSort *pSE2 = (const SEdgeSort*) pSEdge2;
//      Vert *pVMin1, *pVMax1, *pVMin2, *pVMax2;
//      if (pSE1->pV0 < pSE1->pV1) {
//        pVMin1 = pSE1->pV0;
//        pVMax1 = pSE1->pV1;
//      }
//      else {
//        pVMin1 = pSE1->pV1;
//        pVMax1 = pSE1->pV0;
//      }
//      if (pSE2->pV0 < pSE2->pV1) {
//        pVMin2 = pSE2->pV0;
//        pVMax2 = pSE2->pV1;
//      }
//      else {
//        pVMin2 = pSE2->pV1;
//        pVMax2 = pSE2->pV0;
//      }
//      if ( (pVMin1 < pVMin2) ||
//	 ((pVMin1 == pVMin2) && (pVMax1 < pVMax2)) ||
//	 ((pVMin1 == pVMin2) && (pVMax1 == pVMax2)
//	  && (pSE1->iF < pSE2->iF)) )  return (-1);
//      else return (1);
//    }
//  }

static void vSetUpEdgeList(const std::set<Face*> spFHull)
{
  int iNumHullFaces = spFHull.size();
  // 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*3 > iMaxNumEdges) {
    if (aFSEdges != static_cast<FaceSort*>(NULL))
      delete [] aFSEdges;
    aFSEdges = new FaceSort[iNumHullFaces*3];
    iMaxNumEdges = iNumHullFaces*3;
  }
  int i = 0;
  for (std::set<Face*>::iterator iterF = spFHull.begin();
       iterF != spFHull.end(); iterF++, i++) {
    Face *pF = *iterF;
    aFSEdges[3*i  ] = FaceSort(pF->pVVert(0), pF->pVVert(1), i, pF);
    aFSEdges[3*i+1] = FaceSort(pF->pVVert(1), pF->pVVert(2), i, pF);
    aFSEdges[3*i+2] = FaceSort(pF->pVVert(2), pF->pVVert(0), i, pF);
  }
  sort(aFSEdges, aFSEdges + 3*iNumHullFaces);
}

void
vCheckBdryFragmentForSubsegEncroachment(const std::set<BFace*>& spBF,
					std::set<FaceSort>& sFSEncSubseg,
					const double * const adPoint)
{
  sFSEncSubseg.clear();
  std::set<Face*> spFHull;
  std::set<BFace*>::iterator iterBF;
  for (iterBF = spBF.begin(); iterBF != spBF.end(); iterBF++) {
    spFHull.insert((*iterBF)->pFFace());
  }
  SUMAA_LOG_EVENT_BEGIN(SHEWCHUK_HULL);
  vSetUpEdgeList(spFHull);
  int iNumHullFaces = spFHull.size();
  SUMAA_LOG_EVENT_END(SHEWCHUK_HULL);

  SUMAA_LOG_EVENT_BEGIN(SHEWCHUK_SUBSEG);
  for (int iEdge = 0; iEdge < 3*iNumHullFaces; iEdge ++) {
    Vert *pV0 = aFSEdges[iEdge].pVVert(0);
    Vert *pV1 = aFSEdges[iEdge].pVVert(1);
    if (!(iEdge == 3*iNumHullFaces - 1) &&
	((aFSEdges[iEdge].pVVert(0) == aFSEdges[iEdge+1].pVVert(0) &&
	  aFSEdges[iEdge].pVVert(1) == aFSEdges[iEdge+1].pVVert(1))
	 ||
	 (aFSEdges[iEdge].pVVert(1) == aFSEdges[iEdge+1].pVVert(0) &&
	  aFSEdges[iEdge].pVVert(0) == aFSEdges[iEdge+1].pVVert(1)))) {
      // Interior edge in bdry fragment
      assert(aFSEdges[iEdge].iLeftCell()  != aFSEdges[iEdge+1].iLeftCell() );
      iEdge++; // Skip over both halves of the pair.
    }
    vMessage(4, "Checking for bdry subseg encroachments between:\n");
    vMessage(4, "  Vert %p (type %d):  (%12.6g, %12.6g, %12.6g)\n",
	     pV0, pV0->iVertType(), pV0->dX(), pV0->dY(),
	     pV0->dZ());
    vMessage(4, "  Vert %p (type %d):  (%12.6g, %12.6g, %12.6g)\n",
	     pV1, pV1->iVertType(), pV1->dX(), pV1->dY(),
	     pV1->dZ());
    vMessage(4, "Checking whether both fall on bdry segments...");
    if ((pV0->iVertType() != Vert::eBdryApex &&
	 pV0->iVertType() != Vert::eBdryCurve)
	||
	(pV1->iVertType() != Vert::eBdryApex &&
	 pV1->iVertType() != Vert::eBdryCurve)
	) {
      // Can't possibly be a bdry subsegment, because one or both verts
      // don't fall on bdry segments.
      vMessage(4, "no.\n");
      continue;
    }
    vMessage(4, "yes.\n");

    // Now check to see if this really is a bdry subsegment.
    // For now, take the easiest possible way out, by finding out
    // everything about the neighborhood of pV0, and then determining
    // which BFace's also contain pV1.  If there are more than two such
    // BFace's, then pV0-pV1 is the intersection of more than two
    // BdryPatch's.  If there are precisely two, then the critical
    // question is whether the BFace's lie on the same patch.
    std::set<Cell*> spCInc;
    {
      vMessage(4, "  Checking whether this is really a bdry subseg...");
      std::set<Vert*> spVNeigh;
      std::set<BFace*> spBFInc;
      bool qBdryVert;
      vNeighborhood(pV0, spCInc, spVNeigh, &spBFInc, &qBdryVert);
      assert(qBdryVert);
      std::set<BFace*> spBFCommon;
      for (iterBF = spBFInc.begin(); iterBF != spBFInc.end(); iterBF++) {
	BFace *pBF = *iterBF;
	assert(pBF->qHasVert(pV0));
	if (pBF->qHasVert(pV1))
	  spBFCommon.insert(pBF);
      }
      switch (spBFCommon.size()) {
      case 0:
	// No common BFaces; connected across the interior of the mesh.
	vMessage(4, "not even on the bdry.\n");
	continue;
      case 1:
	// Should never happen.
	vFatalError("Exactly one common bdry face between two points?!?",
		    "bdry subseg encroachment check during insertion");
	break;
      case 2:
	// Are the two BFace's on the same or different BdryPatch's?
	{
	  std::set<BFace*>::iterator iter = spBFCommon.begin();
	  BFace *pBF0 = *iter;
	  iter++;
	  BFace *pBF1 = *iter;
	  assert(pBF0->eType() == Cell::eTriBFace ||
		 pBF0->eType() == Cell::eIntTriBFace);
	  assert(pBF1->eType() == Cell::eTriBFace ||
		 pBF1->eType() == Cell::eIntTriBFace);
	  BdryPatch* pPatch0 = pBF0->pPatchPointer();
	  BdryPatch* pPatch1 = pBF1->pPatchPointer();
	  if (pPatch0 == pPatch1) {
	    vMessage(4, "no.\n");
	    continue;
	  }
	}
	// Deliberate fallthrough, because this and all subsequent cases
	// add an encroached subseg to the list.
      default:
	vMessage(4, "yes (%d common BFace's)\n",
		 static_cast<int>(spBFCommon.size()));
      } // Done with switch over # of common BFace's
    } // Done checking whether it's really a bdry subseg

    double dRadius = 0.5 * dDIST3D(pV0->adCoords(), pV1->adCoords());
    double adMid[] = {0.5 * (pV0->dX() + pV1->dX()),
		      0.5 * (pV0->dY() + pV1->dY()),
		      0.5 * (pV0->dZ() + pV1->dZ())};
    if (adPoint) {
      // Check to see whether this particular point encroaches on the
      // subsegment.
      vMessage(4, "  Checking whether the edge is encroached...");
      double dDistance = dDIST3D(adMid, adPoint);
      if (iFuzzyComp(dDistance, dRadius) == -1) {
	vMessage(4, "yes.\n");
	sFSEncSubseg.insert(aFSEdges[iEdge]);
	vMessage(4, "  Added to encroached subseg list\n");
	break;
      } // Done checking this neighboring vert
    } // Done checking for encroachment when a location is given
    else {
      // Check to see if the segment pV0-pV1 is really encroached by
      // some point in one of the tets incident on each when no point is given
      vMessage(4, "  Checking whether the edge is encroached...");
      // First, make a list of all the verts in all the tets incident on
      // both pV0 and pV1.
      std::set<Vert*> spVCommonNeigh;
      std::set<Cell*>::iterator iterC;
      for (iterC = spCInc.begin(); iterC != spCInc.end(); iterC++) {
	Cell *pC = *iterC;
	assert(pC->qHasVert(pV0));
	if (pC->qHasVert(pV1)) {
	  spVCommonNeigh.insert(pC->pVVert(0));
	  spVCommonNeigh.insert(pC->pVVert(1));
	  spVCommonNeigh.insert(pC->pVVert(2));
	  spVCommonNeigh.insert(pC->pVVert(3));
	}
      }
      spVCommonNeigh.erase(pV0);
      spVCommonNeigh.erase(pV1);

      for (std::set<Vert*>::iterator iterV = spVCommonNeigh.begin();
	   iterV != spVCommonNeigh.end(); iterV++) {
	Vert *pV = *iterV;
	double dDistance = dDIST3D(adMid, pV->adCoords());
	if (iFuzzyComp(dDistance, dRadius) == -1) {
	  vMessage(4, "yes.\n");
	  sFSEncSubseg.insert(aFSEdges[iEdge]);
	  vMessage(4, "  Added to encroached subseg list\n");
	  break;
	} // Done checking this neighboring vert
      } // Done with loop over common neighbors
    } // Done checking for actual encroachment
  } // Done looping over all edges in the hull
  SUMAA_LOG_EVENT_END(SHEWCHUK_SUBSEG);
}


//  bool qInsertInThisTet(Cell* const pC)
//  {
//    SUMAA_LOG_EVENT_BEGIN(STEINER_CELL_QUAL);
//    assert(pC->iFullCheck());
//    assert(pC->eType() == Cell::eTet);

//    if (!pC->qWellSized()) {
//      double dLengthScale = 0.25 * (+ pC->pVVert(0)->dLS()
//				  + pC->pVVert(1)->dLS()
//				  + pC->pVVert(2)->dLS()
//				  + pC->pVVert(3)->dLS());
//      double dDiam = 2 * ((TetCell*)pC)->dCircumradius();

//      if (dDiam > dLengthScale) {
//        SUMAA_LOG_EVENT_END(STEINER_CELL_QUAL);
//        return true;
//      }

//      pC->vMarkWellSized();
//    }

//    if (pC->iShapeClass() != Cell::eWellShaped) {
//      int i;
//      double adQual[1];
//      vShortestToRadius(pC, &i, adQual);
//      if (adQual[0] < 0.5) {
//        SUMAA_LOG_EVENT_END(STEINER_CELL_QUAL);
//        return true;
//      }
//      else {
//        pC->vSetShapeClass(Cell::eWellShaped);
//        SUMAA_LOG_EVENT_END(STEINER_CELL_QUAL);
//        return false;
//      }
//    }
//    return false;
//  }

static void vBuildHull(const std::set<Cell*>& spC,
		       std::set<Face*>& spFInterior,
		       std::set<Face*>& spFHull)
{
  spFInterior.clear();
  spFHull.clear();
  std::multiset<Face*> mspFTemp;
  for (std::set<Cell*>::iterator iterC = spC.begin();
       iterC != spC.end(); iterC++) {
    Cell* pC = *iterC;
    for (int iF = 0; iF < 4; iF ++)
      mspFTemp.insert(pC->pFFace(iF));
  }

  for (std::multiset<Face*>::iterator iterF = mspFTemp.begin();
       iterF != mspFTemp.end(); iterF++) {
    Face *pFTmp = *iterF;
    int iCount = mspFTemp.count(pFTmp);
    assert(iCount == 1 || iCount == 2);
    if (iCount == 1) spFHull.insert(pFTmp);
    else             spFInterior.insert(pFTmp);
  }
}


void VolMesh::vGetWatsonData(const double adPoint[3],
			     Cell* const pCSeed,
			     std::set<Cell*>& spCSearch,
			     std::set<Face*>& spFHull,
			     std::set<Face*>& spFToRemove,
			     std::set<BFace*>& spBFEncroached) const
{
  static int iHullsChecked = 0, iTotalHullSize = 0;
  SUMAA_LOG_EVENT_BEGIN(WATSON_NEIGHBORS);
  Vert VDummy;
  VDummy.vSetCoords(3, adPoint);

  spCSearch.clear();
  spFHull.clear();
  spFToRemove.clear();
  // Don't clear LpBFEncroached, as the first entry is set elsewhere and
  // has meaning.  The others may not need to be zeroed either, but this
  // does no harm.

  // This assertion is always true in principle, although some
  // particularly bizarrely-shaped tets cause problems in practice, even
  // if VDummy is at an edge midside.
//    assert(iInsphere(pCSeed->pVVert(0), pCSeed->pVVert(1), pCSeed->pVVert(2),
//		   pCSeed->pVVert(3), &VDummy) == 1);

  // **** Find Watson neighborhood ****
  // Set initial cell search list
  std::queue<Cell*> qpC;
  spCSearch.insert(pCSeed);
  qpC.push(pCSeed);
  int iF;
  for (iF = 0; iF < 4; iF++) {
    Face *pF = pCSeed->pFFace(iF);
    spCSearch.insert(pF->pCCellOpposite(pCSeed));
    qpC.push(pF->pCCellOpposite(pCSeed));
  }

  // Until all cells have been checked...
  while (!qpC.empty()) {
    Cell *pCCand = qpC.front();
    qpC.pop();
    if (pCCand->eType() == Cell::eTet) {
      // For tets:
      if (iInsphere(pCCand->pVVert(0), pCCand->pVVert(1),
		    pCCand->pVVert(2), pCCand->pVVert(3), &VDummy)
	  != -1) {
	// If new point is in circumsphere, add neighbors to search list
	for (iF = 0; iF < 4; iF++) {
	  Face *pF = pCCand->pFFace(iF);
	  Cell *pCTmp = pF->pCCellOpposite(pCCand);
	  if (spCSearch.count(pCTmp) == 0) {
	    spCSearch.insert(pCTmp);
	    qpC.push(pCTmp);
	  }
	}
      }
      // Else delete this cell from search list
      else {
	spCSearch.erase(pCCand);
      }
    }
    else {
      assert(pCCand->eType() == Cell::eTriBFace ||
	     pCCand->eType() == Cell::eIntTriBFace);
      TriBFaceBase* pTBFCand = dynamic_cast<TriBFaceBase*> (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.
      if (pTBFCand->eIsPointEncroaching(adPoint, eEncroachmentType(), true)
	  != eClean) {
	spBFEncroached.insert(pTBFCand);
      }
      // Delete bface from search list
      spCSearch.erase(pTBFCand);
      // Add the cells adjacent to the BFace on each side, if this is an
      // internal BFace.
      if (pTBFCand->eType() == Cell::eIntTriBFace &&
	  !pTBFCand->pPatchPointer()->qDisjointFrom(&VDummy)) {
	Cell *pCTmp = pTBFCand->pFFace(0)->pCCellOpposite(pCCand);
	if (spCSearch.count(pCTmp) == 0) {
	  spCSearch.insert(pCTmp);
	  qpC.push(pCTmp);
	}
	pCTmp = pTBFCand->pFFace(1)->pCCellOpposite(pCCand);
	if (spCSearch.count(pCTmp) == 0) {
	  spCSearch.insert(pCTmp);
	  qpC.push(pCTmp);
	}
      }
    } // Done with BFace case...
  }

  // 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 *pBFCand = *(spBFEncroached.begin());
      if (pBFCand->eType() == Cell::eTriBFace &&
	  pBFCand->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[3];
	pCSeed->vCentroid(adCent);
	if (iOrient3D(pBFCand->pVVert(0)->adCoords(),
		      pBFCand->pVVert(1)->adCoords(),
		      pBFCand->pVVert(2)->adCoords(), adCent) !=
	    iOrient3D(pBFCand->pVVert(0)->adCoords(),
		      pBFCand->pVVert(1)->adCoords(),
		      pBFCand->pVVert(2)->adCoords(), adPoint)) {
	  // Rig things so that this internal bdry face will be split.
	  double adTmp[3];
	  (dynamic_cast<TriBFaceBase*>(pBFCand))->vCircumcenter(adTmp);
	  pBFCand->vSetInsertionLocation(adTmp);
	  pBFCand->vForbidOffsetInsertion();
	}
      }
    }
    break;
  default:
    for (std::set<BFace*>::iterator iterBF = spBFEncroached.begin();
	 iterBF != spBFEncroached.end(); iterBF++) {
      (*iterBF)->vForbidOffsetInsertion();
    }
    break;
  }

  vBuildHull(spCSearch, spFToRemove, spFHull);

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

void VolMesh::vCleanupWatson(WatsonInfo &WI)
  // This is a member function so that a polymorphic call in
  // InsertionQueue::vQualityRefine() will work properly.
{
  Vert VDummy;
  VDummy.vSetCoords(3, WI.adCoords());

  std::set<Cell*>::iterator iterC;
  std::set<Face*>::iterator iterF;
  std::set<BFace*>::iterator iterBF;
    
#ifndef NDEBUG
  for (iterC = WI.spCToRemove.begin(); iterC != WI.spCToRemove.end();
       iterC++) {
    Cell *pC = *iterC;
    assert(iInsphere(pC->pVVert(0), pC->pVVert(1), pC->pVVert(2),
		     pC->pVVert(3), &VDummy) != -1);
  }
  for (iterBF = WI.spBFEncroached.begin(); iterBF != WI.spBFEncroached.end();
       iterBF++) {
    BFace *pBF = *iterBF;
    assert(pBF->eIsPointEncroaching(WI.adCoords(), eEncroachmentType(), true)
	   != eClean);
  }
#endif

  bool qChanged, qRecomputeHull = false;
  do {
    qChanged = false;

    // Are there any faces in the hull that are coplanar with the new
    // vertex?  If there are, they necessarily fall into one of two
    // categories: boundary faces and interior faces.

    // For boundary coplanar hull faces, there are two subcases: faces
    // on patches containing the new vertex, and faces not on patches
    // containing the new vertex.  In the former case, all is well:
    // these boundary faces will be removed from the mesh and new
    // boundary faces substituted.  In the latter case, these boundary
    // faces are separated from the new vertex by a bdry curve, so there
    // is no encroachment on the BFace by this vert (nor on the cell
    // built on the BFace).  So each such must be -DELETED- from the
    // list of BFaces to remove, and the cells opposite them must be
    // -DELETED- from the list of Cells to remove.

    // Also, in some circumstances the angle at the bdry curve can be
    // reflex --- this totally precludes attaching to those BFace's, so
    // the same action should be taken.

    for (iterBF = WI.spBFEncroached.begin(); iterBF != WI.spBFEncroached.end();
	 iterBF++) {
      BFace *pBF = *iterBF;
      // This test is probably wrong for internal boundary faces.
      if (iOrient3D(pBF->pVVert(0), pBF->pVVert(1), pBF->pVVert(2),
		    &VDummy) != 1 &&
	  pBF->pPatchPointer()->qDisjointFrom(&VDummy)) {
	// This is the bad case...
	vMessage(4, "Found a coplanar or reflex BFace in a disjoint patch... removing\n");
	Face *pF = pBF->pFFace();
	Cell *pC = pF->pCCellOpposite(pBF);
	WI.spBFEncroached.erase(pBF);
	WI.spCToRemove.erase(pC);
	qChanged = true;
      }
    }

    // Now remove any "orphaned" BFaces from the encroached list.  Any
    // BFace that is connected to a cell that no longer is slated for
    // removal must be removed from the encroached list, since
    // constraints don't allow us to fix this encroachment.
    for (iterBF = WI.spBFEncroached.begin(); iterBF != WI.spBFEncroached.end();
	 iterBF++) {
      BFace *pBF = *iterBF;
      if (pBF->eType() == Cell::eTriBFace) {
	Face *pF = pBF->pFFace();
	Cell *pC = pF->pCCellOpposite(pBF);
	if (WI.spCToRemove.count(pC) == 0) {
	  WI.spBFEncroached.erase(pBF);
	  qChanged = true;
	}
      }
      else {
	Cell *pC0 = pBF->pFFace(0)->pCCellOpposite(pBF);
	Cell *pC1 = pBF->pFFace(1)->pCCellOpposite(pBF);
	if (WI.spCToRemove.count(pC0) == 0 &&
	    WI.spCToRemove.count(pC1) == 0) {
	  WI.spBFEncroached.erase(pBF);
	  qChanged = true;
	}
      }
    }

    // Now find all the faces that bound the set of tets to be removed,
    // and the faces that are shared by two tets that will be removed.
    // The former will have new tets build on them; the latter will be
    // removed from the mesh.
    if (qChanged || qRecomputeHull) vBuildHull(WI.spCToRemove,
					       WI.spFToRemove, WI.spFHull);

    // After deleting cells adjacent to flat/reflex BFaces (above),
    // there may be orphaned cells, which are still tagged for removal
    // but which have non-removable cells between them and the new
    // vertex.  These will show up as badly-oriented faces in the hull
    // and must be cleaned up; this is accomplished by tagging the
    // adjacent cells to -not- be removed from the mesh and re-building
    // the hull.
    for (iterF = WI.spFHull.begin(); iterF != WI.spFHull.end(); iterF++) {
      Face *pF = *iterF;
      Cell *pCL = pF->pCCellLeft();
      if (pF->qIsBdryFace()) {
	// We shouldn't mess with encroached bdry faces, but we do need
	// to check non-encroached bdry faces.
	BFace *pBF;
	if (pCL->eType() != Cell::eTet) {
	  pBF = dynamic_cast<BFace*>(pCL);
	}
	else {
	  pBF = dynamic_cast<BFace*>(pF->pCCellRight());
	}
	if (WI.spBFEncroached.count(pBF) == 1)
	  continue;
      }
      int iOrient = iOrient3D(pF->pVVert(0), pF->pVVert(1),
			      pF->pVVert(2), &VDummy);
      bool qLeftCellInList = (WI.spCToRemove.count(pCL) != 0);
      if (qLeftCellInList) iOrient *= -1;
      // Now orientation should be 1 for proper hull faces, unless those
      // faces are on the boundary and will be coplanar.
      if (iOrient != 1) {
	// Must be sure to remove the offending cell from the removal
	// list.  Any cells that are now offensive as a result of this
	// will be taken care of on a subsequent pass.
	vMessage(4, "Found a mis-oriented Face in insertion hull... fixing\n");
	qChanged = qRecomputeHull = true;
	if (qLeftCellInList)
	  WI.spCToRemove.erase(pCL);
	else
	  WI.spCToRemove.erase(pF->pCCellRight());
      }
    }
  } while (qChanged);

  vMessage(4, "Watson cleanup: Done sorting faces into hull and interior.\n");
  vMessage(4, "Watson cleanup: # cells to delete: %d\n",
	   static_cast<int>(WI.spCToRemove.size()));
  vMessage(4, "Watson cleanup: # faces to delete: %d\n",
	   static_cast<int>(WI.spFToRemove.size()));
  vMessage(4, "Watson cleanup: # faces in hull: %d\n",
	   static_cast<int>(WI.spFHull.size()));
  assert(4*WI.spCToRemove.size() - 2*WI.spFToRemove.size() ==
	 WI.spFHull.size());
}

void WatsonInfo::vStrip3DHull()
{
  assert(eType() == InsertionQueueEntry::eSubSeg);

  // Retain only cells and bdry faces that contain both of the verts in
  // the subseg.
  for (std::set<Cell*>::iterator iterC = spCToRemove.begin();
       iterC != spCToRemove.end(); ) {
    Cell *pC = *iterC;
    if (pC->qHasVert(IQE.pVVert(0)) && pC->qHasVert(IQE.pVVert(1))) {
      iterC++;
    }
    else {
      std::set<Cell*>::iterator iterTmp = iterC++;
      spCToRemove.erase(iterTmp);
    }
  }

  for (std::set<BFace*>::iterator iterBF = spBFEncroached.begin();
       iterBF != spBFEncroached.end(); ) {
    BFace* pBF = *iterBF;
    if (pBF->qHasVert(IQE.pVVert(0)) && pBF->qHasVert(IQE.pVVert(1))) {
      iterBF++;
    }
    else {
      std::set<BFace*>::iterator iterTmp = iterBF++;
      spBFEncroached.erase(iterTmp);
    }
  }

  vBuildHull(spCToRemove, spFToRemove, spFHull);

  vMessage(4, "Watson Strip3D: Done sorting faces into hull and interior.\n");
  vMessage(4, "Watson Strip3D: # cells to delete: %d\n",
	   static_cast<int>(spCToRemove.size()));
  vMessage(4, "Watson Strip3D: # faces to delete: %d\n",
	   static_cast<int>(spFToRemove.size()));
  vMessage(4, "Watson Strip3D: # faces in hull: %d\n",
	   static_cast<int>(spFHull.size()));
  assert(4*spCToRemove.size() - 2*spFToRemove.size() ==
	 spFHull.size());
}

bool WatsonInfo::qIs3DHullOkay(const enum eEncroachType eET) const
{
  Vert VDummy;
  VDummy.vSetCoords(3, adCoords());

  std::set<Cell*>::iterator iterC;
  std::set<Face*>::iterator iterF;
  std::set<BFace*>::iterator iterBF;
    
#ifndef NDEBUG
  for (iterC = spCToRemove.begin(); iterC != spCToRemove.end();
       iterC++) {
    Cell *pC = *iterC;
    assert(iInsphere(pC->pVVert(0), pC->pVVert(1), pC->pVVert(2),
		     pC->pVVert(3), &VDummy) != -1);
  }
  for (iterBF = spBFEncroached.begin(); iterBF != spBFEncroached.end();
       iterBF++) {
    BFace *pBF = *iterBF;
    assert(pBF->eIsPointEncroaching(adCoords(), eET, true)
	   != eClean);
  }
#endif

  // Are there any faces in the hull that are coplanar with the new
  // vertex?  If there are, they necessarily fall into one of two
  // categories: boundary faces and interior faces.

  // For boundary coplanar hull faces, there are two subcases: faces
  // on patches containing the new vertex, and faces not on patches
  // containing the new vertex.  In the former case, all is well:
  // these boundary faces will be removed from the mesh and new
  // boundary faces substituted.  In the latter case, these boundary
  // faces are separated from the new vertex by a bdry curve, which can
  // cause serious insertion problems.

  // Also, in some circumstances the angle at the bdry curve can be
  // reflex --- this totally precludes attaching to those BFace's, so
  // the same action should be taken.

  for (iterBF = spBFEncroached.begin(); iterBF != spBFEncroached.end();
       iterBF++) {
    BFace *pBF = *iterBF;
    // This test is probably wrong for internal boundary faces.
    if (iOrient3D(pBF->pVVert(0), pBF->pVVert(1), pBF->pVVert(2),
		  &VDummy) != 1 &&
	pBF->pPatchPointer()->qDisjointFrom(&VDummy)) {
      // This is the bad case...
      vMessage(4, "Found a coplanar or reflex BFace in a disjoint patch...\n");
      return false;
    }
  }

  // Now remove any "orphaned" BFaces from the encroached list.  Any
  // BFace that is connected to a cell that no longer is slated for
  // removal must be removed from the encroached list, since
  // constraints don't allow us to fix this encroachment.
  for (iterBF = spBFEncroached.begin(); iterBF != spBFEncroached.end();
       iterBF++) {
    BFace *pBF = *iterBF;
    if (pBF->eType() == Cell::eTriBFace) {
      Cell *pC = pBF->pFFace()->pCCellOpposite(pBF);
      if (spCToRemove.count(pC) == 0) {
	return false;
      }
    }
    else {
      assert(pBF->eType() == Cell::eIntTriBFace);
      Cell *pC0 = pBF->pFFace(0)->pCCellOpposite(pBF);
      Cell *pC1 = pBF->pFFace(1)->pCCellOpposite(pBF);
      if (spCToRemove.count(pC0) == 0 &&
	  spCToRemove.count(pC1) == 0) {
	return false;
      }
    }
  } // Done checking BFace's

  // Are there any badly-oriented faces in the hull?  If so, bail out.
  
  for (iterF = spFHull.begin(); iterF != spFHull.end();
       iterF++) {
    Face *pF = *iterF;
    Cell *pCL = pF->pCCellLeft();
    if (pF->qIsBdryFace()) {
      // We shouldn't mess with encroached bdry faces, but we do need
      // to check non-encroached bdry faces.
      BFace *pBF;
      if (pCL->eType() != Cell::eTet) {
	pBF = dynamic_cast<BFace*>(pCL);
      }
      else {
	pBF = dynamic_cast<BFace*>(pF->pCCellRight());
      }
      if (spBFEncroached.count(pBF) != 0)
	continue;
    }
    int iOrient = iOrient3D(pF->pVVert(0), pF->pVVert(1),
			    pF->pVVert(2), &VDummy);
    bool qLeftCellInList = (spCToRemove.count(pCL) != 0);
    if (qLeftCellInList) iOrient *= -1;
    // Now orientation should be 1 for proper hull faces, unless those
    // faces are on the boundary and will be coplanar.
    if (iOrient != 1) {
      // Must be sure to remove the offending cell from the removal
      // list.  Any cells that are now offensive as a result of this
      // will be taken care of on a subsequent pass.
      vMessage(4, "Found a mis-oriented Face in insertion hull...\n");
      return false;
    }
  }

  return true;
}

bool VolMesh::
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::eTriBFace ||
	 pCBdry->eType() == Cell::eIntTriBFace);
  const Face *pF = pCBdry->pFFace(0);

  std::set<Cell*> spCCheck;
  std::set<Vert*> spVDelete;
  std::set<const Vert*> spVOk;

  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
  const Vert *pV2 = pF->pVVert(2);

  double adTestPt[3]; // Used to calculate

  // These are not considered part of the diametral circle
  spVOk.insert(pV0);
  spVOk.insert(pV1);
  spVOk.insert(pV2);

  pCOpp = pF->pCCellOpposite(pCBdry);

  // Add the first cells to check to the list
  spCCheck.insert(pCOpp);

  // Add its neighbors too...
  for (i = 0; i < 4; 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::eTet) {
	spCCheck.insert(pCOther);
      }
    }
  }

  // spVOk contains three items: pV0, pV1, and pV2.
  // spVDelete is still empty
  assert(spVOk.size() == 3);
  assert(spVDelete.size() == 0);

  // Now go through the list of cells, and check the verts..
  for (std::set<Cell*>::iterator iterC = spCCheck.begin();
       iterC != spCCheck.end(); iterC++) {
    // Get the cell
    pCCheck = *iterC;
    // Check the four verts
    for (j = 0; j < 4; j++) {
      // Get the vert
      pVCheck = pCCheck->pVVert(j);
      if ( (pVCheck->iVertType() == Vert::eInterior) &&
	   (spVDelete.count(pVCheck) == 0) &&
	   (spVOk.count(pVCheck) == 0) ) {
	// If the vert has not been checked yet..
	adTestPt[0] = pVCheck->dX();
	adTestPt[1] = pVCheck->dY();
	adTestPt[2] = pVCheck->dZ();
	if ( iFuzzyComp(dDIST3D(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..
	  spVDelete.insert(pVCheck);
	  for (k = 0; k < 4; k++) {
	    Face *pFCell = pCCheck->pFFace(k);
	    pCOpp = pFCell->pCCellOpposite(pCCheck);
	    if (pCOpp->eType() == Cell::eTet) {
	      spCCheck.insert(pCOpp);
	    }
	  } // Loop over all the faces of the cell
	} // inside diametral circle
	else {
	  // Vert is outside. Add it to the "OK" list
	  spVOk.insert(pVCheck);
	} // if it's inside the diametral circle
      } // if the vert has not been checked yet
    } // loop over the four 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
  bool qSuccess = true;
  if (spVDelete.size() > 0) {
    vMessage(3, "Have %d verts to delete\n",
	     static_cast<int>(spVDelete.size()));
    for (std::set<Vert*>::iterator iterV = spVDelete.begin();
	 iterV != spVDelete.end(); iterV++) {
      // Get the vert
      pVCheck = *iterV;
      // Make sure we can delete it
      assert(pVCheck->iVertType() == Vert::eInterior);
      // Delete it
      int iNSwaps;
      qSuccess = qRemoveVert(pVCheck, iNSwaps);
      iNDeleted ++;
    }
  } // there are verts to delete
  return qSuccess;
}


//@ Generic interface to point insertion via Watson's method.  A virtual
//  function in the Mesh class.
int VolMesh::iInsertPointWatson(WatsonInfo& WI)
{
  assert(qSimplicial());

  Vert *pVNew;
  int iCellsBefore = iNumCells();
  int iBFacesBefore = iNumBdryFaces();
  int iIntBFacesBefore = iNumIntBdryFaces();

  assert(WI.vecpCNew.empty());
  assert(WI.vecpBFNew.empty());

  switch (WI.eType()) {
  case InsertionQueueEntry::eTetCell:
    {
      assert(WI.qIs3DHullOkay(eEncroachmentType()));
      pVNew = pVNewVert();
      pVNew->vSetCoords(3, WI.adNewPoint);

      vInsertWatsonInterior(pVNew, WI.spCToRemove,
			    WI.spFToRemove, WI.spFHull);
      break;
    }
  case InsertionQueueEntry::eSubSeg:
    if (!qAllowBdryChanges) return 0;
    if (!WI.qIs3DHullOkay(eEncroachmentType())) {
      WI.vStrip3DHull();
    }
    {
      pVNew = pVNewVert();
      pVNew->vSetCoords(3, WI.adNewPoint);
      vInsertWatsonBoundary(pVNew, WI.spCToRemove,
			    WI.spFToRemove, WI.spFHull,
			    WI.spBFEncroached);
      break;
    }
  case InsertionQueueEntry::eTriBFace:
    assert(WI.qIs3DHullOkay(eEncroachmentType()));
    {
      BFace *pBF = WI.IQE.pBFBFace();

      if (pBF->qOffsetInsertionRequested()  && pBF->qOffsetInsertionAllowed()) {
	Vert *pV = pVNewVert();
	pV->vSetCoords(3, WI.adCoords());
	vMessage(3, "Attempting offset point insertion at (%f, %f, %f).\n",
		 pV->dX(), pV->dY(), pV->dZ());

	// First, delete all verts inside the ball with center at pV
	// passing through the points of the face.
	int iNDeleted = 0;
	bool qSuccess = qDeleteInteriorPointsInBall
	  (pBF, pV->adCoords(),
	   dDIST3D(pV->adCoords(), pBF->pVVert(0)->adCoords()),
	   iNDeleted);
	if (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, WI.spCToRemove, WI.spFToRemove,
				  WI.spFHull);
	  }
	  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(iInsphere(pCSeed->pVVert(0), pCSeed->pVVert(1),
			     pCSeed->pVVert(2), pCSeed->pVVert(3), pV) == 1);
	    std::set<Cell*> spCToRemove;
	    std::set<Face*> spFToRemove, spFHull;
	    std::set<BFace*> spBFDummy;
	    vGetWatsonData(pV->adCoords(), pCSeed, spCToRemove,
			   spFHull, spFToRemove, spBFDummy);
	    assert(spBFDummy == WI.spBFEncroached);
	    // 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, spCToRemove, spFToRemove, spFHull);
	  }
	}
	else {
	  // Failed to remove all the points that needed removal for
	  // this insertion.
	  return 0;
	}
      }
      else {
	// Split the face instead.

	// If bdry changes aren't allowed, then don't insert.
	if (!qAllowBdryChanges) return 0;

	// Must be sure to get the right insertion point here.
	double adNewLoc[3];
	(dynamic_cast<TriBFaceBase*>(pBF))->vCircumcenter(adNewLoc);

	// First, delete all verts inside the ball with center at pV
	// passing through the points of the face.
	int iNDeleted = 0;
	bool qSuccess = true;
	if (eEncroachmentType() != eBall) {
	  qSuccess = qDeleteInteriorPointsInBall
	    (pBF, adNewLoc, dDIST3D(adNewLoc, pBF->pVVert(0)->adCoords()),
	     iNDeleted);
	}
	if (qSuccess) {
	  vMessage(3, "Deleted %d points to clear ball for bdry insertion.\n",
		   iNDeleted);
	  // If nothing was deleted, just go ahead and insert the point.
	  // Otherwise, recompute the neighborhood.
	  if (iNDeleted != 0) {
	    Face *pFBdry = pBF->pFFace();
	    Cell *pCSeed = pFBdry->pCCellOpposite(pBF);
	    assert(iInsphere(pCSeed->pVVert(0)->adCoords(),
			     pCSeed->pVVert(1)->adCoords(),
			     pCSeed->pVVert(2)->adCoords(),
			     pCSeed->pVVert(3)->adCoords(), adNewLoc) == 1);
	    std::set<BFace*> spBFDummy;
	    vGetWatsonData(adNewLoc, pCSeed, WI.spCToRemove,
			   WI.spFHull, WI.spFToRemove, spBFDummy);
	    assert(spBFDummy == WI.spBFEncroached);
	    // 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.
	  }
	  Face *pF = pBF->pFFace();
	  Cell *pCInit = pF->pCCellOpposite(WI.IQE.pBFBFace());

	  int iFaceIdx;
	  for (iFaceIdx = 0; (iFaceIdx < pCInit->iNumFaces() &&
			      pCInit->pFFace(iFaceIdx) != pF); iFaceIdx++) {}

	  pVNew = pVNewVert();
	  pVNew->vSetCoords(3, adNewLoc);
	  vInsertWatsonBoundary(pVNew, WI.spCToRemove,
				WI.spFToRemove, WI.spFHull,
				WI.spBFEncroached);
	} // Done with successful vertex deletion
	else {
	  // Couldn't remove all the necessary points.
	  return 0;
	}
      } // Done with insertion to split a face.
    }
    break;
  default:
    // Should never have any other cases for VolMesh insertion.
    vFatalError("impossible insertion attempt",
		"VolMesh::qInsertPointWatson(WatsonInfo&)");
  }
  int iCellsAfter = iNumCells();
  int iBFacesAfter = iNumBdryFaces();
  int iIntBFacesAfter = iNumIntBdryFaces();

  vMessage(4, "Started with %d cells,  now have %d.\n", iCellsBefore,
	   iCellsAfter);
  vMessage(4, "Started with %d bfaces, now have %d.\n", iBFacesBefore,
	   iBFacesAfter);
  vMessage(4, "Started with %d interior bfaces, now have %d.\n",
	   iIntBFacesBefore, iIntBFacesAfter);

  // If we inserted on a bdry subseg, then we already know explicitly
  // what entities were changed, so this implicit approach is
  // unnecessary.
//    if (WI.spCNew.size() == 0) {
    for (int iC = iCellsBefore; iC < iCellsAfter; iC++) {
      WI.vecpCNew.push_back(pCCell(iC));
    }

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

    for (int iIBF = iIntBFacesBefore; iIBF < iIntBFacesAfter; iIBF++) {
      WI.vecpBFNew.push_back(pBFIntBFace(iIBF));
    }
//    }
  vUpdateLengthScale(pVNew);
  return 1;
}

//@ The visible interface to point insertion via Watson's method.
// Calls qInsertWatsonInterior, qInsertWatsonBoundary, or
// qInsertWatsonBdryCurve as appropriate, and with appropriate data.
bool VolMesh::qInsertPointWatson(const double adPoint[3],
				 Cell* const pCSeed)
{
  assert(qSimplicial());
  assert(pCSeed->eType() == Cell::eTet);

  Vert *pVNew = pVNewVert();
  pVNew->vSetCoords(3, adPoint);

  return qInsertPointWatson(pVNew, pCSeed);
}

bool VolMesh::qInsertPointWatson(Vert* const pVNew,
				 Cell* const pCSeed)
{
  SUMAA_LOG_EVENT_BEGIN(WATSON);
  const double *adPoint = pVNew->adCoords();

  std::set<Cell*> spCToDelete;
  std::set<BFace*> spBFEncroached, spBFUseOffset;
  std::set<Face*> spFHull, spFToRemove;

  vGetWatsonData(adPoint, pCSeed, spCToDelete, spFHull, spFToRemove,
		 spBFEncroached);

  // Don't pass value judgements on whether this point -should- be
  // inserted, just -do- it.

  // Is the point on the bdry?
  bool qOnBdry = false;
  int iLen = spBFEncroached.size();
  if (iLen > 0) {
    // Can't be on the bdry without encroaching something...
    std::set<TriBFaceBase*> spTBFHasPoint;
    for (std::set<BFace*>::iterator iterBF = spBFEncroached.begin();
	 iterBF != spBFEncroached.end(); iterBF++) {
      TriBFaceBase* pTBF = dynamic_cast<TriBFaceBase*>(*iterBF);
      Face *pF = pTBF->pFFace(0);
      double adPointTemp[] = {adPoint[0], adPoint[1], adPoint[2]};
      pF->vProjOntoFace(adPointTemp);
      if (iFuzzyComp(0, dDIST3D(adPointTemp, adPoint)) == 0 &&
	  pTBF->qPointProjectsInside(adPoint))
	spTBFHasPoint.insert(pTBF);
    }
    switch (spTBFHasPoint.size()) {
    case 0:
      // Not a bdry vert;
      break;
    case 1:
      // A bdry vert, but not on a bdry curve
      qOnBdry = true;
      vMessage(3, "WatsonInsert: On a bdry facet\n");
      break;
    case 2:
      {
	// Have to check whether the two BFace's the point is on are part
	// of the same patch.
	qOnBdry = true;
	std::set<TriBFaceBase*>::iterator iter = spTBFHasPoint.begin();
	// Yes, I'm conservative about sprinkling ++ amonst other code;
	// why do you ask? CFO-G
	BdryPatch *pBP0 = (*iter)->pPatchPointer();
	iter++;
	BdryPatch *pBP1 = (*iter)->pPatchPointer();
	if (pBP0 != pBP1) {
	  vMessage(3, "WatsonInsert: Hit a bdry curve (2 BFace's)\n");
	}
	else {
	  vMessage(3, "WatsonInsert: Two BFace's, but inside a facet\n");
	}
      }
      break;
    default: // More than three BFace's contain the point.
      // Definitely a bdry curve.
      qOnBdry = true;
      vMessage(3, "WatsonInsert: Hit a bdry curve (>2 BFace's)\n");
      break;
    } // Done checking whether it's a bdry vert.
    if (qOnBdry) {
      vInsertWatsonBoundary(pVNew, spCToDelete, spFToRemove, spFHull,
			    spBFEncroached);
      SUMAA_LOG_EVENT_END(WATSON);
      return true;
    }
  } // End cases where one or more bdry faces are encroached.

  // Insert in the interior.
  vInsertWatsonInterior(pVNew, spCToDelete, spFToRemove, spFHull);
  SUMAA_LOG_EVENT_END(WATSON);
  return true;
}

//@ Add a vertex using Watson's method when only interior cells are involved
void VolMesh::vInsertWatsonInterior(Vert* const pVNew,
				    const std::set<Cell*>& spCToRemove,
				    const std::set<Face*>& spFToRemove,
				    const std::set<Face*>& spFHull)
{
  SUMAA_LOG_EVENT_BEGIN(WATSON_INSERT);
  int i, iRegion = (*spCToRemove.begin())->iRegion();

  // Disconnect deleted cells from all their faces and mark for later
  // cleanup.
  std::set<Cell*>::iterator iterC;
  for (iterC = spCToRemove.begin(); iterC != spCToRemove.end(); iterC++) {
    Cell *pC = *iterC;
    for (int iF = 0; iF < 4; iF++) {
      Face *pF = pC->pFFace(iF);
      pF->vRemoveCell(pC);
    }
    pC->vMarkDeleted();
  }

  // Mark deleted faces for later cleanup.
  std::set<Face*>::iterator iterF;
  for (iterF = spFToRemove.begin(); iterF != spFToRemove.end(); iterF++) {
    (*iterF)->vMarkDeleted();
  }

  // Create a new cell for every hull face and connect to hull face.
  // Reset vertex hints for hull face verts.
  int iNumHullFaces = spFHull.size();
  Cell **apCNew = new Cell*[iNumHullFaces];
  i = 0;
  for (iterF = spFHull.begin(); iterF != spFHull.end(); iterF++, i++) {
    Face *pF = *iterF;
    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);
    pF->pVVert(2)->vSetHintFace(pF);
  }

  // Now create a new face for every hull edge, and connect to new
  // cells.
  vSetUpEdgeList(spFHull);
  int iNumNewFaces = (iNumHullFaces * 3) / 2;
  for (i = 0; i < iNumNewFaces; i++) {
    Face *pF = pFNewFace();
    Vert *pV0 = aFSEdges[2*i].pVVert(0);
    Vert *pV1 = aFSEdges[2*i].pVVert(1);
    int iFLeft, iFRight;

    // Now determine which face is on which side of the -new- face
    // pVNew, pV0, pV1.
    {
      int iFA = aFSEdges[2*i + 0].iLeftCell();
      int iFB = aFSEdges[2*i + 1].iLeftCell();
      Face *pFA = aFSEdges[2*i + 0].pFFace();
      Vert *pVOther;
      if (pFA->pVVert(0) != pV0 && pFA->pVVert(0) != pV1)
	pVOther = pFA->pVVert(0);
      else if (pFA->pVVert(1) != pV0 && pFA->pVVert(1) != pV1)
	pVOther = pFA->pVVert(1);
      else {
	assert(pFA->pVVert(2) != pV0 && pFA->pVVert(2) != pV1);
	pVOther = pFA->pVVert(2);
      }
      int iOrient = iOrient3D(pVNew, pV0, pV1, 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("coplanar points (probably cocircular)");
      }
    }

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

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

  // 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(iOrient3D(pC->pVVert(0), pC->pVVert(1), pC->pVVert(2),
		     pC->pVVert(3)) == 1);
  }
  // assert(qIsLocallyDelaunay());
  SUMAA_LOG_EVENT_END(WATSON_CHECK);
#endif
  delete [] apCNew;
  SUMAA_LOG_EVENT_END(WATSON_INSERT);
}

class EdgeToPatch {
private:
  Vert *pV0, *pV1;
  BdryPatch* pBP;
public:
  EdgeToPatch() : pV0(pVInvalidVert), pV1(pVInvalidVert), pBP(NULL) {}
  EdgeToPatch(const EdgeToPatch& ETP) : pV0(ETP.pV0), pV1(ETP.pV1),
    pBP(ETP.pBP) {}
  EdgeToPatch(Vert * const pV0In, Vert * const pV1In,
	      BdryPatch * const pBPIn = NULL) :
    pV0(pV0In), pV1(pV1In), pBP(pBPIn)
    {
      assert(pV0In->qValid() && pV1In->qValid());
      if (pV0 > pV1) {
	Vert *pVTmp = pV1;
	pV1 = pV0;
	pV0 = pVTmp;
      }
    }
  EdgeToPatch& operator=(const EdgeToPatch& ETP)
    {
      if (this != &ETP) {
	pV0 = ETP.pV0;
	pV1 = ETP.pV1;
	pBP = ETP.pBP;
      }
      return (*this);
    }
  Vert *pVVert(const int i) const
    {
      assert(i == 0 || i == 1);
      return (i == 0) ? pV0 : pV1;
    }
  BdryPatch* pPatchPointer() const {return pBP;}
  bool operator!=(const EdgeToPatch& EP) const
    { return !(*this == EP); }
  bool operator==(const EdgeToPatch& EP) const
    {
      return ((pV0 == EP.pV0 && pV1 == EP.pV1) ||
	      (pV0 == EP.pV1 && pV1 == EP.pV0));
    }
  bool operator<(const EdgeToPatch& EP) const
    {
      return ((pV0 < EP.pV0) || (pV0 == EP.pV0 && pV1 < EP.pV1));
    }
};

//@ Add a vertex using Watson's method along the boundary
void VolMesh::vInsertWatsonBoundary(Vert* const pVNew,
				    std::set<Cell*>& spCToRemove,
				    std::set<Face*>& spFToRemove,
				    std::set<Face*>& spFHull,
				    const std::set<BFace*>& spBFEncroached)
{
  // Update for Periodic?
  SUMAA_LOG_EVENT_BEGIN(WATSON_BDRY_INSERT);
  int i;

  std::set<Face*>::iterator iterF;
  std::set<BFace*>::iterator iterBF;
  std::set<Cell*>::iterator iterC;

  bool qIsInternalBdry = false;
  {
    int iReg = (*spCToRemove.begin())->iRegion();
    for (iterC = spCToRemove.begin(); iterC != spCToRemove.end(); iterC++) {
      if ((*iterC)->iRegion() != iReg) {
	qIsInternalBdry = true;
	break;
      }
    }
  }

  // The difference between boundary and interior insertion is that when
  // inserting on the boundary, some of the hull faces are on the
  // boundary and mutually co-planar (or nearly so, for curved
  // surfaces).  These faces are deleted and replaced, as are the
  // BFace's attached to them.

  // To identify these faces, simply find all BFace's in the encroached
  // list that lie on the same BdryPatch as the BFace that triggered the
  // insertion (pBFTarget).  If this point is being inserted on a
  // boundary curve, then there will be two possible BdryPatch's.  These
  // have to be identified, and stored in spBPCand (by adding all
  // encroached BdryPatches, and then weeding out those not containing
  // the new vertex.

  std::set<BdryPatch*> spBPCand;
  for (iterBF = spBFEncroached.begin(); iterBF != spBFEncroached.end();
       iterBF++) {
    spBPCand.insert((*iterBF)->pPatchPointer());
  }
  // This loop is needed even for insertion on a patch interior, because
  // the new vertex may encroach on a triangle on some other patch.
  for (std::set<BdryPatch*>::iterator iterBP = spBPCand.begin();
       iterBP != spBPCand.end(); ) {
    BdryPatch *pBP = *iterBP;
    iterBP++;
    if (pBP->qDisjointFrom(pVNew)) {
      std::set<BdryPatch*>::iterator iterTmp(iterBP);
      spBPCand.erase(--iterTmp);
    }
  }
  assert(spBPCand.size() >= 1);
  vMessage(4, "Number of patches incident on new vert: %d\n",
	   static_cast<int>(spBPCand.size()));
  if (spBPCand.size() == 1) {
    if (qIsInternalBdry) {
      pVNew->vSetType(Vert::eBdryTwoSide);
    }
    else {
      pVNew->vSetType(Vert::eBdry);
    }
  }
  else {
    pVNew->vSetType(Vert::eBdryCurve);
  }

  std::set<BFace*> spBFToRemove;
  for (iterBF = spBFEncroached.begin(); iterBF != spBFEncroached.end();
       iterBF++) {
    BFace* pBF = *iterBF;
    if (spBPCand.count(pBF->pPatchPointer()) != 0) {
      // This bdry face must go.
      spBFToRemove.insert(pBF);
      // The face attached to the bdry face is removed from the hull,
      // because the hull, for this purpose, is the faces to which new
      // cells must be attached.   For internal bdry faces, always do
      // this for the faces on -both- sides.
      for (int ii = 0; ii < pBF->iNumFaces(); ii++) {
	spFHull.erase(pBF->pFFace(ii));
	spFToRemove.insert(pBF->pFFace(ii));
      }
    }
  }

  // For each hull face, identify the region in which the cell to be
  // deleted lies.
  int *aiRegion = new int[spFHull.size()];
  i = 0;
  for (iterF = spFHull.begin(); iterF != spFHull.end(); iterF++, i++) {
    Face *pF = *iterF;
    Cell *pC = pF->pCCellLeft();
    if (pC->eType() == Cell::eTet &&
	spCToRemove.count(pC) != 0) {
      assert(spCToRemove.count(pF->pCCellRight()) == 0);
    }
    else {
      pC = pF->pCCellRight();
      assert(pC->eType() == Cell::eTet);
      assert(spCToRemove.count(pC) != 0);
    }
    aiRegion[i] = pC->iRegion();
  }

  // Disconnect deleted cells from all their faces and mark for later
  // cleanup.
  for (iterC = spCToRemove.begin(); iterC != spCToRemove.end(); iterC++) {
    Cell *pC = *iterC;
    for (int iF = 0; iF < 4; iF++) {
      Face *pF = pC->pFFace(iF);
      pF->vRemoveCell(pC);
    }
    pC->vMarkDeleted();
  }

  // Ditto for the bdry faces that are being deleted.  For each edge,
  // keep track of the patch it lies on, so that the new bdry faces will
  // be attached to the correct patch.  Of course, many of these edges
  // will have a deleted BFace on each side of them, and the duplicates
  // aren't of interest to us, so they'll be removed on the fly.
  //
  // In the end, we have a list of all edges that lie on bdrys (internal
  // or regular bdrys) and also lie on the hull of faces to which new
  // cells will be attached.

  std::set<EdgeToPatch> sEPSingleEdges;
  for (iterBF = spBFToRemove.begin(); iterBF != spBFToRemove.end(); iterBF++) {
    BFace *pBF = *iterBF;
    // This assertion will have to be accompanied by one for the other
    // side for internal boundaries.
    Face *pF = pBF->pFFace(0);
    assert(spFToRemove.count(pF) != 0);
    assert(pBF->iNumFaces() == 1 || spFToRemove.count(pBF->pFFace(1)) != 0);
    pBF->vMarkDeleted();

    EdgeToPatch aEP[3];
    aEP[0] = EdgeToPatch(pF->pVVert(0), pF->pVVert(1),
			 pBF->pPatchPointer());
    aEP[1] = EdgeToPatch(pF->pVVert(1), pF->pVVert(2),
			 pBF->pPatchPointer());
    aEP[2] = EdgeToPatch(pF->pVVert(2), pF->pVVert(0),
			 pBF->pPatchPointer());
    for (int ii = 0; ii < 3; ii++) {
      if (sEPSingleEdges.count(aEP[ii]) != 0) {
	sEPSingleEdges.erase(aEP[ii]);
      }
      else {
	sEPSingleEdges.insert(aEP[ii]);
      }
    }
  }

  // Mark deleted faces for later cleanup.
  for (iterF = spFToRemove.begin(); iterF != spFToRemove.end(); iterF++) {
    (*iterF)->vMarkDeleted();
  }

  // Create a new cell for every hull face and connect to hull face.
  // Reset vertex hints for hull face verts.  Set the region for the
  // cell properly.
  int iNumHullFaces = spFHull.size();
  Cell **apCNew = new Cell*[iNumHullFaces];
  i = 0; 
  for (iterF = spFHull.begin(); iterF != spFHull.end(); iterF++, i++) {
    Face *pF = *iterF;
    Cell *pC = apCNew[i] = pCNewCell();
    pC->vClearFaces();
    pC->vAddFace(pF);
    pC->vSetRegion(aiRegion[i]);
    pF->vAddCell(pC);
    pF->pVVert(0)->vSetHintFace(pF);
    pF->pVVert(1)->vSetHintFace(pF);
    pF->pVVert(2)->vSetHintFace(pF);
  }

  // Now create a new face for every hull edge, and connect to new
  // cells.
  vSetUpEdgeList(spFHull);
  int iOldNumFaces = iNumFaces();
  for (i = 0; i < iNumHullFaces * 3; ) {
    Vert *pV0 = aFSEdges[i].pVVert(0);
    Vert *pV1 = aFSEdges[i].pVVert(1);
    int iFA = aFSEdges[i].iLeftCell();
    Face *pFA = aFSEdges[i].pFFace();
    Vert *pVOther;
    if (pFA->pVVert(0) != pV0 && pFA->pVVert(0) != pV1)
      pVOther = pFA->pVVert(0);
    else if (pFA->pVVert(1) != pV0 && pFA->pVVert(1) != pV1)
      pVOther = pFA->pVVert(1);
    else {
      assert(pFA->pVVert(2) != pV0 && pFA->pVVert(2) != pV1);
      pVOther = pFA->pVVert(2);
    }
    int iOrient = iOrient3D(pVNew, pV0, pV1, pVOther);
    if (iOrient == 0)
      vFoundBug("coplanar points (probably cocircular)");
    assert(abs(iOrient) == 1);

    if ((i < iNumHullFaces*3 - 1) &&
	((aFSEdges[i+1].pVVert(0) == pV0 && aFSEdges[i+1].pVVert(1) == pV1) ||
	 (aFSEdges[i+1].pVVert(0) == pV1 && aFSEdges[i+1].pVVert(1) == pV0))) {
      // This is an edge common to two hull faces, which implies that
      // it's going to have two new cells adjacent to it.

      // Now determine which face is on which side of the -new- face
      // pVNew, pV0, pV1.
      int iFB = aFSEdges[i + 1].iLeftCell();
      int iFRight, iFLeft;
      if (iOrient == 1) {
	// iFA is on the -right- side of the new face.
	iFRight = iFA;
	iFLeft  = iFB;
      }
      else {
	// iFA is on the -left- side of the new face.
	iFRight = iFB;
	iFLeft  = iFA;
      }

      // Are the cells attached to these faces in the same region?
      if (aiRegion[iFRight] == aiRegion[iFLeft]) {
	// This is the case whenever the two faces spawn cells in the
	// same region, regardless of whether -all- new cells are in a
	// single region.
	Face *pFNew = pFNewFace();
	pFNew->vAssign(apCNew[iFLeft], apCNew[iFRight], pVNew, pV0, pV1);
	// pFNew is an interior face by default.
	apCNew[iFLeft ]->vAddFace(pFNew);
	apCNew[iFRight]->vAddFace(pFNew);
	assert(pFNew->iFullCheck());
      }
      else {
	// These cells fall in different regions, so there must be an
	// internal boundary in here.

	// Need to construct two new faces and a new internal bdry face.
	Face *pFNewLeft  = pFNewFace();
	Face *pFNewRight = pFNewFace();
	BFace *pBFNew = pBFNewBFace(3, true);

	// Need to find the patch that this bdry face falls on; the
	// information is stored by edge.
	EdgeToPatch EP(pV0, pV1);
	std::set<EdgeToPatch>::iterator iterEP = sEPSingleEdges.find(EP);
	assert(iterEP != sEPSingleEdges.end());
	BdryPatch3D* pBP3D = dynamic_cast<BdryPatch3D*>
	  ((*iterEP).pPatchPointer());
	pBFNew->vSetPatch(pBP3D);

	pFNewLeft ->vAssign(apCNew[iFLeft ], pBFNew, pVNew, pV0, pV1);
	pFNewRight->vAssign(pBFNew, apCNew[iFRight], pVNew, pV0, pV1);
	pFNewLeft->vSetFaceLoc(Face::eBdryTwoSide);
	pFNewRight->vSetFaceLoc(Face::eBdryTwoSide);

	pBFNew->vAddFace(pFNewLeft);
	pBFNew->vAddFace(pFNewRight);

	apCNew[iFLeft ]->vAddFace(pFNewLeft);
	apCNew[iFRight]->vAddFace(pFNewRight);
      }
      i += 2; // Increment to the next candidate.
    } // Done with the case of an edge common to two new cells
    else {
      // This edge only appears once; a new BFace will have to be
      // created and attached to the mesh here.

      // Now determine which side of the -new- face pVNew, pV0, pV1 the
      // existing hull face is on.
      Face *pFNew = pFNewFace();
      BFace* pBF = pBFNewBFace();
      // Need to find the patch that this bdry face falls on; the
      // information is stored by edge.
      EdgeToPatch EP(pV0, pV1);
      std::set<EdgeToPatch>::iterator iterEP = sEPSingleEdges.find(EP);
      assert(iterEP != sEPSingleEdges.end());
      BdryPatch3D* pBP3D = dynamic_cast<BdryPatch3D*>
	((*iterEP).pPatchPointer());
      pBF->vSetPatch(pBP3D);

      if (iOrient == 1) {
	// iF is on the -right- side of the new face.
	pFNew->vAssign(pBF, apCNew[iFA], pVNew, pV0, pV1);
      }
      else {
	// iF is on the -left- side of the new face.
	pFNew->vAssign(apCNew[iFA], pBF, pVNew, pV0, pV1);
      }

      apCNew[iFA]->vAddFace(pFNew);
      pBF->vAddFace(pFNew);
      pFNew->vSetFaceLoc(Face::eBdryFace);
      assert(pFNew->iFullCheck());
      assert(pBF->iFullCheck());

      i++; // Increment to the next candidate.
    }
  } // Done checking all edges.
  // Set vertex hint for new vertex.
  pVNew->vSetHintFace(pFFace(iNumFaces() - 1));

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

  // Must swap all hull faces, because with boundary insertion, there
  // could be a face that -needs- to be swapped.  Its untouched cell has
  // a visibility problem with the inserted vertex, so the mesh can't be
  // made (even constrained) Delaunay by Watson insertion alone.
  int iNSwaps = 0;
  for (i = iNumFaces() - 1; i >= iOldNumFaces; i--) {
    iNSwaps += iFaceSwap(pFFace(i));
  }
  if (iNSwaps != 0)
    vMessage(4, "Made %d swaps for new faces (bdry insertion).\n",
	     iNSwaps);

  iNSwaps = 0;
  for (iterF = spFHull.begin(); iterF != spFHull.end(); iterF++) {
    iNSwaps += iFaceSwap(*iterF);
  }
  if (iNSwaps != 0)
    vMessage(3, "Made %d swaps for hull faces (bdry insertion).\n",
	     iNSwaps);

  // assert(qIsLocallyDelaunay());

  SUMAA_LOG_EVENT_END(WATSON_BDRY_INSERT);
}

