#include "GR_assert.h"
#include <math.h>
#include <stdlib.h>
#include "GR_misc.h"
#include "GR_Geometry.h"
#include "GR_VolMesh.h"
#include "GR_events.h"
#include "SMsmooth.h"
#include "GR_BFace.h"
#include "GR_Vec.h"

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

static int iOnEdge = 0;
static int iOnBdryEdge = 0;
static int iOnBdryEdge2 = 0;
static int iOnIntBdryEdge = 0;
static int iOnFace = 0;
static int iOnIntBdryFace = 0;
static int iOnBdryFace = 0;
static int iInterior = 0;

static bool qInsertInThisCell(Cell* const pC,
			     double dfLengthScale(const double adCoord[3]))
{
  assert(pC->iFullCheck());
  assert(pC->eType() == Cell::eTet);
  double dRad = dynamic_cast<TetCell*>(pC)->dCircumradius();
//    const double dFactor = M_SQRT2 / 12.;

//    double adCentroid[3];
//    adCentroid[0] = 0.25 * (pC->pVVert(0)->dX() + pC->pVVert(1)->dX() +
//			  pC->pVVert(2)->dX() + pC->pVVert(3)->dX());
//    adCentroid[1] = 0.25 * (pC->pVVert(0)->dY() + pC->pVVert(1)->dY() +
//			  pC->pVVert(2)->dY() + pC->pVVert(3)->dY());
//    adCentroid[2] = 0.25 * (pC->pVVert(0)->dZ() + pC->pVVert(1)->dZ() +
//			  pC->pVVert(2)->dZ() + pC->pVVert(3)->dZ());
  double dLength = 0.25 * (+ dfLengthScale(pC->pVVert(0)->adCoords())
			   + dfLengthScale(pC->pVVert(1)->adCoords())
			   + dfLengthScale(pC->pVVert(2)->adCoords())
			   + dfLengthScale(pC->pVVert(3)->adCoords()));
//    double dCorrectSize = dLength*dLength*dLength * dFactor;

//   double dMinSize = dCorrectSize / 4;
  if (dRad > sqrt(3.)*dLength/2.) return true;
  pC->vMarkWellSized();
  return false;
//   if (dVol < dMinSize) {
//     pC->vMarkWellSized();
//     return false;
//   }

//   double dMaxAng = pC->dMaxDihed();
//   if (dMaxAng > 120.)
//     return true;
//   else {
//     pC->MarkWellSized();
//     return false;
//   }
}

//@ Insert a point on an edge
int VolMesh::iInsertOnEdge(Vert* const pVNew, Vert* const pV0,
			   Vert* const pV1, Face* const pF, const int qSwap)
     // If the edge is actually a boundary edge, iInsertOnBdryEdge will
     // be called when that fact is discovered.
{
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(qSimplicial());
  assert(pF->iFullCheck());
  assert(pF->qHasVert(pV0));
  assert(pF->qHasVert(pV1));

  //@@ Check to see if the face containing this edge is a boundary face
  // If it is, we can cut right the chase, calling iInsertOnBdryEdge
  // immediately.
  if (pF->pCCellLeft ()->eType() == Cell::eTriBFace ||
      pF->pCCellRight()->eType() == Cell::eTriBFace) {
    iOnBdryEdge++;
    return iInsertOnBdryEdge(pVNew, pV0, pV1, pF, qSwap);
  }

  //@@ Likewise, check to see if the given face is an internal bdry face
  bool qIntBdryEdge = false;
  if (pF->iFaceLoc() == Face::eBdryTwoSide)
    qIntBdryEdge = true;
  else {
    assert(pF->pCCellLeft()->iRegion() == pF->pCCellRight()->iRegion());
  }

  //@@ Preamble to walking around the edge
  assert(pF->pCCellLeft()->eType() == Cell::eTet);
  TetCell* pTC = dynamic_cast<TetCell*>(pF->pCCellLeft());

  Vert* pVA = pTC->pVVertOpposite(pF);
  Vert* pVB = pV0;
  Vert* pVC = pV1;
  Vert* pVD;
  if (pF->pVVert(0) != pVB && pF->pVVert(0) != pVC)
    pVD = pF->pVVert(0);
  else if (pF->pVVert(1) != pVB && pF->pVVert(1) != pVC)
    pVD = pF->pVVert(1);
  else {
    assert(pF->pVVert(2) != pVB && pF->pVVert(2) != pVC);
    pVD = pF->pVVert(2);
  }
  if (iOrient3D(pVA, pVB, pVC, pVD) != 1) {
    Vert* pV = pVC; pVC = pVB; pVB = pV;
  }

  //@@ Make lists of cells divided and their faces and verts
  //@@ Create new faces and cell for each original cell
  const int iMaxWrapCells = 25;
  // This is the chain of verts around the perimeter of the shape (Verts
  // B and C are constant).  The extra one at the end will contain a
  // copy of the first one.
  Vert* apVert[iMaxWrapCells+1];
  // Current and new cells
  Cell *apCell[iMaxWrapCells], *apCNew[iMaxWrapCells];
  // apFace contains the original faces separating cells around the edge.
  // Again, a copy is used at the end.
  Face* apFace[iMaxWrapCells+1];
  // These are faces opposite verts B and C, respectively
  Face *apFaceB[iMaxWrapCells], *apFaceC[iMaxWrapCells];
  // New faces formed by dividing apFace members.  Another copy.
  Face* apNewFace[iMaxWrapCells+1];
  // New faces within the original cells
  Face* apIntFace[iMaxWrapCells];

  apVert[0] = pVD;
  apCell[0] = pTC;
  apFace[0] = pF;
  apFaceB[0] = pTC->pFFaceOpposite(pVB);
  apFaceC[0] = pTC->pFFaceOpposite(pVC);
  apCNew[0] = pCNewCell();
  apCNew[0]->vSetRegion(apCell[0]->iRegion());
  apNewFace[0] = pFNewFace();
  apIntFace[0] = pFNewFace();

  int i = 0;
  Cell* pC;
  Face* pFStart = pF;
  //@@@ Walk around the edge, gathering info and setting up new entities
  while (1) {
    i++;
    // These are things that get set up regardless of whether there's
    // another interior cell
    apVert[i] = pTC->pVVertOpposite(apFace[i-1]);
    apFace[i] = pTC->pFFaceOpposite(apVert[i-1]);

    // Continue until returned to starting point...
    if (apFace[i] == pFStart) break;
    if (apFace[i]->iFaceLoc() == Face::eBdryTwoSide) {
      qIntBdryEdge = true;
    }
    // ...or until a boundary is hit
    pC = apFace[i]->pCCellOpposite(pTC);
    if (pC->eType() == Cell::eTriBFace) {
      iOnBdryEdge2++;
      // Clean up the new faces and cells which have been created but
      // which won't be used because this is a boundary edge.  Then
      // insert onto the boundary edge.
      int ii = i;
      while (ii > 0) {
	ii--;
	apCNew[ii]->vMarkDeleted();
	apNewFace[ii]->vMarkDeleted();
	apIntFace[ii]->vMarkDeleted();
      }
      return (iInsertOnBdryEdge(pVNew, pV0, pV1, apFace[i], qSwap));
    }
    assert(pC->eType() == Cell::eTet);
    pTC = dynamic_cast<TetCell*>(pC);

    // Record existing stuff
    apCell[i] = pTC;
    apFaceB[i] = pTC->pFFaceOpposite(pVB);
    apFaceC[i] = pTC->pFFaceOpposite(pVC);

    // Add new stuff
    apCNew[i] = pCNewCell();
    apCNew[i]->vSetRegion(apCell[i]->iRegion());
    apNewFace[i] = pFNewFace();
    apIntFace[i] = pFNewFace();
  }
  if (qIntBdryEdge)
    iOnIntBdryEdge++;
  else
    iOnEdge++;
  // Make a copy of the first new face so that the connectivity will be
  // correct and continuous
  apNewFace[i] = apNewFace[0];
  int iNWrapCells = i;
  assert(iNWrapCells <= iMaxWrapCells);

  //@@ Reconnect everything properly

  //@@@ Loop over cells
  for (i = 0; i < iNWrapCells; i++) {
    // Modify existing faces
    apFace[i]->vReplaceVert(pVB, pVNew);
    apFaceC[i]->vReplaceCell(apCell[i], apCNew[i]);
    // Modify existing cell
    apCell[i]->vReplaceFace(apFaceC[i], apIntFace[i]);
    // Set up two new faces
    if (i == 0)
      apNewFace[i]->vAssign(apCNew[iNWrapCells-1], apCNew[i],
			    pVB, apVert[i], pVNew);
    else
      apNewFace[i]->vAssign(apCNew[i-1], apCNew[i],
			    pVB, apVert[i], pVNew);
    apNewFace[i]->vSetFaceLoc(apFace[i]->iFaceLoc());

    apIntFace[i]->vAssign(apCell[i], apCNew[i],
			  apVert[i], apVert[i+1], pVNew);
    // This new face is correctly tagged as interior by default

    // Set up connectivity for new cell; region has already been set
    apCNew[i]->vAssign(apFaceC[i], apIntFace[i],
			apNewFace[i], apNewFace[i+1]);
  }

  for (i = 0; i < iNWrapCells; i++) {
    // Check cell connectivity
    assert(apCell[i]->iFullCheck());
    // Check cell orientation
    assert(iOrient3D(apVert[i+1], apVert[i], pVNew, pVC) == 1);
    assert(iOrient3D(apVert[i], apVert[i+1], pVNew, pVB) == 1);
    // Check all faces
    assert(apFace[i]->iFullCheck());
    assert(apFaceB[i]->iFullCheck());
    assert(apFaceC[i]->iFullCheck());
    assert(apNewFace[i]->iFullCheck());
    assert(apIntFace[i]->iFullCheck());
  }

  //@@ Set vertex hint for the new vertex and an old one that needs it
  pVNew->vSetHintFace(apNewFace[0]);
  pVB  ->vSetHintFace(apNewFace[0]);

  if (qIntBdryEdge)
    pVNew->vSetType(Vert::eBdryTwoSide);
  else
    pVNew->vSetType(Vert::eInterior);

  //@@ Do face swapping for all the external faces
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);
  SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
  int iRetVal = 0;
  if (qSwap) {
    for (i = 0; i < iNWrapCells; i++) {
      iRetVal += iFaceSwap(apFaceB[i]);
      iRetVal += iFaceSwap(apFaceC[i]);
    }
  }
  SUMAA_LOG_EVENT_END(NET_SWAPPING);
  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  return (iRetVal);
}

//@ Insert a point on a boundary edge
int VolMesh::iInsertOnBdryEdge(Vert* const pVNew, Vert* const pV0,
			       Vert* const pV1, Face* const pF,
			       const int qSwap)
{
  // Update for Periodic?
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(qSimplicial());
  assert(pF->iFullCheck());
  assert(pF->qHasVert(pV0));
  assert(pF->qHasVert(pV1));
  assert(pF->pCCellLeft ()->eType() == Cell::eTriBFace ||
	 pF->pCCellRight()->eType() == Cell::eTriBFace);

  if (!qAllowBdryChanges) return (-1);

  //@@ Identify the boundary face and the interior cell
  BFace* pBF0;
  TetCell* pTC;
  if (pF->pCCellLeft ()->eType() == Cell::eTriBFace) {
    assert(pF->pCCellRight()->eType() == Cell::eTet);
    pBF0 = dynamic_cast<TriBFace*>(pF->pCCellLeft());
    pTC  = dynamic_cast<TetCell*> (pF->pCCellRight());
  }
  else {
    assert(pF->pCCellLeft ()->eType() == Cell::eTet);
    pBF0 = dynamic_cast<TriBFace*>(pF->pCCellRight());
    pTC  = dynamic_cast<TetCell*> (pF->pCCellLeft());
  }

  Vert* pVA = pTC->pVVertOpposite(pF);
  Vert* pVB = pV0;
  Vert* pVC = pV1;
  Vert* pVD;
  if (pF->pVVert(0) != pVB && pF->pVVert(0) != pVC)
    pVD = pF->pVVert(0);
  else if (pF->pVVert(1) != pVB && pF->pVVert(1) != pVC)
    pVD = pF->pVVert(1);
  else {
    assert(pF->pVVert(2) != pVB && pF->pVVert(2) != pVC);
    pVD = pF->pVVert(2);
  }
  if (iOrient3D(pVA, pVB, pVC, pVD) != 1) {
    Vert* pV = pVC; pVC = pVB; pVB = pV;
  }

  //@@ Make lists of cells divided and their faces and verts
  //@@ Create new faces and cell for each original cell
  const int iMaxWrapCells = 25;
  // This is the chain of verts around the perimeter of the shape (Verts
  // B and C are constant)
  Vert* apVert[iMaxWrapCells+1];
  // Current and new cells
  Cell *apCell[iMaxWrapCells], *apCNew[iMaxWrapCells];
  // apFace contains the original faces separating cells around the edge
  Face* apFace[iMaxWrapCells+1];
  // These are faces opposite verts B and C, respectively
  Face *apFaceB[iMaxWrapCells], *apFaceC[iMaxWrapCells];
  // New faces formed by dividing apFace members
  Face* apNewFace[iMaxWrapCells+1];
  // New faces within the original cells
  Face* apIntFace[iMaxWrapCells];

  //@@@ Prolog:  set stuff up for the boundary face at this end
  apCell[0] = pTC;
  apVert[0] = pVD;
  apFace[0] = pF;
  apFaceB[0] = pTC->pFFaceOpposite(pVB);
  apFaceC[0] = pTC->pFFaceOpposite(pVC);
  apCNew[0] = pCNewCell();
  apCNew[0]->vSetRegion(apCell[0]->iRegion());
  apNewFace[0] = pFNewFace();
  apIntFace[0] = pFNewFace();

  int i = 0;
  Cell* pC;
  //@@@ Walk around the edge, gathering info and setting up new entities
  while (1) {
    i++;
    // These are things that get set up regardless of whether there's
    // another interior cell
    apVert[i] = pTC->pVVertOpposite(apFace[i-1]);
    apFace[i] = pTC->pFFaceOpposite(apVert[i-1]);
    apNewFace[i] = pFNewFace();

    // Go on to the next cell, if and only if it's an interior cell
    // rather than a boundary face
    pC = apFace[i]->pCCellOpposite(pTC);
    if (pC->eType() == Cell::eTriBFace) break;
    assert(pC->eType() == Cell::eTet);
    pTC = dynamic_cast<TetCell*>(pC);
    apCell[i] = pTC;
    apFaceB[i] = pTC->pFFaceOpposite(pVB);
    apFaceC[i] = pTC->pFFaceOpposite(pVC);
    apCNew[i] = pCNewCell();
    apCNew[i]->vSetRegion(apCell[i]->iRegion());
    apIntFace[i] = pFNewFace();
  }
  BFace* pBF1 = dynamic_cast<TriBFace*>(pC);
  int iNWrapCells = i;
  assert(iNWrapCells <= iMaxWrapCells);

  //@@ Reconnect everything properly

  //@@@ Prolog:  modify existing bdry face and add data for new one
  // Create a new boundary face
  BFace* pBFNew0 = pBFNewBFace();
  pBFNew0->vSetPatch(pBF0->pPatchPointer());
  pBFNew0->vAssign(apNewFace[0]);

  // Modify existing faces
  apFace[0]->vReplaceVert(pVB, pVNew);
  apFaceC[0]->vReplaceCell(apCell[0], apCNew[0]);
  // Set up two new faces
  apNewFace[0]->vAssign(pBFNew0, apCNew[0],
			pVB, apVert[0], pVNew);
  apNewFace[0]->vSetFaceLoc(Face::eBdryFace);

  apIntFace[0]->vAssign(apCell[0], apCNew[0],
			apVert[0], apVert[1], pVNew);

  // Modify existing cell
  apCell[0]->vReplaceFace(apFaceC[0], apIntFace[0]);

  // Set up connectivity for new cell; the region has already been set.
  apCNew[0]->vAssign(apFaceC[0], apIntFace[0],
		     apNewFace[0], apNewFace[1]);

  assert(apCell[0]->iFullCheck());
  assert(iOrient3D(apVert[1], apVert[0], pVB, pVC) == 1);
  assert(iOrient3D(apVert[1], apVert[0], pVNew, pVC) == 1);
  assert(iOrient3D(apVert[0], apVert[1], pVNew, pVB) == 1);
  // Check all faces
  assert(apFace[0]->iFullCheck());
  assert(apFaceB[0]->iFullCheck());
  assert(apFaceC[0]->iFullCheck());
  assert(apNewFace[0]->iFullCheck());
  assert(apIntFace[0]->iFullCheck());

  //@@@ Loop over cells
  for (i = 1; i < iNWrapCells; i++) {
    // Modify existing faces
    apFace[i]->vReplaceVert(pVB, pVNew);
    apFaceC[i]->vReplaceCell(apCell[i], apCNew[i]);
    // Modify existing cell
    apCell[i]->vReplaceFace(apFaceC[i], apIntFace[i]);
    // Set up two new faces
    apNewFace[i]->vAssign(apCNew[i-1], apCNew[i],
			  pVB, apVert[i], pVNew);
    apNewFace[i]->vSetFaceLoc(apFace[i]->iFaceLoc());

    apIntFace[i]->vAssign(apCell[i], apCNew[i],
			  apVert[i], apVert[i+1], pVNew);
    // This new face is correctly tagged as interior by default

    // Set up new cell, including a region marker.
    apCNew[i]->vAssign(apFaceC[i], apIntFace[i],
			apNewFace[i], apNewFace[i+1]);
    apCNew[i]->vSetRegion(apCell[i]->iRegion());

    assert(apCell[i]->iFullCheck());
    assert(apCNew[i-1]->iFullCheck());
    // Check cell orientation
    assert(iOrient3D(apVert[i+1], apVert[i], pVNew, pVC) == 1);
    assert(iOrient3D(apVert[i], apVert[i+1], pVNew, pVB) == 1);
    // Check all faces
    assert(apFace[i]->iFullCheck());
    assert(apFaceB[i]->iFullCheck());
    assert(apFaceC[i]->iFullCheck());
    assert(apNewFace[i]->iFullCheck());
    assert(apIntFace[i]->iFullCheck());
  }

  //@@ Set vertex hint for the new vertex and an old one that needs it
  pVNew->vSetHintFace(apNewFace[0]);
  pVB->  vSetHintFace(apNewFace[0]);

  // Check to see whether this bdry edge is on a bdry segment, by
  // checking whether the BFace's on either side of the original edge
  // lie in the same bdry patch.

  BdryPatch* pBP0 = pBF0->pPatchPointer();
  BdryPatch* pBP1 = pBF1->pPatchPointer();
  if (pBP0 == pBP1) {
    pVNew->vSetType(Vert::eBdry);
  }
  else {
    pVNew->vSetType(Vert::eBdryCurve);
  }

  //@@@ Epilog: Take care of faces on the second boundary
  // Create a new boundary face
  BFace* pBFNew1 = pBFNewBFace();
  pBFNew1->vSetPatch(pBF1->pPatchPointer());
  pBFNew1->vAssign(apNewFace[iNWrapCells]);

  // Modify existing face
  apFace[iNWrapCells]->vReplaceVert(pVB, pVNew);

  // Set up new face
  apNewFace[iNWrapCells]->vAssign(apCNew[iNWrapCells-1], pBFNew1,
				  pVB, apVert[iNWrapCells], pVNew);
  apNewFace[iNWrapCells]->vSetFaceLoc(Face::eBdryFace);

  // Check the last cell
  assert(apCNew[iNWrapCells-1]->iFullCheck());
  // Check the last faces
  assert(apFace[iNWrapCells]->iFullCheck());
  assert(apNewFace[iNWrapCells]->iFullCheck());

  //@@ Do face swapping for all the external faces
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);
  SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
  int iRetVal = 0;
  if (qSwap) {
    for (i = 0; i < iNWrapCells; i++) {
      iRetVal += iFaceSwap(apFaceB[i]);
      iRetVal += iFaceSwap(apFaceC[i]);
    }
  }
  SUMAA_LOG_EVENT_END(NET_SWAPPING);
  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  return (iRetVal);
}

//@ Insert a point on a face, either interior or boundary
int VolMesh::iInsertOnFace(Vert* const pVNew, Face * const pF,
			   Cell * const pC, const int 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.
{
  assert(pC->iFullCheck());
  assert(pF->iFullCheck());
  assert(pF->pCCellLeft() == pC || pF->pCCellRight() == pC);
  assert(pC->qHasFace(pF));

  if (pF->eType() != Face::eTriFace) return 0;

  Vert * pVA = pF->pVVert(0);
  Vert * pVB = pF->pVVert(1);
  Vert * pVC = pF->pVVert(2);

  Cell* pCD = pF->pCCellLeft ();
  Cell* pCE = pF->pCCellRight();

  //@@ Boundary face case
  if (pCD->eType() == Cell::eTriBFace ||
      pCE->eType() == Cell::eTriBFace) {
    // Update for Periodic?
    // Don't insert on precious boundaries.
    if (!qAllowBdryChanges) return (-1);

    iOnBdryFace++;
    Vert *pVD = pC->pVVertOpposite(pF);
    BFace* pBF;
    if (pCD == pC) {
      assert(pCE->eType() == Cell::eTriBFace);
      pBF = dynamic_cast<TriBFace*>(pCE);
      pVB = pF->pVVert(2), pVC = pF->pVVert(1);
    }
    else { // The verts are already okay.
      assert(pCD->eType() == Cell::eTriBFace);
      pBF = dynamic_cast<TriBFace*>(pCD);
    }
    assert(iOrient3D(pVA, pVB, pVC, pVD) == 1);

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

    assert(iOrient3D(pVA, pVD, pVB, pVNew) == 1);
    assert(iOrient3D(pVB, pVD, pVC, pVNew) == 1);
    assert(iOrient3D(pVC, pVD, pVA, pVNew) == 1);
    Cell* pCA = pC;
    Cell* pCB = pCNewCell();
    Cell* pCC = pCNewCell();

    //@@@ Extract the last info you need from the old bdry face and nuke it
    assert(pBF->eType() == Cell::eTriBFace ||
	   pBF->eType() == Cell::eIntTriBFace);
    // Cast the boundary face to a TriBFaceBase, which will definitely
    // know how to find a bdry patch.
    BdryPatch* pBP = pBF->pPatchPointer();
    pBF->vMarkDeleted();

    //@@@ Create new triangular boundary faces
    BFace* pBFA = pBFNewBFace(3);
    BFace* pBFB = pBFNewBFace(3);
    BFace* pBFC = pBFNewBFace(3);
    pBFA->vSetPatch(pBP);
    pBFB->vSetPatch(pBP);
    pBFC->vSetPatch(pBP);

    //@@@ Create new faces and tag with locations
    Face* pFA1 = pFNewFace();
    pFA1->vAssign(pCB, pCC, pVD, pVA, pVNew);

    Face* pFB1 = pFNewFace();
    pFB1->vAssign(pCC, pCA, pVD, pVB, pVNew);

    Face* pFC1 = pFNewFace();
    pFC1->vAssign(pCA, pCB, pVD, pVC, pVNew);

    Face* pFA2 = pFNewFace();
    pFA2->vAssign(pBFA, pCA, pVB, pVC, pVNew);

    Face* pFB2 = pFNewFace();
    pFB2->vAssign(pBFB, pCB, pVC, pVA, pVNew);

    Face* pFC2 = pF;
    pFC2->vAssign(pBFC, pCC, pVA, pVB, pVNew);

    pFA2->vSetFaceLoc(Face::eBdryFace);
    pFB2->vSetFaceLoc(Face::eBdryFace);
    pFC2->vSetFaceLoc(Face::eBdryFace);

    //@@@ Connect boundary faces to real faces properly
    pBFA->vAssign(pFA2);
    pBFB->vAssign(pFB2);
    pBFC->vAssign(pFC2);

    //@@@ Clean up vertices
    pVNew->vSetHintFace(pFA1);
    pVA->vSetHintFace(pFC2);
    pVB->vSetHintFace(pFA2);
    pVC->vSetHintFace(pFB2);
    pVNew->vSetType(Vert::eBdry); // Can't be a BdryApex or BdryCurve

    //@@@ Fix up new cells
    pCA->vAssign(pFA, pFA2, pFB1, pFC1);
    pCB->vAssign(pFB, pFB2, pFC1, pFA1);
    pCC->vAssign(pFC, pFC2, pFA1, pFB1);
    pCA->vSetRegion(pC->iRegion());
    pCB->vSetRegion(pC->iRegion());
    pCC->vSetRegion(pC->iRegion());

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

    //@@@ Swap locally
    SUMAA_LOG_EVENT_END(INSERT_VERTEX);
    SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
    int iRetVal = 0;
    if (qSwap) {
      iRetVal += iFaceSwap(pFA);
      iRetVal += iFaceSwap(pFB);
      iRetVal += iFaceSwap(pFC);
    }
    SUMAA_LOG_EVENT_END(NET_SWAPPING);
    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
    return(iRetVal);
  } // Done with insertion on boundary face
  //@@ Something other than tets
  else if (pCD->eType() != Cell::eTet ||
	   pCE->eType() != Cell::eTet) return 0;
  //@@ Insertion on an interior face
  else {
    if (pF->iFaceLoc() == Face::eBdryTwoSide)
      iOnIntBdryFace++;
    else
      iOnFace++;
    Vert *pVD = pCD->pVVertOpposite(pF);
    Vert *pVE = pCE->pVVertOpposite(pF);
    assert(iOrient3D(pVA, pVC, pVB, pVD) == 1);
    assert(iOrient3D(pVA, pVB, pVC, pVE) == 1);

    Face* pFA1 = pCE->pFFaceOpposite(pVA);
    Face* pFB1 = pCE->pFFaceOpposite(pVB);
    Face* pFC1 = pCE->pFFaceOpposite(pVC);

    Face* pFA2 = pCD->pFFaceOpposite(pVA);
    Face* pFB2 = pCD->pFFaceOpposite(pVB);
    Face* pFC2 = pCD->pFFaceOpposite(pVC);

    assert(iOrient3D(pVA, pVE, pVB, pVNew) == 1);
    assert(iOrient3D(pVB, pVE, pVC, pVNew) == 1);
    assert(iOrient3D(pVC, pVE, pVA, pVNew) == 1);
    assert(iOrient3D(pVA, pVB, pVD, pVNew) == 1);
    assert(iOrient3D(pVB, pVC, pVD, pVNew) == 1);
    assert(iOrient3D(pVC, pVA, pVD, pVNew) == 1);

    //@@@ Create new cells and tag each with its region
    Cell* pCA1 = pCE;
    Cell* pCB1 = pCNewCell();
    Cell* pCC1 = pCNewCell();
    Cell* pCA2 = pCD;
    Cell* pCB2 = pCNewCell();
    Cell* pCC2 = pCNewCell();

    pCB1->vSetRegion(pCA1->iRegion());
    pCC1->vSetRegion(pCA1->iRegion());
    pCB2->vSetRegion(pCA2->iRegion());
    pCC2->vSetRegion(pCA2->iRegion());

    //@@@ Create new faces. Tag each descendent of pF as the same type
    //as pF.
    Face* pFA  = pF;
    pFA->vAssign(pCA2, pCA1, pVNew, pVB, pVC);
    Face* pFB  = pFNewFace();
    pFB->vAssign(pCB2, pCB1, pVNew, pVC, pVA);
    Face* pFC  = pFNewFace();
    pFC->vAssign(pCC2, pCC1, pVNew, pVA, pVB);
    pFB->vSetFaceLoc(pF->iFaceLoc());
    pFC->vSetFaceLoc(pF->iFaceLoc());

    Face* pFAE = pFNewFace();
    pFAE->vAssign(pCB1, pCC1, pVA, pVNew, pVE);
    Face* pFBE = pFNewFace();
    pFBE->vAssign(pCC1, pCA1, pVB, pVNew, pVE);
    Face* pFCE = pFNewFace();
    pFCE->vAssign(pCA1, pCB1, pVC, pVNew, pVE);

    Face* pFAD = pFNewFace();
    pFAD->vAssign(pCC2, pCB2, pVA, pVNew, pVD);
    Face* pFBD = pFNewFace();
    pFBD->vAssign(pCA2, pCC2, pVB, pVNew, pVD);
    Face* pFCD = pFNewFace();
    pFCD->vAssign(pCB2, pCA2, pVC, pVNew, pVD);

    //@@@ Update some vertex information
    pVNew->vSetHintFace(pFAE);
    // Default is interior; only change if we must.
    if (pF->iFaceLoc() == Face::eBdryTwoSide)
      pVNew->vSetType(Vert::eBdryTwoSide);
    pVA->vSetHintFace(pFC);
    pVB->vSetHintFace(pFA);
    pVC->vSetHintFace(pFB);

    //@@@ Fix up connectivity for cells
    pCA1->vAssign(pFA1, pFBE, pFCE, pFA);
    pCB1->vAssign(pFB1, pFCE, pFAE, pFB);
    pCC1->vAssign(pFC1, pFAE, pFBE, pFC);
    pCA2->vAssign(pFA2, pFBD, pFCD, pFA);
    pCB2->vAssign(pFB2, pFCD, pFAD, pFB);
    pCC2->vAssign(pFC2, pFAD, pFBD, pFC);

    //@@@ Fix up connectivity for faces
    pFA1->vReplaceCell(pCE, pCA1);
    pFB1->vReplaceCell(pCE, pCB1);
    pFC1->vReplaceCell(pCE, pCC1);
    pFA2->vReplaceCell(pCD, pCA2);
    pFB2->vReplaceCell(pCD, pCB2);
    pFC2->vReplaceCell(pCD, pCC2);

    // Check cells
    assert(pCA1->iFullCheck());
    assert(pCB1->iFullCheck());
    assert(pCC1->iFullCheck());
    assert(pCA2->iFullCheck());
    assert(pCB2->iFullCheck());
    assert(pCC2->iFullCheck());

    // Check faces
    assert(pFA1->iFullCheck());
    assert(pFB1->iFullCheck());
    assert(pFC1->iFullCheck());
    assert(pFA2->iFullCheck());
    assert(pFB2->iFullCheck());
    assert(pFC2->iFullCheck());
    assert(pFA ->iFullCheck());
    assert(pFB ->iFullCheck());
    assert(pFC ->iFullCheck());
    assert(pFAD->iFullCheck());
    assert(pFBD->iFullCheck());
    assert(pFCD->iFullCheck());
    assert(pFAE->iFullCheck());
    assert(pFBE->iFullCheck());
    assert(pFCE->iFullCheck());

    //@@@ Swap locally
    SUMAA_LOG_EVENT_END(INSERT_VERTEX);
    SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
    int iRetVal = 0;
    if (qSwap) {
      iRetVal += iFaceSwap(pFA1);
      iRetVal += iFaceSwap(pFB1);
      iRetVal += iFaceSwap(pFC1);
      iRetVal += iFaceSwap(pFA2);
      iRetVal += iFaceSwap(pFB2);
      iRetVal += iFaceSwap(pFC2);
    }
    SUMAA_LOG_EVENT_END(NET_SWAPPING);
    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
    return(iRetVal);
  } // End of interior face insertion
}

//@ Insert a point inside a cell
int VolMesh::iInsertInInterior(Vert* const  pVNew, Cell* const pC,
			       const int qSwap)
{
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(pC->iFullCheck());
  iInterior++;

  // 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);
  Vert* pVD = pC->pVVert(3);

  Face* pFA = pC->pFFaceOpposite(pVA);
  Face* pFB = pC->pFFaceOpposite(pVB);
  Face* pFC = pC->pFFaceOpposite(pVC);
  Face* pFD = pC->pFFaceOpposite(pVD);

  // Get indices for new cells
  Cell* pCNewC0 = pCNewCell();
  Cell* pCNewC1 = pCNewCell();
  Cell* pCNewC2 = pCNewCell();
  pCNewC0->vSetRegion(pC->iRegion());
  pCNewC1->vSetRegion(pC->iRegion());
  pCNewC2->vSetRegion(pC->iRegion());

  // Create new faces.  These all fall in the interior.
  Face* pFNewFAB = pFNewFace();
  pFNewFAB->vAssign(pC, pCNewC0, pVA, pVB, pVNew);

  Face* pFNewFBC = pFNewFace();
  pFNewFBC->vAssign(pC, pCNewC1, pVB, pVC, pVNew);

  Face* pFNewFCA = pFNewFace();
  pFNewFCA->vAssign(pC, pCNewC2, pVC, pVA, pVNew);

  Face* pFNewFAD = pFNewFace();
  pFNewFAD->vAssign(pCNewC0, pCNewC2, pVA, pVD, pVNew);

  Face* pFNewFBD = pFNewFace();
  pFNewFBD->vAssign(pCNewC1, pCNewC0, pVB, pVD, pVNew);

  Face* pFNewFCD = pFNewFace();
  pFNewFCD->vAssign(pCNewC2, pCNewC1, pVC, pVD, pVNew);

  // The vertex falls in the interior, too.
  pVNew->vSetHintFace(pFNewFAB);

  // Set up faces for the new cells
  pC     ->vAssign(pFD, pFNewFAB, pFNewFBC, pFNewFCA);
  pCNewC0->vAssign(pFC, pFNewFAB, pFNewFAD, pFNewFBD);
  pCNewC1->vAssign(pFA, pFNewFBC, pFNewFBD, pFNewFCD);
  pCNewC2->vAssign(pFB, pFNewFCA, pFNewFCD, pFNewFAD);

  // Change the face-cell info for the original hull
  pFA->vReplaceCell(pC, pCNewC1);
  pFB->vReplaceCell(pC, pCNewC2);
  pFC->vReplaceCell(pC, pCNewC0);

  assert(pFA     ->iFullCheck());
  assert(pFB     ->iFullCheck());
  assert(pFC     ->iFullCheck());
  assert(pFD     ->iFullCheck());
  assert(pFNewFAB->iFullCheck());
  assert(pFNewFBC->iFullCheck());
  assert(pFNewFCA->iFullCheck());
  assert(pFNewFAD->iFullCheck());
  assert(pFNewFBD->iFullCheck());
  assert(pFNewFCD->iFullCheck());

  assert(pC     ->iFullCheck());
  assert(pCNewC0->iFullCheck());
  assert(pCNewC1->iFullCheck());
  assert(pCNewC2->iFullCheck());

  assert(iOrient3D(pVA, pVB, pVC, pVNew) == 1);
  assert(iOrient3D(pVA, pVD, pVB, pVNew) == 1);
  assert(iOrient3D(pVB, pVD, pVC, pVNew) == 1);
  assert(iOrient3D(pVC, pVD, pVA, pVNew) == 1);

  // Do some local swapping to clean up
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);
  SUMAA_LOG_EVENT_BEGIN(NET_SWAPPING);
  int iRetVal = 0;
  if (qSwap) {
    iRetVal += iFaceSwap(pFA);
    iRetVal += iFaceSwap(pFB);
    iRetVal += iFaceSwap(pFC);
    iRetVal += iFaceSwap(pFD);
  }
  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 VolMesh::qInsertPoint(const double adPointIn[3], Cell* const pCGuess,
			   int * const piSwaps, const bool qSwap,
			   const bool qForce, Vert* pVNew)
  // pVNew is non-pVInvalidVert if and only if the point being
  // inserted already belongs to the point set of the mesh; in this
  // case, pVNew->adCoords() == adPointIn.
{
#ifdef SUMAA_LOG
  assert(SUMAA_LOG_current == INSERT_VERTEX);
#endif
  assert(qSimplicial());
  assert(qForce || dfLengthScale);
  *piSwaps = 0;
  Cell* pC = pCGuess;
  if (pC->eType() == Cell::eTriBFace) {
    Face *pFTmp = pC->pFFace(0);
    pC = pFTmp->pCCellOpposite(pC);
  }
  double adPoint[] = {adPointIn[0], adPointIn[1], adPointIn[2]};

  if (pVNew->qValid()) {
    assert(qForce); // Can't twiddle the point location if it's an
		    // existing vertex.
    // Check that the points coincide
    assert(dDIST3D(pVNew->adCoords(), adPointIn) < 1.e-12);
  }

  //@@ Compute barycentric coordinates
  double adBary[4];
  dynamic_cast<TetCell*>(pC)->vBarycentrics(adPoint, adBary);
  int iNegBary = ((adBary[0] < -1.e-10 ? 1 : 0) +
		  (adBary[1] < -1.e-10 ? 1 : 0) +
		  (adBary[2] < -1.e-10 ? 1 : 0) +
		  (adBary[3] < -1.e-10 ? 1 : 0));

  if (iNegBary) {
    //@@ Find which cell the site falls within (or near, if outside the mesh)
    assert(pCGuess->iFullCheck());
    Cell *pCNew = pCGuess, *pCOld = pCGuess;
    int iOrient = 0;
    int iSteps = 0;
    bool qDone;
    Face *pF;
    do {
      qDone = true;
      iSteps ++;
      int i;
      double adCent[3];
      pC->vCentroid(adCent);
      for (i = 0; i < 4; i++) {
	pF = pC->pFFace(i);
	// If the face is pointing away from this cell, then signs need to
	// be reversed to walk in the proper direction
	int iSign = (pF->pCCellLeft() == pC) ? -1 : 1;
	iOrient = iOrient3D(pF, adPoint) * iSign;
	// If the vertex lies outside the current cell and behind the
	// current face, walk in that direction.  If the walk would
	// involve going outside the boundary, continue checking to try to
	// find another face to walk across; this improves behavior of
	// walks which try to insert a boundary point.
	if (iOrient == -1) {
	  if (iSteps < 2000) {
	    pCNew = pF->pCCellOpposite(pC);
	    qDone = false;
	    if (pCNew->eType() != Cell::eTriBFace)
	      break;
	  }
	  else {
	    // Assume that we may be entering orbit...
	    int iOrient01 = iOrient3D(adCent, pF->pVVert(0)->adCoords(),
				      pF->pVVert(1)->adCoords(),
				      adPoint);
	    int iOrient12 = iOrient3D(adCent, pF->pVVert(1)->adCoords(),
				      pF->pVVert(2)->adCoords(),
				      adPoint);
	    int iOrient20 = iOrient3D(adCent, pF->pVVert(2)->adCoords(),
				      pF->pVVert(0)->adCoords(),
				      adPoint);
	    // If the target point lies on or inside the infinite
	    // triangular pyramid defined by the cell centroid and the
	    // vertices of pF, then stepping through this face
	    // definitely moves us in the direction of the target.
	    if (iOrient01*iOrient12 != -1 &&
		iOrient12*iOrient20 != -1 &&
		iOrient20*iOrient01 != -1) {
	      pCNew = pF->pCCellOpposite(pC);
	      qDone = false;
	      if (pCNew->eType() != Cell::eTriBFace)
		break;
	    }
	  }
	}
      }
      pCOld = pC;
      pC = pCNew;
    } while (iOrient == -1 && pC->eType() != Cell::eTriBFace &&
	     iSteps < 5000);

    if (!qDone && pC->eType() != Cell::eTriBFace) {
      return 0;
    }

    // Restore identity of last interior cell visited, if needed
    pC = pCOld;
    //@@ Check for near-coincidence with existing verts

    if ((!qForce) && (qInsertInThisCell(pC, dfLengthScale) != 1)) {
      return 0;
    }
    //@@ Re-compute barycentric coordinates
    dynamic_cast<TetCell*>(pC)->vBarycentrics(adPoint, adBary);
  }

  if (!qForce) {
    double dMinDist = 0.6 * dLengthScale(adPoint);
    if (dDIST3D(adPoint, pC->pVVert(0)->adCoords()) < dMinDist ||
	dDIST3D(adPoint, pC->pVVert(1)->adCoords()) < dMinDist ||
	dDIST3D(adPoint, pC->pVVert(2)->adCoords()) < dMinDist ||
	dDIST3D(adPoint, pC->pVVert(3)->adCoords()) < dMinDist) {
      return (0);
    }
  }
  //@@ Snap to face/edge

  // Define the threshold of barycentric coordinate that will induce
  // snap to face/edge
  double dBaryThresh;
  if (qForce) dBaryThresh = 1.e-10;
  else        dBaryThresh = 0.1;

  double dSum = 0.;
  int i, iLowBary = 0;
  for (i = 0; i <= 3; i ++)
    if (adBary[i] < dBaryThresh) {
      adBary[i] = 0;
      iLowBary++;
    }
    else
      dSum += adBary[i];
  assert(iFuzzyComp(dSum, 0) != 0);
  assert(iLowBary <= 3);
  switch (iLowBary) {
  case 0:
    { // Insertion squarely inside a cell
      if (!(pVNew->qValid()))
	pVNew = this->pVNewVert();
      pVNew->vSetCoords(3, adPoint);
      *piSwaps = iInsertInInterior(pVNew, pC, qSwap);
      return (1);
    }
  case 1:
    { // Insertion on a face (possibly a boundary face; let the
      // callee deal with that question)
      // Find the face opposite the zero barycentric coordinate
      const Vert *pVZero;
      if      (adBary[0] < dBaryThresh)
	pVZero = pC->pVVert(0);
      else if (adBary[1] < dBaryThresh)
	pVZero = pC->pVVert(1);
      else if (adBary[2] < dBaryThresh)
	pVZero = pC->pVVert(2);
      else
	pVZero = pC->pVVert(3); // (adBary[3] == 0)
      Face *pF = pC->pFFaceOpposite(pVZero);
      // Project the proposed point onto this face; this will give the
      // face circumcenter
      assert(pF->eType() == Face::eTriFace);
      if (! qForce) dynamic_cast<TriFace*>(pF)->vCircumcenter(adPoint);

      // Check barycentrics again.
      dynamic_cast<TetCell*>(pC)->vBarycentrics(adPoint, adBary);
      double dFaceBaryThresh;
      if (qForce) dFaceBaryThresh = 1.e-10;
      else        dFaceBaryThresh = 0.2;

      iLowBary = ((adBary[0] < dFaceBaryThresh ? 1 : 0) +
		  (adBary[1] < dFaceBaryThresh ? 1 : 0) +
		  (adBary[2] < dFaceBaryThresh ? 1 : 0) +
		  (adBary[3] < dFaceBaryThresh ? 1 : 0));
      if (iLowBary == 1) { // In the interior of the face
	if (!(pVNew->qValid()))
	  pVNew = pVNewVert();
	pVNew->vSetCoords(3, adPoint);
	*piSwaps = iInsertOnFace(pVNew, pF, pC, qSwap);
	if (*piSwaps == -1) {
	  // Tried to insert on a bdry face when bdry's are not allowed
	  // to be changed.
	  pVNew->vMarkDeleted();
	  *piSwaps = 0;
	  return 0;
	}
	else
	  return 1;
      }
      for (int iBary = 0; iBary <= 3; iBary ++)
	if (adBary[iBary] < dFaceBaryThresh) adBary[iBary] = 0;

      assert(iLowBary == 2);
    }
    // Deliberate fallthrough if the face case becomes an edge case
  case 2:
    {
      // Must be insertion on an edge (possibly a boundary edge)
      // As with face insertion, let the callee handle the boundary case

      // Find the two verts w/ non-zero barycentrics; the new vert falls
      // on the edge between them
      Vert *pV0, *pV1;
      Face *pF;

      // Not a very elegant way to identify these things, but so what?
      if (adBary[0] < dBaryThresh) {
	if (adBary[1] < dBaryThresh) {
	  pV0 = pC->pVVert(2);
	  pV1 = pC->pVVert(3);
	}
	else if (adBary[2] < dBaryThresh) {
	  pV0 = pC->pVVert(1);
	  pV1 = pC->pVVert(3);
	}
	else {
	  pV0 = pC->pVVert(2);
	  pV1 = pC->pVVert(1);
	}
	pF  = pC->pFFaceOpposite(pC->pVVert(0));
      }
      else if (adBary[1] < dBaryThresh) {
	if (adBary[2] < dBaryThresh) {
	  pV0 = pC->pVVert(0);
	  pV1 = pC->pVVert(3);
	}
	else {
	  pV0 = pC->pVVert(2);
	  pV1 = pC->pVVert(0);
	}
	pF  = pC->pFFaceOpposite(pC->pVVert(1));
      }
      else {
	pV0 = pC->pVVert(0);
	pV1 = pC->pVVert(1);
	pF  = pC->pFFaceOpposite(pC->pVVert(2));
      }
      assert(pV0->qValid());
      assert(pV1->qValid());
      assert(pF->qValid());
      if (!(pVNew->qValid()))
	pVNew = this->pVNewVert();
      if (! qForce) {
	double adActual[] = {0.5 * (pV0->dX() + pV1->dX()),
			   0.5 * (pV0->dY() + pV1->dY()),
			   0.5 * (pV0->dZ() + pV1->dZ())};
	pVNew->vSetCoords(3, adActual);
	*piSwaps = iInsertOnEdge(pVNew, pV0, pV1, pF, qSwap);
      }
      else {
	pVNew->vSetCoords(3, adPoint);
	*piSwaps = iInsertOnEdge(pVNew, pV0, pV1, pF, qSwap);
      }
      if (*piSwaps == -1) {
	// Tried to insert on a bdry edge when bdry's are not allowed
	// to be changed.
	pVNew->vMarkDeleted();
	*piSwaps = 0;
	return 0;
      }
      return (1);
    } // End edge insertion case
  case 3:
    // Inserting a point that's already in the mesh; this can occur when
    // re-inserting interior structured verts if one wasn't deleted in
    // the first place.  Bail out quietly, after marking the new vertex
    // (which is the attempted duplicate) as deleted.
    if (pVNew->qValid())
      pVNew->vMarkDeleted();
    return (0);
  default:
    assert2(0, "Invalid number of small barycentrics");
    return (0);
  } // End switch on the number of small barycentrics
} // End of point insertion
