//@ Include files, #define's, etc.
#include <math.h>
#include "GR_assert.h"
#include "GR_misc.h"
#include "GR_Geometry.h"
#include "GR_Face.h"
#include "GR_Cell.h"
#include "GR_BFace.h"
#include "GR_Quality.h"
#include "GR_Vertex.h"
#include "GR_VolMesh.h"
#include "GR_Vec.h"

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

static eFaceCat eFC;
static int iSwapDepth = 0;
static const int iMaxSwapDepth = 75;

int aiEdgeReq[10], aiEdgeDone[10];

//@ Decide which type of triangular face this is, anyway.
eFaceCat TriFace::eCategorizeFace(Vert*& pVVertA, Vert*& pVVertB,
				  Vert*& pVVertC, Vert*& pVVertD,
				  Vert*& pVVertE, TetCell* apTCTets[],
				  int& iNTets,
                                  Vert*& pVPivot0, Vert*& pVPivot1,
                                  Vert*& pVOther, const double dMaxAngle) const
{
  //@@ Case: one or both cells is not a tet, including boundaries
  Cell* pCLeft = pCCellLeft();
  Cell* pCRight = pCCellRight();
  if (!pCLeft->qValid() || !pCRight->qValid()) return(eOther);
  if (qIsBdryFace()) return (eBdry);
  if (pCLeft->eType() != Cell::eTet ||
      pCRight->eType() != Cell::eTet) return(eOther);

  //@@ List all tets formed from five verts defining two cells sharing iFace
  iNTets = 0;
  apTCTets[0] = dynamic_cast<TetCell*>(pCLeft);
  apTCTets[1] = dynamic_cast<TetCell*>(pCRight);

  // The (Face*)this is really const, because nothing gets done here.
  // But the Vert*'s really should be non-const, because the data is
  // being returned to a caller who might, potentially, decide to swap
  // the mesh.
  pVVertA = const_cast<Vert*>(pVVert(0));
  pVVertB = const_cast<Vert*>(pVVert(1));
  pVVertC = const_cast<Vert*>(pVVert(2));
  pVVertD = apTCTets[0]->pVVertOpposite(this);
  pVVertE = apTCTets[1]->pVVertOpposite(this);
  iNTets = 2;

  // Any add'l cell must be a neighbor of both of the first two.
  int i;
  for (i = 0; i < 4 ; ++ i) {
    Face* pFCandFace = apTCTets[0]->pFFace(i);
    Cell* pCCandTet = pFCandFace->pCCellOpposite(apTCTets[0]);
    // Don't bother with cells that aren't tets; these are unswappable
    // anyway
    if (pCCandTet->eType() != Cell::eTet) continue;
    if (pCCandTet != apTCTets[1] &&
	pCCandTet->qHasVert(pVVertE))
      apTCTets[iNTets++] = dynamic_cast<TetCell*>(pCCandTet);
  }

  assert(iNTets <= 4);

  switch (iNTets) {
  case 2:
    {
      //@@ Two tets
      assert(iOrient3D(pVVertA, pVVertC, pVVertB, pVVertD) == 1);
      assert(iOrient3D(pVVertA, pVVertB, pVVertC, pVVertE) == 1);
//       if (iOrient3D(pVVertA, pVVertC, pVVertB, pVVertD) != 1 ||
// 	  iOrient3D(pVVertA, pVVertB, pVVertC, pVVertE) != 1)
// 	vFatalError("Orientation error", "FaceSwap");


      // The following code for eN32 face identification uses only as
      // many orientation primitives as required to determine the actual
      // face type.  This is very useful, as this is a critical code
      // path. 

      // Initialize these to ludicrous values.
      int iOrientA = 100, iOrientB = 100, iOrientC = 100;
      bool qApproxFlatA = false, qApproxFlatB = false, qApproxFlatC = false;
      bool qBadOrientation = false;

      // There is always the possibility that the configuration has one
      // or two correctly-oriented but nearly flat (within the tolerance
      // set by VolMesh::dMaxAngleForSurfSwap) faces.  So if the
      // orientation is positive, check whether it should perhaps be
      // tagged as zero instead.

      iOrientA = iOrient3D(pVVertB, pVVertC, pVVertD, pVVertE);
      if (iOrientA == -1) {
        pVPivot0 = pVVertB;
        pVPivot1 = pVVertC;
        pVOther  = pVVertA;
	qBadOrientation = true;
      }
      else {
	if (iOrientA == 1) {
	  // Check whether it's nearly flat
	  double adNorm0[3], adNorm1[3];
	  vNormal3D(pVVertB->adCoords(), pVVertC->adCoords(),
		    pVVertD->adCoords(), adNorm0);
	  vNormal3D(pVVertC->adCoords(), pVVertB->adCoords(),
		    pVVertE->adCoords(), adNorm1);
	  vNORMALIZE3D(adNorm0);
	  vNORMALIZE3D(adNorm1);
	  double dDot = dDOT3D(adNorm0, adNorm1);
	  if (dDot > cos(dMaxAngle * M_PI / 180)) {
	    qApproxFlatA = true;
	  }
	}
	iOrientB = iOrient3D(pVVertC, pVVertA, pVVertD, pVVertE);
	if (iOrientB == -1) {
	  pVPivot0 = pVVertC;
	  pVPivot1 = pVVertA;
	  pVOther  = pVVertB;
	  qBadOrientation = true;
	}
	else if (iOrientB == 1) {
	  // Check whether it's nearly flat
	  double adNorm0[3], adNorm1[3];
	  vNormal3D(pVVertC->adCoords(), pVVertA->adCoords(),
		    pVVertD->adCoords(), adNorm0);
	  vNormal3D(pVVertA->adCoords(), pVVertC->adCoords(),
		    pVVertE->adCoords(), adNorm1);
	  vNORMALIZE3D(adNorm0);
	  vNORMALIZE3D(adNorm1);
	  double dDot = dDOT3D(adNorm0, adNorm1);
	  if (dDot > cos(dMaxAngle * M_PI / 180)) {
	    qApproxFlatB = true;
	  }
	}
	if (iOrientA + iOrientB == 0) return eN20;
	if (!qBadOrientation) {
	  iOrientC = iOrient3D(pVVertA, pVVertB, pVVertD, pVVertE);
	  if (iOrientC == -1) {
	    pVPivot0 = pVVertA;
	    pVPivot1 = pVVertB;
	    pVOther  = pVVertC;
	    qBadOrientation = true;
	  }
	  else if (iOrientC == 1) {
	    // Check whether it's nearly flat
	    double adNorm0[3], adNorm1[3];
	    vNormal3D(pVVertA->adCoords(), pVVertB->adCoords(),
		      pVVertD->adCoords(), adNorm0);
	    vNormal3D(pVVertB->adCoords(), pVVertA->adCoords(),
		      pVVertE->adCoords(), adNorm1);
	    vNORMALIZE3D(adNorm0);
	    vNORMALIZE3D(adNorm1);
	    double dDot = dDOT3D(adNorm0, adNorm1);
	    if (dDot > cos(dMaxAngle * M_PI / 180)) {
	      qApproxFlatC = true;
	    }
	  }
	}
      }
      if (qBadOrientation) {
	// Need to distinguish between eN32 cases and boundary cases.
	// The current code does this incompletely, in that one of the
	// two cells must have a bdry face for the boundary case to be
	// recognized.  Nevertheless, the main point here is to weed out
	// the cases where two bdry tets have a slightly reflex edge
	// between them, and correctly -not- identify these as eN32.
	Face *pF0 = apTCTets[0]->pFFaceOpposite(pVOther);
	Face *pF1 = apTCTets[1]->pFFaceOpposite(pVOther);

	if (pF0->qIsBdryFace() || pF1->qIsBdryFace())
	  return eBdryReflex;
	else
	  return eN32;
      }	    

      switch (iOrientA + iOrientB + iOrientC) {
      case 0: // Impossible, but included for completeness.
	return eOther;
      case 1: // Two orientations are zero (three points co-linear).
	      // Hopelessly unswappable.  Bail out.
	return eN20;
      case 2: {
	//@@@ Four points are coplanar; T22, T44, N44
	// One orientation must be 0; verts are re-labeled to make it
	// iOrientC.  This implies that AB is the coplanar edge.	  
	assert(! (iOrientA && iOrientB && iOrientC));
	if (iOrientA == 0) {
	  Vert* pVTemp;
	  pVTemp = pVVertA;
	  pVVertA = pVVertB;
	  pVVertB = pVVertC;
	  pVVertC = pVTemp;
	}
	else if (iOrientB == 0) {
	  Vert* pVTemp;
	  pVTemp = pVVertA;
	  pVVertA = pVVertC;
	  pVVertC = pVVertB;
	  pVVertB = pVTemp;
	}
	else if (iOrientC == 0) {}
	// The following assertion has been removed deliberately, to
	// accomodate the near-flat cases tested for with dot products
	// above.  Prior to adding support for near-flat swapping, the
	// orientation was in fact re-aligned correctly, so there's no
	// cause for concern.
	// Verify that the relabeling was done correctly
	//	assert( iOrient3D(pVVertA, pVVertB, pVVertD, pVVertE) ==
	//	0 );
	

	// The configuration of these two tets is classified based on the
	// properties of the coplanar faces.
	// 1 If both are BFaces with the same boundary condition,
	//   swappable two tets to two.
	// 2 If both are BFaces with different bdry cond, not swappable.
	// 3 If one is a BFace and the other not, not swappable.
	// 4 If neither is a BFace, both opposite cells are tets, and the
	//   tets have the same fourth vert, swappable four to four.
	// 5 If neither is a BFace, both opposite cells are tets, and the
	//   tets do not have the same fourth vert, not swappable,
	//   although some non-local transformations might make it so.
	// 6 If neither is a BFace and one or both opposite cells is not a
	//   tet, not swappable.
	TetCell *pTCTetA = apTCTets[0], *pTCTetB = apTCTets[1];
	// These are the faces that are coplanar
	Face* pFFaceA = pTCTetA->pFFaceOpposite(pVVertC);
	Face* pFFaceB = pTCTetB->pFFaceOpposite(pVVertC);
	Cell* pCOppCellA = pFFaceA->pCCellOpposite(pTCTetA);
	Cell* pCOppCellB = pFFaceB->pCCellOpposite(pTCTetB);

  	if (pCOppCellA->eType() == Cell::eTriBFace &&
	    pCOppCellB->eType() == Cell::eTriBFace) {
	  // Both faces are on the boundary
	  TriBFace* pTBFA = dynamic_cast<TriBFace*>(pCOppCellA);
	  TriBFace* pTBFB = dynamic_cast<TriBFace*>(pCOppCellB);
	  // These are already -known- to be co-planar, so checking bdry
	  // patches is sufficient.
	  // FIX ME: 0.3  This is a spot where checking patches, not
	  // BC's is the right thing to do.  But until patch
	  // optimization gets written, the BC check is the better
	  // choice, because then you can swap coplanar faces in meshes
	  // that you read in.
	  // 	  if (pTBFA->pPatchPointer() == pTBFB->pPatchPointer())
	  if (pTBFA->iBdryCond() == pTBFB->iBdryCond())
	    return eT22; // case 1
	  else
	    return eOther;  // case 2
	}
	else if (pCOppCellA->eType() == Cell::eTet &&
		 pCOppCellB->eType() == Cell::eTet) {
	  // Both faces are internal to the mesh and separate pairs of tets
	  TetCell* pTCA = dynamic_cast<TetCell*>(pCOppCellA);
	  TetCell* pTCB = dynamic_cast<TetCell*>(pCOppCellB);
	  if (pTCA->pVVertOpposite(pFFaceA) ==
	      pTCB->pVVertOpposite(pFFaceB))
	    return eT44; // case 4
	  else {
            pVPivot0 = pVVertA;
            pVPivot1 = pVVertB;
            pVOther  = pVVertC;
	    return eN44; // case 5
          }
	}
	else
	  // Exactly one face on the boundary or internal faces with
	  // cells other than tets
	  return eOther; // cases 3,6
      } // End two-tet cases with four verts coplanar
      case 3:
	{
	  // This is probably a T23 case, but it might be a perturbed T22
	  // case, which should be tagged so that removal of surface verts
	  // will work better.  (Actually, it could even be a perturbed
	  // N20, which is a disaster during removal if you try to run
	  // with it as a T23! so nuke this case.)
	  if ((qApproxFlatA && qApproxFlatB) ||
	      (qApproxFlatB && qApproxFlatC) ||
	      (qApproxFlatC && qApproxFlatA))
	    return eN20;
	  if (qApproxFlatA || qApproxFlatB || qApproxFlatC) {
	    // It is a perturbed T22
	    if (qApproxFlatA) {
	      Vert* pVTemp;
	      pVTemp = pVVertA;
	      pVVertA = pVVertB;
	      pVVertB = pVVertC;
	      pVVertC = pVTemp;
	    }
	    else if (qApproxFlatB) {
	      Vert* pVTemp;
	      pVTemp = pVVertA;
	      pVVertA = pVVertC;
	      pVVertC = pVVertB;
	      pVVertB = pVTemp;
	    }
	    else if (qApproxFlatC) {}
	    
	    TetCell *pTCTetA = apTCTets[0], *pTCTetB = apTCTets[1];
	    // These are the faces that are coplanar
	    Face* pFFaceA = pTCTetA->pFFaceOpposite(pVVertC);
	    Face* pFFaceB = pTCTetB->pFFaceOpposite(pVVertC);
	    Cell* pCOppCellA = pFFaceA->pCCellOpposite(pTCTetA);
	    Cell* pCOppCellB = pFFaceB->pCCellOpposite(pTCTetB);
	    
	    if (pCOppCellA->eType() == Cell::eTriBFace &&
		pCOppCellB->eType() == Cell::eTriBFace) {
	      // Both faces are on the boundary
	      TriBFace* pTBFA = dynamic_cast<TriBFace*>(pCOppCellA);
	      TriBFace* pTBFB = dynamic_cast<TriBFace*>(pCOppCellB);
	      // These are already -known- to be co-planar, so checking bdry
	      // patches is sufficient.
	      // FIX ME: 0.3  This is a spot where checking patches, not
	      // BC's is the right thing to do.  But until patch
	      // optimization gets written, the BC check is the better
	      // choice, because then you can swap coplanar faces in meshes
	      // that you read in.
	      // 	  if (pTBFA->pPatchPointer() == pTBFB->pPatchPointer())
	      if (pTBFA->iBdryCond() == pTBFB->iBdryCond())
		return eT22; // case 1
	      else
		return eOther;  // case 2
	    }
	  } // If one is nearly flat
	  // Otherwise, the configuration is convex; swappable two for three.
	  return eT23;
	}
      default:
	// No other cases should be possible
	assert2(0, "Impossible convexity case");
	return eOther;
      }
    } // End of two tet cases
  case 3:
    {
      //@@ Three tet cases
      int iOrient1, iOrient2;
      pVVertA = pVVertE;
      pVVertB = pVVertD;
      pVVertC = const_cast<Vert*>(pVVert(0));
      pVVertD = const_cast<Vert*>(pVVert(1));
      pVVertE = const_cast<Vert*>(pVVert(2));
      if (! apTCTets[2]->qHasVert(pVVertC)) {}
      else if (! apTCTets[2]->qHasVert(pVVertE)) {
	Vert* pVTemp = pVVertC;
	pVVertC = pVVertE;
	pVVertE = pVVertD;
	pVVertD = pVTemp;
      }
      else {
	Vert* pVTemp = pVVertC;
	pVVertC = pVVertD;
	pVVertD = pVVertE;
	pVVertE = pVTemp;
      }

      // The first of these returns will save a few orientation
      // primitive evaluations
      iOrient1 = iOrient3D(pVVertA, pVVertC, pVVertB, pVVertD);
      if (iOrient1 == -1) return eOther;
      else if (iOrient1 == 0) return eN30;
      iOrient2 = iOrient3D(pVVertA, pVVertB, pVVertC, pVVertE);
      if (iOrient2 == -1) return eOther;
      else if (iOrient2 == 0) return eN30;
      assert(iOrient1 == 1);
      assert(iOrient2 == 1);
      //@@@@ T32
      assert(! apTCTets[0]->qHasVert(pVVertA) );
      assert(! apTCTets[1]->qHasVert(pVVertB) );
      assert(! apTCTets[2]->qHasVert(pVVertC) );
      return eT32;
    } // End of three tet cases
  case 4:
    //@@@ Four tets => N40
    return eN40; 
    // End of four tet cases
  } // End of classification/weeding
  // Should never reach here, but just in case...
  return eOther;
}

//@ Quality function for purposes of edge swapping. 
// Would also work for face swapping, but there are more efficient
// approaches for that case.  dQualFunc must return a value which is
// larger for more desirable sets of four verts.
// 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];};
// };

double VolMesh::dQualFunc(Vert* const pV0, Vert* const pV1, 
			  Vert* const pV2, Vert* const pV3) const
{
  Vert *apV[] = {pV0, pV1, pV2, pV3};
  TetCellCV FT(apV);
  switch (eSwapMeasure) {
  case eDelaunay:  // Seek to maximize shortest edge to circumradius.
    {
      double adRes[1];
      int i;
      vShortestToRadius(&FT, &i, adRes);
      return adRes[0];
    }
  case eMaxDihed:  // Seek to minimize max dihedral angle.  Return 
                   // negative of max dihedral.
    return -FT.dMaxDihed();
  case eMinSine:   // Seek to maximize min sine of dihedral angle.
                   // Return min difference of dihedral from 0 or 180
		   // degrees. 
    {
      double adDihed[6];
      int iDum;
      vAllDihed(&FT, &iDum, adDihed);
      assert(iDum == 6);
      double dMinDiff = 90;
      for (iDum = 0; iDum < 6; iDum++) {
        double dDiff = 90 - fabs(90-adDihed[iDum]);
        if (dDiff < dMinDiff) dMinDiff = dDiff;
      }
      return dMinDiff;
    }
  case eFromQual:
    {
      double dValue = pQ->dEvaluate(&FT);
      bool qMaximizeMin = pQ->qMaximizeMinValue();
      return qMaximizeMin ? dValue : -dValue;
    }      
  default:
    assert(0);
  }
  return(-1000);
}

// Pull in the >1200 lines of code that initializes data for the final
// configurations after edge swapping.  
#include "InitCanon.cxx"

//@ Locally reconfigure tets incident on an edge to improve mesh quality
// The edge given by pVNorth, pVSouth is a candidate for removal.  If
// the edge doesn't lie on a boundary and a reconfiguration which
// improves mesh quality exists, then a swap is done, replacing 4 tets
// with 4, 5 with 6, 6 with 8, or 7 with 10.  Generally N tets are
// replaced by 2(N-2).  At this point, the largest transformation
// considered is 7 tets to 10.  The number of possible configurations
// for 8 to 12 swapping is somewhere around 150, so it's not
// considered, at least for now.
//
// After the edge swap is complete, all new faces are face swapped.
// The return value is the total number of face swaps done, with edge
// swaps counting as N-2 face swaps.

// FIX ME: Edge swapping could result in botching internal bdrys
int VolMesh::iEdgeSwap3D(Face* const pF, Vert* pVNorth,
			 Vert* pVSouth, Vert* const pVOther,
			 const bool qAllowAnyValid,
			 const Vert* const pVWant0,
			 const Vert* const pVWant1)
{
  // The last three args are optional and used only for forced edge
  // swapping during surface recovery.  qAllowAnyValid (default:0) is
  // set to 1 to indicate that any submesh which removes the specified
  // edge is preferred.  This implies that the initial quality should be
  // set to some artifically low value and initial quality computation
  // skipped.  If there's some edge that we want in the mesh,
  // the verts of that edge are given as pVWant0 and pVWant1 (default:
  // pVInvalidVert).  A quality bonus is artificially given to
  // configurations which have the magic edge.
  //
  assert(pF->qValid());
  assert(pVNorth->qValid());
  assert(pVSouth->qValid());
  assert(pVOther->qValid());
  assert(pF->qHasVert(pVNorth));
  assert(pF->qHasVert(pVSouth));
  assert(pF->qHasVert(pVOther));
  // This call only needs to appear once, but that's okay, because
  // subsequent calls short-circuit and exit immediately.
  vInitCanonConfigs();

  const double dInvalidPenalty = -10000;
  const double dBonus = 100;

  Vert *apVOuter[10];
  Cell *apCOrig [10];
  Face *apFNorth[10];
  Face *apFSouth[10];
  Face *apFOrig [10];
  // Determine the size of the orbit around the pivot edge and list
  // all the verts, faces, and cells involved.
  apFOrig[0] = pF;
  apVOuter[0] = pVOther;
  apCOrig[0] = pF->pCCellLeft();
  if (apCOrig[0]->eType() == Cell::eTriBFace)
    return iBdryEdgeSwap3D(pF, pVNorth, pVSouth, pVOther);
  if (apCOrig[0]->eType() != Cell::eTet) return 0;
  int iReg = apCOrig[0]->iRegion();
  apFNorth[0] = apCOrig[0]->pFFaceOpposite(pVSouth);
  apFSouth[0] = apCOrig[0]->pFFaceOpposite(pVNorth);
  Face *pFNew = apCOrig[0]->pFFaceOpposite(apVOuter[0]);
  int i = 1;
  while ((i < 10) && (pFNew != pF)) {
    // Abort if there is more than one region involved.
    if (pFNew->iFaceLoc() == Face::eBdryTwoSide) return 0;
    apFOrig[i] = pFNew;
    apVOuter[i] = apCOrig[i-1]->pVVertOpposite(apFOrig[i-1]);
    apCOrig[i] = pFNew->pCCellOpposite(apCOrig[i-1]);
    // Hit a boundary or non-tet cell.
    if (apCOrig[i]->eType() == Cell::eTriBFace)
      return iBdryEdgeSwap3D(pFNew, pVNorth, pVSouth, apVOuter[i]);
    if (apCOrig[i]->eType() != Cell::eTet) return 0;
    // Give up if more than one region is incident on the edge. 
    if (iReg != apCOrig[i]->iRegion()) return 0;
    apFNorth[i] = apCOrig[i]->pFFaceOpposite(pVSouth);
    apFSouth[i] = apCOrig[i]->pFFaceOpposite(pVNorth);
    pFNew = apCOrig[i]->pFFaceOpposite(apVOuter[i]);
    i++;
  }
  if (pFNew != pF) return 0; // Ran out of space
  int iNOrig = i;
  aiEdgeReq[iNOrig]++;
  if (iNOrig > 7) {
    vMessage(3, "Too many tets around edge: %d\n", iNOrig);
    return 0;
  }

  double dOrigQual;
  if (qAllowAnyValid)
    dOrigQual = dInvalidPenalty + 1;
  else {
    // Compute the quality of the current configuration.  This is the
    // smallest quality for any of the existing cells.
    dOrigQual = 1000;
    for (i = 0; i < iNOrig; i++) {
      Cell *pC = apCOrig[i];
      // Evaluate quality, depending on the swap criterion in use.
      double dQual = dQualFunc(pC->pVVert(0), pC->pVVert(1),
			       pC->pVVert(2), pC->pVVert(3));
      dOrigQual = min(dQual, dOrigQual);
    }
  }
  
  // Set up info about the possible new configurations.

  // Compute quality for each possible tet pair (one pair per possible
  // triangle). 
  double a3dNewQual[10][10][10];
  int iSign = iOrient3D(apVOuter[0], apVOuter[1], pVSouth, pVNorth);
  int iV0, iV1, iV2;
#ifndef NDEBUG
  for (iV0 = 0; iV0 < 10; iV0++)
    for (iV1 = 0; iV1 < 10; iV1++)
      for (iV2 = 0; iV2 < 10; iV2++)
	a3dNewQual[iV0][iV1][iV2] = 0;
  for (int iDum = 0; iDum < iNOrig; iDum++) 
    assert(iSign == iOrient3D(apVOuter[iDum], apVOuter[(iDum+1)%iNOrig],
			       pVSouth, pVNorth));
#endif
  for (iV0 = 0; iV0 < iNOrig; iV0++)
    for (iV1 = iV0+1; iV1 < iNOrig; iV1++)
      for (iV2 = iV1+1; iV2 < iNOrig; iV2++) {
        if (iOrient3D(apVOuter[iV0], apVOuter[iV1], apVOuter[iV2],
		      pVNorth) == iSign) 
	  // Evaluate quality
	  a3dNewQual[iV0][iV1][iV2] = dQualFunc(apVOuter[iV0], apVOuter[iV1],
						apVOuter[iV2],  pVNorth);
	else
          a3dNewQual[iV0][iV1][iV2] = dInvalidPenalty; 

	// Don't bother evaluating for the other tet if this one is
	// worse that the worst in the original config, because this tet
	// pair can never be in a final config anyway.
	if (a3dNewQual[iV0][iV1][iV2] >= dOrigQual) {
	  double dDummyQual;
	  if (iOrient3D(apVOuter[iV0], apVOuter[iV2], apVOuter[iV1],
			pVSouth) == iSign)
	    dDummyQual = dQualFunc(apVOuter[iV0], apVOuter[iV2],
				   apVOuter[iV1],  pVSouth);
	  else
	    dDummyQual = dInvalidPenalty;
	  a3dNewQual[iV0][iV1][iV2] = min(a3dNewQual[iV0][iV1][iV2],
					  dDummyQual);
	}

	// Copy things around symmetrically
	a3dNewQual[iV2][iV1][iV0] = a3dNewQual[iV2][iV0][iV1] =
	  a3dNewQual[iV1][iV2][iV0] = a3dNewQual[iV1][iV0][iV2] =
	  a3dNewQual[iV0][iV2][iV1] = a3dNewQual[iV0][iV1][iV2];
	
      }
  
  // Find the new configuration with the highest quality.
  double dMaxQual = dOrigQual;
  int iCanon, iPerm;
  int iBestCanon = -1, iBestPerm = -1;
  bool qBonusPossible = (pVWant0 != pVInvalidVert);
  for (iCanon = 0; iCanon < aCS[iNOrig].iNCanon; iCanon++)
    for (iPerm = 0; iPerm < aCS[iNOrig].aConf[iCanon].iNPerm; iPerm++) {
      int iBonusGiven = 0;
      double dConfQual = 1000;
      // Once this config is known to be worse than the best one we know
      // about, we can give up on it.
      for (int iTri = 0;
	   iTri < iNOrig-2 && (dConfQual >= dMaxQual || qBonusPossible);
	   iTri++) {
        iV0 = (aCS[iNOrig].aConf[iCanon].a2iFaceVert[iTri][0] + iPerm)
	  % iNOrig;
        iV1 = (aCS[iNOrig].aConf[iCanon].a2iFaceVert[iTri][1] + iPerm)
	  % iNOrig;
        iV2 = (aCS[iNOrig].aConf[iCanon].a2iFaceVert[iTri][2] + iPerm)
	  % iNOrig;
	// Give a bonus if a particular edge is desired and that edge
	// is included in the configuration.
	if (qBonusPossible && ! (iBonusGiven == 2)) {
	  int iThisBonus = 0;
	  Vert * pV0 = apVOuter[iV0];
	  Vert * pV1 = apVOuter[iV1];
	  Vert * pV2 = apVOuter[iV2];
	  if (pVWant0 == pV0 || pVWant0 == pV1 || pVWant0 == pV2)
	    iThisBonus++;
	  if (pVWant1 == pV0 || pVWant1 == pV1 || pVWant1 == pV2)
	    iThisBonus++;
	  if (iThisBonus > iBonusGiven) iBonusGiven = iThisBonus;
	}
	dConfQual = min(dConfQual, a3dNewQual[iV0][iV1][iV2]);
	// Only give a bonus to valid configurations
	if (qBonusPossible && dConfQual > dOrigQual)
	  dConfQual += iBonusGiven*dBonus;
      }
      if (dConfQual > dMaxQual) {
	dMaxQual = dConfQual;
	iBestCanon = iCanon;
	iBestPerm = iPerm;
      }   
    }
  
  // If the highest-quality alternative is worse than the current
  // configuration, do nothing.  The check for iBestCanon is there so
  // that some stupid roundoff in comparison of dMaxQual and dOrigQual
  // won't cause a reconfiguration when there shouldn't be one.
  if (dMaxQual <= dOrigQual || iBestCanon == -1) {
    return 0;
  }
  else {
    // If the highest-quality alternative is better than the current
    // configuration:
    aiEdgeDone[iNOrig]++;
    
    // Set up arrays to hold all the interior faces and cells
    Face** apFNew = new Face*[3*iNOrig-8];
    Cell** apCNew = new Cell*[2*iNOrig-4];
    
    if (iNOrig == 3) {
      apFNew[0] = apFOrig[0];
      apFOrig[1]->vMarkDeleted();
      apFOrig[2]->vMarkDeleted();
      apCNew[0] = apCOrig[0];
      apCNew[1] = apCOrig[1];
      apCOrig[2]->vMarkDeleted();
    }
    else {
      // List interior faces, including setting up new ones.
      for (i = 0; i < iNOrig; i++) apFNew[i] = apFOrig[i];
      for (i = iNOrig; i < 3*iNOrig-8; i++) 
	apFNew[i] = pFNewFace();

      // List interior cells, including setting up new ones.
      for (i = 0; i < iNOrig; i++) apCNew[i] = apCOrig[i];
      for (i = iNOrig; i < 2*iNOrig-4; i++) {
	apCNew[i] = pCNewCell();
	apCNew[i]->vSetRegion(apCOrig[0]->iRegion());
      }
    }

    // Renumber all the old face, vertex and cell data, offsetting it
    // by iBestPerm to make the desired configuration match the stored
    // canonical configuration.  
    if (iBestPerm != 0) {
      Cell **apCTemp = new Cell*[iNOrig];

      Face **apFNorthTemp = new Face*[iNOrig];
      Face **apFSouthTemp = new Face*[iNOrig];
      Vert **apVOuterTemp = new Vert*[iNOrig];
      // Copy data to temporary arrays.
      for (i = 0; i < iNOrig; i++) {
        apCTemp[i] = apCOrig[i];
        apFNorthTemp[i] = apFNorth[i];
        apFSouthTemp[i] = apFSouth[i];
        apVOuterTemp[i] = apVOuter[i];
      }
      // Copy back with an offset.
      for (i = 0; i < iNOrig; i++) {
        apCOrig[i]  = apCTemp     [(i+iBestPerm)%iNOrig];
        apFNorth[i] = apFNorthTemp[(i+iBestPerm)%iNOrig];
        apFSouth[i] = apFSouthTemp[(i+iBestPerm)%iNOrig];
        apVOuter[i] = apVOuterTemp[(i+iBestPerm)%iNOrig];
      }
      delete [] apCTemp;
      delete [] apFNorthTemp;
      delete [] apFSouthTemp;
      delete [] apVOuterTemp;
    }

    // Disconnect outer faces from cells that are about to be rearranged
    for (i = 0; i < iNOrig; i++) {
      apFNorth[i]->vRemoveCell(apCOrig[i]);
      apFSouth[i]->vRemoveCell(apCOrig[i]);
    }

    // Re-connect the mesh locally
    Config *pConf = &(aCS[iNOrig].aConf[iBestCanon]);
    // Assign data to new faces
    for (i = 0; i < 3*iNOrig-8; i++) {
      Vert *pV0, *pV1, *pV2;
      iV0 = pConf->a2iFaceVert[i][0];
      switch (iV0) {
      case NORTH_POLE: pV0 = pVNorth; break;
      case SOUTH_POLE: pV0 = pVSouth; break;
      default:         pV0 = apVOuter[iV0]; break;
      }
      iV1 = pConf->a2iFaceVert[i][1];
      switch (iV1) {
      case NORTH_POLE: pV1 = pVNorth; break;
      case SOUTH_POLE: pV1 = pVSouth; break;
      default:         pV1 = apVOuter[iV1]; break;
      }
      iV2 = pConf->a2iFaceVert[i][2];
      switch (iV2) {
      case NORTH_POLE: pV2 = pVNorth; break;
      case SOUTH_POLE: pV2 = pVSouth; break;
      default:         pV2 = apVOuter[iV2]; break;
      }
      Cell *pC0 = apCNew[pConf->a2iFaceCell[i][0]];
      Cell *pC1 = apCNew[pConf->a2iFaceCell[i][1]];
      if (iSign == 1) 
	apFNew[i]->vAssign(pC0, pC1, pV0, pV1, pV2);
      else 
	apFNew[i]->vAssign(pC1, pC0, pV0, pV1, pV2);
      
      pV0->vSetHintFace(apFNew[i]);
      pV1->vSetHintFace(apFNew[i]);
      pV2->vSetHintFace(apFNew[i]);
    }

    // Cell setup, including re-connection to outer faces
    for (i = 0; i < 2*iNOrig-4; i++) {
      Face *apF[4];
      for (int ii = 0; ii < 4; ii++) {
	int iFace = pConf->a2iCellFace[i][ii];
	if (iFace >= SOUTH_FACE) {
	  apF[ii] = apFSouth[iFace - SOUTH_FACE];
	  apF[ii]->vAddCell(apCNew[i]);
	}
	else if (iFace >= NORTH_FACE) {
	  apF[ii] = apFNorth[iFace - NORTH_FACE];
	  apF[ii]->vAddCell(apCNew[i]);
	}
	else
	  apF[ii] = apFNew[iFace];
      }
      apCNew[i]->vAssign(apF);
    }

    // Reset vertex hints for polar verts
    pVNorth->vSetHintFace(apFNorth[0]);
    pVSouth->vSetHintFace(apFSouth[0]);

    // Mark all verts for smoothing
    pVNorth->vMarkIllShaped();
    pVSouth->vMarkIllShaped();
    for (i = 0; i < iNOrig; i++)
      apVOuter[i]->vMarkIllShaped();

    // Validate the result
    for (i = 0; i < iNOrig*2 - 4; i++) {
      // Check everything about the tet
      assert(apCNew[i]->iFullCheck());
      // Every tet has to hit one pole or the other
      assert(apCNew[i]->qHasVert(pVNorth) ||
	     apCNew[i]->qHasVert(pVSouth));
      // Assert tetrahedron orientation
      assert(iOrient3D(apCNew[i]->pVVert(0), apCNew[i]->pVVert(1),
		       apCNew[i]->pVVert(2), apCNew[i]->pVVert(3)) == 1);
    }

    for (i = 0; i < 3*iNOrig - 8 ; i++) {
      assert(apFNew[i]->pCCellLeft()->qValid());
      assert(apFNew[i]->pCCellRight()->qValid());
      assert(apFNew[i]->iFullCheck());
    }

    for (i = 0; i < iNOrig; i++) {    
      assert(apFSouth[i]->iFullCheck());
      assert(apFNorth[i]->iFullCheck());
    }

    // If recursion is permitted, check all the new faces for possible
    // swaps, as well as all of the faces on the outside of the orbit.
    int iRetVal = iNOrig-2;
    iSwapDepth++;
    if (iSwapDepth > iMaxSwapDepth)
      vMessage(2, "Swap recursion depth exceeded!\n");
    if (qSwapRecur && iSwapDepth <= iMaxSwapDepth) {
      for (i = 0; i < 3*iNOrig-8; i++) {
	iRetVal += iFaceSwap(apFNew[i]);
      }
      for (i = 0; i < iNOrig; i++) {
	iRetVal += iFaceSwap(apFNorth[i]) + iFaceSwap(apFSouth[i]);
      }
    }
    else {
      // If you aren't going to swap them, mark them for future
      // consideration. 
      for (i = 0; i < 3*iNOrig-8; i++) {
	apFNew[i]->vMarkForSwap();
      }
      for (i = 0; i < iNOrig; i++) {
	apFNorth[i]->vMarkForSwap();
	apFSouth[i]->vMarkForSwap();
      }
    }

    delete [] apCNew;
    delete [] apFNew;
    // Return the total (equivalent) number of swaps performed
    iSwapDepth--;
    return iRetVal;
  }
}

//@ Locally reconfigure tets incident on a bdry edge to improve mesh quality
// The bdry edge given by pVNorth, pVSouth is a candidate for removal.
// pF is required to be a boundary face, and pVOther is the third vertex
// on that face.
// 
// If the bdry edge has two incident bdry faces that are coplanar and
// have the same boundary condition, and if a reconfiguration which
// improves mesh quality exists, then a swap is done, replacing 2 tets
// with 2, 3 with 4, ... or 6 with 10.  Generally N tets are replaced by
// 2(N-1).  At this point, the largest transformation considered is 6
// tets to 10.  The number of possible configurations for 7 to 12
// swapping is somewhere around 150, so it's not considered, at least
// for now. 
//
// All of these reconfigurations are close cousins to interior edge
// swaps.  A 2-to-2 boundary swap is just a 3-to-2 swap with one tet
// flattened into oblivion.  This means that all of the canonical
// configuration info from interior can be re-used here.
//
// After the edge swap is complete, all new interior faces are face
// swapped. The return value is the total number of face swaps done,
// with boundary edge swaps counting as N-3 face swaps. 
int VolMesh::iBdryEdgeSwap3D(Face* const pF, Vert* pVNorth,
			     Vert* pVSouth, Vert* const pVOther)
{
  // Note the optional args used for forced edge swapping during surface
  // recovery aren't used here because they aren't needed.
  //
  assert(pF->eType() == Face::eTriFace);
  assert(pF->qValid());
  assert(pVNorth->qValid());
  assert(pVSouth->qValid());
  assert(pVOther->qValid());
  assert(pF->qHasVert(pVNorth));
  assert(pF->qHasVert(pVSouth));
  assert(pF->qHasVert(pVOther));
  // The face must be on the boundary.
  assert(pF->pCCellLeft()->eType()  == Cell::eTriBFace ||
	 pF->pCCellRight()->eType() == Cell::eTriBFace);

  // This call only needs to appear once, but that's okay, because
  // subsequent calls short-circuit and exit immediately.
  vInitCanonConfigs();
  if (!qOKToChangeSurfaceEdges // Surface edge swapping turned off
      || !qAllowBdryChanges)   // All surface changes turned off
    return 0;

  const double dInvalidPenalty = -10000;

  //@@ Identify all cells and faces incident on the bad edge
  Vert *apVOuter[11];
  Cell *apCOrig [11];
  Face *apFNorth[11];
  Face *apFSouth[11];
  Face *apFOrig [11];
  BFace *pBFBFace0, *pBFBFace1;

  // Determine the size of the orbit around the pivot edge and list
  // all the verts, faces, and cells involved.
  apFOrig[0] = pF;
  apVOuter[0] = pVOther;
  // ID the boundary associated with the original face *pF
  {
    Cell *pCTemp = pF->pCCellLeft();
    if (pCTemp->eType() == Cell::eTriBFace) {
      pBFBFace0 = dynamic_cast<BFace*>(pCTemp);
      apCOrig[0] = pF->pCCellRight();
    }
    else {
      pBFBFace0 = dynamic_cast<BFace*>(pF->pCCellRight());
      apCOrig[0] = pCTemp;
    }
  }
  apFNorth[0] = apCOrig[0]->pFFaceOpposite(pVSouth);
  apFSouth[0] = apCOrig[0]->pFFaceOpposite(pVNorth);
  Face *pFNew = apCOrig[0]->pFFaceOpposite(apVOuter[0]);
  int i = 1;
  while (i <= 10) {
    // Abort if there is more than one region involved.
    if (pFNew->iFaceLoc() == Face::eBdryTwoSide) return 0;
    apFOrig[i] = pFNew;
    apVOuter[i] = apCOrig[i-1]->pVVertOpposite(apFOrig[i-1]);
    apCOrig[i] = pFNew->pCCellOpposite(apCOrig[i-1]);
    // Hit a boundary or non-tet cell.
    if (apCOrig[i]->eType() == Cell::eTriBFace) break;
    else if (apCOrig[i]->eType() != Cell::eTet) return 0;
    // Give up now if there is more than one region incident on the edge.
    if (apCOrig[i]->iRegion() != apCOrig[0]->iRegion()) return 0;
    apFNorth[i] = apCOrig[i]->pFFaceOpposite(pVSouth);
    apFSouth[i] = apCOrig[i]->pFFaceOpposite(pVNorth);
    pFNew = apCOrig[i]->pFFaceOpposite(apVOuter[i]);
    i++;
  }
  // Assign data for the second boundary face
  pBFBFace1 = dynamic_cast<BFace*>(apCOrig[i]);

  int iNOrigCells = i; 
  int iNOrigVerts = iNOrigCells + 1;
  aiEdgeReq[iNOrigCells]++;

  // A set of easy exits for cases where surface edge swapping is not
  // legal.

  // Did we run out of space while identifying the entities around the
  // edge we'd like to swap?
  if (apCOrig[i]->eType() != Cell::eTriBFace) return 0;

  // Do we have more cells incident on the edge than we're prepared to
  // handle?
  if (iNOrigCells > 6) {
    vMessage(4, "Too many tets around bdry edge: %d\n", iNOrigCells);
    return 0;
  }

  // Do the boundary faces have the same BC?
  if (!qBdryFacesOKToSwap(dynamic_cast<TriBFaceBase*>(pBFBFace0),
			  dynamic_cast<TriBFaceBase*>(pBFBFace1)))
    return 0;

  // Are the boundary faces near enough to being coplanar?

  // Strict coplanarity test for surface edge swapping.  Tested by
  // checking whether the four surface points in the configuration form
  // a degenerate tet. 

//   if (iOrient3D(pVNorth, pVSouth, pVOther, apVOuter[iNOrigVerts-1]) != 0)
//     return 0;

  // Allow some slightly non-planar surface edge swapping.
  {
    double adNorm1[3], adNorm2[3];
    vUnitNormal(pVNorth->adCoords(), pVSouth->adCoords(),
		pVOther->adCoords(), adNorm1);
    vUnitNormal(pVNorth->adCoords(), apVOuter[iNOrigVerts-1]->adCoords(),
		pVSouth->adCoords(), adNorm2);
    double dCos = dDOT3D(adNorm1, adNorm2);
    double dMinCos = cos(dMaxAngleForSurfSwap * M_PI / 180.);
    if (iFuzzyComp(dCos, dMinCos) == -1) return 0;
  }

  // Do the two boundary triangles form a convex quad?  If so, the
  // orientation of the tets formed by these triangles and any
  // non-coplanar point (i.e., any other member of apVOuter) will be the
  // same.  This exit is redundant (no swapping would occur for this
  // case anyway), but saves a lot of quality computations.
  if (iOrient3D(pVOther, apVOuter[iNOrigVerts-1], pVNorth, apVOuter[1])
      !=
      iOrient3D(pVOther, pVSouth, apVOuter[iNOrigVerts-1], apVOuter[1])
      ) return 0;

  double dOrigQual;
  //@@ Compute the quality of the current configuration.
  // This is the smallest quality for any of the existing cells.
  dOrigQual = 1000;
  for (i = 0; i < iNOrigCells; i++) {
    Cell *pC = apCOrig[i];
    // Evaluate quality, depending on the swap criterion in use.
    double dQual = dQualFunc(pC->pVVert(0), pC->pVVert(1),
			     pC->pVVert(2), pC->pVVert(3));
    dOrigQual = min(dQual, dOrigQual);
  }
  
  // Set up info about the possible new configurations.

  // Compute quality for each possible tet pair (one pair per possible
  // triangle). 
  double a3dNewQual[10][10][10];
  int iSign = iOrient3D(apVOuter[0], apVOuter[1], pVSouth, pVNorth);
  int iV0, iV1, iV2;
#ifndef NDEBUG
  for (iV0 = 0; iV0 < 10; iV0++)
    for (iV1 = 0; iV1 < 10; iV1++)
      for (iV2 = 0; iV2 < 10; iV2++)
	a3dNewQual[iV0][iV1][iV2] = 0;
  for (int iDum = 0; iDum < iNOrigCells; iDum++) 
    assert(iSign == iOrient3D(apVOuter[iDum], apVOuter[iDum+1],
			       pVSouth, pVNorth));
#endif
  for (iV0 = 0; iV0 < iNOrigVerts; iV0++)
    for (iV1 = iV0+1; iV1 < iNOrigVerts; iV1++)
      for (iV2 = iV1+1; iV2 < iNOrigVerts; iV2++) {
        if (iOrient3D(apVOuter[iV0], apVOuter[iV1], apVOuter[iV2],
		      pVNorth) == iSign) 
	  // Evaluate quality
	  a3dNewQual[iV0][iV1][iV2] = dQualFunc(apVOuter[iV0], apVOuter[iV1],
						apVOuter[iV2],  pVNorth);
	else
          a3dNewQual[iV0][iV1][iV2] = dInvalidPenalty; 

	// Don't bother evaluating for the other tet if this one is
	// worse that the worst in the original config, because this tet
	// pair can never be in a final config anyway.
	if (a3dNewQual[iV0][iV1][iV2] >= dOrigQual) {
	  double dDummyQual;
	  if (iOrient3D(apVOuter[iV0], apVOuter[iV2], apVOuter[iV1],
			pVSouth) == iSign)
	    dDummyQual = dQualFunc(apVOuter[iV0], apVOuter[iV2],
				   apVOuter[iV1],  pVSouth);
	  else
	    dDummyQual = dInvalidPenalty;
	  a3dNewQual[iV0][iV1][iV2] = min(a3dNewQual[iV0][iV1][iV2],
					  dDummyQual);
	}

	// Copy things around symmetrically
	a3dNewQual[iV2][iV1][iV0] = a3dNewQual[iV2][iV0][iV1] =
	  a3dNewQual[iV1][iV2][iV0] = a3dNewQual[iV1][iV0][iV2] =
	  a3dNewQual[iV0][iV2][iV1] = a3dNewQual[iV0][iV1][iV2];
	
      }
  
  // Find the new configuration with the highest quality.
  double dMaxQual = dOrigQual;
  int iCanon, iPerm;
  int iBestCanon = -1, iBestPerm = -1;
  for (iCanon = 0; iCanon < aCS[iNOrigVerts].iNCanon; iCanon++)
    for (iPerm = 0; iPerm < aCS[iNOrigVerts].aConf[iCanon].iNPerm; iPerm++) {
      double dConfQual = 1000;
      // Once this config is known to be worse than the best one we know
      // about, we can give up on it.
      for (int iTri = 0;
	   iTri < iNOrigVerts-2 && (dConfQual >= dMaxQual);
	   iTri++) {
        iV0 = (aCS[iNOrigVerts].aConf[iCanon].a2iFaceVert[iTri][0] + iPerm)
	  % iNOrigVerts;
        iV1 = (aCS[iNOrigVerts].aConf[iCanon].a2iFaceVert[iTri][1] + iPerm)
	  % iNOrigVerts;
        iV2 = (aCS[iNOrigVerts].aConf[iCanon].a2iFaceVert[iTri][2] + iPerm)
	  % iNOrigVerts;
	dConfQual = min(dConfQual, a3dNewQual[iV0][iV1][iV2]);
      }
      if (dConfQual > dMaxQual) {
	dMaxQual = dConfQual;
	iBestCanon = iCanon;
	iBestPerm = iPerm;
      }   
    }
  
  // If the highest-quality alternative is worse than the current
  // configuration, do nothing.  The check for iBestCanon is there so
  // that some stupid roundoff in comparison of dMaxQual and dOrigQual
  // won't cause a reconfiguration when there shouldn't be one.
  if (dMaxQual <= dOrigQual || iBestCanon == -1) {
    return 0;
  }
  else {
    // If the highest-quality alternative is better than the current
    // configuration:
    aiEdgeDone[iNOrigCells]++;

    // Set up arrays to hold all the interior faces and cells
    Face** apFNew = new Face*[3*iNOrigVerts-8];
    Cell** apCNew = new Cell*[2*iNOrigVerts-4];
    
    if (iNOrigVerts == 3) {
      apFNew[0] = apFOrig[0];
      apFOrig[1]->vMarkDeleted();
      apFOrig[2]->vMarkDeleted();
      apCNew[0] = apCOrig[0];
      apCNew[1] = apCOrig[1];
    }
    else {
      // List interior faces, including setting up new ones.
      for (i = 0; i < iNOrigVerts; i++) apFNew[i] = apFOrig[i];
      for (i = iNOrigVerts; i < 3*iNOrigVerts-8; i++) 
	apFNew[i] = pFNewFace();

      // List interior cells, including setting up new ones.
      for (i = 0; i < iNOrigCells; i++) apCNew[i] = apCOrig[i];
      for (i = iNOrigCells; i < 2*iNOrigVerts-4; i++) {
	apCNew[i] = pCNewCell();
	apCNew[i]->vSetRegion(apCOrig[0]->iRegion());
      }
    }

    // Disconnect outer faces from cells that are about to be rearranged
    for (i = 0; i < iNOrigVerts-1; i++) {
      apFNorth[i]->vRemoveCell(apCOrig[i]);
      apFSouth[i]->vRemoveCell(apCOrig[i]);
    }

    // Create a new North face and a new South face to hold the
    // reconfigured boundary faces.  Also, connect them to the existing
    // boundary data.

    apFNorth[iNOrigVerts-1] = pFNewFace();
    apFSouth[iNOrigVerts-1] = pFNewFace(); 
    if (iSign == 1) { // apVOuter[0], apVOuter[1], pVSouth, pVNorth is a
		      // right-handed tet.
      apFNorth[iNOrigVerts-1]->vAssign(pBFBFace0, pCInvalidCell,
				       apVOuter[0],
				       apVOuter[iNOrigVerts-1],
				       pVNorth);
      apFSouth[iNOrigVerts-1]->vAssign(pBFBFace1, pCInvalidCell,
				       apVOuter[iNOrigVerts-1],
				       apVOuter[0],
				       pVSouth);
    }
    else {
      apFNorth[iNOrigVerts-1]->vAssign(pCInvalidCell, pBFBFace0,
				       apVOuter[0],
				       apVOuter[iNOrigVerts-1],
				       pVNorth);
      apFSouth[iNOrigVerts-1]->vAssign(pCInvalidCell, pBFBFace1,
				       apVOuter[iNOrigVerts-1],
				       apVOuter[0],
				       pVSouth);
    }
    pBFBFace0->vAssign(apFNorth[iNOrigVerts-1]);
    pBFBFace1->vAssign(apFSouth[iNOrigVerts-1]);
    
    // Renumber all the old face, vertex and cell data, offsetting it
    // by iBestPerm to make the desired configuration match the stored
    // canonical configuration.  
    if (iBestPerm != 0) {
      Face **apFNorthTemp = new Face*[iNOrigVerts];
      Face **apFSouthTemp = new Face*[iNOrigVerts];
      Vert **apVOuterTemp = new Vert*[iNOrigVerts];
      // Copy data to temporary arrays.
      for (i = 0; i < iNOrigVerts; i++) {
        apFNorthTemp[i] = apFNorth[i];
        apFSouthTemp[i] = apFSouth[i];
        apVOuterTemp[i] = apVOuter[i];
      }
      // Copy back with an offset.
      for (i = 0; i < iNOrigVerts; i++) {
        apFNorth[i] = apFNorthTemp[(i+iBestPerm)%iNOrigVerts];
        apFSouth[i] = apFSouthTemp[(i+iBestPerm)%iNOrigVerts];
        apVOuter[i] = apVOuterTemp[(i+iBestPerm)%iNOrigVerts];
      }
      delete [] apFNorthTemp;
      delete [] apFSouthTemp;
      delete [] apVOuterTemp;
    }

    // Re-connect the mesh locally
    Config *pConf = &(aCS[iNOrigVerts].aConf[iBestCanon]);
    // Assign data to new faces
    for (i = 0; i < 3*iNOrigVerts-8; i++) {
      Vert *pV0, *pV1, *pV2;
      iV0 = pConf->a2iFaceVert[i][0];
      switch (iV0) {
      case NORTH_POLE: pV0 = pVNorth; break;
      case SOUTH_POLE: pV0 = pVSouth; break;
      default:         pV0 = apVOuter[iV0]; break;
      }
      iV1 = pConf->a2iFaceVert[i][1];
      switch (iV1) {
      case NORTH_POLE: pV1 = pVNorth; break;
      case SOUTH_POLE: pV1 = pVSouth; break;
      default:         pV1 = apVOuter[iV1]; break;
      }
      iV2 = pConf->a2iFaceVert[i][2];
      switch (iV2) {
      case NORTH_POLE: pV2 = pVNorth; break;
      case SOUTH_POLE: pV2 = pVSouth; break;
      default:         pV2 = apVOuter[iV2]; break;
      }
      Cell *pC0 = apCNew[pConf->a2iFaceCell[i][0]];
      Cell *pC1 = apCNew[pConf->a2iFaceCell[i][1]];
      if (iSign == 1) 
	apFNew[i]->vAssign(pC0, pC1, pV0, pV1, pV2);
      else 
	apFNew[i]->vAssign(pC1, pC0, pV0, pV1, pV2);
      
      pV0->vSetHintFace(apFNew[i]);
      pV1->vSetHintFace(apFNew[i]);
      pV2->vSetHintFace(apFNew[i]);
    }

    // Cell setup, including re-connection to outer faces
    for (i = 0; i < 2*iNOrigVerts-4; i++) {
      Face *apF[4];
      for (int ii = 0; ii < 4; ii++) {
	int iFace = pConf->a2iCellFace[i][ii];
	if (iFace >= SOUTH_FACE) {
	  apF[ii] = apFSouth[iFace - SOUTH_FACE];
	  apF[ii]->vAddCell(apCNew[i]);
	}
	else if (iFace >= NORTH_FACE) {
	  apF[ii] = apFNorth[iFace - NORTH_FACE];
	  apF[ii]->vAddCell(apCNew[i]);
	}
	else
	  apF[ii] = apFNew[iFace];
      }
      apCNew[i]->vAssign(apF);
    }

    // Reset vertex hints for polar verts
    pVNorth->vSetHintFace(apFNorth[0]);
    pVSouth->vSetHintFace(apFSouth[0]);

    // Mark all verts for smoothing
    pVNorth->vMarkIllShaped();
    pVSouth->vMarkIllShaped();
    for (i = 0; i < iNOrigVerts; i++)
      apVOuter[i]->vMarkIllShaped();

    // Validate the result
    for (i = 0; i < iNOrigVerts*2 - 4; i++) {
      // Check everything about the tet
      assert(apCNew[i]->iFullCheck());
      // Every tet has to hit one pole or the other
      assert(apCNew[i]->qHasVert(pVNorth) ||
	     apCNew[i]->qHasVert(pVSouth));
      // Assert tetrahedron orientation
      assert(iOrient3D(apCNew[i]->pVVert(0), apCNew[i]->pVVert(1),
		       apCNew[i]->pVVert(2), apCNew[i]->pVVert(3)) == 1);
    }

    for (i = 0; i < 3*iNOrigVerts - 8 ; i++)
      assert(apFNew[i]->iFullCheck());

    for (i = 0; i < iNOrigVerts; i++) {
      assert(apFSouth[i]->iFullCheck());
      assert(apFNorth[i]->iFullCheck());
    }

    // If recursion is permitted, check all the new faces for possible
    // swaps, as well as all of the faces on the outside of the orbit.
    int iRetVal = iNOrigCells-1;
    iSwapDepth++;
    if (qSwapRecur && iSwapDepth <= iMaxSwapDepth) {
      for (i = 0; i < 3*iNOrigVerts-8; i++) {
	iRetVal += iFaceSwap(apFNew[i]);
      }
      for (i = 0; i < iNOrigVerts; i++) {
	iRetVal += iFaceSwap(apFNorth[i]) + iFaceSwap(apFSouth[i]);
      }
    }
    else {
      // If you aren't going to swap them, mark them for future
      // consideration. 
      for (i = 0; i < 3*iNOrigVerts-8; i++) {
	apFNew[i]->vMarkForSwap();
      }
      for (i = 0; i < iNOrigVerts; i++) {
	apFNorth[i]->vMarkForSwap();
	apFSouth[i]->vMarkForSwap();
      }
    }

    delete [] apCNew;
    delete [] apFNew;
    // Return the total (equivalent) number of swaps performed
    iSwapDepth--;
    return iRetVal;
  }
}

//@ Swap away face if this is legal and improves the mesh quality measure
// Swapping typically propagates through the mesh, and the return value from
// iFaceSwap is the total number of swaps done during this invocation.
int VolMesh::iFaceSwap(Face* const pF)
{
  // If the face has been removed, don't bother 
  if (!pF->iFullCheck()) return(0);

  // If this is an internal boundary face, we can't swap it.
  if (pF->iFaceLoc() == Face::eBdryTwoSide) return(0);

  if (iSwapDepth > iMaxSwapDepth) return(0);
  iSwapDepth++;

  Vert *pVVertA, *pVVertB, *pVVertC, *pVVertD, *pVVertE;
  Vert *pVPivot0, *pVPivot1, *pVOther;
  TetCell* apTCTets[4];
  int iNTets;
  TriFace* pTF = (dynamic_cast<TriFace*>(pF));
  if (!pTF->qValid()) return 0;
  eFC = pTF->eCategorizeFace
    (pVVertA, pVVertB, pVVertC, pVVertD, pVVertE, apTCTets, iNTets,
     pVPivot0, pVPivot1, pVOther); 
  
  //@@ Determine the preferable configuration and swap if necessary
  switch (eFC) {
  default: // Catch-all for bad cases
    vWarning("Shouldn't be here: face wasn't categorized");
    pF->vClearSwapFlag();
    iSwapDepth--;
    return 0;
  // These cases are handled by edge swapping.
  case eN44:
  case eN32:
    {
      int iRetVal = 0;
      if (eSwapType() != eDelaunay && qEdgeSwappingAllowed())
	iRetVal = iEdgeSwap3D(pF, pVPivot0, pVPivot1, pVOther);
      iSwapDepth--;
      return iRetVal;
    }
    // These cases are definitely unswappable
  case eBdry:
  case eN40:
  case eN30:
  case eN20:
  case eBdryReflex:
  case eOther:
    pF->vClearSwapFlag();
    iSwapDepth--;
    return(0);
  case eT44:
  case eT22:
  case eT23:
  case eT32:
    if (qDoSwap(pVVertA, pVVertB, pVVertC, pVVertD, pVVertE)) {
      int iRes =  iReconfigure(apTCTets, dynamic_cast<TriFace*>(pF), pVVertA,
			       pVVertB, pVVertC, pVVertD, pVVertE);
      iSwapDepth--;
      return iRes;
    }
    else {
      pF->vClearSwapFlag();
      iSwapDepth--;
      return(0);
    }
  }
}

//@@ Determines whether swapping is needed for local Delaunay-ness.
// Returns true to swap, false to remain the same.
static bool
qDoSwapDelaunay(const Vert* pVVertA, const Vert* pVVertB,
		const Vert* pVVertC, const Vert* pVVertD,
		const Vert* pVVertE)
{
  int iInSphereRes = iInsphere(pVVertA, pVVertB, pVVertC, pVVertD, pVVertE);

  switch (eFC) {
  case eT32:
    // Current tets are ABDE, BCDE, CADE.  Swap if pVVertE is on (0) or
    // outside (-1) the circumsphere of ABCD.  This is an arbitrary
    // tie-break; configurations which are co-spherical go to the two
    // tet configuration.
    return (iInSphereRes != 1);
  case eT23:
  case eT22:
  case eT44: // This one is the same as eT22 for all situations.
    // Current tets are ABCD, CBAE.  Swap if pVVertE is inside (1) the
    // circumsphere of ABCD.  Again, configurations which are
    // co-spherical go to the two-tet configuration.
    return (iInSphereRes == 1);
  default:
    return(false);
  }
}

//@@ Determines whether swapping will improve the maximum face angle.
// Returns true to swap, false to remain the same.
static bool
qDoSwapMaxDihed(const Vert* pVVertA, const Vert* pVVertB,
		const Vert* pVVertC, const Vert* pVVertD,
		const Vert* pVVertE)
{
  // This algorithm finds face angles between adjacent faces by dotting
  // their unit normals.  The largest magnitude loses.
  //
  // To prevent pairs from flopping back and forth, a swap is only
  // authorized if the inequality is bigger than eps.

  double dEps = 1.e-10;
  switch (eFC) {
  case eT23:
  case eT32:
    {
      //@@@ The case where the five pVVerts form a genuinely convex set
      double adNormABD[3], adNormABE[3], adNormADE[3];
      double adNormBCD[3], adNormBCE[3], adNormBDE[3];
      double adNormCAD[3], adNormCAE[3], adNormCDE[3];
      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertD->adCoords(), adNormABD);
      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertE->adCoords(), adNormABE);
      vUnitNormal(pVVertA->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormADE);

      vUnitNormal(pVVertB->adCoords(), pVVertC->adCoords(),
		  pVVertD->adCoords(), adNormBCD);
      vUnitNormal(pVVertB->adCoords(), pVVertC->adCoords(),
		  pVVertE->adCoords(), adNormBCE);
      vUnitNormal(pVVertB->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormBDE);

      vUnitNormal(pVVertC->adCoords(), pVVertA->adCoords(),
		  pVVertD->adCoords(), adNormCAD);
      vUnitNormal(pVVertC->adCoords(), pVVertA->adCoords(),
		  pVVertE->adCoords(), adNormCAE);
      vUnitNormal(pVVertC->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormCDE);

      double dLastDot, dWorst2 = -1., dWorst3 = -1.;

      // Care must be taken to be certain that I'm taking the dot
      // product of two vectors that both point into the enclosed angle
      // or both point out, or correct the sign.

      dLastDot = dDOT3D(adNormABD, adNormBCD);   /* Both in */
      dWorst2 = max(dWorst2, dLastDot);
      dLastDot = dDOT3D(adNormBCD, adNormCAD);   /* Both in */
      dWorst2 = max(dWorst2, dLastDot);
      dLastDot = dDOT3D(adNormCAD, adNormABD);   /* Both in */
      dWorst2 = max(dWorst2, dLastDot);

      dLastDot = dDOT3D(adNormABE, adNormBCE);   /* Both out */
      dWorst2 = max(dWorst2, dLastDot);
      dLastDot = dDOT3D(adNormBCE, adNormCAE);   /* Both out */
      dWorst2 = max(dWorst2, dLastDot);
      dLastDot = dDOT3D(adNormCAE, adNormABE);   /* Both out */
      dWorst2 = max(dWorst2, dLastDot);

      dLastDot = - dDOT3D(adNormABD, adNormABE); /* One in, one out */
      dWorst3 = max(dWorst3, dLastDot);
      dLastDot = - dDOT3D(adNormBCD, adNormBCE); /* One in, one out */
      dWorst3 = max(dWorst3, dLastDot);
      dLastDot = - dDOT3D(adNormCAD, adNormCAE); /* One in, one out */
      dWorst3 = max(dWorst3, dLastDot);

      dLastDot = - dDOT3D(adNormADE, adNormBDE); /* One in, one out */
      dWorst3 = max(dWorst3, dLastDot);
      dLastDot = - dDOT3D(adNormBDE, adNormCDE); /* One in, one out */
      dWorst3 = max(dWorst3, dLastDot);
      dLastDot = - dDOT3D(adNormCDE, adNormADE); /* One in, one out */
      dWorst3 = max(dWorst3, dLastDot);

      if      (dWorst3 > dWorst2 + dEps) return(eFC == eT32);
      else if (dWorst2 > dWorst3 + dEps) return(eFC == eT23);
      else                               return(false);
    }
  case eT22:
    {
      double adNormACE[3], adNormECB[3], adNormBCD[3], adNormDCA[3];
      double adNormABC[3], adNormCDE[3], adNormABE[3];
      vUnitNormal(pVVertA->adCoords(), pVVertC->adCoords(),
		  pVVertE->adCoords(), adNormACE);
      vUnitNormal(pVVertE->adCoords(), pVVertC->adCoords(),
		  pVVertB->adCoords(), adNormECB);
      vUnitNormal(pVVertB->adCoords(), pVVertC->adCoords(),
		  pVVertD->adCoords(), adNormBCD);
      vUnitNormal(pVVertD->adCoords(), pVVertC->adCoords(),
		  pVVertA->adCoords(), adNormDCA);

      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertC->adCoords(), adNormABC);
      vUnitNormal(pVVertC->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormCDE);
      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertE->adCoords(), adNormABE);

      double dLastDot, dWorstNow = -1., dWorstSwapped = -1.;

      dLastDot = dDOT3D(adNormACE, adNormECB);
      dWorstNow = max(dWorstNow, dLastDot);
      dLastDot = dDOT3D(adNormBCD, adNormDCA);
      dWorstNow = max(dWorstNow, dLastDot);
      dLastDot = fabs(dDOT3D(adNormABC, adNormABE));
      dWorstNow = max(dWorstNow, dLastDot);

      dLastDot = dDOT3D(adNormACE, adNormDCA);
      dWorstSwapped = max(dWorstSwapped, dLastDot);
      dLastDot = dDOT3D(adNormBCD, adNormECB);
      dWorstSwapped = max(dWorstSwapped, dLastDot);
      dLastDot = fabs(dDOT3D(adNormCDE, adNormABE));
      dWorstSwapped = max(dWorstSwapped, dLastDot);

      if   (dWorstNow > dWorstSwapped + dEps) return(true);
      else                                    return(false);
    }
  case eT44:
    return(false);
  default:
    assert2(0, "Bad face type in finding max dihedral angle");
    return(false);
  }
}

//@@ Determines whether swapping will improve the maximum face angle.
// Returns true to swap, false to remain the same.
static bool
qDoSwapMinSine(const Vert* pVVertA, const Vert* pVVertB,
	       const Vert* pVVertC, const Vert* pVVertD,
	       const Vert* pVVertE)
{
  // This algorithm finds sine of face angles between adjacent faces by
  // crossing their unit normals.  The idea is to pick the configuration
  // which maximizes this value.
  //
  // To prevent pairs from flopping back and forth, a swap is only
  // authorized if the inequality is bigger than eps.

  double dEps = 1.e-10;
  switch (eFC) {
  case eT23:
  case eT32:
    {
      double adNormABD[3], adNormABE[3], adNormADE[3], adNormABC[3];
      double adNormBCD[3], adNormBCE[3], adNormBDE[3];
      double adNormCAD[3], adNormCAE[3], adNormCDE[3];
      //@@@ The case where the five pVVerts form a genuinely convex set
      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertD->adCoords(), adNormABD);
      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertE->adCoords(), adNormABE);
      vUnitNormal(pVVertA->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormADE);

      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertC->adCoords(), adNormABC); // Points to E

      vUnitNormal(pVVertB->adCoords(), pVVertC->adCoords(),
		  pVVertD->adCoords(), adNormBCD);
      vUnitNormal(pVVertB->adCoords(), pVVertC->adCoords(),
		  pVVertE->adCoords(), adNormBCE);
      vUnitNormal(pVVertB->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormBDE);

      vUnitNormal(pVVertC->adCoords(), pVVertA->adCoords(),
		  pVVertD->adCoords(), adNormCAD);
      vUnitNormal(pVVertC->adCoords(), pVVertA->adCoords(),
		  pVVertE->adCoords(), adNormCAE);
      vUnitNormal(pVVertC->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormCDE);

      double dWorst2 = 1., dWorst3 = 1.;

      // Because I'm taking the MAGNITUDE of the cross product, the
      // actual sign of the result doesn't matter (vector can be
      // reversed).  This means that normal orientation doesn't matter.

      // Sines of dihedral angles of tet ABCD
      double adTemp[3], dMag;

      vCROSS3D(adNormABD, adNormBCD, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormBCD, adNormCAD, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormCAD, adNormABD, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormABD, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormBCD, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormCAD, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      // Sines of dihedral angles of tet ABCE
      vCROSS3D(adNormABE, adNormBCE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormBCE, adNormCAE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormCAE, adNormABE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormABE, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormBCE, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      vCROSS3D(adNormCAE, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst2 = min(dWorst2, dMag);

      // Now have the smallest sine of dihedral angle in tet ABCD and ABCE

      //Sines of dihedral angles of tet ABDE
      vCROSS3D(adNormABD, adNormABE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormABD, adNormADE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormABD, adNormBDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormABE, adNormADE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormABE, adNormBDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormADE, adNormBDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      //Sines of dihedral angles of tet BCDE
      vCROSS3D(adNormBCD, adNormBCE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormBCD, adNormBDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormBCD, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormBCE, adNormBDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormBCE, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormBDE, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      //Sines of dihedral angles of tet CADE
      vCROSS3D(adNormCAD, adNormCAE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormCAD, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormCAD, adNormADE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormCAE, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormCAE, adNormADE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);

      vCROSS3D(adNormCDE, adNormADE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorst3 = min(dWorst3, dMag);
      // Now have the smallest sine of dihedral in tets ABDE, BCDE, CADE

      if      (dWorst3 > dWorst2 + dEps) return(eFC == eT23);
      else if (dWorst2 > dWorst3 + dEps) return(eFC == eT32);
      else                               return(false);
    }
  case eT22:
    {
      double adNormACE[3], adNormECB[3], adNormBCD[3], adNormDCA[3];
      double adNormABC[3], adNormCDE[3], adNormABE[3];
      vUnitNormal(pVVertA->adCoords(), pVVertC->adCoords(),
		  pVVertE->adCoords(), adNormACE);
      vUnitNormal(pVVertE->adCoords(), pVVertC->adCoords(),
		  pVVertB->adCoords(), adNormECB);
      vUnitNormal(pVVertB->adCoords(), pVVertC->adCoords(),
		  pVVertD->adCoords(), adNormBCD);
      vUnitNormal(pVVertD->adCoords(), pVVertC->adCoords(),
		  pVVertA->adCoords(), adNormDCA);

      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertC->adCoords(), adNormABC);
      vUnitNormal(pVVertC->adCoords(), pVVertD->adCoords(),
		  pVVertE->adCoords(), adNormCDE);
      vUnitNormal(pVVertA->adCoords(), pVVertB->adCoords(),
		  pVVertE->adCoords(), adNormABE);

      double dWorstNow = 1., dWorstSwap = 1.;
      double adTemp[3], dMag;

      // Sines of dihedral angles of tet ABCD, except those that are the
      // same both ways
      vCROSS3D(adNormBCD, adNormDCA, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);

      vCROSS3D(adNormBCD, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);

      vCROSS3D(adNormDCA, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);


      // Sines of dihedral angles of tet ABCE, except those that are the
      // same both ways
      vCROSS3D(adNormECB, adNormACE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);

      vCROSS3D(adNormECB, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);

      vCROSS3D(adNormACE, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);


      // Only check the center face once; the sine is the same on both
      // sides
      vCROSS3D(adNormABE, adNormABC, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstNow = min(dWorstNow, dMag);

      // Now have the smallest sine of dihedral angle in tet ABCD and ABCE

      // Sines of dihedral angles of tet DECA, except those that are the
      // same both ways
      vCROSS3D(adNormDCA, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);

      vCROSS3D(adNormACE, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);

      vCROSS3D(adNormDCA, adNormACE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);


      // Sines of dihedral angles of tet DECB, except those that are the
      // same both ways
      vCROSS3D(adNormBCD, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);

      vCROSS3D(adNormECB, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);

      vCROSS3D(adNormBCD, adNormECB, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);


      // Only check the center face once; the sine is the same on both
      // sides
      vCROSS3D(adNormABE, adNormCDE, adTemp);
      dMag = dMAG3D(adTemp);
      dWorstSwap = min(dWorstSwap, dMag);

      // Now have the smallest sine of dihedral angle in tet DECA and DECB

      if   (dWorstNow + dEps < dWorstSwap) return(true);
      else                                 return(false);
    }
  case eT44:
    return(false);
  default:
    assert2(0, "Bad face type in finding max dihedral angle");
    return(false);
  }
}

bool VolMesh::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, pVVertE);
  case eMaxDihed:
    return qDoSwapMaxDihed(pVVertA, pVVertB, pVVertC, pVVertD, pVVertE);
  case eMinSine:
    return qDoSwapMinSine(pVVertA, pVVertB, pVVertC, pVVertD, pVVertE);
  case eFromQual:
    {
      if (eFC != eT23 && eFC != eT32) return false;
      assert(pQ != NULL);
      TetCellCV TCCV11(const_cast<Vert*>(pVVertA),
		       const_cast<Vert*>(pVVertB),
		       const_cast<Vert*>(pVVertC),
		       const_cast<Vert*>(pVVertD));
      TetCellCV TCCV12(const_cast<Vert*>(pVVertA),
		       const_cast<Vert*>(pVVertC),
		       const_cast<Vert*>(pVVertB),
		       const_cast<Vert*>(pVVertE));
      double dQual11 = pQ->dEvaluate(&TCCV11);
      double dQual12 = pQ->dEvaluate(&TCCV12);

      TetCellCV TCCV21(const_cast<Vert*>(pVVertD),
		       const_cast<Vert*>(pVVertE),
		       const_cast<Vert*>(pVVertB),
		       const_cast<Vert*>(pVVertA));
      TetCellCV TCCV22(const_cast<Vert*>(pVVertD),
		       const_cast<Vert*>(pVVertE),
		       const_cast<Vert*>(pVVertC),
		       const_cast<Vert*>(pVVertB));
      TetCellCV TCCV23(const_cast<Vert*>(pVVertD),
		       const_cast<Vert*>(pVVertE),
		       const_cast<Vert*>(pVVertA),
		       const_cast<Vert*>(pVVertC));
      double dQual21 = pQ->dEvaluate(&TCCV21);
      double dQual22 = pQ->dEvaluate(&TCCV22);
      double dQual23 = pQ->dEvaluate(&TCCV23);

      bool qMaximize = pQ->qMaximizeMinValue();
      static const double dEps = 1.e-10;
      if (qMaximize) {
	double dQual1 = min(dQual11, dQual12);
	double dQual2 = MIN3(dQual21, dQual22, dQual23);
	if      (dQual1 > dQual2 + dEps) return (eFC == eT32);
	else if (dQual2 > dQual1 + dEps) return (eFC == eT23);
	else                             return false;
      }
      else {
	double dQual1 = max(dQual11, dQual12);
	double dQual2 = MAX3(dQual21, dQual22, dQual23);
	if      (dQual1 < dQual2 - dEps) return (eFC == eT32);
	else if (dQual2 < dQual1 - dEps) return (eFC == eT23);
	else                             return false;
      }
    }
  default:
    assert2(0, "Bad swap measure specified");
    return(false);
  }
}

int VolMesh::iReconfigure(TetCell* apTCTets[], TriFace* const pF,
			  Vert* const pVVertA, Vert* const pVVertB,
			  Vert* const pVVertC, Vert* const pVVertD,
			  Vert* const pVVertE)
{
  assert(pF->iFullCheck());
  assert(pVVertA->qValid());
  assert(pVVertB->qValid());
  assert(pVVertC->qValid());
  assert(pVVertD->qValid());
  assert(pVVertE->qValid());

  switch (eFC) {
  case eT22:
    return (iReconfTet22(apTCTets, pF,
			 pVVertA, pVVertB, pVVertC, pVVertD, pVVertE));
  case eT23:
    return (iReconfTet23(apTCTets, pF,
			 pVVertA, pVVertB, pVVertC, pVVertD, pVVertE));
  case eT32:
    return (iReconfTet32(apTCTets, pF,
			 pVVertA, pVVertB, pVVertC, pVVertD, pVVertE));
  case eT44:
    return (iReconfTet44(apTCTets, pF,
			 pVVertA, pVVertB, pVVertC, pVVertD, pVVertE));
  default:
    assert2(0, "Bad face type");
    return(0);
  }
}

//@@ Swap two tets for three.
int
VolMesh::iReconfTet23(TetCell* apTCTets[], TriFace* const pF,
		      Vert* const pVVertA, Vert* const pVVertB,
		      Vert* const pVVertC, Vert* const pVVertD,
		      Vert* const pVVertE)
{
  // Identify both cells and all seven faces in the two-tet config.
  TetCell* pTCTetD = apTCTets[0];
  TetCell* pTCTetE = apTCTets[1];
  if ( pTCTetD->iRegion() != pTCTetE->iRegion() ) return 0;
  assert(pTCTetD->qHasVert(pVVertD));
  assert(pTCTetE->qHasVert(pVVertE));
  Face* pFFaceD = pF;

  // Top half
  Face* pFFaceA = pTCTetE->pFFaceOpposite(pVVertA);
  Face* pFFaceB = pTCTetE->pFFaceOpposite(pVVertB);
  Face* pFFaceC = pTCTetE->pFFaceOpposite(pVVertC);

  // Bottom half
  Face* pFFaceA1 = pTCTetD->pFFaceOpposite(pVVertA);
  Face* pFFaceB1 = pTCTetD->pFFaceOpposite(pVVertB);
  Face* pFFaceC1 = pTCTetD->pFFaceOpposite(pVVertC);

/* Now switch configurations.  This involves adding one cell and two
 * faces, net.  First assign indices for this new stuff, preparing to
 * write over pTCTetD, pTCTetE, and faceD, which are all going away.
 */
  Cell* pTCTetA = pTCTetD;
  Cell* pTCTetB = pTCTetE;
  Cell* pTCTetC = pCNewCell();
  pTCTetC->vSetRegion(pTCTetA->iRegion());

  Face* pFFaceA2 = pFNewFace();
  Face* pFFaceB2 = pFNewFace();
  Face* pFFaceC2 = pFFaceD;
  // Note that code in the Delaunay vertex removal routine in
  // Remove3D.cxx requires that pFFaceC2 is the one that ends up in the
  // place of the original pF.  (As of Feb, 2004).

  // All nine faces get touched; the internal ones are known to be okay
  // at this point, so they get their swap flags switched back. 

  /* Fix up the cell neighbors of the six faces on the hull. */
  pFFaceA->vReplaceCell(pTCTetE, pTCTetA);
  pFFaceB->vReplaceCell(pTCTetE, pTCTetB);
  pFFaceC->vReplaceCell(pTCTetE, pTCTetC);

  pFFaceA1->vReplaceCell(pTCTetD, pTCTetA);
  pFFaceB1->vReplaceCell(pTCTetD, pTCTetB);
  pFFaceC1->vReplaceCell(pTCTetD, pTCTetC);

  /* Set up the face_vert and face_cell lists for the three new faces. */
  pFFaceA2->vAssign(pTCTetB, pTCTetC, pVVertA, pVVertD, pVVertE);
  pFFaceB2->vAssign(pTCTetC, pTCTetA, pVVertB, pVVertD, pVVertE);
  pFFaceC2->vAssign(pTCTetA, pTCTetB, pVVertC, pVVertD, pVVertE);

  pFFaceA2->vClearSwapFlag();
  pFFaceB2->vClearSwapFlag();
  pFFaceC2->vClearSwapFlag();

  /* Reset vertex hints for the "equator" verts */
  pVVertA->vSetHintFace(pFFaceA2);
  pVVertB->vSetHintFace(pFFaceB2);
  pVVertC->vSetHintFace(pFFaceC2);

  // Mark all verts to be smoothed
  pVVertA->vMarkIllShaped();
  pVVertB->vMarkIllShaped();
  pVVertC->vMarkIllShaped();
  pVVertD->vMarkIllShaped();
  pVVertE->vMarkIllShaped();
  
  /* Set up the cell_vert and cell_face lists for the three new cells. */
  pTCTetA->vAssign(pFFaceA, pFFaceA1, pFFaceB2, pFFaceC2);
  pTCTetB->vAssign(pFFaceB, pFFaceB1, pFFaceC2, pFFaceA2);
  pTCTetC->vAssign(pFFaceC, pFFaceC1, pFFaceA2, pFFaceB2);

  assert(!pTCTetA->qHasVert(pVVertA));
  assert(!pTCTetB->qHasVert(pVVertB));
  assert(!pTCTetC->qHasVert(pVVertC));
  assert(iOrient3D(pVVertA, pVVertB, pVVertD, pVVertE) == 1);
  assert(iOrient3D(pVVertB, pVVertC, pVVertD, pVVertE) == 1);
  assert(iOrient3D(pVVertC, pVVertA, pVVertD, pVVertE) == 1);

  assert(pTCTetA->iFullCheck());
  assert(pTCTetB->iFullCheck());
  assert(pTCTetC->iFullCheck());

  assert(pFFaceA ->iFullCheck());
  assert(pFFaceB ->iFullCheck());
  assert(pFFaceC ->iFullCheck());
  assert(pFFaceA1->iFullCheck());
  assert(pFFaceB1->iFullCheck());
  assert(pFFaceC1->iFullCheck());
  assert(pFFaceA2->iFullCheck());
  assert(pFFaceB2->iFullCheck());
  assert(pFFaceC2->iFullCheck());

  int iRetVal = 1;
  if (qSwapRecur && iSwapDepth <= iMaxSwapDepth) {
    iRetVal += iFaceSwap(pFFaceA );
    iRetVal += iFaceSwap(pFFaceB );
    iRetVal += iFaceSwap(pFFaceC );
    iRetVal += iFaceSwap(pFFaceA1);
    iRetVal += iFaceSwap(pFFaceB1);
    iRetVal += iFaceSwap(pFFaceC1);
  }
  else {
    pFFaceA ->vMarkForSwap();
    pFFaceB ->vMarkForSwap();
    pFFaceC ->vMarkForSwap();
    pFFaceA1->vMarkForSwap();
    pFFaceB1->vMarkForSwap();
    pFFaceC1->vMarkForSwap();
  }
  return(iRetVal);
}
//@@ Swap three tets for two.
int
VolMesh::iReconfTet32(TetCell* apTCTets[], TriFace* const pF,
		      Vert* const pVVertA, Vert* const pVVertB,
		      Vert* const pVVertC, Vert* const pVVertD,
		      Vert* const pVVertE)
{
  // Identify the three cells of the three-tet config.  The cell ID's
  // require that the cell NOT containing pVVertA is pTCTetA, etc.
  TetCell* pTCTetA = apTCTets[0];
  TetCell* pTCTetB = apTCTets[1];
  TetCell* pTCTetC = apTCTets[2];
  assert(! pTCTetA->qHasVert(pVVertA));
  assert(! pTCTetB->qHasVert(pVVertB));
  assert(! pTCTetC->qHasVert(pVVertC));
  if ( (pTCTetA->iRegion() != pTCTetB->iRegion()) ||
       (pTCTetA->iRegion() != pTCTetC->iRegion()) ) return 0;

  assert(pF->qHasVert(pVVertA) || pF->qHasVert(pVVertB)
	 || pF->qHasVert(pVVertC));
  assert(pF->qHasVert(pVVertD));
  assert(pF->qHasVert(pVVertE));

  // Identify the nine faces in the three-tet configuration. pFFaceA and
  // pFFaceA1 are faces in the shell not incident on pVVertA. pFFaceA2 IS
  // incident on pVVertA, because that's the easiest way to define it
  // uniquely.
  Face* pFFaceA  = pTCTetA->pFFaceOpposite( pVVertE );
  Face* pFFaceB  = pTCTetB->pFFaceOpposite( pVVertE );
  Face* pFFaceC  = pTCTetC->pFFaceOpposite( pVVertE );

  Face* pFFaceA1 = pTCTetA->pFFaceOpposite( pVVertD );
  Face* pFFaceB1 = pTCTetB->pFFaceOpposite( pVVertD );
  Face* pFFaceC1 = pTCTetC->pFFaceOpposite( pVVertD );

  Face* pFFaceA2 = pTCTetC->pFFaceOpposite( pVVertB );
  Face* pFFaceB2 = pTCTetA->pFFaceOpposite( pVVertC );
  Face* pFFaceC2 = pTCTetB->pFFaceOpposite( pVVertA );

  // Set the indices for the interior face and the two cells that will
  // result from all this.

  // Code in Delaunay vertex removal in Remove3D.cxx requires that pF be
  // the face that survives.
  Face* pFFaceD = pF;
  Cell* pTCTetD = pTCTetA;
  Cell* pTCTetE = pTCTetB;

  // Fix up the cell neighbors of the six faces on the hull.
  pFFaceA->vReplaceCell(pTCTetA, pTCTetD);
  pFFaceB->vReplaceCell(pTCTetB, pTCTetD);
  pFFaceC->vReplaceCell(pTCTetC, pTCTetD);

  pFFaceA1->vReplaceCell(pTCTetA, pTCTetE);
  pFFaceB1->vReplaceCell(pTCTetB, pTCTetE);
  pFFaceC1->vReplaceCell(pTCTetC, pTCTetE);

  // Set stuff up for the interior face.
  pFFaceD->vAssign(pTCTetD, pTCTetE, pVVertA, pVVertB, pVVertC);

  // All seven faces have been touched; clear the swap request for the
  // interior face.
  pFFaceD->vClearSwapFlag();

  // Set up the two new cells.
  pTCTetD->vAssign(pFFaceA , pFFaceB , pFFaceC , pFFaceD);
  pTCTetE->vAssign(pFFaceA1, pFFaceB1, pFFaceC1, pFFaceD);

  // Reset all the vertex hints
  pVVertA->vSetHintFace(pFFaceD);
  pVVertB->vSetHintFace(pFFaceD);
  pVVertC->vSetHintFace(pFFaceD);
  pVVertD->vSetHintFace(pFFaceA);
  pVVertE->vSetHintFace(pFFaceA1);

  // Mark all verts to be smoothed
  pVVertA->vMarkIllShaped();
  pVVertB->vMarkIllShaped();
  pVVertC->vMarkIllShaped();
  pVVertD->vMarkIllShaped();
  pVVertE->vMarkIllShaped();

  // Tag the extra faces and the extra cell for later deletion.  It's
  // easier to mark all three interior faces for deletion and then
  // retrieve the one you're keeping than to do the logic.  
  pFFaceA2->vMarkDeleted();
  pFFaceB2->vMarkDeleted();
  pFFaceC2->vMarkDeleted();
  pF->vMarkNotDeleted();
  pTCTetC->vMarkDeleted();

  assert(pTCTetD->qHasVert(pVVertD));
  assert(pTCTetE->qHasVert(pVVertE));
  assert(iOrient3D(pVVertA, pVVertC, pVVertB, pVVertD) == 1);
  assert(iOrient3D(pVVertA, pVVertB, pVVertC, pVVertE) == 1);

  assert(pTCTetD->iFullCheck());
  assert(pTCTetE->iFullCheck());

  assert(pFFaceA ->iFullCheck());
  assert(pFFaceB ->iFullCheck());
  assert(pFFaceC ->iFullCheck());
  assert(pFFaceD ->iFullCheck());
  assert(pFFaceA1->iFullCheck());
  assert(pFFaceB1->iFullCheck());
  assert(pFFaceC1->iFullCheck());

  int iRetVal = 1;
  if (qSwapRecur && iSwapDepth <= iMaxSwapDepth) {
    iRetVal += iFaceSwap(pFFaceA );
    iRetVal += iFaceSwap(pFFaceB );
    iRetVal += iFaceSwap(pFFaceC );
    iRetVal += iFaceSwap(pFFaceA1);
    iRetVal += iFaceSwap(pFFaceB1);
    iRetVal += iFaceSwap(pFFaceC1);
  }
  else {
    pFFaceA ->vMarkForSwap();
    pFFaceB ->vMarkForSwap();
    pFFaceC ->vMarkForSwap();
    pFFaceA1->vMarkForSwap();
    pFFaceB1->vMarkForSwap();
    pFFaceC1->vMarkForSwap();
  }
  return(iRetVal);
}

//@@ Swap two tets for two, in the case where two faces are coplanar.
int
VolMesh::iReconfTet22(TetCell* apTCTets[], TriFace* const pF,
		      Vert* const pVVertA, Vert* const pVVertB,
		      Vert* const pVVertC, Vert* const pVVertD,
		      Vert* const pVVertE)
{
  // This routine does only the 2-to-2 case.  4-to-4 is done elsewhere.
  // The cells and verts are assumed to have been set up properly already.

/*      	        C
        	       /\
        	      /||\
        	     / |  \                 A,D are on the base;
       	       	    /  | | \                B,E are front, C,F are back.
  		   /   |    \               Cell A is on the left, B on
        	  /    |  |  \	            the right.
        	 /     |      \
  		/      |   |   \
  	       /     C |     F  \
  	      /	       |  E |    \
  	     /   B     | _ - B    \
  	    /  _ _ - - |    /  - _ \
  	   D -         |           -E
  	    -__	    A  |  /  D  __--
  	       -__     |     __-
  	     	  -___ |/__--
  	     	      -A-

*/
  if (!qSurfaceEdgeChangesAllowed() || !qBdryChangesAllowed()) return 0;
  TetCell* pTCTetA = apTCTets[0];
  TetCell* pTCTetB = apTCTets[1];
  assert(pTCTetA->qHasVert(pVVertD));
  assert(pTCTetB->qHasVert(pVVertE));
  assert(pTCTetA->iRegion() == pTCTetB->iRegion());

  Face* pFFaceA = pTCTetA->pFFaceOpposite(pVVertC);
  Face* pFFaceB = pTCTetA->pFFaceOpposite(pVVertB);
  Face* pFFaceC = pTCTetA->pFFaceOpposite(pVVertA);

  Face* pFFaceD = pTCTetB->pFFaceOpposite(pVVertC);
  Face* pFFaceE = pTCTetB->pFFaceOpposite(pVVertB);
  Face* pFFaceF = pTCTetB->pFFaceOpposite(pVVertA);

  // Check that faces A and D are bdry faces.  Then, if they belong to
  // different patches, bail out, so that we won't end up losing a
  // subsegment (from A to B, in this case) when we swap.
  {
    Cell* pCOppA = pFFaceA->pCCellOpposite(pTCTetA);
    Cell* pCOppD = pFFaceD->pCCellOpposite(pTCTetB);
    assert(pCOppA->eType() == Cell::eTriBFace);
    assert(pCOppD->eType() == Cell::eTriBFace);
    TriBFace* pTBFA = dynamic_cast<TriBFace*>(pCOppA);
    TriBFace* pTBFD = dynamic_cast<TriBFace*>(pCOppD);
    if (!qBdryFacesOKToSwap(pTBFA, pTBFD)) return 0;
  }
  
  // Do the actual two-for-two swap here.  The faces all stay in the
  // same places except for pFFace, the dividing face, which gets turned.
  // Cell A ends up in front, cell B in back.  The only vert hints
  // affected are for verts 1 and 2, which can be set explicitly.

  // Fix up the two faces which are changing cell neighbors:
  pFFaceC->vReplaceCell(pTCTetA, pTCTetB);
  pFFaceE->vReplaceCell(pTCTetB, pTCTetA);
  // Mark the other two external faces as needing a swap check
  pFFaceB->vMarkForSwap();
  pFFaceF->vMarkForSwap();

  // Fix up the co-planar faces to indicate their new vertex
  // affiliations
  pFFaceA->vReplaceVert(pVVertB, pVVertE);
  pFFaceD->vReplaceVert(pVVertA, pVVertD);

  // Update the cell-face table for cells A and B.
  pTCTetA->vReplaceFace(pFFaceC, pFFaceE);
  pTCTetB->vReplaceFace(pFFaceE, pFFaceC);

  // Reset two vertex hints
  pVVertA->vSetHintFace(pFFaceB);
  pVVertB->vSetHintFace(pFFaceC);
  
  // Mark all verts to be smoothed
  pVVertA->vMarkIllShaped();
  pVVertB->vMarkIllShaped();
  pVVertC->vMarkIllShaped();
  pVVertD->vMarkIllShaped();
  pVVertE->vMarkIllShaped();

  // Finally, reset the dividing face
  pF->vAssign(pTCTetB, pTCTetA, pVVertE, pVVertC, pVVertD);
  //  pF->vClearSwapFlag();

  assert(pTCTetA->iFullCheck());
  assert(pTCTetB->iFullCheck());

  assert(pF ->iFullCheck());
  assert(pFFaceA->iFullCheck());
  assert(pFFaceB->iFullCheck());
  assert(pFFaceC->iFullCheck());
  assert(pFFaceD->iFullCheck());
  assert(pFFaceE->iFullCheck());
  assert(pFFaceF->iFullCheck());

  assert(pTCTetB->pVVertOpposite(pF) == pVVertB);
  assert(iOrient3D(pVVertC, pVVertD, pVVertE, pVVertB) == -1);
  assert(pTCTetA->pVVertOpposite(pF) == pVVertA);
  assert(iOrient3D(pVVertC, pVVertD, pVVertE, pVVertA) ==  1);

  int iRetVal = 1;
  if (qSwapRecur) {
    iRetVal += iFaceSwap(pFFaceB);
    iRetVal += iFaceSwap(pFFaceC);
    iRetVal += iFaceSwap(pFFaceE);
    iRetVal += iFaceSwap(pFFaceF);
  }
  else {
    pFFaceB ->vMarkForSwap();
    pFFaceC ->vMarkForSwap();
    pFFaceE ->vMarkForSwap();
    pFFaceF ->vMarkForSwap();
  }
  return(iRetVal);
}

//@@ Swap four tets for four, in the case where two faces are coplanar.
int
VolMesh::iReconfTet44(TetCell* apTCTets[], TriFace* const pF,
		      Vert* const pVVertA, Vert* const pVVertB,
		      Vert* const pVVertC, Vert* const pVVertD,
		      Vert* const pVVertE)
{
  // The cells and verts are assumed to have been set up properly already.

  // Here's the before picture:
/*      	        C
        	       /\
        	      /||\
        	     / |  \                 A,D are on the base;
       	       	    /  | | \                B,E are front, C,F are back.
  		   /   |    \               Cell A is on the left, B on
        	  /    |  |  \	            the right.
        	 /     |      \
  		/      |   |   \
  	       /     C |     F  \
  	      /	       |  E |    \
  	     /   B     | _ - B    \
  	    /  _ _ - - |    /  - _ \
  	   D -         |           -E
  	    -__	    A  |  /  D  __--
  	       -__     |     __-
  	     	  -___ |/__--
  	     	      -A-

  	               ___---B_
  	       ____----	    /  --__
  	   D---            /       -E
  	    - _	    A     / |D  _ -
  	    \  - _       /   _ -   /
             \     -_   / _-      /
  	      \       -A-  |     /
  	       \     C1|     F1 /
  		\      |  |    /
        	 \ B1  |  E1  /
        	  \    |     /	            A,D are on the base;
  		   \   | |  /               B1,E1 are front, C1,F1 are back.
       	       	    \  |   /                Cell A1 is on the left, B1 on
        	     \ || /                 the right.
        	      \| /                  Face G divides these two cells.
        	       \/
                       C1

*/
  TetCell* pTCTetA = apTCTets[0];
  TetCell* pTCTetB = apTCTets[1];
  assert(pTCTetA->qHasVert(pVVertD));
  assert(pTCTetB->qHasVert(pVVertE));

  Face* pFFaceA = pTCTetA->pFFaceOpposite(pVVertC);
  Face* pFFaceB = pTCTetA->pFFaceOpposite(pVVertB);
  Face* pFFaceC = pTCTetA->pFFaceOpposite(pVVertA);

  Face* pFFaceD = pTCTetB->pFFaceOpposite(pVVertC);
  Face* pFFaceE = pTCTetB->pFFaceOpposite(pVVertB);
  Face* pFFaceF = pTCTetB->pFFaceOpposite(pVVertA);

  // Now identify the cells, faces and vert on the opposite side of the plane
  TetCell* pTCTetA1 = dynamic_cast<TetCell*>(pFFaceA->pCCellOpposite(pTCTetA));
  TetCell* pTCTetB1 = dynamic_cast<TetCell*>(pFFaceD->pCCellOpposite(pTCTetB));

  // Check to make sure that there isn't something fishy going on with
  // region boundaries.
  if ( (pTCTetA->iRegion() != pTCTetB->iRegion()) ||
       (pTCTetA1->iRegion() != pTCTetB1->iRegion()) ) return 0;
  
  assert(pTCTetA1->qHasVert(pVVertD));
  assert(pTCTetB1->qHasVert(pVVertE));

  Face* pFFaceB1 = pTCTetA1->pFFaceOpposite(pVVertB);
  Face* pFFaceC1 = pTCTetA1->pFFaceOpposite(pVVertA);

  Face* pFFaceE1 = pTCTetB1->pFFaceOpposite(pVVertB);
  Face* pFFaceF1 = pTCTetB1->pFFaceOpposite(pVVertA);
  Face* pFFaceG  = pTCTetB1->pFFaceOpposite(pVVertE);

  Vert* pVVertC1 = pTCTetB1->pVVertOpposite(pFFaceD);

  // Do the actual four-for-four swap here.  The faces all stay in the
  // same places except for pF, the dividing face, which gets turned.
  // Cell A ends up in front, cell B in back.  The only vert hints
  // affected are for verts A and B, which are set explicitly.
  // Below the plane, pFFaceG gets moves, cell A1 ends up in front, cell
  // B1 in back. 

  // Fix up the two faces which are changing cell neighbors:
  pFFaceC->vReplaceCell(pTCTetA, pTCTetB);
  pFFaceE->vReplaceCell(pTCTetB, pTCTetA);
  pFFaceC1->vReplaceCell(pTCTetA1, pTCTetB1);
  pFFaceE1->vReplaceCell(pTCTetB1, pTCTetA1);
  pFFaceB->vMarkForSwap();
  pFFaceF->vMarkForSwap();
  pFFaceB1->vMarkForSwap();
  pFFaceF1->vMarkForSwap();

  // Fix up the co-planar faces to indicate their new vertex
  // affiliations
  pFFaceA->vReplaceVert(pVVertB, pVVertE);
  pFFaceD->vReplaceVert(pVVertA, pVVertD);

  // Update the cell-face table for cells A and B.
  pTCTetA->vReplaceFace(pFFaceC, pFFaceE);
  pTCTetB->vReplaceFace(pFFaceE, pFFaceC);
  pTCTetA1->vReplaceFace(pFFaceC1, pFFaceE1);
  pTCTetB1->vReplaceFace(pFFaceE1, pFFaceC1);

  // Reset two vertex hints
  pVVertA->vSetHintFace(pFFaceB);
  pVVertB->vSetHintFace(pFFaceC);
  
  // Mark all verts to be smoothed
  pVVertA->vMarkIllShaped();
  pVVertB->vMarkIllShaped();
  pVVertC->vMarkIllShaped();
  pVVertC1->vMarkIllShaped();
  pVVertD->vMarkIllShaped();
  pVVertE->vMarkIllShaped();

  // Finally, reset the dividing faces
  pF->vAssign(pTCTetB, pTCTetA, pVVertE, pVVertC, pVVertD);
  //  pF->vClearSwapFlag();
  pFFaceG->vAssign(pTCTetB1, pTCTetA1, pVVertD, pVVertC1, pVVertE);
  //  pFFaceG->vClearSwapFlag();

  assert(pTCTetA->iFullCheck());
  assert(pTCTetB->iFullCheck());
  assert(pTCTetA1->iFullCheck());
  assert(pTCTetB1->iFullCheck());

  assert(pF ->iFullCheck());
  assert(pFFaceA->iFullCheck());
  assert(pFFaceB->iFullCheck());
  assert(pFFaceC->iFullCheck());
  assert(pFFaceD->iFullCheck());
  assert(pFFaceE->iFullCheck());
  assert(pFFaceF->iFullCheck());

  assert(pFFaceB1->iFullCheck());
  assert(pFFaceC1->iFullCheck());
  assert(pFFaceE1->iFullCheck());
  assert(pFFaceF1->iFullCheck());
  assert(pFFaceG->iFullCheck());

  assert(pTCTetB->pVVertOpposite(pF) == pVVertB);
  assert(iOrient3D(pVVertC, pVVertD, pVVertE, pVVertB) == -1);
  assert(pTCTetA->pVVertOpposite(pF) == pVVertA);
  assert(iOrient3D(pVVertC, pVVertD, pVVertE, pVVertA) ==  1);

  assert(pTCTetB1->pVVertOpposite(pFFaceG) == pVVertB);
  assert(iOrient3D(pVVertC1, pVVertD, pVVertE, pVVertB) ==  1);
  assert(pTCTetA1->pVVertOpposite(pFFaceG) == pVVertA);
  assert(iOrient3D(pVVertC1, pVVertD, pVVertE, pVVertA) == -1);

  int iRetVal = 1;
  if (qSwapRecur) {
    iRetVal += iFaceSwap(pFFaceB);
    iRetVal += iFaceSwap(pFFaceC);
    iRetVal += iFaceSwap(pFFaceE);
    iRetVal += iFaceSwap(pFFaceF);
    iRetVal += iFaceSwap(pFFaceB1);
    iRetVal += iFaceSwap(pFFaceC1);
    iRetVal += iFaceSwap(pFFaceE1);
    iRetVal += iFaceSwap(pFFaceF1);
  }
  else {
    pFFaceB ->vMarkForSwap();
    pFFaceC ->vMarkForSwap();
    pFFaceE ->vMarkForSwap();
    pFFaceF ->vMarkForSwap();
    pFFaceB1->vMarkForSwap();
    pFFaceC1->vMarkForSwap();
    pFFaceE1->vMarkForSwap();
    pFFaceF1->vMarkForSwap();
  }
  return(iRetVal);
}

