#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <sys/stat.h>
#include "GR_config.h"
#include "GR_events.h"
#include "GR_Geometry.h"
#include "GR_SurfMesh.h"
#include "GR_VertTree.h"
#include "GR_VolMesh.h"
#include "GR_Constrain.h"

static const int iMaxSplittings = 12;

void VolMesh::vSplitEdge(const Vert* const pV0, const Vert* const pV1,
			 SurfMesh& SM)
{
  // Force a split of the given surface edge within both the surface and
  // volume meshes.
  int iV0 = iVertIndex(pV0);
  int iV1 = iVertIndex(pV1);
  Vert *pVS0 = SM.pVVert(iV0);
  Vert *pVS1 = SM.pVVert(iV1);
  assert(dDIST3D(pVS0->adCoords(), pV0->adCoords()) < 1.e-12);
  assert(dDIST3D(pVS1->adCoords(), pV1->adCoords()) < 1.e-12);
  Face *pFS = pFInvalidFace;
  {
    std::set<Vert*> spVNear;
    std::set<Face*> spFNear;
    std::set<Cell*> spCTmp;
    std::set<BFace*> spBFTmp;
    bool qTmp;
    vNeighborhood(pVS0, spCTmp, spVNear, &spBFTmp, &qTmp, &spFNear);
    if (spVNear.count(pVS1) == 0) return;
    std::set<Face*>::iterator iterF;
    for (iterF = spFNear.begin(); iterF != spFNear.end(); iterF++) {
      pFS = *iterF;
      if (pFS->qHasVert(pVS0) && pFS->qHasVert(pVS1)) break;
    }
    assert(pFS->qHasVert(pVS0));
    assert(pFS->qHasVert(pVS1));
  }

  // Edge swap failed; have to solve this by splitting.
  double adNewLoc[3];
  bool qSmall0 = pVS0->qIsSmallAngleCorner();
  bool qSmall1 = pVS1->qIsSmallAngleCorner();
  if ((qSmall0 && qSmall1) ||
      (!qSmall0 && !qSmall1)) {
    // Bisect
    adNewLoc[0] = 0.5*(pV0->dX() + pV1->dX());
    adNewLoc[1] = 0.5*(pV0->dY() + pV1->dY());
    adNewLoc[2] = 0.5*(pV0->dZ() + pV1->dZ());
    vMessage(3, "Splitting edge at midpoint: %12.6f %12.6f %12.6f\n",
	     adNewLoc[0], adNewLoc[1], adNewLoc[2]);
  }
  else if (qSmall0) {
    // Split near pVS0.
    double dSplitDist = pVS0->dLS();
    double adVec[] = adDIFF3D(pV1->adCoords(), pV0->adCoords());
    double dEdgeLen = dMAG3D(adVec);
    while (dSplitDist > dEdgeLen/1.5) dSplitDist /= 2;
    assert(dSplitDist < 2./3. * dEdgeLen);

    vNORMALIZE3D(adVec);

    adNewLoc[0] = pV0->dX() + adVec[0] * dSplitDist;
    adNewLoc[1] = pV0->dY() + adVec[1] * dSplitDist;
    adNewLoc[2] = pV0->dZ() + adVec[2] * dSplitDist;
    vMessage(3, "Splitting edge near vert %4d: %12.6f %12.6f %12.6f\n",
	     iV0, adNewLoc[0], adNewLoc[1], adNewLoc[2]);
  }
  else {
    assert(qSmall1);
    // Split near pVS1.
    double dSplitDist = pVS1->dLS();
    double adVec[] = adDIFF3D(pV0->adCoords(), pV1->adCoords());
    double dEdgeLen = dMAG3D(adVec);
    while (dSplitDist > dEdgeLen/1.5) dSplitDist /= 2;
    assert(dSplitDist < 2./3. * dEdgeLen);

    vNORMALIZE3D(adVec);

    adNewLoc[0] = pV1->dX() + adVec[0] * dSplitDist;
    adNewLoc[1] = pV1->dY() + adVec[1] * dSplitDist;
    adNewLoc[2] = pV1->dZ() + adVec[2] * dSplitDist;
    vMessage(3, "Splitting edge near vert %4d: %12.6f %12.6f %12.6f\n",
	     iV1, adNewLoc[0], adNewLoc[1], adNewLoc[2]);
  }

  Vert *pVNewSurf = SM.pVNewVert();
  pVNewSurf->vSetCoords(3, adNewLoc);
  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  int iSwaps =  SM.iInsertOnFace(pVNewSurf, pFS, pFS->pCCellLeft(), true);
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);

  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  qInsertPoint(adNewLoc, pCBeginPipe(pV0, pV1), &iSwaps, 1, 1);
  Vert *pVNew = pVVert(iNumVerts() - 1);
  // This vertex was added during bdry recovery, along a curve.  Mark it
  // as a bdry curve point, and tag for possible deletion, in the case
  // that we're interested in removing as many constrained verts as possible.
  pVNew->vSetType(Vert::eBdryCurve);
  pVNew->vMarkToDelete();
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);
}

static void vRecoverAllEdges(VolMesh * const pVM,
			     SurfMesh * const pSM,
			     std::set<ConstrainEdge>& sCESubSegs,
			     std::set<ConstrainEdge>& sCEInFacets,
			     std::set<ConstrainFace>& sCF)
{
  std::set<ConstrainEdge> sCE(sCESubSegs);
  if (sCE.empty()) sCE.insert(sCEInFacets.begin(), sCEInFacets.end());

  int iSplittings = 0;
  if (!sCE.empty()) {
    //@@ Recover edges
    // Verts added to the volume mesh to constrain it; must be added to
    // surface mesh on return for constraint checking to work properly.

    int iChange;
    bool qDoInsertion = false;
    int iTotalPasses;
    do { // Begin recovering edges
      iTotalPasses = 0;
      int iNFruitlessPasses = 0;
      do { // Begin a pass to check all unconstrained edges
	vMessage(3, "Recovery pass %d (%d fruitless); insertion %s\n",
		 iTotalPasses, iNFruitlessPasses,
		 (qDoInsertion ? "allowed" : "not allowed"));
	iTotalPasses ++;
	iNFruitlessPasses ++;
	iChange = 0;
	int iVertsBefore = pVM->iNumVerts();
	std::set<ConstrainEdge>::iterator iter; 
	for (iter = sCE.begin(); iter != sCE.end(); ) {
	  Vert *pV0 = pVM->pVVert(iter->iV0);
	  Vert *pV1 = pVM->pVVert(iter->iV1);
	  int iRecovered = pVM->iRecoverEdge(pV0, pV1, *pSM, qDoInsertion);
	  switch (iRecovered) {
	  case 1: // Edge recovered
	    {
	      std::set<ConstrainEdge>::iterator iterTmp = iter++;
	      sCE.erase(iterTmp);
	      iNFruitlessPasses = 0;
	      iChange ++;
	      break;
	    }
	  case -1:
	    iChange ++;
	    // Deliberate fallthrough
	  default:
	    iter++;
	    break;
	  }
	} // Done checking non-constrained segments
	int iVertsAfter = pVM->iNumVerts();
	vMessage(2, "Status:%2d passes (%sinserting; last %d fruitless). Changes: %3d Left: %3d\n",
		 iTotalPasses, qDoInsertion ? "" : "not ",
		 iNFruitlessPasses, iChange, static_cast<int>(sCE.size()));
	if (sCE.empty()) {
	  qIsBdryConstrained(pVM, pSM, sCESubSegs, sCEInFacets, sCF);
	  // Always recover subsegments first, then work on edges
	  // within subfacets.
	  assert(sCE.empty());
	  if (!sCESubSegs.empty()) {
	    sCE.insert(sCESubSegs.begin(), sCESubSegs.end());
	  }
	  else {
	    sCE.insert(sCEInFacets.begin(), sCEInFacets.end());
	  }
	  // If swapping fixed all the edges, don't turn off insertion,
	  // or you could loop forever trying to do surface recovery.
	  if (qDoInsertion && (iVertsAfter != iVertsBefore)) {
	    qDoInsertion = false;
	    assert(iNFruitlessPasses == 0);
	    iTotalPasses = 0;
	  }
	}
	for (iter = sCE.begin(); iter != sCE.end(); iter++) {
	  ConstrainEdge CE = *iter;
	  Vert *pV0 = pVM->pVVert(CE.iV0);
	  Vert *pV1 = pVM->pVVert(CE.iV1);
	  vMessage(4, "Bad edge:\n");
	  vMessage(4, "%d %12.6f %12.6f %12.6f\n",
		   CE.iV0, pV0->dX(), pV0->dY(),
		   pV0->dZ());
	  vMessage(4, "%d %12.6f %12.6f %12.6f\n\n",
		   CE.iV1, pV1->dX(), pV1->dY(),
		   pV1->dZ());
	}
      } while (iNFruitlessPasses < 2 && iChange && (!sCE.empty())
	       && iTotalPasses < 10); // Done with edge recovery, at
				      // least for now

      // We might even be totally done...
      if (sCE.empty()) return;
      // Don't insert verts onto boundaries that are immutable.
      if (pVM->qBdryChangesAllowed()) {
	if (qDoInsertion) {
	  iSplittings++;
	  if (iSplittings <= iMaxSplittings) {
	    vMessage(2, "Splitting all unconstrained edges (%dth time).\n",
		     iSplittings);
	    vMessage(2, "  Constraint stats before:\n");
	    qIsBdryConstrained(pVM, pSM, sCESubSegs, sCEInFacets, sCF);
	    // Always recover subsegments first, then work on edges
	    // within subfacets.
	    sCE.clear();
	    if (!sCESubSegs.empty()) {
	      sCE.insert(sCESubSegs.begin(), sCESubSegs.end());
	    }
	    else {
	      sCE.insert(sCEInFacets.begin(), sCEInFacets.end());
	    }
	    // Since allowing insertion wasn't enough, force a split of
	    // all edges not yet constrained.  Enable swap recursion in
	    // the volume mesh while splitting.
	    pVM->vAllowSwapRecursion();
	    // ENCROACHMENT CHECK BEGIN
	    VertTree VT(*(pVM->pVAGetVerts()));
	    int iNumEncroached = 0;
	    // ENCROACHMENT CHECK END
	    int iSmallVerts = 0, iSmallEdges = 0, iTotalEdges = sCE.size();
	    std::set<ConstrainEdge>::iterator iter;
	    int iE = 0;
	    for (iter = sCE.begin(); iter != sCE.end(); iter++, iE++) {
	      ConstrainEdge CE = *iter;
	      Vert *pV0 = pVM->pVVert(CE.iV0);
	      Vert *pV1 = pVM->pVVert(CE.iV1);
	      vMessage(3, "   Edge %d.  Vert 0: %s.  Vert 1: %s.\n",
		       iE, (pV0->qIsSmallAngleCorner() ? "small" : " okay"),
		       (pV1->qIsSmallAngleCorner() ? "small" : " okay"));
	      iSmallVerts += (pV0->qIsSmallAngleCorner() ? 1 : 0) +
		(pV1->qIsSmallAngleCorner() ? 1 : 0);
	      iSmallEdges += (pV0->qIsSmallAngleCorner() ||
			      pV1->qIsSmallAngleCorner() ? 1 : 0);
	      // ENCROACHMENT CHECK BEGIN
	      // Check edges for encroachment; if the edge isn't
	      // encroached, we should be able to avoid getting here.
	      {
		double adMid[] = {(pV0->dX() + pV1->dX()) * 0.5,
				  (pV0->dY() + pV1->dY()) * 0.5,
				  (pV0->dZ() + pV1->dZ()) * 0.5};
		double dRad = dDIST3D(pV0->adCoords(), pV1->adCoords()) * 0.5001;
		std::vector<int> veciTmp = VT.veciBallQuery(adMid, dRad);
		assert(veciTmp.size() >= 2);
		// Will confirm later that the endpoint are included.
		bool qGot0 = false, qGot1 = false;
		if (veciTmp.size() > 2) {
		  iNumEncroached++;
		  // What's the actual radius ratio?
		  vMessage(4, "  Edge from %d to %d is encroached %d times.\n",
			   CE.iV0, CE.iV1,
			   static_cast<int>(veciTmp.size()) - 2);
		  double adEdge[] = adDIFF3D(pV0->adCoords(), pV1->adCoords());
		  for (int ii = veciTmp.size() - 1; ii >= 0; ii--) {
		    int iVert = veciTmp[ii];
		    if (iVert == CE.iV0) {
		      qGot0 = true;
		    }
		    else if (iVert == CE.iV1) {
		      qGot1 = true;
		    }
		    else {
		      double adDiff[] =
			adDIFF3D(pVM->pVVert(iVert)->adCoords(), adMid);
		      vMessage(4, "    Vert %4d is %5.2f%%r; angle wrt edge is %5.2f deg.\n",
			       iVert, dMAG3D(adDiff) / dRad * 100,
			       GR_acos(dDOT3D(adEdge, adDiff)/
				       (dMAG3D(adDiff)*dMAG3D(adEdge)))*180/M_PI);
		    }
		  }
		  assert(qGot0 && qGot1);
		}
	      }
	      // ENCROACHMENT CHECK END
	      pVM->vSplitEdge(pV0, pV1, *pSM);
	    }
	    vMessage(2, "  Out of %d edges, %d had a small angle vert (%d verts)\n",
		     iTotalEdges, iSmallEdges, iSmallVerts);
	    // ENCROACHMENT CHECK BEGIN
	    vMessage(2, "  %d edges were encroached.\n", iNumEncroached);
	    // ENCROACHMENT CHECK END
	    pVM->vDisallowSwapRecursion();
	    qDoInsertion = false;

	    vMessage(2, "  Constraint stats after:\n");
	    // Reset the list of stuff to recover
	    qIsBdryConstrained(pVM, pSM, sCESubSegs, sCEInFacets, sCF);
	    // Always recover subsegments first, then work on edges
	    // within subfacets.
	    sCE.clear();
	    if (!sCESubSegs.empty()) {
	      sCE.insert(sCESubSegs.begin(), sCESubSegs.end());
	    }
	    else {
	      sCE.insert(sCEInFacets.begin(), sCEInFacets.end());
	    }
	  }
	} // Done splitting all missing edges
	else {
	  // Try insertion if you haven't recently.
	  qDoInsertion = true;
	}
      }
    } while (iSplittings <= iMaxSplittings && (!sCE.empty()));
    // Done recovering edges
  } // Edges only need to be recovered if they are any missing in the first
    // place.
}

//@ Triangulation of a generalized boundary geometry
VolMesh::VolMesh(SurfMesh& SM, const int iQualMeas,
		 const double dRes, const double dGrade,
		 const double dMaxAngle)
  : Mesh(*(SM.pVAGetVerts())),
    qOKToChangeSurfaceEdges(true), qDoEdgeSwapping(true),
    qStrictPatchChecking(true),
    dMaxAngleForSurfSwap(dMaxAngle), VATriF(), VAQuadF(), VATriBF(),
    VAQuadBF(), VATet(), VAPyr(), VAPrism(), VAHex(), VAIntTriBF(),
    VAIntQuadBF(), pB3D()
{
  SUMAA_LOG_EVENT_BEGIN(INIT_TETRA);
  dMeshResolution = dRes;
  dLipschitzAlpha = dGrade;
  pB3D = SM.pBoundaryData();
  qBdryFill = false;
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(3, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);

  // Set up an outer triangulation of a big brick
  double dXMin, dXMax, dYMin, dYMax, dZMin, dZMax;
  dXMin = dYMin = dZMin =  1.e100;
  dXMax = dYMax = dZMax = -1.e100;

  int iVert, iNVerts = iNumVerts();
  for (iVert = 0; iVert < iNVerts; iVert++) {
    double dX = pVVert(iVert)->dX();
    double dY = pVVert(iVert)->dY();
    double dZ = pVVert(iVert)->dZ();
    dXMin = min(dXMin, dX);
    dXMax = max(dXMax, dX);
    dYMin = min(dYMin, dY);
    dYMax = max(dYMax, dY);
    dZMin = min(dZMin, dZ);
    dZMax = max(dZMax, dZ);
  }

  double dXMean = 0.5 * (dXMin + dXMax);
  double dYMean = 0.5 * (dYMin + dYMax);
  double dZMean = 0.5 * (dZMin + dZMax);
  double dDX = dXMax - dXMin;
  double dDY = dYMax - dYMin;
  double dDZ = dZMax - dZMin;

  // Make the region roughly cubical
  double dMaxAspect = 0.1;
  if (dDX < dMaxAspect * dDY) dDX = dMaxAspect*dDY;
  if (dDX < dMaxAspect * dDZ) dDX = dMaxAspect*dDZ;

  if (dDY < dMaxAspect * dDX) dDY = dMaxAspect*dDX;
  if (dDY < dMaxAspect * dDZ) dDY = dMaxAspect*dDZ;

  if (dDZ < dMaxAspect * dDX) dDZ = dMaxAspect*dDX;
  if (dDZ < dMaxAspect * dDY) dDZ = dMaxAspect*dDY;

  double dXLo = dXMean - dDX;
  double dXHi = dXMean + dDX;
  double dYLo = dYMean - dDY;
  double dYHi = dYMean + dDY;
  double dZLo = dZMean - dDZ;
  double dZHi = dZMean + dDZ;

  int iNumOrigVerts = iNumVerts();
  vMeshInBrick(dXLo, dXHi, dYLo, dYHi, dZLo, dZHi);

  vMessage(1, "Inserting surface points to create unconstrained volume mesh.\n");
  // Do Delaunay insertion of all the existing verts in the mesh.  This
  // is a little weird, because the vert data is already -stored- in the
  // mesh, it just isn't currently connected to anything.
  vSetSwapType(eDelaunay);
  int iNumSwaps = 0;
  int iStatusInterval = max(25, iNumOrigVerts / 10);
  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  std::set<Vert*> spVInsertLater;
  for (iVert = 0; iVert < iNumOrigVerts; iVert++) {
    Vert *pV = pVVert(iVert);
    // Don't insert this point yet if it isn't connected to anything in
    // the surface mesh.
    Vert *pVS = SM.pVVert(iVert);
    if (!pVS->pFHintFace()->qValid()) {
      spVInsertLater.insert(pVS);
      continue;
    }

    // Make sure that the small angle flag is copied over from the
    // surface mesh.
    pV->vSetSmallAngleCorner(pVS->qIsSmallAngleCorner());

    // Insert, then swap.  Force insertion at the specified point and
    // use the existing vertex location.
    int iNewSwaps;
    int iC = 0;

    while (pCCell(iC)->qDeleted()) iC++;

    const double * adCoords = pV->adCoords();

#ifndef NDEBUG
    bool qSuccess = qInsertPoint(adCoords, pCCell(iC), &iNewSwaps, true,
				 true, pV);
    assert(qSuccess);
#else
    (void) qInsertPoint(adCoords, pCCell(iC), &iNewSwaps, true,
			true, pV);
#endif
    iNumSwaps += iNewSwaps;
    if ((iVert+1) % iStatusInterval == 0) {
      vMessage(2, "Inserted %5d verts out of %5d; made %5d swaps.\n",
	       iVert+1, iNumOrigVerts, iNumSwaps);
    }
  }
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);
  vSetAllHintFaces();
  assert(qValid());

  bool qRealRecur = qSwapRecursionAllowed();
  vDisallowSwapRecursion();

  vMessage(1, "Recovering surface mesh...\n");
  // Add some dummy verts to surface mesh so that identical surface and
  // volume mesh indices are identical points.  Mark them for cleanup
  // before we exit.
  int ii;
  for (ii = 0; ii < 8; ii++) {
    Vert *pV = SM.pVNewVert();
    pV->vMarkDeleted();
  }

  //@@ Take a census to see which edges and faces need work to recover
  std::set<ConstrainEdge> sCESubSegs;
  std::set<ConstrainEdge> sCEInFacets;
  std::set<ConstrainFace> sCFSubFacets;

  // Find out how many bdry (sub)segments are not currently in the mesh.
  bool qOK = qIsBdryConstrained(this, &SM, sCESubSegs,
				sCEInFacets, sCFSubFacets);

  int iOuterPasses = 0;
  while (!qOK && iOuterPasses < 3) {
    iOuterPasses++;
    vMessage(2, "Recovering surface edges...\n");
    vRecoverAllEdges(this, &SM, sCESubSegs, sCEInFacets, sCFSubFacets);

    if (!sCESubSegs.empty()) {
      vWriteFile_Surface(SM, "failed", ".surf");
      vMessage(0, "Failed to recover an input bdry edge.\n");
      vFatalError("Increasing iMaxSplittings in " __FILE__ " -may- help.",
		  "3D constrained Delaunay tetrahedralization.");
    }

    if (!sCEInFacets.empty()) {
      vMessage(2, "%d %s\n", static_cast<int>(sCEInFacets.size()),
	       "generated surface edges within input facets are not present.");
      vMessage(2, "If all goes well, facet recovery will recover these.\n");
    }

    if (!sCFSubFacets.empty()) {
      vMessage(2, "Recovering surface facets...\n");
      // If all edges weren't recovered, there's a definite problem
      assert(sCESubSegs.empty());
      //@@ Recover faces
      qDoEdgeSwapping = false;
      int iTotalPasses = 0, iChange;
#ifndef NDEBUG
      vMessage(4, "Faces before attempt to recover: \n");
      for (int i_ = 0; i_ < iNumFaces(); i_++) {
	Face *pF = pFFace(i_);
	int iVA = iVertIndex(pF->pVVert(0));
	int iVB = iVertIndex(pF->pVVert(1));
	int iVC = iVertIndex(pF->pVVert(2));
	vMessage(4, "%4d %4d %4d %4d\n", i_, iVA, iVB, iVC);
      }
#endif

      do {
	iTotalPasses++;
	do {
	  iChange = 0;
	  std::set<ConstrainFace>::iterator iter;
	  for (iter = sCFSubFacets.begin(); iter != sCFSubFacets.end(); ) {
	    ConstrainFace CF = *iter;
	    int aiVerts[] = {CF.iV0, CF.iV1, CF.iV2};
	    ii = iRecoverFace(&SM, aiVerts);
	    if (ii) {
	      iChange += 1;
	      std::set<ConstrainFace>::iterator iterTmp = iter++;
	      sCFSubFacets.erase(iterTmp);
	    }
	    else iter++;
	  }
	} while (iChange && !sCFSubFacets.empty());
	if (iChange == 0) qDoEdgeSwapping = true;
	qOK = qIsBdryConstrained(this, &SM, sCESubSegs, sCEInFacets,
				 sCFSubFacets);
      } while ((!sCFSubFacets.empty()) &&
	       (sCESubSegs.empty() && sCEInFacets.empty()) &&
	       iTotalPasses < 10);
    } // Had to recover faces as a separate pass
  }
  if (!sCFSubFacets.empty() || !sCEInFacets.empty() || !sCESubSegs.empty())
    vFoundBug("creation of constrained 3D mesh from boundary data");

  //@@ Strip stuff outside the domain
  vMessage(1, "Removing extraneous tets and cleaning up...\n");

  // Mark eight verts at the brick corners for deletion
  for (iVert = iNumOrigVerts; iVert < iNumOrigVerts + 8; iVert++) {
    pVVert(iVert)->vMarkDeleted();
  }

  vDeleteExternal(&SM);
  if (!qValid())
    vFoundBug("creation of constrained 3D mesh from boundary data");

  // Mark vertices on internal boundaries.
  for (int iFace = 0; iFace < iNumFaces(); iFace++) {
    Face *pF = pFFace(iFace);
    if (pF->iFaceLoc() == Face::eBdryTwoSide) {
      pF->pVVert(0)->vSetType(Vert::eBdryTwoSide);
      pF->pVVert(1)->vSetType(Vert::eBdryTwoSide);
      pF->pVVert(2)->vSetType(Vert::eBdryTwoSide);
    }
  }

  // Mark vertices on external boundaries.  Must be done after internal
  // boundaries.
  for (int iBFace = 0; iBFace < iNumBdryFaces(); iBFace++) {
    BFace *pBF = pBFBFace(iBFace);
    if (pBF->eType() == Cell::eTriBFace ||
	pBF->eType() == Cell::eQuadBFace) {
      Face *pF = pBF->pFFace();
      pF->vSetFaceLoc(Face::eBdryFace);
      for (int iV = 0; iV < 3; iV++) {
	Vert *pV = pF->pVVert(iV);
	if (pV->iVertType() == Vert::eBdryTwoSide) {
	  pV->vSetType(Vert::eBdryCurve);
	}
	int iVertInd = iVertIndex(pV);
	if (iVertInd < iNumOrigVerts) {
	  // Original verts can be tagged with the type that they have
	  // in the surface mesh.
	  pV->vSetType(SM.pVVert(iVertInd)->iVertType());
	}
	assert(pV->iVertType() != Vert::eInterior &&
	       pV->iVertType() != Vert::eUnknown &&
	       pV->iVertType() != Vert::eInvalid);
      } // end loop over this face
    } // Done dealing with external BFace's
  } // end loop over BFace's

  // Swap to ensure (?) a Delaunay mesh!
  vMessage(3, "Swapping towards a Delaunay mesh.\n");
  vAllowSwapRecursion();
  vSetSwapType(eDelaunay);
  iSwap();

  // Now insert all the vertices that were removed from the bdry rep
  // before meshing, as well as verts that were extras to begin with.

  iStatusInterval = max(25, int(spVInsertLater.size()/10));
  int iStatusValue = iStatusInterval;
  int iVertsInserted = 0;
  iNumSwaps = 0;
  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  std::set<Vert*>::iterator iter = spVInsertLater.begin();
  std::set<Vert*>::iterator iterEnd = spVInsertLater.end();

  while (iter != iterEnd) {
    Vert *pVS = *iter;
    iter++;

    assert(!pVS->pFHintFace()->qValid() && !pVS->qDeleted());

    const double * adCoords = pVS->adCoords();
    vMessage(4, "Attempting to insert a vertex at (%10G %10G %10G)\n",
	     adCoords[0], adCoords[1], adCoords[2]);

    // This is a slow way to do the re-insertion, but it -will- work.
    bool qHit = false;
    Cell *pC;
    double adBary[4];
    for (int iC = 0; !qHit && iC < iNumTetCells(); iC++) {
      pC = pCCell(iC);
      if (pC->qDeleted()) continue;
      TetCell* pTC = dynamic_cast<TetCell*>(pC);
      assert(pTC != NULL);
      pTC->vBarycentrics(adCoords, adBary);
      double dMinBary = min(min(adBary[0], adBary[1]),
			     min(adBary[2], adBary[3]));
      qHit = (iFuzzyComp(dMinBary, 0) != -1);
    }
    if (!qHit) {
      vFatalError("Tried to insert a point outside the mesh",
		  "VolMesh constructor");
    }

    // This vertex may already have been inserted during surface
    // recovery, in which case we don't need to worry about it.
    double dMaxBary = max(max(adBary[0], adBary[1]),
			  max(adBary[2], adBary[3]));
    if (iFuzzyComp(dMaxBary, 1) == 0) {
      continue;
    }

    int iNewSwaps = 0;

#ifndef NDEBUG
    bool qSuccess =
#endif
      qInsertPoint(adCoords, pC, &iNewSwaps, true, true);
#ifndef NDEBUG
    assert(qSuccess);
#endif
    iVertsInserted++;
    iNumSwaps += iNewSwaps;
    if (iVertsInserted > iStatusValue) {
      vMessage(2, "Inserted %5d extra verts from the surface mesh, and made %5d swaps.\n",
	       iVertsInserted, iNumSwaps);
      iStatusValue += iStatusInterval;
    }
  }
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);
  if (iVertsInserted > 0) iSwap();

  vSetAllHintFaces();
  assert(qValid());
  SUMAA_LOG_EVENT_END(INIT_TETRA);

  if (qRealRecur) {
    vAllowSwapRecursion();
  }
  else {
    vDisallowSwapRecursion();
  }

  // Set up length scale
  vInitLengthScale(dMeshResolution, dLipschitzAlpha);

  // Clean up
  vPurge();
  SM.vPurge();
}
