#include "GR_assert.h"
#include "GR_Geometry.h"
#include "GR_SurfMesh.h"
#include "GR_events.h"
#include "GR_Vec.h"
#include "GR_Bdry3D.h"

//@ Insert a point on an interior face
int SurfMesh::iInsertOnFace(Vert* const pVNew, Face* const pF,
			    Cell* const pC, const bool qSwap)
  // Note: Geometrically, pVNew may not fall exactly on pF because
  // of round-off.  Once things are done, though, the connectivity will
  // be the same as if it had.
  //
  // Before:                 After:
  //        B                      B
  //       /|\		          /|\       .
  //	  / | \		         / F \      .
  //  FA1/  |  \FA2	     FA1/  M  \FA2
  //    /   |	\	       /C3 B C4\    .
  //   /    |	 \	      /    |    \   .
  //  C	 C1 | C2  D	     C-FN1---FN2-D
  //   \    |	 /	      \    |    /
  //    \   |	/	       \C1 F C2/
  //  FB1\  |  /FB2	     FB1\  M  /FB2
  //	  \ | /		         \ A /
  //       \|/		          \|/
  //        A	           A
  //
{
  // Update for Periodic?
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(pC->iFullCheck());
  assert(pF->iFullCheck());

  if (pF->eType() == Face::eMultiEdge) {
    return iInsertOnMultiEdge(pVNew, dynamic_cast<MultiEdge*>(pF), qSwap);
  }

  Cell* pC1 = pF->pCCellLeft ();
  Cell* pC2 = pF->pCCellRight();
  assert(pC1->iFullCheck());
  assert(pC2->iFullCheck());
  assert(pC1 == pC || pC2 == pC);
  BdryPatch3D *pBP1 = pBdryPatch(iCellIndex(pC1));
  BdryPatch3D *pBP2 = pBdryPatch(iCellIndex(pC2));

  iOnEdge++;

  assert(pC1->eType() == Cell::eTriCell);
  assert(pC2->eType() == Cell::eTriCell);

  TriCell* const pTC1 = dynamic_cast<TriCell*>(pC1);
  TriCell* const pTC2 = dynamic_cast<TriCell*>(pC2);

  Vert *pVA = pF->pVVert(0);
  Vert *pVB = pF->pVVert(1);
  Vert *pVC = pTC1->pVVertOpposite(pF);
  Vert *pVD = pTC2->pVVertOpposite(pF);

  // Orientation for surface is a bit tricky
//   assert(iOrient2D(pVA, pVB, pVC) == 1);
//   assert(iOrient2D(pVB, pVA, pVD) == 1);

  Face* pFA1 = pTC1->pFFaceOpposite(pVA);
  Face* pFB1 = pTC1->pFFaceOpposite(pVB);
  Face* pFA2 = pTC2->pFFaceOpposite(pVA);
  Face* pFB2 = pTC2->pFFaceOpposite(pVB);
  Face* pFM  = pF;

  // Orientation for surface is a bit tricky
//   assert(iOrient2D(pVA, pVD, pVNew) == 1);
//   assert(iOrient2D(pVD, pVB, pVNew) == 1);
//   assert(iOrient2D(pVB, pVC, pVNew) == 1);
//   assert(iOrient2D(pVC, pVA, pVNew) == 1);

  //@@ Create new cells
  Cell* pC3 = pCNewCell();
  Cell* pC4 = pCNewCell();

  //@@ Create and set up new faces
  Face* pFMA = pFNewFace();
  Face* pFMB = pFM;
  Face* pFN1 = pFNewFace();
  Face* pFN2 = pFNewFace();

  pFMA->vAssign(pC3, pC4,   pVA, pVNew);
  pFMB->vAssign(pC1, pC2, pVNew,   pVB);
  pFN1->vAssign(pC1, pC3,   pVC, pVNew);
  pFN2->vAssign(pC4, pC2,   pVD, pVNew);

  //@@ Fix two outer faces which have new cells
  pFB1->vReplaceCell(pC1, pC3);
  pFB2->vReplaceCell(pC2, pC4);

  //@@ Set up the new cells
  pC1->vAssign(pFA1, pFN1, pFMB);
  pC2->vAssign(pFA2, pFN2, pFMB);
  pC3->vAssign(pFB1, pFN1, pFMA);
  pC4->vAssign(pFB2, pFN2, pFMA);
  vSetBdryPatch(pC1, pBP1);
  vSetBdryPatch(pC2, pBP2);
  vSetBdryPatch(pC3, pBP1);
  vSetBdryPatch(pC4, pBP2);

  //@@ Add a vertex hint for the new vert
  pVNew->vSetHintFace(pFN1);
  //@@ Update hints for the old verts that need it
  pVA->vSetHintFace(pFMA);
  pVB->vSetHintFace(pFMB);

  //@@ Swap locally
  // Swapping is not yet defined for surface meshes
  int iRetVal = 0;
  if (qSwap) {
    SUMAA_LOG_EVENT_END(INSERT_VERTEX);
    SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
    iRetVal += iFaceSwap(pFA1);
    iRetVal += iFaceSwap(pFA2);
    iRetVal += iFaceSwap(pFB1);
    iRetVal += iFaceSwap(pFB2);
    SUMAA_LOG_EVENT_END(NET_SWAPPING);
    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  }
  return(iRetVal);
}

//@ Insert a point on a MultiEdge
int SurfMesh::iInsertOnMultiEdge(Vert* const pVNew, MultiEdge* const pME,
				 const bool qSwap)
  // Note: Geometrically, pVNew may not fall exactly on pF because
  // of round-off.  Once things are done, though, the connectivity will
  // be the same as if it had. Each old cell incident on the MultiEdge
  // gets split in the same way.
  //
  // Before:                 After:
  //        B                      B
  //       /|		          /|
  //	  / |		         / M
  //  FA1/  |	             FA1/  E
  //    /   |		       /C2 2
  //   /    |		      /    |
  //  C	 C  M	             C-FN---
  //   \    E		      \    |
  //    \   |		       \C1 M
  //  FB1\  |	             FB1\  E
  //	  \ |		         \ 1
  //       \|		          \|
  //        A	           A
  //
{
  // Update for Periodic?
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(pME->iFullCheck());

  int iNCells = pME->iNumCells();
  // A place to store the new cells; the other set of new cells
  // overwrites the old ones.
  Cell **apC2 = new Cell*[iNCells];

  // All the faces that will need to be checked for swapping.
  Face **apFCheck = new Face*[iNCells * 2];

  MultiEdge *pME1 = pME;
  MultiEdge *pME2 = pMENewMultiEdge();

  Vert *pVA = pME->pVVert(0);
  Vert *pVB = pME->pVVert(1);

  iOnMultiEdge++;

  int iC;
  for (iC = 0; iC < iNCells; iC++) {
    Cell *pC1 = pME->pCCell(iC);
    assert(pC1->iFullCheck());
    BdryPatch3D *pBP = pBdryPatch(iCellIndex(pC1));

    assert(pC1->eType() == Cell::eTriCell);
    TriCell* const pTC = dynamic_cast<TriCell*>(pC1);

    Vert *pVC = pTC->pVVertOpposite(pME);

    Face* pFA1 = pTC->pFFaceOpposite(pVA);
    Face* pFB1 = pTC->pFFaceOpposite(pVB);

    apFCheck[iC*2    ] = pFA1;
    apFCheck[iC*2 + 1] = pFB1;

    //@@ Create new cell
    Cell * pC2 = apC2[iC] = pCNewCell();

    //@@ Create and set up new faces
    Face* pFN = pFNewFace();

    pFN->vAssign(pC1, pC2, pVC, pVNew);

    //@@ Fix outer face which has a new cell
    pFA1->vReplaceCell(pC1, pC2);

    //@@ Set up the new cells
    if (pFB1->eType() != Face::eMultiEdge) {
      pC1->vAssign(pFB1, pFN, pME1);
    }
    else {
      // Go ahead and assign anyway, and hope for the best...
      pC1->vAssign(pFN, pME1, pFB1);
    }

    if (pFA1->eType() != Face::eMultiEdge) {
      pC2->vAssign(pFA1, pFN, pME2);
    }
    else {
      // Go ahead and assign anyway, and hope for the best...
      pC2->vAssign(pFN, pME2, pFA1);
    }

    vSetBdryPatch(pC1, pBP);
    vSetBdryPatch(pC2, pBP);

    assert(pFN->iFullCheck() == 1);
    assert(pFA1->iFullCheck() == 1);
    assert(pFB1->iFullCheck() == 1);
    assert(pC1->iFullCheck() == 1);
  }
  // Update the two MultiEdge's
  pME1->vReplaceVert(pVB, pVNew);
  pME2->vAssign(pVNew, pVB, apC2, iNCells);
  assert(pME1->iNumCells() == pME2->iNumCells());

#ifndef NDEBUG
  for (iC = 0; iC < iNCells; iC++) {
    assert(apC2[iC]->iFullCheck() == 1);
  }
#endif

  //@@ Add a vertex hint for the new vert
  pVNew->vSetHintFace(pME1);

  //@@ Update hints for the old verts that need it
  pVA->vSetHintFace(pME1);
  pVB->vSetHintFace(pME2);

  //@@ Swap locally
  // Swapping is not yet defined for surface meshes
  int iRetVal = 0;
  if (qSwap) {
    SUMAA_LOG_EVENT_END(INSERT_VERTEX);
    SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
    for (int i = 0; i < iNCells*2; i++)
      iRetVal += iFaceSwap(apFCheck[i]);
    SUMAA_LOG_EVENT_END(NET_SWAPPING);
    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  }
  delete [] apC2;
  delete [] apFCheck;
  return(iRetVal);
}

//@ Insert a point inside a cell
int SurfMesh::iInsertInInterior(Vert* const pVNew, Cell* const pC,
				const bool qSwap)
{
  // Update for Periodic?
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(pC->iFullCheck());
  assert(pC->eType() == Cell::eTriCell);
  TriCell* pTC = dynamic_cast<TriCell*>(pC);

  // Verts ABCD are ordered in a right-handed way.  The faces are
  // defined relative to them so that there is no possible ambiguity.
  // Note that the face definitions here DO NOT correspond to
  // iFace[ABCD](iCell).
  Vert* pVA = pC->pVVert(0);
  Vert* pVB = pC->pVVert(1);
  Vert* pVC = pC->pVVert(2);

  Face* pFA = pTC->pFFaceOpposite(pVA);
  Face* pFB = pTC->pFFaceOpposite(pVB);
  Face* pFC = pTC->pFFaceOpposite(pVC);

  // Get new cells
  Cell* pCA = pC;
  Cell* pCB = pCNewCell();
  Cell* pCC = pCNewCell();

  // Identify the bdry patch
  BdryPatch3D *pBP = pBdryPatch(iCellIndex(pC));

  // Create new faces and set them up
  Face* pFA1 = pFNewFace();
  Face* pFB1 = pFNewFace();
  Face* pFC1 = pFNewFace();

  pFA1->vAssign(pCB, pCC, pVA, pVNew);
  pFB1->vAssign(pCC, pCA, pVB, pVNew);
  pFC1->vAssign(pCA, pCB, pVC, pVNew);

  // Fix old faces
  pFA->vReplaceCell(pC, pCA);
  pFB->vReplaceCell(pC, pCB);
  pFC->vReplaceCell(pC, pCC);

  // Set up new cells
  pCA->vAssign(pFA, pFB1, pFC1);
  pCB->vAssign(pFB, pFC1, pFA1);
  pCC->vAssign(pFC, pFA1, pFB1);

  assert(pCA ->iFullCheck());
  assert(pCB ->iFullCheck());
  assert(pFC ->iFullCheck());
  assert(pFA ->iFullCheck());
  assert(pFB ->iFullCheck());
  assert(pFC ->iFullCheck());
  assert(pFA1->iFullCheck());
  assert(pFB1->iFullCheck());
  assert(pFC1->iFullCheck());

  // Add a vertex hint for the new vert
  pVNew->vSetHintFace(pFC1);

  vSetBdryPatch(pCA, pBP);
  vSetBdryPatch(pCB, pBP);
  vSetBdryPatch(pCC, pBP);

  // Do some local swapping to clean up
  // Swapping not yet defined for surface meshes
  int iRetVal = 0;
  if (qSwap) {
    SUMAA_LOG_EVENT_END(INSERT_VERTEX);
    SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
    // Swapping for bdry faces not yet implemented
    iRetVal += iFaceSwap(pFA);
    iRetVal += iFaceSwap(pFB);
    iRetVal += iFaceSwap(pFC);
    SUMAA_LOG_EVENT_END(NET_SWAPPING);
    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  }
  return (iRetVal);
}

//@ Add a site (or some variant if it lies outside the mesh) to the mesh
// Return value indicates success
bool SurfMesh::qInsertPoint(const double adPoint[3], Cell* const pC,
			    int * const piSwaps, const bool qSwap,
			    const bool qForce, Vert* pVNew)
{
  *piSwaps = 0;
  //@@ The given cell had better contain the site
  // Project the cell and site onto a plane using the cell normal.
  double adNorm[3];
  pC->vVecSize(adNorm);
  vNORMALIZE3D(adNorm);
#ifndef NDEBUG
  {
    double dFaceLoc = dDOT3D(adNorm, pC->pVVert(0)->adCoords());
    double dSiteLoc = dDOT3D(adNorm, adPoint);
    assert(iFuzzyComp(dFaceLoc, dSiteLoc) == 0);
  }
#endif
  if (pVNew->qValid()) {
    assert(qForce); // Can't twiddle the point location if it's an
		    // existing vertex.
    assert(1.e-12 > dDIST3D(pVNew->adCoords(), adPoint));
  }

  // Coordinates in the plane of the triangular face
  Vert *pV0 = pC->pVVert(0);
  Vert *pV1 = pC->pVVert(1);
  Vert *pV2 = pC->pVVert(2);
  double adX[] = adDIFF3D(pV0->adCoords(), pV1->adCoords());
  vNORMALIZE3D(adX);
  double adY[3];
  vCROSS3D(adNorm, adX, adY);
  vNORMALIZE3D(adY);

  // Now project onto a plane
  double adLoc0[] = {dDOT3D(pV0->adCoords(), adX),
		     dDOT3D(pV0->adCoords(), adY)};
  double adLoc1[] = {dDOT3D(pV1->adCoords(), adX),
		     dDOT3D(pV1->adCoords(), adY)};
  double adLoc2[] = {dDOT3D(pV2->adCoords(), adX),
		     dDOT3D(pV2->adCoords(), adY)};
  double adLocSite[] = {dDOT3D(adPoint, adX), dDOT3D(adPoint, adY)};

  int iOrient01 = iOrient2D(adLoc0, adLoc1, adLocSite);
  int iOrient12 = iOrient2D(adLoc1, adLoc2, adLocSite);
  int iOrient20 = iOrient2D(adLoc2, adLoc0, adLocSite);
  // Make sure the inserted site lies inside or on the boundary of the
  // cell
  if (iOrient01*iOrient12 == -1 ||
      iOrient01*iOrient20 == -1 ||
      iOrient12*iOrient20 == -1) {
    vWarning("Point not inserted in surface mesh.  Lies outside cell.");
    return false;
  }

  // If the site lies strictly inside, insert it inside
  if (iOrient01*iOrient12*iOrient20 != 0) {
    iInterior++;
    if (!pVNew->qValid())
      pVNew = pVNewVert();
    pVNew->vSetCoords(3, adPoint);
    *piSwaps = iInsertInInterior(pVNew, pC, qSwap);
    return true;
  }
  else {
    // Find the face on which the site lies
    Face *pF = pFInvalidFace;
    for (int iF = 0; iF < 3; iF++) {
      pF = pC->pFFace(iF);
      int qHas0 = pF->qHasVert(pV0);
      int qHas1 = pF->qHasVert(pV1);
      if (iOrient01 == 0 && qHas0 && qHas1) break;
      int qHas2 = pF->qHasVert(pV2);
      if (iOrient12 == 0 && qHas1 && qHas2) break;
      if (iOrient20 == 0 && qHas2 && qHas0) break;
    }
    if (!pVNew->qValid())
      pVNew = pVNewVert();
    pVNew->vSetCoords(3, adPoint);
    *piSwaps = iInsertOnFace(pVNew, pF, pC, qSwap);
    return true;
  }
}

//@@ Determines whether swapping is needed for local Delaunay-ness.
static bool
qDoSwapDelaunay (const Vert * pVVertA, const Vert * pVVertB,
		 const Vert * pVVertC, const Vert * pVVertD)
  // Is pVVertD inside the diametral sphere of pVVertA, pVVertB, pVVertC?
{
  // Find the circumcenter (in-plane) of ABC by finding intersection of
  // the perpendicular bisectors of AB and BC with the plane ABC.
  double adVecAB[] = adDIFF3D(pVVertB->adCoords(), pVVertA->adCoords());
  double adVecBC[] = adDIFF3D(pVVertC->adCoords(), pVVertB->adCoords());

  vNORMALIZE3D(adVecAB);
  vNORMALIZE3D(adVecBC);

  double adNorm[3];
  vCROSS3D(adVecAB, adVecBC, adNorm);
  vNORMALIZE3D(adNorm);

  // The components of each of these vectors, dotted with coordinate
  // vectors from the origin, give level sets of distance parallel to
  // the vector.  Calibrate for the level set we need.
  double adRHS[3];
  {
    double adMidPtAB[] = {0.5 * (pVVertA->dX() + pVVertB->dX()),
			  0.5 * (pVVertA->dY() + pVVertB->dY()),
			  0.5 * (pVVertA->dZ() + pVVertB->dZ())};
    double adMidPtBC[] = {0.5 * (pVVertB->dX() + pVVertC->dX()),
			  0.5 * (pVVertB->dY() + pVVertC->dY()),
			  0.5 * (pVVertB->dZ() + pVVertC->dZ())};
    adRHS[0] = dDOT3D(adVecAB, adMidPtAB);
    adRHS[1] = dDOT3D(adVecBC, adMidPtBC);
    adRHS[2] = dDOT3D(adNorm,  pVVertA->adCoords());
  }
  double adCircCent[3];
  vSolve3By3(adVecAB, adVecBC, adNorm, adRHS[0], adRHS[1], adRHS[2],
	     adCircCent);
#ifndef NDEBUG
  {
    // Make sure we found the right point.
    double dDistA = dDIST3D(adCircCent, pVVertA->adCoords());
    double dDistB = dDIST3D(adCircCent, pVVertB->adCoords());
    double dDistC = dDIST3D(adCircCent, pVVertC->adCoords());
    assert(iFuzzyComp(dDistA, dDistB) == 0);
    assert(iFuzzyComp(dDistA, dDistC) == 0);
  }
#endif
  double dRadius = dDIST3D(adCircCent, pVVertA->adCoords());
  double dDistance = dDIST3D(adCircCent, pVVertD->adCoords());

  switch (iFuzzyComp(dRadius, dDistance)) {
  case -1:
    return false;		// pVVertD is outside
  case 1:
    return true;		// pVVertD is inside
  case 0:
    // Gives a consistent answer, regardless of which order the
    // configuration is presented in.  (ABCD, BADC, CDAB, DCBA all
    // result in the same triangulation.)
    return (min(pVVertA, pVVertB) > min(pVVertC, pVVertD));
  default:
    assert(0);
    return false;
  }
}

//@@ Determines whether swapping will improve the maximum face angle.
static bool
qDoSwapMaxDihed (const Vert * pVVertA, const Vert * pVVertB,
		 const Vert * pVVertC, const Vert * pVVertD)
{
  // This algorithm finds face angles between adjacent faces by dotting
  // unit vectors along the faces.  Small magnitude (close to -1)
  // indicates edges pointing in opposite directions, and this is bad,
  // so we choose the configuration that gives the largest min dot
  // product.
  //
  // To prevent pairs from flopping back and forth, the tie-breaker is
  // invoked if the inequality is small.

  // All of these normals are intended to point inwards.
  double adVecDA[] = adDIFF3D(pVVertD->adCoords(), pVVertA->adCoords());
  double adVecBD[] = adDIFF3D(pVVertB->adCoords(), pVVertD->adCoords());
  double adVecCB[] = adDIFF3D(pVVertC->adCoords(), pVVertB->adCoords());
  double adVecAC[] = adDIFF3D(pVVertA->adCoords(), pVVertC->adCoords());

  vNORMALIZE3D (adVecDA);
  vNORMALIZE3D (adVecBD);
  vNORMALIZE3D (adVecCB);
  vNORMALIZE3D (adVecAC);

  double dDotBDA = - dDOT3D (adVecBD, adVecDA);
  double dDotACB = - dDOT3D (adVecAC, adVecCB);
  double dDotDAC = - dDOT3D (adVecDA, adVecAC);
  double dDotCBD = - dDOT3D (adVecCB, adVecBD);

  double dMinThis  = min (dDotBDA, dDotACB);
  double dMinOther = min (dDotDAC, dDotCBD);

  switch (iFuzzyComp(dMinThis, dMinOther)) {
  case -1:
    return true;
  case 0:
    // Gives a consistent answer, regardless of which order the
    // configuration is presented in.
    return (min(pVVertA, pVVertB) > min(pVVertC, pVVertD));
  case 1:
    return false;
  default:
    assert(0);
    return false;
  }
}

bool
SurfMesh::qDoSwap (const Vert * const pVVertA, const Vert * const pVVertB,
		   const Vert * const pVVertC, const Vert * const pVVertD,
		   const Vert * const pVVertE) const
{
  assert (pVVertA->qValid ());
  assert (pVVertB->qValid ());
  assert (pVVertC->qValid ());
  assert (pVVertD->qValid ());
  assert(!pVVertE->qValid ());

  switch (eSwapType()) {
  case eDelaunay:
    return qDoSwapDelaunay (pVVertA, pVVertB, pVVertC, pVVertD);
  case eMaxDihed:
    return qDoSwapMaxDihed (pVVertA, pVVertB, pVVertC, pVVertD);
  default:
    return (false);
  }
}

int SurfMesh::iFaceSwap(Face* const pF)
{
  assert (eSwapType () != eNone);
  if (!pF->iFullCheck ())
    return (0);

  // Never swap a face that divides two regions in a multi-region mesh.
  if (pF->iFaceLoc () == Face::eBdryTwoSide)
    return (0);

  Cell *pCLeft = pF->pCCellLeft ();
  Cell *pCRight = pF->pCCellRight ();
  if (!pCLeft->qValid () || !pCRight->qValid () ||
      pCLeft->eType () != Cell::eTriCell ||
      pCRight->eType () != Cell::eTriCell ||
      (pBdryPatch(iCellIndex(pCLeft)) !=
       pBdryPatch(iCellIndex(pCRight)))
      )
    {
      return (0);
    }

  Vert *pVVertA = pF->pVVert (0);
  Vert *pVVertB = pF->pVVert (1);
  Vert *pVVertC = pCLeft->pVVertOpposite (pF);
  Vert *pVVertD = pCRight->pVVertOpposite (pF);

  // The orientation check that the 2D code does, can't be done here.

  if (qDoSwap (pVVertA, pVVertB, pVVertC, pVVertD)) {
    return iReconfigure(pF);
  }
  else
    return (0);
}

int SurfMesh::iReconfigure(Face * const pF)
{
  int iRetVal = Mesh2D::iReconfigure (pF);
  if (iRetVal) {
    // Check to make sure that neither cell incident on pF now has a
    // MultiEdge as its first face, as this can cause orientation
    // problems.
    Cell *pCL = pF->pCCellLeft();
    Face *pF0 = pCL->pFFace(0);
    if (pF0->eType() == Face::eMultiEdge) {
      vMessage(2, "Cleaning up after swap (verts %d %d)\n",
	       iVertIndex(pF0->pVVert(0)),
	       iVertIndex(pF0->pVVert(1)));
      Face *pF1 = pCL->pFFace(1);
      Face *pF2 = pCL->pFFace(2);
      if (pF1->eType() == Face::eEdgeFace)
	pCL->vAssign(pF1, pF2, pF0);
      else {
	assert(pF2->eType() == Face::eEdgeFace);
	pCL->vAssign(pF2, pF0, pF1);
      }
    }
    Cell *pCR = pF->pCCellRight();
    pF0 = pCR->pFFace(0);
    if (pF0->eType() == Face::eMultiEdge) {
      vMessage(2, "Cleaning up after swap (verts %d %d)\n",
	       iVertIndex(pF0->pVVert(0)),
	       iVertIndex(pF0->pVVert(1)));
      Face *pF1 = pCR->pFFace(1);
      Face *pF2 = pCR->pFFace(2);
      if (pF1->eType() == Face::eEdgeFace)
	pCR->vAssign(pF1, pF2, pF0);
      else {
	assert(pF2->eType() == Face::eEdgeFace);
	pCR->vAssign(pF2, pF0, pF1);
      }
    }
    return (iRetVal);
  }
  else {
    assert(iRetVal == 0);
    return 0;
  }
}
