#include <set>
#include <map>

#include <stdlib.h>
#include "GR_ADT.h"
#include "GR_BFace.h"
#include "GR_Geometry.h"
#include "GR_Vec.h"
#include "GR_VolMesh.h"
#include "GR_events.h"
#include "SMsmooth.h"

#define iMaxVertNeigh MAX_NUM_PTS

// This is an ugly way to deal with these variables, but why create and
// destroy them every time?
static double **a2dNeighPts = NULL;
static int **a2iFaceConn = NULL;
static std::set<Cell*>  spCIncident;
static std::set<Vert*>  spVNeighbors;
static std::set<BFace*> spBFIncident;
static std::set<Face*>  spFNearby; // All faces in tets incident on pV.

extern int aiEdgeReq[], aiEdgeDone[];

int VolMesh::iSmoothVertex(Vert * const pV)
{
  SUMAA_LOG_EVENT_BEGIN(NET_SMOOTHING);
  bool qBdryVert;
  // Set up coordinates of the current vertex
  double adCurrPt[3];
  adCurrPt[0] = pV->dX();
  adCurrPt[1] = pV->dY();
  adCurrPt[2] = pV->dZ();

  if (pV->qDeleted() || pV->qWellShaped() || !pV->qSmoothable()) {
    SUMAA_LOG_EVENT_END(NET_SMOOTHING);
    return 0;
  }

  // If the vertex is moved, this marking is changed.
  pV->vMarkWellShaped();

  if (a2dNeighPts == NULL) {
    assert(a2iFaceConn == NULL);
    a2dNeighPts = new double*[iMaxVertNeigh];
    a2iFaceConn = new int*[2*iMaxVertNeigh];
    for (int i = 0; i < iMaxVertNeigh; i++) {
      a2dNeighPts[i] = new double[3];
      a2iFaceConn[i] = new int[3];
      a2iFaceConn[i+iMaxVertNeigh] = new int[3];
    }
  }

  // Set up coordinates of the current vertex
  adCurrPt[0] = pV->dX();
  adCurrPt[1] = pV->dY();
  adCurrPt[2] = pV->dZ();

  // Reset lists of nearby entities
  spCIncident.clear();
  spVNeighbors.clear();
  spBFIncident.clear();
  spFNearby.clear();

  vNeighborhood(pV, spCIncident, spVNeighbors, &spBFIncident,
		&qBdryVert, &spFNearby);

  // Verify that the right number of each kind of entity exists.
  assert(!qSimplicial() || ((!qBdryVert &&
			     spCIncident.size() == spVNeighbors.size()*2 - 4)
			    ||
			    (qBdryVert &&
			     spCIncident.size() == (spVNeighbors.size()*2 -
						    spBFIncident.size() - 2))
			    )
	 );

  if (qBdryVert) {
    // Can only smooth a bdry vertex if it lies in a coplanar patch
    // of the surface and if all its surrounding bdry faces have the
    // same boundary condition.  If these conditions are met, then
    // we need to set up the surface normal prior to smoothing.
    bool qSurfSmoothOK = true;
    std::set<BFace*>::iterator iterBF = spBFIncident.begin();
    TriBFaceBase* pTBF = dynamic_cast<TriBFaceBase*>(*iterBF);
    if (!pTBF->qValid()) {
      // This is a quad bface instead
      return 0;
    }
    int iBC = pTBF->iBdryCond();
    double adNorm[3];
    pTBF->pFFace()->vUnitNormal(adNorm);
    for (iterBF++ ; iterBF != spBFIncident.end() && qSurfSmoothOK; iterBF++) {
      TriBFace* pTBFTest = dynamic_cast<TriBFace*>(*iterBF);
      if (!pTBFTest->qValid()) {
	// This is a quad bface instead
	return 0;
      }
      if (pTBFTest->iBdryCond() != iBC)
	qSurfSmoothOK = false;
      else {
	double adNormTest[3];
	pTBFTest->pFFace()->vUnitNormal(adNormTest);
	double dDot = fabs(dDOT3D(adNorm, adNormTest));
	if (iFuzzyComp(dDot, 1.) != 0)
	  qSurfSmoothOK = false;
      }
    }
    if (qSurfSmoothOK && qAllowBdryChanges)
      SMsetNormal(adNorm);
    else {
      SUMAA_LOG_EVENT_END(NET_SMOOTHING);
      return 0;
    }
  }
  else if (pV->iVertType() == Vert::eBdryTwoSide) {
    // Can only smooth an internal bdry vertex if it lies in a coplanar
    // patch of the internal surface and if its incident cells fall into
    // only two regions.
    bool qSurfSmoothOK = true;

    // First check for region membership.
    int aiIncRegs[] = {iInvalidRegion, iInvalidRegion};
    std::set<Cell*>::iterator iterC = spCIncident.begin();
    aiIncRegs[0] = (*iterC)->iRegion();
    for (iterC++ ; iterC != spCIncident.end(); iterC++) {
      int iReg = (*iterC)->iRegion();
      if (iReg != aiIncRegs[0]) {
	if (aiIncRegs[1] == iInvalidRegion)
	  aiIncRegs[1] = iReg;
	else if (iReg != aiIncRegs[1]) {
	  // More than two regions incident; bail out now.
	  SUMAA_LOG_EVENT_END(NET_SMOOTHING);
	  return 0;
	}
      }
    }
    // Make sure there really are two regions in this case.
    assert(aiIncRegs[1] != iInvalidRegion);

    // Identify all faces that separate the two regions.
    std::set<TriFace*> spTFRegionBdry;
    std::set<Face*>::iterator iterF;
    for (iterF = spFNearby.begin(); iterF != spFNearby.end(); iterF++) {
      Face *pF = *iterF;
      if (pF->qHasVert(pV) &&
	  pF->iFaceLoc() == Face::eBdryTwoSide) {
	assert(pF->pCCellLeft()->iRegion() !=
	       pF->pCCellRight()->iRegion());
	if (pF->eType() != Face::eTriFace) return 0;
	spTFRegionBdry.insert(dynamic_cast<TriFace*>(pF));
      }
    }

    // Now check normals of all faces that separate the two regions.

    std::set<TriFace*>::iterator iterTF = spTFRegionBdry.begin();
    TriFace* pTF = *iterTF;
    double adNorm[3];
    pTF->vUnitNormal(adNorm);
    for (iterTF++; iterTF != spTFRegionBdry.end(); iterTF++) {
      TriFace* pTFTest = *iterTF;
      double adNormTest[3];
      pTFTest->vUnitNormal(adNormTest);
      double dDot = fabs(dDOT3D(adNorm, adNormTest));
      if (iFuzzyComp(dDot, 1.) != 0)
	qSurfSmoothOK = false;
    }
    if (qSurfSmoothOK && qAllowBdryChanges) {
      // Set stuff up so that the internal bdry point remains on the
      // internal bdry.
      SMsetNormal(adNorm);
      qBdryVert = true;
    }
    else {
      SUMAA_LOG_EVENT_END(NET_SMOOTHING);
      return 0;
    }
  }

  // If there are too many neighbors, punt
  if (spVNeighbors.size() > iMaxVertNeigh ||
      spCIncident.size() > 2*iMaxVertNeigh) {
    SUMAA_LOG_EVENT_END(NET_SMOOTHING);
    return 0;
  }

  // List the coordinates of all vertices
  int ii = 0;
  std::set<Vert*>::iterator iterV;
  std::map<Vert*, int> mVertIndex;
  for (iterV = spVNeighbors.begin(); iterV != spVNeighbors.end();
       iterV++, ii++) {
    Vert *pVTemp = *iterV;
    mVertIndex.insert(std::pair<Vert*, int>(pVTemp, ii));
    a2dNeighPts[ii][0] = pVTemp->dX();
    a2dNeighPts[ii][1] = pVTemp->dY();
    a2dNeighPts[ii][2] = pVTemp->dZ();
  }

  // List all opposite faces (w/ correct orientation) using indices
  // into the location array
  std::set<Cell*>::iterator iterC;
  for (iterC = spCIncident.begin(), ii = 0;
       iterC != spCIncident.end(); iterC++, ii++) {
    Cell *pC = *iterC;
    if (pC->eType() != Cell::eTet) return 0;
    Face *pF = dynamic_cast<TetCell*>(pC)->pFFaceOpposite(pV);
    if (pF->pCCellLeft() == pC) {
      a2iFaceConn[ii][0] = (mVertIndex.find(pF->pVVert(0)))->second;
      a2iFaceConn[ii][1] = (mVertIndex.find(pF->pVVert(2)))->second;
      a2iFaceConn[ii][2] = (mVertIndex.find(pF->pVVert(1)))->second;
    }
    else {
      a2iFaceConn[ii][0] = (mVertIndex.find(pF->pVVert(0)))->second;
      a2iFaceConn[ii][1] = (mVertIndex.find(pF->pVVert(1)))->second;
      a2iFaceConn[ii][2] = (mVertIndex.find(pF->pVVert(2)))->second;
    }
  }

  double adPrevPt[] = {adCurrPt[0], adCurrPt[1], adCurrPt[2]};
  SMsmooth(spVNeighbors.size(),spCIncident.size(),
	   adCurrPt,a2dNeighPts,a2iFaceConn,
	   pvSmoothData, (qBdryVert ? 1 : 0) );
  double dDistMoved = dDIST3D(adCurrPt, adPrevPt);

  // Only mark this vertex and its neighborhood for future
  // improvement if the worst angle in the submesh is less than 30
  // degrees or more than 150 degrees.
  // FIX ME: 0.3.0  Someday when (if?) you can routinely hit 30/150
  // degrees, make this value a floating threshold.
  double dQualValue =
    (static_cast<SMsmooth_data*>(pvSmoothData))
    ->local_mesh->current_active_value;
  SMconvertToDegrees(static_cast<SMsmooth_data*>(pvSmoothData)->smooth_param
		     ->function_id, &dQualValue);
  if (dQualValue < 30) {
    pV->vMarkIllShaped();
    // Also, mark all faces in incident tets as dirty.
    std::set<Face*>::iterator iterF;
    for (iterF = spFNearby.begin(); iterF != spFNearby.end(); iterF++)
      (*iterF)->vMarkForSwap();
  }

  SUMAA_LOG_EVENT_END(NET_SMOOTHING);
  if ( dDistMoved > 1.e-8 ) {
    // Make sure that the neighbor is optimally placed.
    for (iterV = spVNeighbors.begin(); iterV != spVNeighbors.end(); iterV++)
      (*iterV)->vMarkIllShaped();
    pV->vSetCoords(3, adCurrPt);
    return 1;
  }
  else
    return 0;
}

//@ Optimization driver for 3D mesh smoothing
int VolMesh::iSmooth(const int iMaxPasses)
{
  SUMAA_LOG_EVENT_BEGIN(SMOOTHING_PASSES);
  int iSmooths = 1;
  int iTotal = 0;
  int iPasses = 0;

  int i;

  while (iSmooths && iPasses < iMaxPasses) {
    SMinitSmoothStats(pvSmoothData);
    ++iPasses;
    iSmooths = 0;
    for (i = 0; i < iNumVerts(); i++) {
      Vert *pV = pVVert(i);
      iSmooths += iSmoothVertex(pV);
    }

    SMprintSmoothStats(pvSmoothData);
    vResetThreshold();

    iTotal += iSmooths;
    vMessage(1, "   Pass %d:  %5d vertices moved by smoothing, %6d total.\n",
	     iPasses, iSmooths, iTotal);
  }

  vMessage(2, "Vertex smoothing complete.\n");
  SUMAA_LOG_EVENT_END(SMOOTHING_PASSES);
  return iSmooths;
}

//@ Optimization driver for face swapping.  Swap measure is set elsewhere.
int VolMesh::iSwap(const int iMaxPasses, const bool qAlreadyMarked)
{
  SUMAA_LOG_EVENT_BEGIN(SWAP_PASSES);
  SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
  int iSwaps = 1;
  assert(eSwapMeasure != eNone);

  if (!qAlreadyMarked) {
    for (int iFace = 0; iFace < iNumFaces(); iFace++)
      pFFace(iFace)->vMarkForSwap();
  }

  int iTotal = 0;
  int iPasses = 0;

  while (iSwaps && iPasses < iMaxPasses) {
    ++ iPasses;
    if (iMaxPasses > 1)
      vMessage(2, "   Pass %d:", iPasses);
    iSwaps = 0;
    for (int i = 0; i < iNumTriFaces(); i++) {
      Face *pF = pFFace(i);
      // Don't swap if the face no longer exists or if it doesn't need
      // to be swapped.
      if (!pF->qDeleted() && pF->qCheckForSwap()) {
	// Always swap recursively
	iSwaps += iFaceSwap(pF);
      }
    }
    iTotal += iSwaps;
    vMessage(2, "%6d swaps", iSwaps);
    if (iMaxPasses > 1)
      vMessage(2, ",%6d total", iTotal);
    vMessage(2, "\n");
  }
  if (iMaxPasses > 1) vMessage(2, "Face swapping complete.\n");
  for (int ii = 2; ii < 10; ii++) {
    vMessage(3, "    Edge swapping from %d tets to %2d: %5d done out of %5d requests\n",
	     ii, 2*ii-4, aiEdgeDone[ii], aiEdgeReq[ii]);
  }

  SUMAA_LOG_EVENT_END(NET_SWAPPING);

  if (!(qValid()))
    vFoundBug("3D mesh reconnection");
  SUMAA_LOG_EVENT_END(SWAP_PASSES);
  vPurge();
  return iTotal;
}

static bool qEdgeHitsBFace(const Vert * const pV0, const Vert * const pV1,
			   const ADT& BFaceTree, const VolMesh& VM)
{
  double adBBox[] = {min(pV0->dX(), pV1->dX()),
		     max(pV0->dX(), pV1->dX()),
		     min(pV0->dY(), pV1->dY()),
		     max(pV0->dY(), pV1->dY()),
		     min(pV0->dZ(), pV1->dZ()),
		     max(pV0->dZ(), pV1->dZ())};
  std::vector<int> veciQueryResult = BFaceTree.veciRangeQuery(adBBox);
  if (veciQueryResult.size() >= 20) {
    vMessage(3, "Warning: query returned %d possible bdry faces.\n",
	     static_cast<int>(veciQueryResult.size()));
  }

  for (int ii = veciQueryResult.size() - 1; ii >= 0; ii--) {
    Face *pF = VM.pBFBFace(veciQueryResult[ii])->pFFace();
    // Now check whether the edge actually penetrates the face.  If so,
    // return true.  (Actually, return this even if the edge hits an
    // edge or corner of the face.)  For the face to be okay (edge
    // outside) then there must be one orientation of each sign (one
    // plus and one minus).  This is so regardless of which sign is the
    // inside orientation and which is the outside orientation.
    bool qHavePlus = false, qHaveMinus = false;
    switch (iOrient3D(pF->pVVert(0), pF->pVVert(1), pV0, pV1)) {
      case 1: qHavePlus = true; break;
      case -1: qHaveMinus = true; break;
      default: break;
    }
    switch (iOrient3D(pF->pVVert(1), pF->pVVert(2), pV0, pV1)) {
      case 1: qHavePlus = true; break;
      case -1: qHaveMinus = true; break;
      default: break;
    }
    switch (iOrient3D(pF->pVVert(2), pF->pVVert(0), pV0, pV1)) {
      case 1: qHavePlus = true; break;
      case -1: qHaveMinus = true; break;
      default: break;
    }
    if (! (qHavePlus && qHaveMinus)) return true;
  }
  return false;
}

static bool qIsVisible(const Cell * const pC, const Vert * const pV,
		       const ADT& BFaceTree, const VolMesh& VM)
{
  assert(pC->eType() == Cell::eTet);
  // Visibility is defined here as each vertex of the tet being visible
  // to the vertex inside the circumsphere.  Technically, one could have
  // situations where part of an edge was obscured, for example, but not
  // any vertex.  The general case would be, I think, very difficult to
  // test correctly.
  for (int iV = 0; iV < 4; iV++) {
    const Vert *pVSource = pC->pVVert(iV);
    if (qEdgeHitsBFace(pV, pVSource, BFaceTree, VM)) return false;
  }
  return true;
}

static int iIdentifyNonDelaunayCells(const VolMesh& VM,
				     bool aqDelaunay[])
{
  int iRetVal = 0;

  // Create a vertex tree so that insphere searching can be done
  // efficiently.
  ADT VT(VM.pVAGetVerts());

  // Also create a tree of bdry faces to make visibility checking
  // easier.
  double (*a2dBFace)[6] = new double[VM.iNumBdryFaces()][6];
  for (int iBFace = 0; iBFace < VM.iNumBdryFaces(); iBFace++) {
    BFace *pBF = VM.pBFBFace(iBFace);

    a2dBFace[iBFace][0] = MIN3(pBF->pVVert(0)->dX(), pBF->pVVert(1)->dX(),
			       pBF->pVVert(2)->dX());
    a2dBFace[iBFace][1] = MAX3(pBF->pVVert(0)->dX(), pBF->pVVert(1)->dX(),
			       pBF->pVVert(2)->dX());
    a2dBFace[iBFace][2] = MIN3(pBF->pVVert(0)->dY(), pBF->pVVert(1)->dY(),
			       pBF->pVVert(2)->dY());
    a2dBFace[iBFace][3] = MAX3(pBF->pVVert(0)->dY(), pBF->pVVert(1)->dY(),
			       pBF->pVVert(2)->dY());
    a2dBFace[iBFace][4] = MIN3(pBF->pVVert(0)->dZ(), pBF->pVVert(1)->dZ(),
			       pBF->pVVert(2)->dZ());
    a2dBFace[iBFace][5] = MAX3(pBF->pVVert(0)->dZ(), pBF->pVVert(1)->dZ(),
			       pBF->pVVert(2)->dZ());
  }
  ADT BFaceTree(VM.iNumBdryFaces(), 3, ADT::eBBoxes,
		reinterpret_cast<double*>(a2dBFace));
  delete [] a2dBFace;

  // Now check each cell to see how many verts fall into its
  // circumsphere.  Keep track of the max as well.
  int iMaxEncroach = 0;
  int *aiEncroach = new int[VM.iNumCells()];
  std::vector<int> veciQueryResult;
  double dEps = 1.e-10;
  int iCell;
  for (iCell = 0; iCell < VM.iNumCells(); iCell++) {
    Cell *pCCell = VM.pCCell(iCell);
    assert(pCCell->eType() == Cell::eTet);
    TetCell *pTC = dynamic_cast<TetCell*>(pCCell);
    double adCircCent[3], dCircRad;
    dCircRad = pTC->dCircumradius();
    pTC->vCircumcenter(adCircCent);
    // Get all verts in the bounding box of the circumsphere
    double adRange[] = {adCircCent[0] - dCircRad - dEps,
			adCircCent[0] + dCircRad + dEps,
			adCircCent[1] - dCircRad - dEps,
			adCircCent[1] + dCircRad + dEps,
			adCircCent[2] - dCircRad - dEps,
			adCircCent[2] + dCircRad + dEps};
    veciQueryResult = VT.veciRangeQuery(adRange);
    assert(veciQueryResult.size() >= 4);
    Vert *pV0 = pTC->pVVert(0);
    Vert *pV1 = pTC->pVVert(1);
    Vert *pV2 = pTC->pVVert(2);
    Vert *pV3 = pTC->pVVert(3);
    int iNumEnc = 0;
    for (int ii = veciQueryResult.size() - 1; ii >= 0; ii --) {
      Vert *pV = VM.pVVert(veciQueryResult[ii]);
      if ((pV == pV0) || (pV == pV1) || (pV == pV2) || (pV == pV3))
	continue;
      int iIsEncroached = iInsphere(pTC->pVVert(0), pTC->pVVert(1),
				    pTC->pVVert(2), pTC->pVVert(3), pV);
#ifndef NDEBUG
      if (iIsEncroached != 1) {
	// On or outside the sphere
	const double *adLoc = pV->adCoords();
	double dTmpDist = dDIST3D(adLoc, adCircCent);
	assert (iFuzzyComp(dTmpDist, dCircRad) >= 0);
      }
#endif
      if (iIsEncroached == 1 && qIsVisible(pTC, pV, BFaceTree, VM)) {
	iNumEnc ++;
      }
    }
    aiEncroach[iCell] = iNumEnc;
    if (iNumEnc == 0) {
      aqDelaunay[iCell] = true;
    }
    else {
      iRetVal ++;
    }
    iMaxEncroach = max(iNumEnc, iMaxEncroach);
  }
  vMessage(1, "Maximum number of extra verts in a circumsphere is %d\n",
	   iMaxEncroach);
  vMessage(1, "Num verts encroaching      Num tets\n");
  int *aiCount = new int[iMaxEncroach+1];
  for (int ii = 0; ii < iMaxEncroach+1; aiCount[ii++] = 0) {}
  for (iCell = 0; iCell < VM.iNumCells(); iCell++)
    aiCount[aiEncroach[iCell]]++;
  for (int ii = 0; ii <= iMaxEncroach; ii++)
    if (aiCount[ii] != 0)
      vMessage(1, "%12d%20d\n", ii, aiCount[ii]);
  delete [] aiEncroach;
  delete [] aiCount;
  return (iRetVal);
}

bool VolMesh::qDelaunayize()
{
  assert(qSimplicial());
  // First, do insphere face swapping; this often gives a Delaunay mesh
  // right away.
  vAllowSwapRecursion();
  vSetSwapType(eDelaunay);
  iSwap();

  int iOldNumNonDelaunay = 0, iNumNonDelaunay;
  // Identify non-Delaunay tets.
  bool *aqDelaunay = new bool[iNumCells()];
  iNumNonDelaunay = iIdentifyNonDelaunayCells(*this, aqDelaunay);

  while (iNumNonDelaunay != 0 && iNumNonDelaunay != iOldNumNonDelaunay) {
    // Tag all verts as not in need of smoothing; for those in
    // non-Delaunay tets, this will change.
    for (int iVert = 0; iVert < iNumVerts(); iVert++) {
      pVVert(iVert)->vMarkWellShaped();
    }

    // Tag (only) the verts of non-Delaunay tets for smoothing.
    for (int iCell = 0; iCell < iNumCells(); iCell++) {
      if (!aqDelaunay[iCell]) {
	Cell *pC = pCCell(iCell);
	pC->pVVert(0)->vMarkIllShaped();
	pC->pVVert(1)->vMarkIllShaped();
	pC->pVVert(2)->vMarkIllShaped();
	pC->pVVert(3)->vMarkIllShaped();
      }
    }

    // Smooth, then swap only the affected faces.
    vMessage(2, "Smoothing...\n");
    iSmooth(1);
    vMessage(2, "Swapping...\n");
    iSwap(5, true);

    // Identify non-Delaunay tets.
    iOldNumNonDelaunay = iNumNonDelaunay;
    iNumNonDelaunay = iIdentifyNonDelaunayCells(*this, aqDelaunay);
    vMessage(2, "Had %d tets non-Delaunay last time; now %d\n",
	     iOldNumNonDelaunay, iNumNonDelaunay);

    // While still making progress, repeat.
  }

  int iFace;
  Vert *pVVertA, *pVVertB, *pVVertC, *pVVertD, *pVVertE;
  Vert *pVPivot0, *pVPivot1, *pVOther;
  TetCell* apTCTets[4];
  int iNTets;
  int iN44Cases, iAttempts, iSuccesses, iPasses = 0;
  do {
    iN44Cases = iAttempts = iSuccesses = 0;
    iPasses++;
    for (iFace = 0; iFace < iNumFaces(); iFace++) {
      Face *pF = pFFace(iFace);
      if (pF->qDeleted() || pF->qIsBdryFace()) continue;
      int iCL = iCellIndex(pF->pCCellLeft());
      int iCR = iCellIndex(pF->pCCellRight());
      if (aqDelaunay[iCL] || aqDelaunay[iCR]) continue;
      assert(pF->eType() == Face::eTriFace);
      eFaceCat eFC =
	dynamic_cast<TriFace*>(pF)->
	eCategorizeFace(pVVertA, pVVertB, pVVertC,
			pVVertD, pVVertE, apTCTets, iNTets,
			pVPivot0, pVPivot1, pVOther);
      if (eFC == eOther || eFC == eBdry) continue;
      int iPointInCircsphere = iInsphere(pVVertA, pVVertB, pVVertC,
					 pVVertD, pVVertE);
      switch (eFC) {
      default:
	break;
      case eN44:
	iN44Cases++;
      case eN32:
	if (iPointInCircsphere == 1) {
	  assert(pVVertD != pVPivot0 && pVVertD != pVPivot1 &&
		 pVVertD != pVOther);
	  assert(pVVertE != pVPivot0 && pVVertE != pVPivot1 &&
		 pVVertE != pVOther);
	  vMessage(2, "Attempting to eliminate edge (%4d %4d) in favor of (%4d %4d)\n",
		   iVertIndex(pVPivot0)+1, iVertIndex(pVPivot1)+1,
		   iVertIndex(pVVertD)+1, iVertIndex(pVVertE)+1);
	  int iRes = iEdgeSwap3D(pF, pVPivot0, pVPivot1, pVOther, true,
				 pVVertD, pVVertE);
	  iAttempts++;
	  if (iRes > 0) {
	    iSuccesses++;
	  }
	  break;
	}
      }
    }
    vMessage(2, "Number of N44 cases attempted: %d\n", iN44Cases);
    vMessage(2, "Succeeded on %d of %d attempts to remove pivot edges in N32 and N44 faces\n", iSuccesses, iAttempts);
  } while (iSuccesses > 0 && iPasses < 10);
  vPurge();

  delete [] aqDelaunay;
  return (iNumNonDelaunay == 0);
}
