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

static int iDelReq, iBdryDelReq, iBdryApexDelReq;
static int iDelDone, iBdryDelDone;

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) {}
};

// 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 Mesh2D::
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.
{
  iNewSwaps = 0;
  iDelReq++;

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

  if (pV->iVertType () == Vert::eBdryApex) {
    vMessage (2, "Deletion of boundary apexes is forbidden.\n");
    iBdryApexDelReq++;
    return false;
  }
  vMessage (4, "Deleting vertex at (%10f, %10f)\n", pV->dX (), pV->dY ());

  if (qRemoveVertCheaply(pV, iNewSwaps))
    return true;

  // 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<Cell*>::iterator iterC;
  std::set<Vert*>::iterator iterV;
  std::set<BFace*>::iterator iterBF;
  std::set<Face*>::iterator iterF;

  // Put the preferred verts at the head of the list.
  for (int ii = 0; ii < iNPref; spVCand.insert (apVPreferred[ii++])) {}
  bool qBdryVert = false;

  vNeighborhood (pV, spCIncident, spVCand, &spBFIncident, &qBdryVert,
		 &spFUpdate);

  if (qBdryVert) {
    // Abort if boundaries are precious
    if (!qBdryChangesAllowed ())
      return false;

    spVCand.clear ();
    for (iterBF = spBFIncident.begin(); iterBF != spBFIncident.end();
	 iterBF++) {
      BFace *pBF = *iterBF;
      spVCand.insert (pBF->pVVert (0));
      spVCand.insert (pBF->pVVert (1));
    }
    spVCand.erase (pV);
    assert (spVCand.size() == spBFIncident.size());
    iBdryDelReq++;
  }
  else {
    assert (spVCand.size() == spCIncident.size());
    assert (2 * spVCand.size() == spFUpdate.size());
  }

  // 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;
    if (pC->eType() == Cell::eQuadCell) return false;
    Face *pFOpp = dynamic_cast<TriCell *>(pC)->pFFaceOpposite (pV);
    vecpFOpposite.push_back(pFOpp);
    spFUpdate.erase(pFOpp);
  }
  // Used only for a pair of assertions just before 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;
  Vert *pV1 = pVInvalidVert, *pVBest = pVInvalidVert;
  for (iterV = spVCand.begin(); iterV != spVCand.end(); iterV++) {
    qVertOK = true;
    dMinDihed = 180;
    pV1 = *iterV;
    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->pCCellLeft ())->qHasVert (pV)) {
	  // For this case, the right-hand rule applies for orientation.
	  assert (iOrient2D (pF->pVVert (0), pF->pVVert (1), pV) != -1);
	  qResult = (iOrient2D (pF->pVVert (0), pF->pVVert (1),
				pV1) == 1);
	  if (qResult) {
	    Vert* apV[] = {pF->pVVert (0), pF->pVVert (1), pV1};
	    TriCellCV 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 (iOrient2D (pF->pVVert (1), pF->pVVert (0),
			     pV) != -1);
	  qResult = (iOrient2D (pF->pVVert (1), pF->pVVert (0),
				pV1) == 1);
	  if (qResult) {
	    Vert *apV[] = {pF->pVVert (1), pF->pVVert (0), pV1};
	    TriCellCV 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 (qVertOK) {
      qAnyVertOK = true;
      if (dMinDihed > dMaxMinDihed) {
	vMessage (4, "Best min dihedral yet.  Old: %10.6f  New: %10.6f\n",
		  dMaxMinDihed, dMinDihed);
	dMaxMinDihed = dMinDihed;
	pVBest = pV1;
      }
    }
  }

  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 will
  // need careful modification for structured meshes. (Version 0.1.8 or
  // later.) 
  // if (dMaxMinDihed < 3.) return false;
  vMessage (4, "%s %d, with min dihedral of %5.2f degrees after.\n",
	    "Contracting to remove vertex", iVertIndex (pV), dMaxMinDihed);

  // 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.
  // In 2D, there are precisely two such cells.
  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() == 1);	// In 2D, there is precisely
    // one of these. 

  }

  // 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.
  // In 2D, there is exactly one such face, which could be found by
  // determining the face common to the two cells to be deleted.
  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::eTriCell);
    TriCell *pTC = dynamic_cast<TriCell *>(*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.
  // These assertions are very different in 2D as compared to 3D.
  if (qBdryVert) {
    assert (iNInc == 1);	// Different in 2D and 3D.
    assert (int(spFUpdate.size() + spFLateral.size() + iNInc)
	    == iTotalIncFaces);
  }
  else {
    assert (iNInc == 2);	// Different in 2D and 3D.
    assert (int(spFUpdate.size() + spFLateral.size() + iNInc)
	    == 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.
  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);
    }
    pC->vReplaceFace (pF, pF1);
    pF1->vReplaceCell (pCInc, pC);
#ifndef NDEBUG
    assert (pC->iFullCheck () == 1);
    assert (pF1->iFullCheck () == 1);
    {
      Cell *pCTest = pF1->pCCellLeft ();
      assert (pCTest->iFullCheck () == 1);
      pCTest = pF1->pCCellRight ();
      assert (pCTest->iFullCheck () == 1);
    }
#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));
  }

#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::eBdryEdge ||
	    iOrient2D (pC->pVVert (0), pC->pVVert (1), pC->pVVert (2)) == 1);
    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 ()));
  }
  // Mark the cells that were changed as bad cells
  for (iterC = spCCheck.begin(); iterC != spCCheck.end(); iterC++) {
    Cell *pC = *iterC;
    pC->vSetShapeClass(Cell::eUnknown);
    pC->vMarkIllSized();		
  }
#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;
    for (int iV = 0; iV < 2; iV++) {
      Vert *pVTmp = pF1->pVVert (iV);
      pVTmp->vSetHintFace (pF1);
      pVTmp->vMarkIllShaped ();
    }
  }

  int iNumSwaps = 0;
  // Swap every face that was incident to the deleted vertex that survived
  for (iterF = spFUpdate.begin(); iterF != spFUpdate.end(); iterF++) {
    Face *pF = *iterF;
    if (!pF->qDeleted()) {
      iNumSwaps += iFaceSwap(pF);
    }
  }
  
  // Success!
  iDelDone++;
  if (qBdryVert)
    iBdryDelDone++;

  vUpdateLengthScale();

  return true;
}

// Yet another removal method: this one deleted everything in the
// neighborhood of the vertex being deleted, then remeshes the resulting
// star-shaped region.  (For 2D, this is guaranteed to work, so there's
// no problem with this.)  The remeshing uses a front (implemented as a
// set of faces), with each face in the front used to create a Delaunay
// triangle.  When creating triangles, new faces are created iff they
// don't already exist in the front.

// This code fails for boundary verts.

class FrontFace2D {
  Vert *pV0, *pV1;
  Face *pF;
public:
  FrontFace2D(Vert * const pVA, Vert * const pVB,
	    Face * const pFIn = pFInvalidFace) 
    : pV0(pVA), pV1(pVB), pF(pFIn) {
    assert(pV0 != pV1);
  }
  Vert *pVVert(const int i) const {
    assert(i == 0 || i == 1);
    return (i == 0 ? pV0 : pV1);
  }
  Face *pFFace() const {return pF;}
  bool operator==(const FrontFace2D& FF) const {
    return ((pV0 == FF.pV0 && pV1 == FF.pV1) ||
	    (pV0 == FF.pV1 && pV1 == FF.pV0));
  }
  bool operator<(const FrontFace2D& FF) const {
    Vert *pVMinThis, *pVMaxThis, *pVMinOther, *pVMaxOther;
    if (pV0 < pV1) {
      pVMinThis = pV0; pVMaxThis = pV1;
    }
    else {
      pVMinThis = pV1; pVMaxThis = pV0;
    }
    if (FF.pV0 < FF.pV1) {
      pVMinOther = FF.pV0; pVMaxOther = FF.pV1;
    }
    else {
      pVMinOther = FF.pV1; pVMaxOther = FF.pV0;
    }
    return ((pVMinThis < pVMinOther) ||
	    ((pVMinThis == pVMinOther) && (pVMaxThis < pVMaxOther)));
  }
};

bool Mesh2D::qRemoveVert(Vert * const pV, int& iNSwaps)
{
  iNSwaps = 0;
  // Find the neighborhood of pV.

  std::set<Cell*> spCInc;
  std::set<Vert*> spVNeigh;
  std::set<Face*> spFNearby;

  {
    std::set<BFace*> spBFDummy;
    bool qBdryVert;
    vNeighborhood (pV, spCInc, spVNeigh, &spBFDummy, &qBdryVert,
		   &spFNearby);
    
    
    if (qBdryVert) {
      assert(pV->qIsBdryVert());
      assert(spBFDummy.size() >= 2); // Could be greater for a corner
      // where regions meet.
      return false;
    }
    else {
      assert(!pV->qIsBdryVert());
      assert(spBFDummy.empty());
    }
  }
  
  // Rip out all the cells in the neighborhood...
  for (std::set<Cell*>::iterator iterC = spCInc.begin();
       iterC != spCInc.end(); iterC++) {
    Cell *pC = *iterC;
    pC->vMarkDeleted();
  }
  
  // For faces, rip out those incident, and add the other to the set of
  // front faces.
  std::set<FrontFace2D> sFF;
  for (std::set<Face*>::iterator iterF = spFNearby.begin();
       iterF != spFNearby.end(); ) {
    Face *pF = *iterF;
    std::set<Face*>::iterator iterTmp = iterF++;
    if (pF->qHasVert(pV)) {
      pF->vMarkDeleted();
      spFNearby.erase(iterTmp);
    }
    else {
      Vert *pV0 = pF->pVVert(0);
      Vert *pV1 = pF->pVVert(1);
      // Update the vertex hints for all vertices in the hull, since
      // their hints may currently be faces that are being deleted.
      pV0->vSetHintFace(pF);
      pV1->vSetHintFace(pF);

      if (pF->pCCellLeft()->qDeleted()) {
	pF->vRemoveCell(pF->pCCellLeft());
	sFF.insert(FrontFace2D(pV0, pV1, pF));
      }
      else {
	assert(pF->pCCellRight()->qDeleted());
	pF->vRemoveCell(pF->pCCellRight());
	sFF.insert(FrontFace2D(pV1, pV0, pF));
      }	
    }
  } // Done with loop over faces.

  // Mark the vertex as deleted.
  pV->vMarkDeleted();
  int iNumNewCells = 0, iNumNewFaces = 0;

  // For each face in the hull, build the Delaunay triangle that belongs
  // on it.
  do {
    assert(sFF.size() >= 3);
    // Maintaining this list correctly is easier if you just rebuild it
    // every time.
    std::set<Vert*> spVCand;
    {
      std::set<FrontFace2D>::iterator iter;
      for (iter = sFF.begin(); iter != sFF.end(); iter++) {
	spVCand.insert(iter->pVVert(0));
	spVCand.insert(iter->pVVert(1));
      }
      assert(spVCand.size() <= sFF.size());
    }

    FrontFace2D FF = *sFF.begin();
    sFF.erase(sFF.begin());
    Vert *pV0 = FF.pVVert(0);
    Vert *pV1 = FF.pVVert(1);

    // Find the right vertex: the one on the correct side of the face
    // with no point in the resulting circumcircle.  In case of ties
    // between candidate verts, pick the one nearest the face midpoint.
    std::set<Vert*>::iterator iterV;
    for (iterV = spVCand.begin();
	 iterV != spVCand.end() && iOrient2D(pV0, pV1, *iterV) != 1;
	 iterV++) {}
    assert(iterV != spVCand.end());
    Vert *pVBest = *iterV++;
    for ( ; iterV != spVCand.end(); iterV++) {
      Vert *pVCand = *iterV;
      if (pVCand == pV0 || pVCand == pV1) continue;
      if (iOrient2D(pV0, pV1, pVCand) == 1) {
	int iIncircleResult = iIncircle(pV0, pV1, pVBest, pVCand);
	if (iIncircleResult == 1) {
	  pVBest = pVCand;
	}
	else if (iIncircleResult == 0) {
	  double adMid[] = {0.5 * (pV0->dX() + pV1->dX()),
			    0.5 * (pV0->dY() + pV1->dY())};
	  if (dDIST2D(adMid, pVCand->adCoords()) <
	      dDIST2D(adMid, pVBest->adCoords()))
	    pVBest = pVCand;
	}
      }
    }
    // At this point, pVBest is definitely the one.
    
    // Create the new cell and connect it to the known face.
    iNumNewCells++;
    Cell *pCNew = pCNewCell();
    pCNew->vClearFaces();
    pCNew->vAddFace(FF.pFFace());
    FF.pFFace()->vAddCell(pCNew);

    // Does the first "new" face already exist?
    Face *pF0 = pFInvalidFace, *pF1 = pFInvalidFace;
    std::set<FrontFace2D>::iterator iter = sFF.find(FrontFace2D(pV0, pVBest));
    if (iter != sFF.end()) {
      // Got a hit on this face.  Identify it, and remove its entry from
      // the front.
      pF0 = iter->pFFace();
      assert(pF0->qHasVert(pV0));
      assert(pF0->qHasVert(pVBest));
      pF0->vAddCell(pCNew);
      // Check the orientation
      assert((pF0->pVVert(0) == pV0 && pF0->pCCellRight() == pCNew) ||
	     (pF0->pVVert(1) == pV0 && pF0->pCCellLeft () == pCNew));
      sFF.erase(iter);
    }
    else {
      // This face is new.  Create a face, and add an entry to the front.
      pF0 = pFNewFace();
      iNumNewFaces++;
      pF0->vAssign(pCInvalidCell, pCNew, pV0, pVBest);
      sFF.insert(FrontFace2D(pV0, pVBest, pF0));
    }
    pCNew->vAddFace(pF0);

    // Does the second "new" face already exist?
    iter = sFF.find(FrontFace2D(pV1, pVBest));
    if (iter != sFF.end()) {
      // Got a hit on this face.  Identify it, and remove its entry from
      // the front.
      pF1 = iter->pFFace();
      assert(pF1->qHasVert(pV1));
      assert(pF1->qHasVert(pVBest));
      pF1->vAddCell(pCNew);
      // Check the orientation
      assert((pF1->pVVert(1) == pV1 && pF1->pCCellRight() == pCNew) ||
	     (pF1->pVVert(0) == pV1 && pF1->pCCellLeft () == pCNew));
      sFF.erase(iter);
    }
    else {
      // This face is new.  Create a face, and add an entry to the front.
      pF1 = pFNewFace();
      iNumNewFaces++;
      pF1->vAssign(pCNew, pCInvalidCell, pV1, pVBest);
      sFF.insert(FrontFace2D(pVBest, pV1, pF1));
    }
    pCNew->vAddFace(pF1);
  } while (!sFF.empty());

  assert(iNumNewCells == int(spFNearby.size()) - 2);
  assert(iNumNewFaces == int(spFNearby.size()) - 3);
  
  // Always succeeds.
  return true;
}

