#include "GR_assert.h"
#include "GR_events.h"
#include "GR_BFace.h"
#include "GR_Geometry.h"
#include "GR_Vec.h"
#include "GR_VolMesh.h"

#include <map>
#include <vector>
#include <algorithm>
using std::multimap;
using std::pair;
using std::set;

struct RemCell {
  Cell *pC;
  Face *pFOnlyV, *pFOnlyV1;
  RemCell() : pC(pCInvalidCell), pFOnlyV(pFInvalidFace),
	      pFOnlyV1(pFInvalidFace) {}
  RemCell(Cell * const pCIn, Face * const pF0, Face * const pF1)
    : pC(pCIn), pFOnlyV(pF0), pFOnlyV1(pF1) {}
};

static int iDelReq = 0, iBdryDelReq = 0, iBdryCurveDelReq = 0;
static int iBdryApexDelReq = 0;
static int iDelDone = 0, iBdryDelDone = 0, iBdryCurveDelDone = 0;

static int iStarRemeshFailed = 0, iStarRemeshSucceeded = 0,
  iStarRestoreFailed = 0;

// class FakeTet : public TetCell {
// private:
//   /// Copy construction disallowed.
//   FakeTet(const FakeTet&) : TetCell() {assert(0);};
//   Vert *apV[4];
// public:
//   eCellType eType() const {return eFakeTet;};
//   FakeTet(Vert* const pV0, Vert* const pV1,
//           Vert* const pV2, Vert* const pV3) : TetCell(NULL)
//     {
//        apV[0] = pV0; apV[1] = pV1; apV[2] = pV2; apV[3] = pV3;
//     };
//   Vert * pVVert(const int i) const {assert(i>=0 && i<4); return apV[i];};
// };

// To remove a given vertex V from a mesh, there must be a vertex V1
// adjacent to V which is visible from all other vertices adjacent to
// V.  (Note that this specifically does NOT require that the
// neighborhood be convex.)  If this requirement is met, then all faces
// and all cells incident on both V and V1 can be removed from the mesh.
// The faces in the deleted cells which are incident on only one of V
// and V1 occur in pairs and can be matched up and reconnected to ensure
// a valid mesh on output.

bool VolMesh::qRemoveVertByContraction(Vert * const pV, int& iNewSwaps,
				       Vert * apVPreferred[], const int iNPref)
  // Remove vertex *pV from the mesh.  Decisions about whether *pV
  // SHOULD be removed are made elsewhere.  The vertices listed in
  // apVPreferred are given a weak preference in determining to which
  // vertex pV should be moved.  This preference takes the form of
  // checking these vertices first in making removal choices.
{
  // Only works for tetrahedral meshes at this point.
  assert(qSimplicial());
  iNewSwaps = 0;
  iDelReq++;

  if (pV->iVertType() == Vert::eBdryApex) {
    vMessage(3, "Deletion of a boundary apex requested! Denied.\n");
    iBdryApexDelReq++;
    return false;
  }

  if (pV->iVertType() == Vert::eBdryTwoSide) {
    vMessage(4, "Deletion of verts on internal boundaries not yet supported.\n");
    return false;
  }

  // List all incident cells on *pV and all vertices in those cells.  As
  // a side effect, find out whether this is a boundary vertex.
  std::set<Cell*> spCIncident;
  std::set<Vert*> spVCand;
  std::set<BFace*> spBFIncident;
  std::set<Face*> spFUpdate;

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

  bool qBdryVert = false;
  bool qBdryCurveVert = false;

  vNeighborhood(pV, spCIncident, spVCand, &spBFIncident, &qBdryVert,
		&spFUpdate);
  
  if (iNPref > 0) {
    // If there are any preferred verts, they will be the only ones
    // considered as contraction targets.  If there are none, proceed as
    // usual.
    spVCand.clear();
    spVCand.insert(apVPreferred, apVPreferred + iNPref);
  }
  
  std::set<Vert*> spVAllBdryCands;
  if (qBdryVert) {
    // BdryApex calls have already been bounced.
    assert(pV->iVertType() == Vert::eBdryCurve ||
	   pV->iVertType() == Vert::eBdryTwoSide ||
	   pV->iVertType() == Vert::eBdry);
    
    // Abort if boundaries are precious
    if (!qBdryChangesAllowed()) return false;
    
    // Permanently disqualify all vertices which don't meet boundary rules
    // for appropriate contraction.  Those rules are:
    //   1.  BdryApex verts never get removed.
    //   2.  BdryCurve verts only go onto BdryCurve and BdryApex verts
    //       that fall on the same curve (i.e., the connecting edge must
    //       lie on a curve, so the two BFaces shared by the two verts
    //       must lie on different patches.
    //   3.  Bdry verts that are structured only go onto BdryCurve or
    //       BdryApex verts, or onto kept, structured Bdry verts.
    //   4.  Other Bdry verts can go onto any flavor of bdry vert.
    
    // First, find all of the faces incident on pV (faces in spFUpdate)
    // that are boundary faces.
    spVCand.clear();
    for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();
	 iterBF++) {
      BFace *pBF = *iterBF;
      Face *pFBdry = pBF->pFFace();
      assert(pFBdry->qHasVert(pV));
      spVCand.insert(pFBdry->pVVert(0));
      spVCand.insert(pFBdry->pVVert(1));
      spVCand.insert(pFBdry->pVVert(2));
    }
    spVCand.erase(pV);
    assert(spVCand.size() == spBFIncident.size());
    spVAllBdryCands.insert(spVCand.begin(), spVCand.end());
    
    switch (pV->iVertType()) {
    case Vert::eBdryCurve:
      {
	qBdryCurveVert = true;
	for (iterV = spVCand.begin(); iterV != spVCand.end(); ) {
	  Vert *pVCand = *iterV;
	  std::set<Vert*>::iterator iterTmp = iterV++;
	  // Eliminate all verts that aren't bdry curves or bdry apexes.
	  if (pVCand->iVertType() != Vert::eBdryApex &&
	      pVCand->iVertType() != Vert::eBdryCurve) {
	    spVCand.erase(iterTmp);
	  }
	  // Now for each remaining candidate, find the pair of bdry faces
	  // that connect it to the vertex to be removed, and eliminate
	  // any candidates that aren't on the same bdry curve.
	  else {
	    BFace* apBFTemp[2];
	    int ii = 0;
	    for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();
		 iterBF++) {
	      BFace *pBFCand = *iterBF;
	      assert(pBFCand->qHasVert(pV));
	      if (pBFCand->qHasVert(pVCand))
		apBFTemp[ii++] = pBFCand;
	    }
	    assert(ii == 2);
	    // The following code will be useful only when patch
	    // optimization is implemented.  Until then, adjacent
	    // triangles in a planar region are officially in different
	    // patches.
	    //	    if (apBFTemp[0]->pPatchPointer() ==
	    //		apBFTemp[1]->pPatchPointer()) {
	    //	      spVCand.vDeleteByIndex(iV);
	    //	    }
	    // The following code will be obsolete when patch
	    // optimization is implemented, because the definitive patch
	    // check (above) will work correctly then.
	    BFace *pBF0 = apBFTemp[0];
	    BFace *pBF1 = apBFTemp[1];
	    double adNorm0[3], adNorm1[3];
	    pBF0->vUnitNormal(adNorm0);
	    pBF1->vUnitNormal(adNorm1);
	    double dDot = dDOT3D(adNorm0, adNorm1);
	    if (dDot > cos(dMaxAngleForSurfSwap * M_PI/180)) {
	      spVCand.erase(iterTmp);
	    }
	  } // Done with the case where the candidate vertex is on a
	    // bdry curve somewhere
	} // Done checking all candidate verts when the vert to be
	  // removed lies on a bdry
      } // end case
      iBdryCurveDelReq++;
      break;
    case Vert::eBdryTwoSide:
    case Vert::eBdry:
      {
	// Any candidate verts must be connected to the vert to be
	// removed by an edge that actually lies on the boundary.
	for (iterV = spVCand.begin(); iterV != spVCand.end(); ) {
	  Vert *pVCand = *iterV;
	  std::set<Vert*>::iterator iterTmp = iterV++;
	  // For each candidate, find the pair of bdry faces that
	  // connect it to the vertex to be removed, and eliminate any
	  // candidates that aren't on the same bdry curve.
	  int ii = 0;
	  for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();
	       iterBF++) {
	    BFace *pBFCand = *iterBF;
	    assert(pBFCand->qHasVert(pV));
	    if (pBFCand->qHasVert(pVCand))
	      ii++;
	  }
	  if (ii != 2) {
	    vMessage(2, "Blocking possible attempt to remove a bdry vert along a bad edge.\n");
	    vMessage(2, "# bdry faces shared by verts: %d\n",
		     ii);
	    spVCand.erase(iterTmp);
	  }
	} // Done checking all candidate verts when the vert to be
	  // removed lies on a bdry
	
	  // If the vertex is structured, then it should only be slid onto
	  // another structured vertex or a vert on a curve or apex.
	if (pV->qIsStructured()) {
	  for (iterV = spVCand.begin(); iterV != spVCand.end(); ) {
	    Vert *pVCand = *iterV;
	    std::set<Vert*>::iterator iterTmp = iterV++;
	    if (! (pVCand->qIsStructured() ||
		   pVCand->iVertType() == Vert::eBdryApex ||
		   pVCand->iVertType() == Vert::eBdryCurve) ) {
	      spVCand.erase(iterTmp);
	    }
	  }
	}
      }
      iBdryDelReq++;
      break;
    case Vert::eBdryApex:
    default:
      // Should never get here.
      assert(0);
    }
  }

  // List all faces in incident cells that are opposite *pV
  std::vector<Face*> vecpFOpposite;
  for (iterC = spCIncident.begin(); iterC != spCIncident.end(); iterC++) {
    Cell *pC = *iterC;
    assert(pC->eType() == Cell::eTet);
    Face *pFOpp = dynamic_cast<TetCell*>(pC)->pFFaceOpposite(pV);
    vecpFOpposite.push_back(pFOpp);
    spFUpdate.erase(pFOpp);
  }
  
  // Used only in assertions just before the actual deletion.
#ifndef NDEBUG
  int iTotalIncFaces = spFUpdate.size();
#endif
  
  // For each neighboring vertex, determine whether it is a valid target
  // for edge contraction.  Keep track of which faces affect which
  // vertices.  If no vertices are okay, start with the one with the
  // fewest interfering faces and try to remove that face from the hull
  // by using edge removal to break the connection from *pV to any of
  // the vertices of the face.  If that succeeds, call this routine
  // recursively to do the contraction; otherwise continue until an
  // attempt has been made to remove each problematic face.
  
  bool qVertOK;
  bool qAnyVertOK = false;
  double dMinDihed;
  double dMaxMinDihed = 0.; // Force a minimum dihedral angle before
  // attempting to remove a vertex.
  Vert *pV1 = pVInvalidVert, *pVBest = pVInvalidVert;
  for (iterV = spVCand.begin(); iterV != spVCand.end(); iterV++) {
    pV1 = *iterV;
    qVertOK = true;
    dMinDihed = 180;
    for (unsigned int iF = 0; iF < vecpFOpposite.size() && qVertOK; iF++) {
      Face *pF = vecpFOpposite[iF];
      if (!pF->qHasVert(pV1)) {
	// Determine whether a given face has the correct orientation
	// for the candidate contraction target.  If the orientation
	// test fails, pV1 can't see the correct side of the face pF.
	bool qResult;
	if ( (pF->pCCellRight())->qHasVert(pV)) {
	  // For this case, the right-hand rule applies for orientation.
	  assert(iOrient3D(pF, pV) == 1);
	  qResult = (iOrient3D(pF, pV1) == 1);
	  if (qResult) {
	    Vert *apV[] = {pF->pVVert(0), pF->pVVert(1), pF->pVVert(2), pV1};
	    TetCellCV FT(apV);
	    double dMin = FT.dMinDihed();
	    double dMax = FT.dMaxDihed();
	    dMinDihed = min(dMinDihed, min(dMin, 180-dMax));
	    vMessage(4, "Min dihed (face): %10.6f  (vert): %10.6f\n",
		     FT.dMinDihed(), dMinDihed);
	  }
	}
	else {
	  // Otherwise, the left-hand rule applies for orientation.
	  assert((pF->pCCellLeft())->qHasVert(pV));
	  assert(iOrient3D(pF, pV) == -1);
	  qResult = (iOrient3D(pF, pV1) == -1);
	  if (qResult) {
	    Vert *apV[] = {pF->pVVert(1), pF->pVVert(0), pF->pVVert(2), pV1};
	    TetCellCV FT(apV);
	    double dMin = FT.dMinDihed();
	    double dMax = FT.dMaxDihed();
	    dMinDihed = min(dMinDihed, min(dMin, 180-dMax));
	    vMessage(4, "Min dihed (face): %10.6f  (vert): %10.6f\n",
		     FT.dMinDihed(), dMinDihed);
	  }
	}
	if (!qResult) {
	  qVertOK = false;
	}
	if (dMinDihed < dMaxMinDihed) {
	  qVertOK = false;
	}
      } // Done checking a face
    } // Done checking a vertex

    if (qVertOK && qBdryVert) {
      // For bdry verts, make certain that there aren't three (or
      // more) boundary verts shared by pV and pV1.  If there are,
      // then contracting this edge will collapse two boundary faces
      // onto each other.  While having BFaces collapse
      // "front-to-front" (interior sides together) is recoverable,
      // the "back-to-back" collapse is not recoverable.  And even the
      // front-to-front collapse may potentially delete a vertex that
      // wasn't tagged for deletion, so it isn't a great idea anyway.

      // This check is made easier because spVCand already holds a
      // list of all bdry verts incident on pV.

      std::set<Vert*> spVNeigh1;
      std::set<Cell*> spCJunk;
      vNeighborhood(pV1, spCJunk, spVNeigh1);
      
      assert(spVNeigh1.count(pV) == 1);
      assert(spVCand.count(pV1) == 1);
      assert(spVCand.count(pV) == 0);
      
      int iCommon = 0;
      for (std::set<Vert*>::iterator iterV2 = spVNeigh1.begin();
	   iterV2 != spVNeigh1.end(); iterV2++) {
	Vert *pVTest = *iterV2;
	if (pVTest->qIsBdryVert() && spVAllBdryCands.count(pVTest) == 1)
	  iCommon++;
      }
      
      switch (iCommon) {
      case 0:
      case 1:
	assert(0);
	vFoundBug("bizarre result for bdry vertex removal");
      case 2:
	// Everything's okay.
	break;
      default:
	// Too many bdry verts incident!
	vMessage(4, "%d common bdry vert neighbors in bdry edge contraction; 2 is max allowed\n", iCommon);
	qVertOK = false;
      } // Done with switch
    } // Done with check for weird bdry circumstances.

    if (qVertOK) {
      if (dMinDihed > dMaxMinDihed) {
	// The following line has been moved inside the 'if' so that the
	// oddball case of adMinDihed[] == 0 can't cause a problem.
	// Hey, it happened once!
	qAnyVertOK = true;
	vMessage(4, "Best min dihedral yet.  Old: %10.6f  New: %10.6f\n",
		 dMaxMinDihed, dMinDihed);
	dMaxMinDihed = dMinDihed;
	pVBest = pV1;
      }
      if (dMinDihed > 0) {
	// Check to see whether this is a preferred vertex
	for (int iPref = 0; iPref < iNPref; iPref++) {
	  if (apVPreferred[iPref] == pV1) {
	    pVBest = pV1;
	    break;
	  }
	}
      } // Done with preferred vert check
    } // Done updating best vertex
  } // Done checking all candidate targets for edge contraction

  if (!qAnyVertOK) {
    return false;
  }

  // If we get to here, we know we CAN remove *pV by contracting edge
  // (*pV, *pVBest), and that this gives the largest minimum dihedral
  // angle of any possible contraction.
  assert(qAnyVertOK);
  assert(pVBest->qValid());
  // If the best you can do is really lousy, don't do it.  This doesn't
  // apply for structured verts, where traditional measures of mesh
  // quality are (at best) difficult to interpret.
  // Change angle threshhold here!
  if (dMaxMinDihed < 3. && !pV->qIsStructured()) return false;
  vMessage(3, "Contracting vert %5d onto %5d; min dihed after: %5.2f deg.\n",
	   iVertIndex(pV), iVertIndex(pVBest), dMaxMinDihed);
  vMessage(3, "Contracting vert %5d at (%10g, %10g, %10g)\n",
	   iVertIndex(pV), pV->dX(), pV->dY(), pV->dZ());
  vMessage(3, "       onto vert %5d at (%10g, %10g, %10g)\n",
	   iVertIndex(pVBest), pVBest->dX(), pVBest->dY(), pVBest->dZ());

  // Remove from the list of cells incident on *pV those not also
  // incident on *pVBest, to get a list of cells which will be deleted.
  for (iterC = spCIncident.begin(); iterC != spCIncident.end(); ) {
    Cell *pC = *iterC;
    std::set<Cell*>::iterator iterTmp = iterC++;
    if ((!pC->qHasVert(pV)) || (!pC->qHasVert(pVBest)))
      spCIncident.erase(iterTmp);
  }

  if (qBdryVert) {
    // Remove from the list of bfaces incident on *pV those not also
    // incident on *pVBest, to get a list of bfaces which will be deleted.
    for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();) {
      BFace *pBF = *iterBF;
      std::set<BFace*>::iterator iterTmp = iterBF++;
      if ((!pBF->qHasVert(pV)) || (!pBF->qHasVert(pVBest))) {
	spBFIncident.erase(iterTmp);
      }
    }
    assert(spBFIncident.size() == 2);
  }

  // List all faces incident on both *pV and *pVBest.  First list all faces
  // in the incident cells, then remove those which aren't incident on
  // both *pV and *pVBest.
  std::set<Face*> spFLateral;
  for (iterC = spCIncident.begin(); iterC != spCIncident.end(); iterC++) {
    Cell *pC = *iterC;
    for (int iF = 0; iF < pC->iNumFaces(); iF++)
      spFLateral.insert(pC->pFFace(iF));
  }
  for (iterF = spFLateral.begin(); iterF != spFLateral.end(); ) {
    Face *pF = *iterF;
    std::set<Face*>::iterator iterTmp = iterF++;
    if ((!pF->qHasVert(pV)) || (!pF->qHasVert(pVBest)))
      spFLateral.erase(iterTmp);
    else
      // These faces are going to be deleted from the mesh, so we won't
      // update them to reflect replacing *pV by *pVBest.
      spFUpdate.erase(pF);
  }

  // Produce matched lists of the faces of the incident cells which
  // contain only one of *pV and *pVBest, and a coordinated list of
  // vertices in these faces such that the vertices match up with faces
  // that contain them.
  int iNInc = spCIncident.size();
  std::vector<struct RemCell> vecRC(iNInc);
  int iC;
  for (iterC = spCIncident.begin(), iC = 0;
       iterC != spCIncident.end(); iterC++, iC++) {
    assert ((*iterC)->eType () == Cell::eTet);
    TetCell *pTC = dynamic_cast<TetCell *>(*iterC);
    vecRC[iC].pC = pTC;
    Face *pFOpp = vecRC[iC].pFOnlyV  = pTC->pFFaceOpposite(pVBest);
    vecRC[iC].pFOnlyV1 = pTC->pFFaceOpposite(pV);
    // These faces are going to be deleted from the mesh, so we won't
    // update them to reflect replacing *pV by *pVBest.
    spFUpdate.erase(pFOpp);
  }

  // Verify that for interior contraction the lists are all the same
  // length or that for boundary contraction, there is one more lateral
  // face than there is of anything else.
  if (qBdryVert) {
    assert(spCIncident.size() + 1 == spFLateral.size());
    assert(spFUpdate.size() + spFLateral.size() +
	   spCIncident.size() == unsigned(iTotalIncFaces));
  }
  else {
    assert(spCIncident.size()     == spFLateral.size());
    assert(spFUpdate.size() + spFLateral.size() +
	   spCIncident.size() == unsigned(iTotalIncFaces));
  }

  // Delete the faces incident on both *pV and *pVBest.
  for (iterF = spFLateral.begin(); iterF != spFLateral.end(); iterF++) {
    Face *pF = *iterF;
    assert(pF->qHasVert(pV));
    assert(pF->qHasVert(pVBest));
    pF->vMarkDeleted();
  }

  // Change connectivity so that cells remaining in the mesh pointing to
  // faces incident only on *pV will point to faces incident on *pVBest.
  // Then delete faces incident only on *pV.
  std::set<Vert*> spVResetHints;
  for (int iF = 0; iF < iNInc; iF++) {
    Face *pF = vecRC[iF].pFOnlyV;
    Face *pF1 = vecRC[iF].pFOnlyV1;
    assert (pF != pF1);
    assert (vecRC[iF].pC->qHasFace (pF));
    assert (vecRC[iF].pC->qHasFace (pF1));
    // Grab the cell that is still connected to the mesh.
    Cell *pCInc = vecRC[iF].pC;
    Cell *pC = pF->pCCellLeft ();
    if (pCInc == pC) {
      pC = pF->pCCellRight ();
      assert (pCInc != pC);
    }
    else {
      // Never assume anything!
      assert(pF->pCCellRight() == pCInc);
    }

    pC->vReplaceFace (pF, pF1);
    pF1->vReplaceCell (pCInc, pC);

    if (pF1->pCCellLeft()->eType() == Cell::eTriBFace &&
	pF1->pCCellRight()->eType() == Cell::eTriBFace) {
      vMessage(2, "Collapsed two BFace's onto each other!\n");
      vMessage(2, "Vert %d (type %d) is sliding onto vert %d (type %d).\n",
	       iVertIndex(pV), pV->iVertType(), iVertIndex(pVBest),
	       pVBest->iVertType());
      vMessage(2, "Verts in the collapsed bdry faces:\n");
      vMessage(2, "  Left  BFace: %6d, %6d, %6d\n",
	       iVertIndex(pF1->pCCellLeft()->pVVert(0)),
	       iVertIndex(pF1->pCCellLeft()->pVVert(1)),
	       iVertIndex(pF1->pCCellLeft()->pVVert(2)));
      vMessage(2, "  Right BFace: %6d, %6d, %6d\n",
	       iVertIndex(pF1->pCCellRight()->pVVert(0)),
	       iVertIndex(pF1->pCCellRight()->pVVert(1)),
	       iVertIndex(pF1->pCCellRight()->pVVert(2)));
      vMessage(2, "There were originally %d BFace's incident on both verts.\n",
	       static_cast<int>(spBFIncident.size()));
      vMessage(2, "Trying to fix the resulting connectivity mess.\n");

      // Remove pF1 and the two BFaces incident on it.  Queue all verts
      // in question for possible hint updates.
      for (int ii = 0; ii < 3; ii++) {
	Vert *pVTmp = pF1->pVVert(ii);
	pVTmp->vClearHintFace();
	spVResetHints.insert(pVTmp);
      }

      pF1->vMarkDeleted();
      pF1->pCCellRight()->vMarkDeleted();
      pF1->pCCellLeft()->vMarkDeleted();

    } // Done with face collapse cleanup.
#ifndef NDEBUG
    else {
      assert(pC->iFullCheck() == 1);
      assert(pF1->iFullCheck() == 1);
      Cell *pCTest = pF1->pCCellLeft();
      assert(pCTest->iFullCheck() == 1);
      pCTest = pF1->pCCellRight();
      assert(pCTest->iFullCheck() == 1);
      assert(pF1->pCCellLeft() == pC || pF1->pCCellRight() == pC);
      assert(pC->qHasFace(pF1));
    }
#endif
    // Delete old faces incident only on *pV.
    pF->vMarkDeleted();
  }

  // Change face data so that faces which did point to *pV now point to
  // *pVBest.
#ifndef NDEBUG
  std::set<Cell*> spCCheck;
#endif
  for (iterF = spFUpdate.begin(); iterF != spFUpdate.end(); iterF++) {
    Face *pF = *iterF;
    assert( pF->qHasVert(pV ));
    assert(!pF->qHasVert(pVBest));
    pF->vReplaceVert(pV, pVBest);
#ifndef NDEBUG
    spCCheck.insert(pF->pCCellLeft());
    spCCheck.insert(pF->pCCellRight());
#endif
    assert( pF->qHasVert(pVBest));
    assert(!pF->qHasVert(pV ));
    assert(pF->pVVert(0) != pF->pVVert(1));
    assert(pF->pVVert(1) != pF->pVVert(2));
    assert(pF->pVVert(2) != pF->pVVert(0));
    assert(pF->iFullCheck());
  }
  
#ifndef NDEBUG
  for (iterC = spCCheck.begin(); iterC != spCCheck.end(); iterC++) {
    Cell *pC = *iterC;
    assert(pC->qHasVert(pVBest));
    assert(!pC->qHasVert(pV));
    assert(pC->iFullCheck() == 1);
    assert(pC->eType() == Cell::eTriBFace ||
	   (iOrient3D(pC->pVVert(0), pC->pVVert(1), pC->pVVert(2),
		      pC->pVVert(3)) == 1
	    && pC->qIsClosed()));
    int iTest;
    for (iTest = 0; iTest < pC->iNumVerts(); iTest++)
      assert(pC->pVVert(iTest) != pC->pVVert((iTest+1) % pC->iNumVerts()));
    if (pC->iNumFaces() > 1)
      for (iTest = 0; iTest < pC->iNumFaces(); iTest++)
	assert(pC->pFFace(iTest) != pC->pFFace((iTest+1) % pC->iNumFaces()));
  }
#endif

  // Delete the cells that were incident on both *pV and *pVBest.
  for (iterC = spCIncident.begin(); iterC != spCIncident.end(); iterC++) {
    (*iterC)->vMarkDeleted();
  }

  if (qBdryVert) {
    // Delete the bdry faces that were incident on both *pV and *pVBest.
    for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();
	 iterBF++) {
      (*iterBF)->vMarkDeleted();
    }
  }

  // Delete *pV.
  pV->vMarkDeleted();

  // Reset the vertex hints for *pVBest and all the peripheral vertices.
  // Some of these hints won't have been broken, but setting them all is
  // easier to code than checking whether it needs to be done.
  for (int iF = 0; iF < iNInc; iF++) {
    Face *pF1 = vecRC[iF].pFOnlyV1;
    if (pF1->qDeleted()) continue;
    for (int iV = 0; iV < 3; iV++) {
      Vert *pVTmp = pF1->pVVert(iV);
      pVTmp->vSetHintFace(pF1);
      pVTmp->vMarkIllShaped();
    }
  }

  for (iterV = spVResetHints.begin(); iterV != spVResetHints.end(); iterV++) {
    Vert *pVTmp = *iterV;
    Face *pFNewHint = pVTmp->pFSetHintFace(this);
    if (!pFNewHint->qValid()) {
      vMessage(2, "Vertex %d left hanging; deleting.\n", iVertIndex(pVTmp));
      pVTmp->vMarkDeleted();
      iDelDone ++;
    }
  }

#ifndef NDEBUG
  // Here's a post-check to be sure that bdry faces haven't collapsed
  // onto each other.  For now, it just spews output.
  if (qBdryVert) {
    std::set<Cell*> spCJunk;
    std::set<Vert*> spVJunk;
    std::set<BFace*> spBFInc;
    bool qBdryChk = false;
    vNeighborhood(pVBest, spCJunk, spVJunk, &spBFInc, &qBdryChk);

    assert(qBdryChk);
    spVJunk.clear();
    vMessage(3, "pVBest has %d incident BFaces.\n",
	     static_cast<int>(spBFInc.size()));
    for (iterBF = spBFInc.begin(); iterBF != spBFInc.end(); iterBF++) {
      BFace *pBF = *iterBF;
      vMessage(3, "  Bdryface:%6d.  Verts:%7d%7d%7d\n",
	       iBFaceIndex(pBF), iVertIndex(pBF->pVVert(0)),
	       iVertIndex(pBF->pVVert(1)), iVertIndex(pBF->pVVert(2)));
      spVJunk.insert(pBF->pVVert(0));
      spVJunk.insert(pBF->pVVert(1));
      spVJunk.insert(pBF->pVVert(2));
    }
    spVJunk.erase(pVBest);
    assert(spVJunk.size() == spBFInc.size());
  }
#endif

  // Try swapping all the faces that were changed.
  // Helpful for 3D isotropic mesh quality; makes no difference for 3D
  // anisotropic removal success, which is a little surprising.
  for (iterF = spFUpdate.begin(); iterF != spFUpdate.end(); iterF++) {
    iNewSwaps += iFaceSwap(*iterF);
  }

  // Success!
  iDelDone ++;
  if (qBdryCurveVert) {
    iBdryCurveDelDone ++;
  }
  else if (qBdryVert) {
    iBdryDelDone ++;
  }
  return true;
}

// class FrontFace : public FaceSort {
//   // Since FrontFace is publicly a FaceSort object, storage and
//   // retrieval of vertex data is trivial.  Also, the FaceSort operator<
//   // is convenient.
// private:
//   mutable Face *pF;
// public:
//   // Gotta have this default constructor to keep picky compilers happy.
//   FrontFace() : FaceSort(), pF(pFInvalidFace) {}
//   FrontFace(Face *pFIn) :
//     FaceSort(pFIn->pVVert(0), pFIn->pVVert(1), pFIn->pVVert(2), -1),
//     pF(pFIn)
//     {
//       assert(pF != NULL);
//       assert(pF->eType() == Face::eTriFace);
//     }
//   FrontFace(Vert * const pV0In, Vert * const pV1In, Vert * const pV2In) :
//     FaceSort(pV0In, pV1In, pV2In, -1), pF(pFInvalidFace) {}
//   ~FrontFace() {}
//   Face *pFFace() const {return pF;}
//   void vStoreFace(Face * const pFIn) const
//     {
//       assert(pFIn->qValid());
//       pF = pFIn;
//       assert(pF->qHasVert(pVVert(0)));
//       assert(pF->qHasVert(pVVert(1)));
//       assert(pF->qHasVert(pVVert(2)));
//     }
// //   bool operator<(const FrontFace& FF) const {
// //     return (this->FaceSort::operator<(FF));
// //   }
//   bool qVertInFrontOfFace(const Vert* const pV) const
//     {return iOrient3D(pVVert(0), pVVert(1), pVVert(2), pV) == 1;}
// };

class StarEdge {
private:
  Vert *pV0, *pV1;
  Face *pFL, *pFR, *pFMid;
public:
  StarEdge(Vert *pV0In = pVInvalidVert, Vert *pV1In = pVInvalidVert,
	   Face* pFLIn = pFInvalidFace, Face* pFRIn = pFInvalidFace,
	   Face* pFMidIn = pFInvalidFace) :
    pV0(pV0In), pV1(pV1In), pFL(pFLIn), pFR(pFRIn), pFMid(pFMidIn) {}
  Vert *pVVert(int i) const {
    assert(i == 0 || i == 1);
    assert(pV0->qValid());
    assert(pV1->qValid());
    return (i == 0 ? pV0 : pV1);
  }
  Face *pFFaceLeft() const {
    assert(pFL->qValid());
    return pFL;
  }
  Face *pFFaceRight() const {
    assert(pFR->qValid());
    return pFR;
  }
  bool qHasFace(Face *pF) const {
    return (pFR == pF || pFL == pF);
  }
  void vReplaceFace(Face *pFOld, Face *pFNew) {
    assert(pFOld->qValid());
    assert(pFNew->qValid());
    if (pFOld == pFL)
      pFL = pFNew;
    else {
      assert(pFOld == pFR);
      pFR = pFNew;
    }
  }
  Face *pFFaceMiddle() const {
    assert(pFMid->qValid());
    return pFMid;
  }
  double dPower(const double adDel[3]) const {
    assert(pFL->qValid());
    assert(pFR->qValid());
    assert(pV0->qValid());
    assert(pV1->qValid());
    // First, ID all the verts of the hypothetical tet
    Vert *pV2 = pFL->pVVert(0);
    if (pV2 == pV0 || pV2 == pV1) {
      pV2 = pFL->pVVert(1);
      if (pV2 == pV0 || pV2 == pV1) {
	pV2 = pFL->pVVert(2);
      }
    }
    assert(pV2 != pV0 && pV2 != pV1);

    Vert *pV3 = pFR->pVVert(0);
    if (pV3 == pV0 || pV3 == pV1) {
      pV3 = pFR->pVVert(1);
      if (pV3 == pV0 || pV3 == pV1) {
	pV3 = pFR->pVVert(2);
      }
    }
    assert(pV3 != pV0 && pV3 != pV1 && pV3 != pV2);

    // Now verify orientation.  If this test fails, this ridge in the
    // star surface is reflex.
    if (iOrient3D(pV0, pV1, pV2, pV3) !=
	iOrient3D(pV0->adCoords(), pV1->adCoords(), pV2->adCoords(), adDel))
      return -LARGE_DBL;

    // Find the circumcenter of the candidate tet and calculate the
    // square of the distance between adDel and adCent, minus the square
    // of the circumradius.  This is the power measure for the candidate
    // tet.
    Vert *apV[] = {pV0, pV1, pV2, pV3};
    TetCellCV FT(apV);
    double adCent[3];
    FT.vCircumcenter(adCent);
    double dDist2 = dDIST_SQ_3D(adCent, adDel);
    double dRad2  = dDIST_SQ_3D(adCent, pV0->adCoords());
    return (dDist2 - dRad2);
  }
  bool operator<(const StarEdge& SE) const {
    Vert *pVMinThis, *pVMaxThis, *pVMinThat, *pVMaxThat;
    if (pV0 < pV1) {
      pVMinThis = pV0; pVMaxThis = pV1;
    }
    else {
      pVMinThis = pV1; pVMaxThis = pV0;
    }
    if (SE.pV0 < SE.pV1) {
      pVMinThat = SE.pV0; pVMaxThat = SE.pV1;
    }
    else {
      pVMinThat = SE.pV1; pVMaxThat = SE.pV0;
    }
    return ((pVMinThis < pVMinThat) ||
	    ((pVMinThis == pVMinThat) && (pVMaxThis < pVMaxThat)));
  }
  bool operator==(const StarEdge& SE) const {
    return ( (((pV0 == SE.pV0) && (pV1 == SE.pV1) &&
	       (pFL == SE.pFL) && (pFR == SE.pFR)) ||
	      ((pV0 == SE.pV1) && (pV1 == SE.pV0) &&
	       (pFL == SE.pFR) && (pFR == SE.pFL)))
	     && pFMid == SE.pFMid );
  }
};

class BdryEdge {
private:
  Vert *pV0, *pV1;
  Face *pFStar;
  BdryPatch *pBP;
public:
  BdryEdge(Vert *pV0In = pVInvalidVert, Vert *pV1In = pVInvalidVert,
	   Face* pFIn = pFInvalidFace, BdryPatch* pBPIn = NULL) :
    pV0(pV0In), pV1(pV1In), pFStar(pFIn), pBP(pBPIn) {}
  Vert *pVVert(int i) const {
    assert(i == 0 || i == 1);
    assert(pV0->qValid());
    assert(pV1->qValid());
    return (i == 0 ? pV0 : pV1);
  }
  Face *pFStarFace() const {
    assert(pFStar->qValid());
    return pFStar;
  }
  BdryPatch *pBPPatchPointer() const {
    assert(pBP->qValid());
    return pBP;
  }
};

static void vUpdateDataStructs(Vert* pV0, Vert* pV1,
			       Face* pFOld, Face* pFNew,
			       const double adLoc[], set<StarEdge>& sSE,
			       multimap<double, StarEdge>& mRM)
{
  set<StarEdge>::iterator iter;
  double dPower;
  iter = sSE.find(StarEdge(pV0, pV1));
  // For bdry verts, it's possible for the alleged star edge to not be a
  // star edge because it lies on the bdry.
  if (iter != sSE.end()) {
    StarEdge SETmp = *iter;
    sSE.erase(iter);
    SETmp.vReplaceFace(pFOld, pFNew);
    sSE.insert(SETmp);
    dPower = SETmp.dPower(adLoc);
    if (dPower > -0.9*LARGE_DBL)
      mRM.insert(pair<double, StarEdge> (-dPower, SETmp));
  }
}

bool VolMesh::qRemoveVert(Vert * const pVDelete, int& iNumSwaps)
{
  assert(qSimplicial());

  iNumSwaps = 0;
  if (pVDelete->iVertType() == Vert::eBdryApex) return false;

  // Find out everything about the neighborhood of the vertex.
  std::set<Cell*>::iterator iterC;
  std::set<Face*>::iterator iterF;
  std::set<BFace*>::iterator iterBF;

  std::set<Face*> spFUpdate;
  std::set<BFace*> spBFIncident;
  bool qBdryVert;
  unsigned int uNNeighV, uNIncCell, uNIncBF;
  {
    std::set<Vert*> spVNeigh;
    std::set<Cell*> spCIncident;
    vNeighborhood(pVDelete, spCIncident, spVNeigh, &spBFIncident, &qBdryVert,
		  &spFUpdate);
    uNNeighV = spVNeigh.size();
    uNIncCell = spCIncident.size();
    uNIncBF = spBFIncident.size();
  }

  int iFinalFaceCount = 0;
  switch (pVDelete->iVertType()) {
  case Vert::eBdryCurve:
    assert(uNIncCell >= 2);
    iFinalFaceCount = 2;
    break;
  case Vert::eBdry:
    assert(uNIncCell >= 3);
    iFinalFaceCount = 3;
    break;
  case Vert::eInterior:
    assert(uNIncCell >= 4);
    iFinalFaceCount = 4;
    break;
  default:
  case Vert::eBdryApex:
    // This catches eUnknown, eBdryTwoSide, ePseudoSurface, and eInvalid
    // as well as the well-advertised eBdryApex.  None of these cases
    // can be coarsened, and two of them should never even occur in a
    // mesh.  One (ePseudoSurface) has never been used, although vague
    // plans do exist for it in the future context of anisotropic meshing.
    assert(0);
    break;
  }

  bool qSuccess = true;
  bool qOldRecursionSetting = qSwapRecursionAllowed();
  vDisallowSwapRecursion();

  // This set contains all the faces on the surface of the star
  set<Face*> spFStarFaces;

  // This map makes it trivial to identify the cell adjacent to a face
  // of the original star that lies outside the star; essential for easy
  // recovery if vertex removal fails.
  std::map<Face*, Cell*> mStarFaceRecovery;

  // This set contains connectivity info for the star (edge -> face)
  set<StarEdge> sSEStarEdges;

  // This map holds the priority queue for tets to create; not all star
  // edges will be present in the queue (only those that aren't reflex
  // are queued), and some entries will correspond to star edges that
  // have been modified (queue entries are cleaned lazily).
  multimap<double, StarEdge> mRidgeMap;

  // Go through all the incident faces.  This makes it easy to find all
  // the star edges and their adjacent star faces in one go.
  for (iterF = spFUpdate.begin(); iterF != spFUpdate.end(); iterF++) {
    Face *pF = *iterF;
    if (pF->qHasVert(pVDelete) && !pF->qIsBdryFace()) {
      // This face has a star edge associated with it.
      Vert *pV0 = pF->pVVert(0);
      Vert *pV1 = pF->pVVert(1);
      Vert *pV2 = pF->pVVert(2);
      Vert *pVA, *pVB;
      if (pV0 == pVDelete) {
	pVA = pV1; pVB = pV2;
      }
      else if (pV1 == pVDelete) {
	pVA = pV2; pVB = pV0;
      }
      else {
	assert(pV2 == pVDelete);
	pVA = pV0; pVB = pV1;
      }
      Cell *pCL = pF->pCCellLeft();
      Cell *pCR = pF->pCCellRight();
      Face *pFL = dynamic_cast<TetCell*>(pCL)->pFFaceOpposite(pVDelete);
      Face *pFR = dynamic_cast<TetCell*>(pCR)->pFFaceOpposite(pVDelete);

      // These faces are in the star surface
      spFStarFaces.insert(pFL);
      spFStarFaces.insert(pFR);

      // Create the star edge
      StarEdge SE(pVA, pVB, pFL, pFR, pF);
      sSEStarEdges.insert(SE);

      // Create a durable copy of the original star face data.  Every
      // StarFace gets added at least once,
      assert(!pFL->pCCellOpposite(pCL)->qHasVert(pVDelete));
      assert(!pFR->pCCellOpposite(pCR)->qHasVert(pVDelete));
      mStarFaceRecovery.insert(pair<Face*, Cell*>
			       (pFL, pFL->pCCellOpposite(pCL)));
      mStarFaceRecovery.insert(pair<Face*, Cell*>
			       (pFR, pFR->pCCellOpposite(pCR)));

      // Find its power value and enqueue
      double dPower = SE.dPower(pVDelete->adCoords());
      if (dPower > -0.9 * LARGE_DBL)
	mRidgeMap.insert(pair<double, StarEdge>(-dPower, SE));
    }
  }
  assert(spFStarFaces.size() == uNIncCell);
  assert(spFStarFaces.size() == 2*uNNeighV - 4 ||
	 (qBdryVert &&
	  spFStarFaces.size() == (2*uNNeighV - 2 - uNIncBF)));
  assert(sSEStarEdges.size() == 3*uNNeighV - 6 ||
	 (qBdryVert &&
	  sSEStarEdges.size() == (3*uNNeighV - 3 - 2*uNIncBF)));

  // This array holds edge recovery info; no associative access is
  // needed, and the size is known in advance, so use an array.
  std::vector<StarEdge> vSERecovery(3*uNNeighV - 6);
  std::vector<StarEdge>::iterator SEiter =
    std::copy(sSEStarEdges.begin(), sSEStarEdges.end(), vSERecovery.begin());
  assert(int(SEiter - vSERecovery.begin()) <= int(3*uNNeighV - 6));

  std::vector<BdryEdge> aBERecovery(uNIncBF);
  int iBF = 0;
  for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();
       iterBF++, iBF++) {
    BFace *pBF = *iterBF;
    Face *pF = pBF->pFFace();
    TetCell *pTC = dynamic_cast<TetCell*>(pF->pCCellOpposite(pBF));
    Face *pFStar = pTC->pFFaceOpposite(pVDelete);
    BdryPatch* pBP = pBF->pPatchPointer();

    Vert *pV0, *pV1, *pV2;
    pV0 = pBF->pVVert(0);
    pV1 = pBF->pVVert(1);
    pV2 = pBF->pVVert(2);

    Vert *pVBegin, *pVEnd;
    if (pV0 == pVDelete) {
      pVBegin = pV1; pVEnd = pV2;
    }
    else if (pV1 == pVDelete) {
      pVBegin = pV2; pVEnd = pV0;
    }
    else {
      assert(pV2 == pVDelete);
      pVBegin = pV0; pVEnd = pV1;
    }
    aBERecovery[iBF] = BdryEdge(pVBegin, pVEnd, pFStar, pBP);
  }

  // Until the star shrinks to the point where we can do it all in one go...
  while (qSuccess && (spFStarFaces.size() > unsigned(iFinalFaceCount))
	 && !mRidgeMap.empty()) {
    // Grab the next ridge
    StarEdge SE = (*mRidgeMap.begin()).second;
    vMessage(4, "Ridge map has a size of %d; best score = %f\n",
	     static_cast<int>(mRidgeMap.size()), (*mRidgeMap.begin()).first);
    mRidgeMap.erase(mRidgeMap.begin());

    // If this ridge doesn't actually match an edge remaining in the
    // star surface, then its neighborhood has been modified since this
    // entry was enqueued, which means that the swap originally implied
    // by this entry can't be made.  So skip it and go on.  (The reason
    // for the SE != SECurrent check is that operator< is less
    // discriminating in checking for matches than operator==
    StarEdge SECurrent = *sSEStarEdges.find(SE);
    if (!(SE == SECurrent)) continue;

    // Identify the face defined by the ridge and the vertex to be
    // removed.
    TriFace *pFMid = dynamic_cast<TriFace*>(SE.pFFaceMiddle());
    Vert *pVA, *pVB, *pVC, *pVD, *pVE, *pVPivot0, *pVPivot1, *pVOther;
    TetCell* apTCTets[4];
    int iNTets;

    eFaceCat eFC = pFMid->
      eCategorizeFace(pVA, pVB, pVC, pVD, pVE,
		      apTCTets, iNTets, pVPivot0, pVPivot1, pVOther,
		      dMaxAngleForSurfSwap);
    // Depending on classification, process the face
    switch (eFC) {
    case eT23:
      // Same number of everything in the star afterwards, although
      // volume inside the star shrinks.
      {
	// In this case, verts pVA, pVB, and pVC are the verts of face
	// pFMid, although there's no guarantee about order or even
	// orientation.  Rotate things if necessary to get pVDelete ==
	// pVC.  The new edge in the star surface will be pVD-pVE, and
	// pVA-pVD, pVD-pVB, pVB-pVE, pVE-pVA will need to be updated
	// (each has a new star face next to it).  Also, all but pVD-pVE
	// will need to have their queue entries updated / replaced.
	iNumSwaps ++;
	// Rotate so that pVDelete is pVC.
	if (pVDelete == pVA) {
	  pVA = pVB;
	  pVB = pVC;
	  pVC = pVDelete;
	}
	else if (pVDelete == pVB) {
	  pVB = pVA;
	  pVA = pVC;
	  pVC = pVDelete;
	}
	else {}
	assert(pVC == pVDelete);

	// The old star faces that you'll be deleting are associated
	// with this ridge in the StarEdge.
	Face *pFFaceD, *pFFaceE;
	if (SE.pFFaceLeft()->qHasVert(pVD)) {
	  pFFaceD = SE.pFFaceLeft();
	  pFFaceE = SE.pFFaceRight();
	}
	else {
	  pFFaceE = SE.pFFaceLeft();
	  pFFaceD = SE.pFFaceRight();
	}
	assert(pFFaceD->qHasVert(pVD));
	assert(pFFaceE->qHasVert(pVE));
	spFStarFaces.erase(pFFaceD);
	spFStarFaces.erase(pFFaceE);

	// Now remove this StarEdge from the master list
	sSEStarEdges.erase(SE);

	// Force the reconfiguration
	int iRes = iReconfTet23(apTCTets, pFMid, pVA, pVB, pVC, pVD, pVE);
	assert(iRes == 1);

	// Now find the new stuff.  pFMid should be incident on pVDelete
	// still, so start there and work out.
	assert(pFMid->qValid());
	assert(pFMid->qHasVert(pVDelete));
	TetCell *pTCLeft  = dynamic_cast<TetCell*>(pFMid->pCCellLeft());
	TetCell *pTCRight = dynamic_cast<TetCell*>(pFMid->pCCellRight());
	Face *pFLeft  = pTCLeft ->pFFaceOpposite(pVDelete);
	Face *pFRight = pTCRight->pFFaceOpposite(pVDelete);
	assert(pFLeft->qHasVert(pVB));
	assert(pFLeft->qHasVert(pVD));
	assert(pFLeft->qHasVert(pVE));
	assert(pFRight->qHasVert(pVA));
	assert(pFRight->qHasVert(pVD));
	assert(pFRight->qHasVert(pVE));

	// Add two new faces to star face set
	pair<set<Face*>::iterator, bool> pRes =
	  spFStarFaces.insert(pFLeft);
	assert(pRes.second);
	pRes =
	  spFStarFaces.insert(pFRight);
	assert(pRes.second);

	// Add edge  to the star surface
	StarEdge SENew(pVD, pVE, pFLeft, pFRight, pFMid);
	sSEStarEdges.insert(SENew);

	// Update four other star edges, including new queue entries
	vUpdateDataStructs(pVA, pVD, pFFaceD, pFRight,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVD, pVB, pFFaceD, pFLeft,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVB, pVE, pFFaceE, pFLeft,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVE, pVA, pFFaceE, pFRight,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
      }
      break;
    case eT32:
      // Two fewer incident tets, two fewer star faces, three fewer star edges.
      {
	// In this case, verts pVA, pVB, and pVC are the verts on the
	// "equator" of the three-tet configuration.  No advance
	// guarantees on which end is vertex pVDelete, so check on that.
	// No new edges in the star surface.  Edges pVA-pVB, pVB-pVC,
	// and pVC-pVA will need to be updated (each has a new star face
	// next to it), including their queue entries.
	iNumSwaps ++;

	assert(pVDelete == pVD || pVDelete == pVE);
	bool qDelIsD = (pVDelete == pVD);
	Vert *pVOtherPole = qDelIsD ? pVE : pVD;

	// The old star faces that you'll be deleting are incident on
	// pVOtherPole.  More specifically, they're the faces opposite
	// pVDelete in the three tets in the list.
	Face *pFFaceBC = apTCTets[0]->pFFaceOpposite(pVDelete);
	Face *pFFaceCA = apTCTets[1]->pFFaceOpposite(pVDelete);
	Face *pFFaceAB = apTCTets[2]->pFFaceOpposite(pVDelete);
	assert(pFFaceBC->qHasVert(pVOtherPole));
	assert(pFFaceBC->qHasVert(pVB));
	assert(pFFaceBC->qHasVert(pVC));
	assert(pFFaceCA->qHasVert(pVOtherPole));
	assert(pFFaceCA->qHasVert(pVC));
	assert(pFFaceCA->qHasVert(pVA));
	assert(pFFaceAB->qHasVert(pVOtherPole));
	assert(pFFaceAB->qHasVert(pVA));
	assert(pFFaceAB->qHasVert(pVB));
	spFStarFaces.erase(pFFaceBC);
	spFStarFaces.erase(pFFaceCA);
	spFStarFaces.erase(pFFaceAB);

	// Now remove three StarEdges from the master list
#ifndef NDEBUG
	int iRes;
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVA, pVOtherPole));
	assert(iRes == 1);
#ifndef NDEBUG
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVB, pVOtherPole));
	assert(iRes == 1);
#ifndef NDEBUG
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVC, pVOtherPole));
	assert(iRes == 1);

	// Force the reconfiguration
#ifndef NDEBUG
	iRes =
#endif
	  iReconfTet32(apTCTets, pFMid, pVA, pVB, pVC, pVD, pVE);
	assert(iRes == 1);

	// Now find the new stuff.  pFMid should be the new face in the
	// interior of the shell.  Not a lot of finding to do here.
	assert(pFMid->qValid());
	assert(pFMid->qHasVert(pVA));
	assert(pFMid->qHasVert(pVB));
	assert(pFMid->qHasVert(pVC));

	// Add one new face to star face set
	pair<set<Face*>::iterator, bool> pRes =
	  spFStarFaces.insert(pFMid);
	assert(pRes.second);

	// Update three star edges, including new queue entries
	vUpdateDataStructs(pVA, pVB, pFFaceAB, pFMid,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVB, pVC, pFFaceBC, pFMid,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVC, pVA, pFFaceCA, pFMid,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
      }
      break;
    case eT22:
      // One fewer incident tet, one fewer star face, one fewer star edge.
      {
	// In this case, verts pVA, pVB, and pVC are the face on one
	// side of the coplanar quad.  Edge AB is a diagonal of the quad.
	iNumSwaps ++;
	// This assertion was removed when support for near-flat T22
	// swapping was added.  This assertion is uniformly satisfied by
	// the output from eCategorizeFace, which has tracked true for a
	// zillion calls, though, so there's no reason to be worried.
	// assert(iOrient3D(pVA, pVB, pVD, pVE) == 0);
	if (pVDelete == pVC) {
	  break;
	}

	bool qDelIsB = (pVDelete == pVB);
	Vert *pVOtherPole = qDelIsB ? pVA : pVB;

	// The old star faces that you'll be deleting are incident on
	// pVOtherPole.  More specifically, they're the faces opposite
	// pVDelete in the four tets in the list.  Must explicitly
	// identify these faces for later use.
	Face *pFFaceCD, *pFFaceEC;

	TetCell* pTCD = apTCTets[0];
	TetCell* pTCE = apTCTets[1];
	assert(pTCD->qHasVert(pVA));
	assert(pTCD->qHasVert(pVB));
	assert(pTCD->qHasVert(pVC));
	assert(pTCD->qHasVert(pVD));

	assert(pTCE->qHasVert(pVA));
	assert(pTCE->qHasVert(pVB));
	assert(pTCE->qHasVert(pVC));
	assert(pTCE->qHasVert(pVE));

	pFFaceCD = pTCD->pFFaceOpposite(pVDelete);
	pFFaceEC = pTCE->pFFaceOpposite(pVDelete);

	assert(pFFaceCD->qHasVert(pVOtherPole));
	assert(pFFaceCD->qHasVert(pVC));
	assert(pFFaceCD->qHasVert(pVD));
	assert(pFFaceEC->qHasVert(pVOtherPole));
	assert(pFFaceEC->qHasVert(pVE));
	assert(pFFaceEC->qHasVert(pVC));
	spFStarFaces.erase(pFFaceCD);
	spFStarFaces.erase(pFFaceEC);

	// Now remove one StarEdge from the master list
#ifndef NDEBUG
	int iRes;
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVC, pVOtherPole));
	assert(iRes == 1);

	// Force the reconfiguration
#ifndef NDEBUG
	iRes =
#endif
	  iReconfTet22(apTCTets, pFMid, pVA, pVB, pVC, pVD, pVE);
	assert(iRes == 1);

	// Now find the new stuff.  This code depends on the current
	// internal behavior of iReconfTet22!
	assert(pFMid->qHasVert(pVC));
	assert(pFMid->qHasVert(pVD));
	assert(pFMid->qHasVert(pVE));

	// Add the new star face
	pair<set<Face*>::iterator, bool> pRes =
	  spFStarFaces.insert(pFMid);
	assert(pRes.second);

	// No new star edge

	// Update two star edges, including new queue entries
	vUpdateDataStructs(pVC, pVD, pFFaceCD, pFMid,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVE, pVC, pFFaceEC, pFMid,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
      }
      break;
    case eT44:
      // Two fewer incident tets, two fewer star faces, three fewer star edges.
      {
	// It's possible to get a T44 case where the axis of the four
	// tets does NOT include the vertex to be deleted.  In that
	// case, do nothing.
	if (pVDelete == pVC) {
	  break;
	}

	// In this case, verts pVA, pVB, and pVC are the face on one
	// side of the coplanar quad.  Edge AB is a diagonal of the quad.
	iNumSwaps ++;
	// The following assertion was removed to allow near-flat bdry
	// edges to be swapped as flat.  Before that, it had always been
	// satisfied, so there's no real danger here.
	// assert(iOrient3D(pVA, pVB, pVD, pVE) == 0);

	bool qDelIsB = (pVDelete == pVB);
	Vert *pVOtherPole;
	if (qDelIsB) {
	  pVOtherPole = pVA;
	}
	else {
	  assert(pVDelete == pVA);
	  pVOtherPole = pVB;
	}


	// The old star faces that you'll be deleting are incident on
	// pVOtherPole.  More specifically, they're the faces opposite
	// pVDelete in the four tets in the list.  Must explicitly
	// identify these faces for later use.
	Face *pFFaceCD, *pFFaceDF, *pFFaceFE, *pFFaceEC;

	TetCell* pTCD = apTCTets[0];
	TetCell* pTCE = apTCTets[1];
	assert(pTCD->qHasVert(pVA));
	assert(pTCD->qHasVert(pVB));
	assert(pTCD->qHasVert(pVC));
	assert(pTCD->qHasVert(pVD));

	assert(pTCE->qHasVert(pVA));
	assert(pTCE->qHasVert(pVB));
	assert(pTCE->qHasVert(pVC));
	assert(pTCE->qHasVert(pVE));

	pFFaceCD = pTCD->pFFaceOpposite(pVDelete);
	pFFaceEC = pTCE->pFFaceOpposite(pVDelete);

	// Identify the vert on the back side of the coplanar quad.
	// It's on the common face of the two back tets, and it isn't
	// pVA or pVB.
	Face *pFFaceABD = pTCD->pFFaceOpposite(pVC);
	Face *pFFaceABE = pTCE->pFFaceOpposite(pVC);

	TetCell *pTCCellDF = dynamic_cast<TetCell*>
	  (pFFaceABD->pCCellOpposite(pTCD));
	TetCell *pTCCellEF = dynamic_cast<TetCell*>
	  (pFFaceABE->pCCellOpposite(pTCE));

	Face *pFCommon = pFCommonFace(pTCCellDF, pTCCellEF);
	Vert *pVF = pFCommon->pVVert(0);
	if (pVF == pVA || pVF == pVB) {
	  pVF = pFCommon->pVVert(1);
	  if (pVF == pVA || pVF == pVB) {
	    pVF = pFCommon->pVVert(2);
	  }
	}
	assert(pVF != pVA && pVF != pVB);

	assert(pTCCellDF->qHasVert(pVA));
	assert(pTCCellDF->qHasVert(pVB));
	assert(pTCCellDF->qHasVert(pVD));
	assert(pTCCellDF->qHasVert(pVF));
	assert(pTCCellEF->qHasVert(pVA));
	assert(pTCCellEF->qHasVert(pVB));
	assert(pTCCellEF->qHasVert(pVE));
	assert(pTCCellEF->qHasVert(pVF));

	pFFaceDF = pTCCellDF->pFFaceOpposite(pVDelete);
	pFFaceFE = pTCCellEF->pFFaceOpposite(pVDelete);

	assert(pFFaceCD->qHasVert(pVOtherPole));
	assert(pFFaceCD->qHasVert(pVC));
	assert(pFFaceCD->qHasVert(pVD));
	assert(pFFaceDF->qHasVert(pVOtherPole));
	assert(pFFaceDF->qHasVert(pVD));
	assert(pFFaceDF->qHasVert(pVF));
	assert(pFFaceFE->qHasVert(pVOtherPole));
	assert(pFFaceFE->qHasVert(pVF));
	assert(pFFaceFE->qHasVert(pVE));
	assert(pFFaceEC->qHasVert(pVOtherPole));
	assert(pFFaceEC->qHasVert(pVE));
	assert(pFFaceEC->qHasVert(pVC));
	spFStarFaces.erase(pFFaceCD);
	spFStarFaces.erase(pFFaceDF);
	spFStarFaces.erase(pFFaceFE);
	spFStarFaces.erase(pFFaceEC);

	// Now remove four StarEdges from the master list
#ifndef NDEBUG
	int iRes;
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVC, pVOtherPole));
	assert(iRes == 1);
#ifndef NDEBUG
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVD, pVOtherPole));
	assert(iRes == 1);
#ifndef NDEBUG
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVE, pVOtherPole));
	assert(iRes == 1);
#ifndef NDEBUG
	iRes =
#endif
	  sSEStarEdges.erase(StarEdge(pVF, pVOtherPole));
	assert(iRes == 1);

	// Force the reconfiguration
#ifndef NDEBUG
	iRes =
#endif
	  iReconfTet44(apTCTets, pFMid, pVA, pVB, pVC, pVD, pVE);
	assert(iRes == 1);

	// Now find the new stuff.  This code depends on the current
	// internal behavior of iReconfTet44!
	TetCell *pTCC = apTCTets[0];
	if (!pTCC->qHasVert(pVDelete)) {
	  pTCC = apTCTets[1];
	}
	Face *pFNewMid = pTCC->pFFaceOpposite(pVC);
	TetCell *pTCF = dynamic_cast<TetCell*>(pFNewMid->pCCellOpposite(pTCC));
	assert(pTCC->qHasVert(pVC));
	assert(pTCC->qHasVert(pVD));
	assert(pTCC->qHasVert(pVE));
	assert(pTCC->qHasVert(pVDelete));
	assert(pTCF->qHasVert(pVF));
	assert(pTCF->qHasVert(pVD));
	assert(pTCF->qHasVert(pVE));
	assert(pTCF->qHasVert(pVDelete));

	// Got the new tets incident on pVDelete; now find the new star
	// faces.
	Face *pFNewStarC = pTCC->pFFaceOpposite(pVDelete);
	Face *pFNewStarF = pTCF->pFFaceOpposite(pVDelete);
	assert(pFNewStarC->qHasVert(pVC));
	assert(pFNewStarC->qHasVert(pVD));
	assert(pFNewStarC->qHasVert(pVE));
	assert(pFNewStarF->qHasVert(pVF));
	assert(pFNewStarF->qHasVert(pVD));
	assert(pFNewStarF->qHasVert(pVE));

	// Add the new star faces
	pair<set<Face*>::iterator, bool> pRes =
	  spFStarFaces.insert(pFNewStarC);
	assert(pRes.second);
	pRes =
	  spFStarFaces.insert(pFNewStarF);
	assert(pRes.second);

	// Add new star edge
	StarEdge SENew(pVD, pVE, pFNewStarC, pFNewStarF, pFNewMid);
	sSEStarEdges.insert(SENew);

	// Update four star edges, including new queue entries
	vUpdateDataStructs(pVC, pVD, pFFaceCD, pFNewStarC,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVD, pVF, pFFaceDF, pFNewStarF,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVF, pVE, pFFaceFE, pFNewStarF,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
	vUpdateDataStructs(pVE, pVC, pFFaceEC, pFNewStarC,
			   pVDelete->adCoords(), sSEStarEdges, mRidgeMap);
      }
      break;
    case eN32:
      if (pVOther == pVDelete) {
	break;
      }
      // else fall through
    default:
      {
	// Others:  I think the others are all signs that it's impossible to
	//       remove the vertex.  At least, I hope so, because that'll
	//       make the code logic much, much simpler.
	// Output face classification and a warning message.
	vMessage(3, "Hit a face categorized as type %d; aborting removal\n", eFC);
	qSuccess = false;
      }
    }
  }

  if (qSuccess) {
    // Did all the required swapping, now it's just time to clean up the
    // small remaining part of the star (four tets for interior points,
    // three for bdry, two for bdry curve).  Credit one swap for this
    // final reconfiguration.
    iNumSwaps ++;
    switch (iFinalFaceCount) {
    case 4:
      // Done with an interior point removal
      {
	assert(spFStarFaces.size() == 4);
	// Identify all four remaining faces.
	Face *apF[4];
	std::copy(spFStarFaces.begin(), spFStarFaces.end(), apF);

	// Create new tet
	Cell *pCNew = pCNewCell();
	pCNew->vAssign(apF);

	// Find the cells and faces incident on pVDelete
	Cell *apC[4];
	set<Face*> spFIncidentFaces;
	for (int iF = 0; iF < 4; iF++) {
	  Face *pF = apF[iF];
	  Cell *pCCand = pF->pCCellLeft();
	  if (!pCCand->qHasVert(pVDelete)) {
	    pCCand = pF->pCCellRight();
	  }
	  assert(pCCand->qHasVert(pVDelete));
	  pF->vReplaceCell(pCCand, pCNew);
	  apC[iF] = pCCand;
	  for (int iF2 = 0; iF2 < 4; iF2++) {
	    Face *pFTmp = pCCand->pFFace(iF2);
	    if (pFTmp != pF) {
	      assert(pFTmp->qHasVert(pVDelete));
	      spFIncidentFaces.insert(pFTmp);
	    }
	  } // Loop adding incident faces
	  // Yes, every hint will get set three times, but at least
	  // they'll all be right.
	  pF->pVVert(0)->vSetHintFace(pF);
	  pF->pVVert(1)->vSetHintFace(pF);
	  pF->pVVert(2)->vSetHintFace(pF);
	} // Loop over remaining star faces

	// Mark all incident faces as deleted
	assert(spFIncidentFaces.size() == 6);
	set<Face*>::iterator iter;
	for (iter = spFIncidentFaces.begin(); iter != spFIncidentFaces.end();
	     iter++) {
	  (*iter)->vMarkDeleted();
	}

	apC[0]->vMarkDeleted();
	apC[1]->vMarkDeleted();
	apC[2]->vMarkDeleted();
	apC[3]->vMarkDeleted();
	assert(pCNew->iFullCheck() == 1);
	assert(apF[0]->iFullCheck() == 1);
	assert(apF[1]->iFullCheck() == 1);
	assert(apF[2]->iFullCheck() == 1);
	assert(apF[3]->iFullCheck() == 1);
      }
      break;
    case 3:
      // Done with a boundary point removal (on a facet, not a segment)
      {
	assert(spFStarFaces.size() == 3);
	// Identify all three remaining faces.
	Face *apF[3];
	std::copy(spFStarFaces.begin(), spFStarFaces.end(), apF);

	// Create new tet
	Cell *pCNew = pCNewCell();
	pCNew->vClearFaces();
	pCNew->vAddFace(apF[0]);
	pCNew->vAddFace(apF[1]);
	pCNew->vAddFace(apF[2]);

	// Find the cells and faces incident on pVDelete
	Cell *apC[3];
	set<Face*> spFIncidentFaces;
	for (int iF = 0; iF < 3; iF++) {
	  Face *pF = apF[iF];
	  Cell *pCCand = pF->pCCellLeft();
	  if (!pCCand->qHasVert(pVDelete)) {
	    pCCand = pF->pCCellRight();
	  }
	  assert(pCCand->qHasVert(pVDelete));
	  pF->vReplaceCell(pCCand, pCNew);
	  apC[iF] = pCCand;
	  for (int iF2 = 0; iF2 < 4; iF2++) {
	    Face *pFTmp = pCCand->pFFace(iF2);
	    if (pFTmp != pF) {
	      assert(pFTmp->qHasVert(pVDelete));
	      spFIncidentFaces.insert(pFTmp);
	    }
	  } // Loop adding incident faces
	} // Loop over remaining star faces

	// Mark all incident faces as deleted
	assert(spFIncidentFaces.size() == 6);
	set<Face*>::iterator iter;
	// The sizes of these arrays are known precisely in advance.
	BFace* apBFBFaces[3];
	Vert* apVAllVerts[4]; // Not including the one to be deleted.
	Vert* apVBdryVerts[3]; // Not including the one to be deleted.
	int iNumBF = 0, iNumVAll = 0, iNumVBdry = 0;
	for (iter = spFIncidentFaces.begin(); iter != spFIncidentFaces.end();
	     iter++) {
	  Face *pFTmp = *iter;
	  if (pFTmp->qIsBdryFace()) {
	    Cell *pCCand = pFTmp->pCCellLeft();
	    if (pCCand->eType() != Cell::eTriBFace) {
	      pCCand = pFTmp->pCCellRight();
	    }
	    assert(pCCand->eType() == Cell::eTriBFace);
	    if (pCCand->qHasVert(pVDelete)) {
	      // This way, I don't inadventantly get other nearby bdry
	      // faces mucking things up.
	      bool qHit;
	      for (int iNumV = 0; iNumV < 3; iNumV++) {
		qHit = false;
		Vert *pVTmp = pFTmp->pVVert(iNumV);
		for (int ii = 0; ii < iNumVBdry; ii++) {
		  if (apVBdryVerts[ii] == pVTmp) qHit = true;
		}
		if (!qHit) apVBdryVerts[iNumVBdry++] = pVTmp;
	      }
	      qHit = false;
	      BFace *pBF = dynamic_cast<BFace*>(pCCand);
	      for (int ii = 0; ii < iNumBF; ii++) {
		if (apBFBFaces[ii] == pBF) qHit = true;
	      }
	      if (!qHit) apBFBFaces[iNumBF++] = pBF;
	    }
	  }
	  for (int iV = 0; iV < 3; iV++) {
	    bool qHit = false;
	    Vert *pVTmp = pFTmp->pVVert(iV);
	    for (int ii = 0; ii < iNumVAll; ii++) {
	      if (apVAllVerts[ii] == pVTmp) qHit = true;
	    }
	    if (!qHit) apVAllVerts[iNumVAll++] = pVTmp;
	  }
	  (*iter)->vMarkDeleted();
	}
	assert(iNumVAll == 4);
	assert(iNumVBdry == 3);
	Vert *pVPole;
	for (int ii = 0; ii <= 3; ii++) { // We know how many there are...
	  Vert *pVCand = apVAllVerts[ii];
	  if (apVBdryVerts[0] != pVCand &&
	      apVBdryVerts[1] != pVCand &&
	      apVBdryVerts[2] != pVCand) {
	    pVPole = pVCand;
	    break;
	  }
	}

	assert(iNumBF == 3);
	// Can't check that they all lie on the same patch, because this
	// might be a mesh that was read in, with qStrictPatchChecking
	// not in effect.  In such a case, every surface triangle is its
	// own patch, even though there may be vast co-planar sections
	// of bdry with the same BC.
	// FIX ME: Once patch optimization is written, assert that the
	// bdry faces lie on the same patch.
	assert(apBFBFaces[0]->iBdryCond() ==
	       apBFBFaces[1]->iBdryCond());
	assert(apBFBFaces[0]->iBdryCond() ==
	       apBFBFaces[2]->iBdryCond());

	// Create a new bdry face; it will be attached to a face later.
	BFace *pBFNew = pBFNewBFace();
	pBFNew->vSetPatch(apBFBFaces[0]->pPatchPointer());
        apBFBFaces[0]->vMarkDeleted();
	apBFBFaces[1]->vMarkDeleted();
	apBFBFaces[2]->vMarkDeleted();

	apC[0]->vMarkDeleted();
	apC[1]->vMarkDeleted();
	apC[2]->vMarkDeleted();

	// Need to get the bdry verts (verts of the new face) into
	// in cyclic order.
	Vert *pVNew0 = apVBdryVerts[0];
	Vert *pVNew1 = apVBdryVerts[1];
	Vert *pVNew2 = apVBdryVerts[2];
	Face *pFNew = pFNewFace();
	if (iOrient3D(pVNew0, pVNew1, pVNew2, pVPole) == 1) {
	  pFNew->vAssign(pBFNew, pCNew, pVNew0, pVNew1, pVNew2);
	}
	else {
	  pFNew->vAssign(pBFNew, pCNew, pVNew0, pVNew2, pVNew1);
	}

	pVNew0->vSetHintFace(pFNew);
	pVNew1->vSetHintFace(pFNew);
	pVNew2->vSetHintFace(pFNew);
	assert(apF[0]->qHasVert(pVPole));
	pVPole->vSetHintFace(apF[0]);

	pBFNew->vAddFace(pFNew);
	pCNew->vAddFace(pFNew);

	assert(pCNew->iFullCheck() == 1);
	assert(pBFNew->iFullCheck() == 1);
	assert(pFNew->iFullCheck() == 1);
	assert(apF[0]->iFullCheck() == 1);
	assert(apF[1]->iFullCheck() == 1);
	assert(apF[2]->iFullCheck() == 1);
      }
      break;
    case 2:
      // Done with removal of a vert on a bdry curve
      {
	// Verts 1, d, and 2 lie along a common line.  Verts 3 and 4 are
	// also on the surface.  Tet 1 has verts 134d; tet 2 has 234d.
	assert(spFStarFaces.size() == 2);
	// Identify all remaining faces.
	Face *apF[2];
	std::copy(spFStarFaces.begin(), spFStarFaces.end(), apF);

	Cell *pCTmp = apF[0]->pCCellLeft();
	if (!pCTmp->qHasVert(pVDelete)) {
	  pCTmp = apF[0]->pCCellRight();
	}
	TetCell *pTC1 = dynamic_cast<TetCell*>(pCTmp);
	assert(pTC1->qValid());

	pCTmp = apF[1]->pCCellLeft();
	if (!pCTmp->qHasVert(pVDelete)) {
	  pCTmp = apF[1]->pCCellRight();
	}
	TetCell *pTC2 = dynamic_cast<TetCell*>(pCTmp);
	assert(pTC2->qValid());

	Face *pFMid = pFCommonFace(pTC1, pTC2);
	assert(pFMid->qValid());
	assert(pFMid->qHasVert(pVDelete));

	Vert *pV1 = pTC1->pVVertOpposite(pFMid);
	Vert *pV2 = pTC2->pVVertOpposite(pFMid);
	Vert *pV3, *pV4;
	if (pFMid->pVVert(0) == pVDelete) {
	  pV3 = pFMid->pVVert(1);
	  pV4 = pFMid->pVVert(2);
	}
	else if (pFMid->pVVert(1) == pVDelete) {
	  pV3 = pFMid->pVVert(2);
	  pV4 = pFMid->pVVert(0);
	}
	else {
	  assert(pFMid->pVVert(2) == pVDelete);
	  pV3 = pFMid->pVVert(0);
	  pV4 = pFMid->pVVert(1);
	}
	assert(pV1 != pV3);
	assert(pV1 != pV4);
	assert(pV2 != pV3);
	assert(pV2 != pV4);
	assert(pTC1->qHasVert(pV1));

	Face *pFFace13d = pTC1->pFFaceOpposite(pV4);
	Face *pFFace14d = pTC1->pFFaceOpposite(pV3);
	Face *pFFace234 = pTC2->pFFaceOpposite(pVDelete);
	Face *pFFace23d = pTC2->pFFaceOpposite(pV4);
	Face *pFFace24d = pTC2->pFFaceOpposite(pV3);

	// Remove BFaces
	Face *pF23d = pTC2->pFFaceOpposite(pV4);
	Cell *pCTemp1 = pF23d->pCCellOpposite(pTC2);
	assert(pCTemp1->eType() == Cell::eTriBFace);
	pCTemp1->vMarkDeleted();

	Face *pF24d = pTC2->pFFaceOpposite(pV3);
	Cell *pCTemp2 = pF24d->pCCellOpposite(pTC2);
	assert(pCTemp2->eType() == Cell::eTriBFace);
	pCTemp2->vMarkDeleted();

	// Update faces and the remaining cell.  Must do this -after-
	// all uses of pTC2 above, to avoid an assertion in
	// pFFaceOpposite.
	pFFace13d->vReplaceVert(pVDelete, pV2);
	pFFace14d->vReplaceVert(pVDelete, pV2);
	pFFace234->vReplaceCell(pTC2, pTC1);
	pTC1->vReplaceFace(pFMid, pFFace234);

	// Now delete extra faces and one cell
	pFMid->vMarkDeleted();
	pFFace23d->vMarkDeleted();
	pFFace24d->vMarkDeleted();
	pTC2->vMarkDeleted();

	// Reset hints that could have been broken.
	pV2->vSetHintFace(pFFace234);
	pV3->vSetHintFace(pFFace234);
	pV4->vSetHintFace(pFFace234);
	assert(pFFace13d->iFullCheck() == 1);
	assert(pFFace14d->iFullCheck() == 1);
	assert(pFFace234->iFullCheck() == 1);
	assert(pTC1->iFullCheck() == 1);
      }
      break;
    default:
      assert(0); // Should never get here.
      break;
    } // Done with end game.

    // Mark vertex as deleted
    pVDelete->vMarkDeleted();
    iStarRemeshSucceeded++;
  }
  else {
    iStarRemeshFailed++;
    // Back up to checkpointed mesh connectivity, and verify correctness.

    if (qBdryVert) {
      iStarRestoreFailed++;
    }
    else {
      // Identify and delete all cells and faces inside the original
      // star.  To do this, take advantage of knowing the original star
      // faces and the cells adjacent to them -outside- the star; this
      // info was saved up front.
      std::map<Face*, Cell*>::iterator mIter = mStarFaceRecovery.begin();
      std::set<Cell*> spCToDelete;
      std::set<Face*> spFToDelete;
      for ( ; mIter != mStarFaceRecovery.end(); mIter++) {
	Cell *pCBad = (mIter->first)->pCCellOpposite(mIter->second);
	pCBad->vMarkDeleted();
	spCToDelete.insert(pCBad);
      }
    
      for (iterC = spCToDelete.begin(); iterC != spCToDelete.end(); iterC++) {
	Cell *pC = *iterC;
	for (int iF = 0; iF < 4; iF++) {
	  Face *pF = pC->pFFace(iF);
	  mIter = mStarFaceRecovery.find(pF);
	  if (mIter == mStarFaceRecovery.end()) {
	    // The face isn't in the original star bdry, so it's correct
	    // to traverse across the face.
	    Cell *pCOpp = pF->pCCellOpposite(pC);
	    spCToDelete.insert(pCOpp);
	    spFToDelete.insert(pF);
	  }
	}
      }
    
      assert(spCToDelete.size()*4 ==
	     spFToDelete.size()*2 + mStarFaceRecovery.size());
      for (iterC = spCToDelete.begin(); iterC != spCToDelete.end(); iterC++) {
	(*iterC)->vMarkDeleted();
      }
      for (iterF = spFToDelete.begin(); iterF != spFToDelete.end(); iterF++) {
	(*iterF)->vMarkDeleted();
      }

      // Now build a new tet for each of the original star faces; other
      // faces will be assigned on the fly.
      for (mIter = mStarFaceRecovery.begin();
	   mIter != mStarFaceRecovery.end();
	   mIter++) {
	Face *pF = mIter->first;
	Cell *pCOut = mIter->second;
	Cell *pCIn = pF->pCCellOpposite(pCOut);
	Cell *pCNew = pCNewCell();
	pF->vReplaceCell(pCIn, pCNew);
	pCNew->vClearFaces();
	pCNew->vAddFace(pF);
      }

      // For each edge in the original star, create a new face and connect
      // it (bi-directionally) to the cells on either side of it.
      for (unsigned int iE = 0; iE < 3*uNNeighV - 6; iE++) {
	StarEdge SE = vSERecovery[iE];
	Vert *pV0 = SE.pVVert(0);
	Vert *pV1 = SE.pVVert(1);
	Face *pFL = SE.pFFaceLeft();
	Face *pFR = SE.pFFaceRight();

	mIter = mStarFaceRecovery.find(pFL);
	assert(mIter->first == pFL);
	Cell *pCLOut = mIter->second;

	mIter = mStarFaceRecovery.find(pFR);
	assert(mIter->first == pFR);
	Cell *pCROut = mIter->second;

	Cell *pCLIn = pFL->pCCellOpposite(pCLOut);
	Cell *pCRIn = pFR->pCCellOpposite(pCROut);
	Face *pFNew = pFNewFace();
	pFNew->vAssign(pCLIn, pCRIn, pVDelete, pV0, pV1);
	pCLIn->vAddFace(pFNew);
	pCRIn->vAddFace(pFNew);
	pVDelete->vSetHintFace(pFNew);
	pV0->vSetHintFace(pFNew);
	pV1->vSetHintFace(pFNew);
      }

      // Verify that all new faces and cells are properly connected, and
      // that cells are right-handed.
#ifndef NDEBUG
      {
	std::set<Cell*> spCInc;
	std::set<Face*> spFNearby;
	std::set<BFace*> spBFDummy;
	std::set<Vert*> spVDummy;
	bool qBVert = false;
	vNeighborhood(pVDelete, spCInc, spVDummy, &spBFDummy, &qBVert,
		      &spFNearby);
	assert(!qBVert && spBFDummy.size() == 0);
	for ( iterC = spCInc.begin();
	     iterC != spCInc.end(); iterC++) {
	  Cell *pC = *iterC;
	  assert(pC->iFullCheck());
	}
	for ( iterF = spFNearby.begin();
	     iterF != spFNearby.end(); iterF++) {
	  Face *pF = *iterF;
	  assert(pF->iFullCheck());
	  Vert *pV0 = pF->pVVert(0);
	  Vert *pV1 = pF->pVVert(1);
	  Vert *pV2 = pF->pVVert(2);
	  TetCell *pTCL = dynamic_cast<TetCell*>(pF->pCCellLeft());
	  if (pTCL) {
	    // The left side isn't a bdry.
	    Vert *pVL = pTCL->pVVertOpposite(pF);
	    assert(iOrient3D(pV0, pV1, pV2, pVL) == -1);
	  }
	  TetCell *pTCR = dynamic_cast<TetCell*>(pF->pCCellRight());
	  if (pTCR) {
	    // The right side isn't a bdry.
	    Vert *pVR = pTCR->pVVertOpposite(pF);
	    assert(iOrient3D(pV0, pV1, pV2, pVR) == 1);
	  }
	}
      }
#endif

    } // Done with restoration for interior verts.
  }

  // Restore the old swap recursion setting
  if (qOldRecursionSetting) {
    vAllowSwapRecursion();
  }
  // else do nothing, because recursion is already off.

  return qSuccess;
}
  
// TO DO:
//   o Fix for co-circular points (see hex test case).
//   o Bdry re-triangulation
//   o Checking for encroached bdry faces
//   o Backtracking for failures
