#include <set>
#include "GR_Bdry3D.h"
#include "GR_config.h"
#include "GR_events.h"
#include "GR_misc.h"
#include "GR_Sort.h"
#include "GR_SurfMesh.h"
#include "GR_Bdry3D.h"

// This is here so that configure can find the library.
extern "C"
{
  void vGRUMMP_Surf_Lib (void) {}
}

//@ Copy constructor
SurfMesh::SurfMesh(const SurfMesh& SM)
     // Attach to the point set and connectivity, and compute the vertex
     // hint data.  NOTE:  THE CONNECTIVITY MUST NOT BE DESTROYED BY THE
     // ROUTINE WHICH CALLS THIS CONSTRUCTOR.  The assumption is that
     // the connectivity is created with 'new' and not 'delete'-ed until
     // the SurfMesh dies.
  : Mesh2D(SM), VAME(), iOnMultiEdge(0), pB3D(), apBPatch(NULL),
    iBPatchSize(0)
{
  assert(0);
}

SurfMesh& SurfMesh::operator=(const SurfMesh& SM)
{
  assert(0);
  // Do nothing for self-assignment
  if (this != &SM) {
    VAVerts  = SM.VAVerts;
    pQ       = SM.pQ;
    VAEdgeF  = SM.VAEdgeF;
    VAEdgeBF = SM.VAEdgeBF;
    VATri    = SM.VATri;
    VAQuad   = SM.VAQuad;
    pB3D     = SM.pB3D;
    VAME     = SM.VAME;

    vSetAllHintFaces();
  }
  if (!(qValid()))
    vFoundBug("surface mesh copy");
  return (*this);
}

//@ Create an initial surface triangulation from a Bdry3D object.

// For now, this code is set up strictly for getting initial
// triangulations from polyhedral input data.  As such, each surface
// patch is triangulated (constrained Delaunay) individually, and the
// SurfMesh constructor combines the topology info into a mesh
// (probably by using a cell-vertex to face-cell data rep conversion).

// In the longer term, this constructor will have to negotiate
// consistent discretizations of patch bounding curves and pass
// appropriate subsets of this info to the patches so they can
// discretize themselves properly.  Each patch (or possibly composite of
// patches, a la the 2D interpolated spline stuff) should triangulate
// itself on demand, returning a set of points and a connectivity.  In
// this scenario, verts may well be added by patch triangulation, and
// the constructor will have to handle this.
SurfMesh::SurfMesh(Bdry3D * const pB3DIn, const int iQualMeas) :
  Mesh2D(iQualMeas), VAME(), iOnMultiEdge(0), pB3D(pB3DIn),
  apBPatch(NULL), iBPatchSize(0)
{
  SUMAA_LOG_EVENT_BEGIN (INIT_SURF_TRI);
  vSetVerts(pB3D->pVAGetVerts());

  // Count up, in advance, how many cells there will be in the final
  // surface mesh.  Can't just use the Euler formula, because we don't
  // know how many closed surfaces and/or handles there are in the
  // input.
  int iP, iTotalCells = 0;
  for (iP = 0; iP < pB3D->iNumBdryPatches(); iP++) {
    int iNumSegs = pB3D->pBdry(iP)->iNumSegments();
    int iNumKids = pB3D->pBdry(iP)->iNumChildren();
    iTotalCells += iNumSegs + 2 * (iNumKids - 1);
  }
  apBPatch = new BdryPatch3D*[iTotalCells];
  iBPatchSize = iTotalCells;

  // Set up a place to put surface triangulation data as it comes
  // through.
  VATriCellCV VATriCV;

  vMessage(2, "Setting up surface mesh based on geometric data...\n");
  // Turn messaging down a notch during surface triangulation.
  iMessageStdoutLevel++;
  iMessageFileLevel++;
  // For every patch, triangulate it, and add its triangles to the list
  // for future conversion.
  int iCumCell = 0;
  for (iP = 0; iP < pB3D->iNumBdryPatches(); iP++) {
    BdryPatch3D *pBPatch = pB3D->pBdry(iP);
    // Can only handle polygons without negotiating curve discretization
    // in advance, so check for that.
    assert(pBPatch->eType() == BdryPatch3D::ePolygon);
    BdryPolygon *pBP = dynamic_cast<BdryPolygon*>(pBPatch);
    if (pBP->iNumSegments() == 3 && pBP->iNumChildren() == 0) {
      // We have a triangle with nothing inside it.  This is an easy case.
      Vert *pV0 = pVVert(pB3D->iVertIndex(pBP->pVVert(0)));
      Vert *pV1 = pVVert(pB3D->iVertIndex(pBP->pVVert(1)));
      Vert *pV2 = pVVert(pB3D->iVertIndex(pBP->pVVert(2)));
      TriCellCV *pCCV = VATriCV.pTNewEntry();
      if (pBPatch->iLeftRegion() == iOutsideRegion) {
	pCCV->vAssign(pV0, pV2, pV1);
      }
      else {
	pCCV->vAssign(pV0, pV1, pV2);
      }
#ifndef NDEBUG
      double adNormPatch[3], adNormVerts[3];
      vUnitNormal(pV0->adCoords(), pV1->adCoords(), pV2->adCoords(),
		  adNormVerts);
      // First arg is a dummy anyway for polygons...
      pBP->vUnitNormal(adNormVerts, adNormPatch);
      double dDot = dDOT3D(adNormPatch, adNormVerts);
      if (pBPatch->iLeftRegion() == iOutsideRegion)
	vSCALE3D(adNormVerts, -1);
      vMessage(4, "Cell: %3d  Verts: %3d, %3d, %3d.\n\t(%7f %7f %7f) Dot: %6f  Regs: %3d %3d\n",
	       iCumCell, iVertIndex(pV0), iVertIndex(pV1), iVertIndex(pV2),
	       adNormVerts[0], adNormVerts[1], adNormVerts[2],
	       dDot, pBPatch->iLeftRegion(), pBPatch->iRightRegion());
      assert((dDot > 0 && pBPatch->iRightRegion() == iOutsideRegion) ||
	     (dDot < 0 && pBPatch->iLeftRegion()  == iOutsideRegion) ||
	     (pBPatch->iRightRegion() != iOutsideRegion &&
	      pBPatch->iLeftRegion()  != iOutsideRegion));
#endif
      vSetBdryPatch(iCumCell, pBPatch);
      iCumCell++;
    }

    else {
      // Something that we actually have to triangulate...
      Mesh2D M2D(pBP->B2DProjection(*pB3D), 2, false);
      double adOffset[3], adBasis1[3], adBasis2[3];
      pBP->vFindProjectionBasis(adOffset, adBasis1, adBasis2);

      // Now the Mesh2D has the triangulation of the surface polygon, but
      // the indices of the verts do not match indices in the surface
      // mesh!  So we have to find the real vertices
      // So copy that data over.
      for (int iC = 0; iC < M2D.iNumCells(); iC++) {
	Cell *pC = M2D.pCCell(iC);
	Vert *pV2D0 = pC->pVVert(0);
	Vert *pV2D1 = pC->pVVert(1);
	Vert *pV2D2 = pC->pVVert(2);

	double adLoc3D0[3], adLoc3D1[3], adLoc3D2[3];
	adLoc3D0[0] = (+ pV2D0->dX()*adBasis1[0]
		       + pV2D0->dY()*adBasis2[0]
		       + adOffset[0]);
	adLoc3D0[1] = (+ pV2D0->dX()*adBasis1[1]
		       + pV2D0->dY()*adBasis2[1]
		       + adOffset[1]);
	adLoc3D0[2] = (+ pV2D0->dX()*adBasis1[2]
		       + pV2D0->dY()*adBasis2[2]
		       + adOffset[2]);

	adLoc3D1[0] = (+ pV2D1->dX()*adBasis1[0]
		       + pV2D1->dY()*adBasis2[0]
		       + adOffset[0]);
	adLoc3D1[1] = (+ pV2D1->dX()*adBasis1[1]
		       + pV2D1->dY()*adBasis2[1]
		       + adOffset[1]);
	adLoc3D1[2] = (+ pV2D1->dX()*adBasis1[2]
		       + pV2D1->dY()*adBasis2[2]
		       + adOffset[2]);

	adLoc3D2[0] = (+ pV2D2->dX()*adBasis1[0]
		       + pV2D2->dY()*adBasis2[0]
		       + adOffset[0]);
	adLoc3D2[1] = (+ pV2D2->dX()*adBasis1[1]
		       + pV2D2->dY()*adBasis2[1]
		       + adOffset[1]);
	adLoc3D2[2] = (+ pV2D2->dX()*adBasis1[2]
		       + pV2D2->dY()*adBasis2[2]
		       + adOffset[2]);

	// These points may only be very close approximations to
	// actually points in the patch if the patch is not quite
	// precisely coplanar.  So find the nearest point.

	int iLocalInd0 = pBP->iLocalVertIndex(adLoc3D0);
	int iLocalInd1 = pBP->iLocalVertIndex(adLoc3D1);
	int iLocalInd2 = pBP->iLocalVertIndex(adLoc3D2);

	assert(iLocalInd0 != -1 && iLocalInd1 != -1 && iLocalInd2 != -1);

	Vert *pV0 = pVVert(pB3D->iVertIndex(iP, iLocalInd0));
	Vert *pV1 = pVVert(pB3D->iVertIndex(iP, iLocalInd1));
	Vert *pV2 = pVVert(pB3D->iVertIndex(iP, iLocalInd2));

	TriCellCV *pCCV = VATriCV.pTNewEntry();
	if (pBPatch->iLeftRegion() != iOutsideRegion) {
	  pCCV->vAssign(pV0, pV2, pV1);
	}
	else {
	  pCCV->vAssign(pV0, pV1, pV2);
	}
	vSetBdryPatch(iCumCell, pBPatch);
	iCumCell++;
      }
    }
  }
  assert(iCumCell == iTotalCells);
  // Turn messaging back up a notch.
  iMessageStdoutLevel--;
  iMessageFileLevel--;

  // Convert from cell-vert to face-cell data structure and store here
  // in the SurfMesh.
  vConvertFromCellVert(VATriCV);
  vSetAllHintFaces();

  // Is the thing right side out?
  double dSize = dInteriorSize();
  if (dSize < 0) {
    vMessage(3, "Turning surface mesh right side out.\n");
    // Need to flip the surface orientation, both for the SurfMesh and
    // the Bdry3D.
    for (int iF = 0; iF < iNumFaces(); iF++) {
      pFFace(iF)->vInterchangeCells();
    }
    for (int iBP = 0; iBP < pB3D->iNumPatches(); iBP++) {
      pB3D->pBdry(iBP)->vSwitchOrientation();
      pB3D->pBdry(iBP)->vFlipNormal();
    }
    dSize = dInteriorSize();
  }
  assert(dSize > 0);
  vMessage(1, "Total volume enclosed in surface mesh: %10.5g\n",
	   dSize);

  // Now "fix" all small angles in the input.
  //
  vMarkSmallAngles();

  if (!(qValid ()))
    vFoundBug ("surface mesh creation from Bdry3D data");
  SUMAA_LOG_EVENT_END (INIT_SURF_TRI);
}

MultiEdge* SurfMesh::pMENewMultiEdge()
{
  MultiEdge *pME = VAME.pTNewEntry ();
  pME->vSetDefaultFlags ();
  return pME;
}

double SurfMesh::dBoundarySize() const
{
  double dRetVal = 0;
  for (int iC = iNumCells() - 1; iC >= 0; iC--)
    dRetVal += pCCell(iC)->dSize();
  return dRetVal;
}

double SurfMesh::dInteriorSize() const
{
  double dVolX = 0;
  // The checks for size consistency can fail when there are faces that
  // should be coplanar but just fail to be.  So don't do the checks.
// #ifndef NDEBUG
//   double dVolY = 0, dVolZ = 0;
// #endif
  double adNorm[3], adCent[3];
  for (int iC = iNumCells() - 1; iC >= 0; iC--) {
    if (pCCell(iC)->qDeleted()) continue;
    // If this cell of the SurfMesh lies on an interior patch, it
    // shouldn't contribute to the integrals for interior size.
    BdryPatch *pBP = pBdryPatch(iC);
    if (pBP->qIsBoundary()) {
      Cell *pC = pCCell(iC);
      pC->vCentroid(adCent);
      double dSize = pC->dSize();
      pBP->vUnitNormal(adCent, adNorm);
      if (pBP->iLeftRegion() != iOutsideRegion) {
	vSCALE3D(adNorm, -1);
      }

      dVolX += adCent[0]*adNorm[0]*dSize;
      // See comment above.
// #ifndef NDEBUG
//       dVolY += adCent[1]*adNorm[1]*dSize;
//       dVolZ += adCent[2]*adNorm[2]*dSize;
//       double adNorm2[3];
//       vNormal(pC->pVVert(0)->adCoords(), pC->pVVert(1)->adCoords(),
//	      pC->pVVert(2)->adCoords(), adNorm2);
//       vMessage(2, "Cell: %3d  Norm: (%7f %7f %7f)  Cent: (%7f %7f %7f)\n",
//	       iC, adNorm[0], adNorm[1], adNorm[2],
//	       adCent[0], adCent[1], adCent[2]);
//       vMessage(2, "   Scaled cell norm: (%7f %7f %7f)\n",
//	       adNorm2[0]/dSize, adNorm2[1]/dSize, adNorm2[2]/dSize);
// #endif
    }
  }
  // See comment above.
// #ifndef NDEBUG
//   assert(iFuzzyComp(dVolX, dVolY) == 0);
//   assert(iFuzzyComp(dVolX, dVolZ) == 0);
// #endif
  return dVolX;
}

bool SurfMesh::qWatertight() const
{
  double dAreaX = 0, dAreaY = 0, dAreaZ = 0;
  double adNorm[3];
  for (int iC = iNumCells() - 1; iC >= 0; iC--) {
    pCCell(iC)->vVecSize(adNorm);
    dAreaX += adNorm[0];
    dAreaY += adNorm[1];
    dAreaZ += adNorm[2];
  }
  double dCompSize = dBoundarySize();
  double dInvCompSize = 0.001/dCompSize;
  dAreaX *= dInvCompSize;
  dAreaY *= dInvCompSize;
  dAreaZ *= dInvCompSize;

  return ( (iFuzzyComp(dAreaX, 0.) == 0) &&
	   (iFuzzyComp(dAreaY, 0.) == 0) &&
	   (iFuzzyComp(dAreaZ, 0.) == 0) );
}

void SurfMesh::vMarkSmallAngles()
  // Find all small angles between bdry segments in the surface mesh and
  // check whether they could lead to infinite refinement because of
  // mutual encroachment.  Because of the possibility of uneven edge
  // splits, the edge lengths aren't considered as a factor here.
  // Instead, only the angle is used (it must be greater than 45 degrees
  // to prevent problems).  Vertices where we could get infinite
  // refinement (let's call this implosion) are marked as having a small
  // angle; this marking gets carried over into the volume mesh.  Once
  // the initial length scale is established, each edge incident on a
  // small angle apex is split at a small distance from the corner.
  // Later on, tets incident on the small angle corner are protected
  // from quality-based refinement.
{
  vMessage(2, "Indentifying small angles in the surface mesh.\n");
  int iNumToProtect = 0;

  for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
    // Find all incident edges.
    Vert *pV = pVVert(iV);
    if (!pV->pFHintFace()->qValid()) continue;
    std::set<Face*> spFNearby;
    {
      std::set<Vert*> spVNeigh;
      std::set<Cell*> spCTmp;
      std::set<BFace*> spBFTmp;
      bool qBdryVert;
      vNeighborhood(pV, spCTmp, spVNeigh, &spBFTmp, &qBdryVert,
		    &spFNearby);
    }

    // Eliminate faces that aren't incident on pV or don't lie along
    // bdry segments.  These are easily identifiable, because they have
    // incident cells from the same bdry patch.
    std::set<Face*>::iterator iterF;
    for (iterF = spFNearby.begin(); iterF != spFNearby.end(); ) {
      Face *pF = *iterF;
      if (!pF->qHasVert(pV)) {
	std::set<Face*>::iterator iterFTmp = iterF++;
	spFNearby.erase(iterFTmp);
	continue;
      }

      int iCLeft  = iCellIndex(pF->pCCellLeft() );
      int iCRight = iCellIndex(pF->pCCellRight());
      BdryPatch3D* pBPLeft  = pBdryPatch(iCLeft );
      BdryPatch3D* pBPRight = pBdryPatch(iCRight);
      if (pBPLeft == pBPRight) {
	std::set<Face*>::iterator iterFTmp = iterF++;
	spFNearby.erase(iterFTmp);
      }
      else iterF++;
    }

    // For each vert, find whether any pair of edges can implode.
    const int iNumEdges = spFNearby.size();
    Vert **apVNeigh = new Vert*[iNumEdges];
    double *adEdgeLengths = new double[iNumEdges];
    bool *aqImplodable = new bool[iNumEdges];
    bool qAnyImplodable = false;

    int iEdge = 0;
    for (iterF = spFNearby.begin(); iterF != spFNearby.end();
	 iterF++, iEdge++) {
      aqImplodable[iEdge] = false;
      Face *pF = *iterF;
      apVNeigh[iEdge] = (pF->pVVert(0) == pV) ? pF->pVVert(1) :
	pF->pVVert(0);
    }

    for (int iOuter = 0; iOuter < iNumEdges; iOuter++) {
      double adVecOut[] = adDIFF3D(apVNeigh[iOuter]->adCoords(),
				   pV->adCoords());
      double dLenOut = dMAG3D(adVecOut);
      adEdgeLengths[iOuter] = dLenOut;
      for (int iInner = iOuter+1; iInner < iNumEdges; iInner++) {
	double adVecIn [] = adDIFF3D(apVNeigh[iInner]->adCoords(),
				     pV->adCoords());
	double dLenIn  = dMAG3D(adVecIn );
	double dCosTheta = dDOT3D(adVecOut, adVecIn) / (dLenOut*dLenIn);
	// If the angle between the edges is less than 45 degrees, then
	// there is a chance that this pair could implode.
	if (dCosTheta >= M_SQRT1_2) {
	  double dRatio = dLenOut / dLenIn;
	  // Now adjust dRatio to be in the range [0.5,1].  This works
	  // whether dRatio > 1 or < 1.
	  double dK = floor(-log(dRatio) / M_LN2);
	  double dRatio_Star = dRatio * pow(2.,dK);
	  vMessage(4, " Cos(theta) is %g.  Unscaled ratio is %g; scaled is %g.  This is",
		   dCosTheta, dRatio, dRatio_Star);
	  // Using iFuzzyComp here to avoid possible roundoff weirdness
	  // with lengths that differ by powers of two.
	  assert(iFuzzyComp(dRatio_Star, 0.5) >= 0 &&
		 iFuzzyComp(dRatio_Star, 1) <= 0);
	  if (iFuzzyComp(dRatio_Star, 0.5/dCosTheta) >= 0 &&
	      iFuzzyComp(dRatio_Star, dCosTheta) <= 0) {
	    // So we're in the evil range.  This pair of edges could
	    // implode.
	    vMessage(4, " evil\n");
	    aqImplodable[iOuter] = aqImplodable[iInner] = true;
	    qAnyImplodable = true;
	  }
	  else {
	    vMessage(4, " okay\n");
	  }
	}
      } // Done with inner loop for edge-to-edge comparison
    } // Done with outer loop for edge-to-edge comparison

    if (qAnyImplodable) {
      pV->vSetSmallAngleCorner(true);
      iNumToProtect++;

      // Find the shortest and 2nd shortest edges in the bundle.  If the
      // shortest and 2nd shortest have a length ratio > 1.5, split all
      // edges at the shortest edge length from the common vert.  Otherwise,
      // split at half the shortest edge length.  This will maximize the
      // minimum edge length after insertion.
      double dShortestLen = LARGE_DBL, dSecondLen = LARGE_DBL;
//       int iShortest = -1;
//       int iSecond = -1;
      for (iEdge = 0; iEdge < iNumEdges; iEdge++) {
	if (aqImplodable[iEdge]) {
	  double dThisLen = adEdgeLengths[iEdge];
	  if (dThisLen < dShortestLen) {
//	    iSecond = iShortest;
	    dSecondLen = dShortestLen;
//	    iShortest = iEdge;
	    dShortestLen = dThisLen;
	  }
	  else if (dThisLen < dSecondLen) {
//	    iSecond = iEdge;
	    dSecondLen = dThisLen;
	  }
	}
      }

      assert(dShortestLen <= dSecondLen);
//       assert(iSecond != -1 && iShortest != -1);

      double dInsertionLen = dShortestLen / 2;
//       bool qSplitShortest = true;
      if (dSecondLen / dShortestLen > 1.5) {
	dInsertionLen = dShortestLen;
//	qSplitShortest = false;
      }
      pV->vSetLS(dInsertionLen);
    }
    delete [] aqImplodable;
    delete [] apVNeigh;
    delete [] adEdgeLengths;
  } // Done with this vertex
  vMessage(2, "Out of %d vertices, found %d that may need small angle protection.\n",
	   iNumVerts(), iNumToProtect);
}

// Following is the body of SurfMesh::vProtectSmallAngles, which was
// designed to split edges so that infinite encroachment could not
// occur.  That attempt to prevent implosion failed, but I'm not yet
// ready to write it off completely by eliminating the code.

void SurfMesh::vProtectSmallAngles()
 {
   vMessage(2, "Protecting small angles in the surface mesh.\n");

   bool qChanged = false;
   // Loop from the last initial vertex down, so that we don't end up
   // working on verts inserted during the process.  Mind you,
   // newly-inserted verts should be no problem, but this is more
   // efficient, regardless.
   do {
     int iNumChanges = 0;
     qChanged = false;
     for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
       // Find all incident edges.
       Vert *pV = pVVert(iV);
       if (!pV->pFHintFace()->qValid()) continue;
       std::set<Face*> spFNearby;
       {
	 std::set<Vert*> spVNeigh;
	 std::set<Cell*> spCTmp;
	 std::set<BFace*> spBFTmp;
	 bool qBdryVert;
	 vNeighborhood(pV, spCTmp, spVNeigh, &spBFTmp, &qBdryVert,
		       &spFNearby);
       }

       // Eliminate faces that aren't incident on pV or don't lie along
       // bdry segments.  These are easily identifiable, because they have
       // incident cells from the same bdry patch.
       std::set<Face*>::iterator iterF;
       for (iterF = spFNearby.begin(); iterF != spFNearby.end(); iterF++) {
	 Face *pF = *iterF;
	 if (!pF->qHasVert(pV)) {
	   std::set<Face*>::iterator iterFTmp = iterF--;
	   spFNearby.erase(iterFTmp);
	   continue;
	 }
	 int iCLeft  = iCellIndex(pF->pCCellLeft() );
	 int iCRight = iCellIndex(pF->pCCellRight());
	 BdryPatch3D* pBPLeft  = pBdryPatch(iCLeft );
	 BdryPatch3D* pBPRight = pBdryPatch(iCRight);
	 if (pBPLeft == pBPRight) {
	   std::set<Face*>::iterator iterFTmp = iterF--;
	   spFNearby.erase(iterFTmp);
	 }
       }

       // For each vert, find all pairs of edges that can implode.  For
       // now, assume that all such pairs are going to spiral down as a
       // single group.  This saves trying to seperate them into mutual
       // encroaching groups.
       const int iNumEdges = spFNearby.size();
       bool *aqImplodable = new bool[iNumEdges];
       bool qAnyImplodable = false;
       double *adEdgeLengths = new double[iNumEdges];
       Vert **apVNeigh = new Vert*[iNumEdges];

       int iEdge = 0;
       for (iterF = spFNearby.begin(); iterF != spFNearby.end();
	    iterF++, iEdge++) {
	 Face *pF = *iterF;
	 aqImplodable[iEdge] = false;
	 apVNeigh[iEdge] = (pF->pVVert(0) == pV) ? pF->pVVert(1) :
	   pF->pVVert(0);
       }

       // Outer loop goes all the way to iNumEdges-1 to set edge length for
       // [iNumEdges-1].
       for (int iOuter = 0; iOuter < iNumEdges; iOuter++) {
	double adVecOut[] = adDIFF3D(apVNeigh[iOuter]->adCoords(),
				     pV->adCoords());
	double dLenOut = dMAG3D(adVecOut);
	adEdgeLengths[iOuter] = dLenOut;
	for (int iInner = iOuter+1; iInner < iNumEdges; iInner++) {
	  double adVecIn [] = adDIFF3D(apVNeigh[iInner]->adCoords(),
				       pV->adCoords());
	  double dLenIn  = dMAG3D(adVecIn );
	  double dCosTheta = dDOT3D(adVecOut, adVecIn) / (dLenOut*dLenIn);
	  // If the angle between the edges exceeds 45 degrees, then there
	  // is no way that this pair can implode.
	  if (dCosTheta >= M_SQRT1_2) {
	    double dRatio = dLenOut / dLenIn;
	    // Now adjust dRatio to be in the range [0.5,1].  This works
	    // whether dRatio > 1 or < 1.
	    double dK = floor(-log(dRatio) / M_LN2);
	    double dRatio_Star = dRatio * pow(2.,dK);
	    vMessage(4, " Cos(theta) is %g.  Unscaled ratio is %g; scaled is %g.  This is",
		     dCosTheta, dRatio, dRatio_Star);
	    // Using iFuzzyComp here to avoid possible roundoff weirdness
	    // with lengths that differ by powers of two.
	    assert(iFuzzyComp(dRatio_Star, 0.5) >= 0 &&
		   iFuzzyComp(dRatio_Star, 1) <= 0);
	    if (iFuzzyComp(dRatio_Star, 0.5/dCosTheta) >= 0 &&
		iFuzzyComp(dRatio_Star, dCosTheta) <= 0) {
	      // So we're in the evil range.  This pair of edges could
	      // implode.
	      vMessage(4, " evil\n");
	      aqImplodable[iOuter] = aqImplodable[iInner] = true;
	      qAnyImplodable = true;
	    }
	    else {
	      vMessage(4, " okay\n");
	    }
	  }
	} // Done with inner loop for edge-to-edge comparison
       } // Done with outer loop for edge-to-edge comparison


       if (qAnyImplodable) {
	 pV->vSetSmallAngleCorner(true);
	// Find the shortest and 2nd shortest edges in the bundle.  If the
	// shortest and 2nd shortest have a length ratio > 1.5, split all
	// edges at the shortest edge length from the common vert.  Otherwise,
	// split at half the shortest edge length.  This will maximize the
	// minimum edge length after insertion.
	double dShortestLen = LARGE_DBL, dSecondLen = LARGE_DBL;
	int iShortest = -1;
//	int iSecond = -1;
	for (iEdge = 0; iEdge < iNumEdges; iEdge++) {
	  if (// true ||
	      aqImplodable[iEdge]) {
	    double dThisLen = adEdgeLengths[iEdge];
	    if (dThisLen < dShortestLen) {
//	      iSecond = iShortest;
	      dSecondLen = dShortestLen;
	      iShortest = iEdge;
	      dShortestLen = dThisLen;
	    }
	    else if (dThisLen < dSecondLen) {
//	      iSecond = iEdge;
	      dSecondLen = dThisLen;
	    }
	  }
	}

	assert(dShortestLen <= dSecondLen);
//	assert(iSecond != -1 && iShortest != -1);

	// Now split these edges.
	double dInsertionLen = dShortestLen / 2;
//	bool qSplitShortest = true;
	if (dSecondLen / dShortestLen > 1.5) {
	  dInsertionLen = dShortestLen;
//	  qSplitShortest = false;
	  vMessage(3, "Splitting at a distance of %f from the common vert.\n",
		   dInsertionLen);
	}

	for (iEdge = 0; iEdge < iNumEdges; iEdge++) {
	  if ((// true ||
	       aqImplodable[iEdge]) &&
	      ((iEdge != iShortest)
//	       || (qSplitShortest)
	       )) {
	    // Split this edge.
	    Vert *pVNeigh = apVNeigh[iEdge];
	    double adVec[] = adDIFF3D(pVNeigh->adCoords(),
				      pV->adCoords());
	    vNORMALIZE3D(adVec);
	    adVec[XDIR] = pV->dX() + adVec[XDIR] * dInsertionLen;
	    adVec[YDIR] = pV->dY() + adVec[YDIR] * dInsertionLen;
	    adVec[ZDIR] = pV->dZ() + adVec[ZDIR] * dInsertionLen;

	    // Since cell ID's are changing constantly with insertion, find a
	    // place to insert the new vertex the hard way.
	    Face *pFCand = pV->pFHintFace();
	    assert(pFCand->iFullCheck() && pFCand->qHasVert(pV));
	    Face *pFStart = pFCand;
	    Cell *pCCand = pFCand->pCCellLeft();
	    while (!pFCand->qHasVert(pVNeigh)) {
	      for (int i = 0; i < pCCand->iNumFaces(); i++) {
		Face *pF = pCCand->pFFace(i);
		if (pF != pFCand && pF->qHasVert(pV)) {
		  pFCand = pF;
		  break;
		}
	      }
	      pCCand = pFCand->pCCellOpposite(pCCand);
	      if (pFCand == pFStart) break; // A way out in case something
	      // goes wrong and the loop goes infinite.
	    }
	    assert(pFCand->qHasVert(pV));
	    assert(pFCand->qHasVert(pVNeigh));
	    Vert *pVNew = pVNewVert();
	    pVNew->vSetCoords(3, adVec);
	    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
	    iInsertOnFace(pVNew, pFCand, pCCand, true);
	    qChanged = true;
	    iNumChanges++;
	    pVNew->vSetType(Vert::eBdryCurve);
	    SUMAA_LOG_EVENT_END(INSERT_VERTEX);
	    vMessage(4, "Inserted between verts at (%12g, %12g, %12g)",
		     pV->dX(), pV->dY(), pV->dZ());
	    vMessage(4, "and (%12g, %12g, %12g)\n",
		     pVNeigh->dX(), pVNeigh->dY(), pVNeigh->dZ());
	    vMessage(4, "  New point located at (%12g, %12g, %12g)\n",
		     pVNew->dX(), pVNew->dY(), pVNew->dZ());
	  } // Done splitting an implodable edge
	} // Done checking all edges for splitting
       } // Done with if (any edges are implodable)
       delete [] aqImplodable;
       delete [] adEdgeLengths;
       delete [] apVNeigh;
     }
     vMessage(2, "Done with this pass. Made %d changes.\n", iNumChanges);
   } while (qChanged);
 }

static void
vBuildSortedEdgeList (const VATriCellCV & VATriCV,
		      const VAQuadCellCV & VAQuadCV,
		      const VAEdgeBFaceCV & VABFCV,
		      const int iNEdge,
		      FaceSort *&aFSEdges)
{
  if (aFSEdges != NULL) delete [] aFSEdges;
  aFSEdges = new FaceSort[iNEdge];
  int iEdge = 0, iCell = 0, i;
  // Add edges of triangular cells

  for (i = 0; i < VATriCV.iSize (); i++) {
    TriCellCV *pTCCV = VATriCV.pTEntry (i);
    aFSEdges[iEdge++] = FaceSort(pTCCV->pVVert (0), pTCCV->pVVert (1), iCell);
    aFSEdges[iEdge++] = FaceSort(pTCCV->pVVert (1), pTCCV->pVVert (2), iCell);
    aFSEdges[iEdge++] = FaceSort(pTCCV->pVVert (2), pTCCV->pVVert (0), iCell);
    iCell++;
  }

  for (i = 0; i < VAQuadCV.iSize (); i++) {
    QuadCellCV *pQCCV = VAQuadCV.pTEntry (i);
    aFSEdges[iEdge++] = FaceSort(pQCCV->pVVert (0), pQCCV->pVVert (1), iCell);
    aFSEdges[iEdge++] = FaceSort(pQCCV->pVVert (1), pQCCV->pVVert (2), iCell);
    aFSEdges[iEdge++] = FaceSort(pQCCV->pVVert (2), pQCCV->pVVert (3), iCell);
    aFSEdges[iEdge++] = FaceSort(pQCCV->pVVert (3), pQCCV->pVVert (0), iCell);
    iCell++;
  }

  for (i = 0; i < VABFCV.iSize (); i++) {
    EdgeBFaceCV *pEBFCV = VABFCV.pTEntry (i);
    assert (pEBFCV->iBdryCond () > 0);
    aFSEdges[iEdge++] = FaceSort(pEBFCV->pVVert (0), pEBFCV->pVVert (1),
				 -(pEBFCV->iBdryCond ()));
  }

  assert (iEdge == iNEdge);
  std::sort(aFSEdges, aFSEdges + iEdge);
}

//@@ This is a general conversion, for [mixed] 2D and surface meshes
// Meshes with only triangles will simply have no quads.
// Meshes with no boundary edges (surface meshes) have an empty VABFCV.
void SurfMesh::
vConvertFromCellVert (const VATriCellCV & VATriCV,
		      const VAQuadCellCV & VAQuadCV,
		      const VAEdgeBFaceCV & VABFCV)
{
  int iNTri = VATriCV.iSize ();
  int iNQuad = VAQuadCV.iSize ();
  int iNCell = iNTri + iNQuad;
  int iNBEdge = VABFCV.iSize ();
  assert(iNBEdge == 0); // Code for fixing orientation requires this.
  int iNEdge = 3 * iNTri + iNBEdge + 4 * iNQuad;

  // Set up enough blank entries in the connectivity tables to hold all
  // the data.  Don't set up edges yet, because for a non-manifold
  // surface mesh we can't tell in advance how many there will be.
  VATri.vSetup (iNTri);
  VAQuad.vSetup (iNQuad);
  VAEdgeBF.vSetup (iNBEdge);

  int i;
  // This code is needed if the SurfMesh already had data before.
  for (i = 0; i < iNTri; i++) {
    VATri.pTEntry(i)->vClearFaces();
  }
  for (i = 0; i < iNQuad; i++) {
    VAQuad.pTEntry(i)->vClearFaces();
  }
  for (i = 0; i < iNBEdge; i++) {
    VAEdgeBF.pTEntry(i)->vClearFaces();
  }

  // If there are no quads in the mesh, tag it as simplicial so that
  // swapping, insertion, etc, will be globally legal on the mesh.
  qSimplex = (iNQuad == 0);

  // This struct is used to hold edge descriptions for sorting.

  // Set up the individual edges and sort them
  FaceSort *aFSEdges = 0;

  vBuildSortedEdgeList(VATriCV, VAQuadCV, VABFCV, iNEdge, aFSEdges);

  int *aiOrientErrors = new int[iNCell];
  for (i = 0; i < iNCell; aiOrientErrors[i++] = 0) {}
  int iNumErrors = 0, iPasses = 0;

  // Attempt to fix orientation errors.

  std::set<BdryPatch3D*> sBPFlip;
  do {
    iNumErrors = 0;
    iPasses++;
    for (i = 0; i < iNCell; aiOrientErrors[i++] = 0) {}
    for (i = 0; i < iNEdge; ) {
      Vert *pV0 = aFSEdges[i].pVVert(0);
      Vert *pV1 = aFSEdges[i].pVVert(1);

      // How many faces are incident on this edge?
      int iNumInc = 1;
      while (i + iNumInc < iNEdge &&
	     qIdenticalVerts(aFSEdges[i], aFSEdges[i + iNumInc]))
	iNumInc++;
      assert(i + iNumInc <= iNEdge);

      if (iNumInc == 2) {
	if ((pV0 != aFSEdges[i + 1].pVVert(1)) ||
	    (pV1 != aFSEdges[i + 1].pVVert(0))) {
	  iNumErrors++;
	  aiOrientErrors[aFSEdges[i  ].iLeftCell()]++;
	  aiOrientErrors[aFSEdges[i+1].iLeftCell()]++;
	}
      }
      i += iNumInc;
    }

    // Now try to flip all the faces with more than one badly-oriented edge.
    if (iNumErrors > 0) {
      int aiErrorFreq[] = {0,0,0,0};
      for (i = 0; i < iNCell; i++) {
	aiErrorFreq[aiOrientErrors[i]]++;
      }
      vMessage(2, "Out of %d cells, %d had no orientation errors.\n",
	       iNCell, aiErrorFreq[0]);
      vMessage(2, "%d had one bad edge, %d had two, %d had three.\n",
	       aiErrorFreq[1], aiErrorFreq[2], aiErrorFreq[3]);
      vMessage(1, "Surface mesh has inconsistent orientation!\n");
      vMessage(1, "You're in real trouble here, but I'll try to fix it.\n");
      vMessage(2, "Flipping all cells with 2 or 3 bad edges.\n");
      int iNTriCV = VATriCV.iSize();
      for (i = 0; i < iNCell; i++) {
	if (aiOrientErrors[i] >= 2) {
	  if (i < iNTriCV) {
	    VATriCV.pTEntry(i)->vSwitchOrientation();
	    sBPFlip.insert(pBdryPatch(i));
	  }
	  else {
	    VAQuadCV.pTEntry(i-iNTriCV)->vSwitchOrientation();
	    sBPFlip.insert(pBdryPatch(i));
	  }
	}
      }
      vMessage(2, "Rebuilding edge list.\n");
      vBuildSortedEdgeList(VATriCV, VAQuadCV, VABFCV, iNEdge, aFSEdges);
    }
  } while (iNumErrors != 0 && iPasses < 25);
  // Now flip all the patches that had a cell flipped.
  std::set<BdryPatch3D*>::iterator iterBP;
  for (iterBP = sBPFlip.begin(); iterBP != sBPFlip.end(); iterBP++) {
    (*iterBP)->vSwitchOrientation();
    // Flipping normals here fixes tet2 at the expense of breaking a
    // bunch of other cases.
  }

  // Read the data off
  int iBFace = 0, iCell;
  bool qOrientationOK = true;
  for (i = 0; i < iNEdge; ) {
    Vert *pV0 = aFSEdges[i].pVVert(0);
    Vert *pV1 = aFSEdges[i].pVVert(1);

    // How many faces are incident on this edge?
    int iNumInc = 1;
    while (i + iNumInc < iNEdge &&
	   qIdenticalVerts(aFSEdges[i], aFSEdges[i + iNumInc]))
      iNumInc++;
    assert(i + iNumInc <= iNEdge);

    switch (iNumInc) {
    case 0:
    case 1:
      // These cases are impossible.
      vMessage(0, "Warning: Orphan edge (no mate) along edge\n");
      switch (pV0->iSpaceDimen()) {
      case 2:
	vMessage(0, "         between vert %5d (%10g, %10g)\n",
		 iVertIndex(pV0), pV0->dX(), pV0->dY());
	vMessage(0, "         and          %5d (%10g, %10g)\n",
		 iVertIndex(pV1), pV1->dX(), pV1->dY());
	break;
      case 3:
	vMessage(0, "         between vert %5d (%10g, %10g, %10g)\n",
		 iVertIndex(pV0), pV0->dX(), pV0->dY(), pV0->dZ());
	vMessage(0, "         and          %5d (%10g, %10g, %10g)\n",
		 iVertIndex(pV1), pV1->dX(), pV1->dY(), pV1->dZ());
	break;
      default:
	assert(0);
      }
      break;
    case 2:
      // The standard case.  This always happens in 2D, and for surface
      // meshes only fails at edges where the mesh is locally non-manifold.
      {
	if ((pV0 != aFSEdges[i + 1].pVVert(1)) ||
	    (pV1 != aFSEdges[i + 1].pVVert(0))) {
	  qOrientationOK = false;
	  iNumErrors++;
	  aiOrientErrors[aFSEdges[i  ].iLeftCell()]++;
	  aiOrientErrors[aFSEdges[i+1].iLeftCell()]++;
	  vMessage(2, "Warning: Inconsistent orientation of input polygons along edge\n");
	  switch (pV0->iSpaceDimen()) {
	  case 2:
	    vMessage(2, "         between vert %5d (%10g, %10g)\n",
		     iVertIndex(pV0), pV0->dX(), pV0->dY());
	    vMessage(2, "         and          %5d (%10g, %10g)\n",
		     iVertIndex(pV1), pV1->dX(), pV1->dY());
	    break;
	  case 3:
	    vMessage(2, "         between vert %5d (%10g, %10g, %10g)\n",
		     iVertIndex(pV0), pV0->dX(), pV0->dY(), pV0->dZ());
	    vMessage(2, "         and          %5d (%10g, %10g, %10g)\n",
		     iVertIndex(pV1), pV1->dX(), pV1->dY(), pV1->dZ());
	    break;
	  default:
	    assert(0);
	  }
	}
	if (qOrientationOK) {
	  // Only do this if orientation is good, to prevent potential
	  // problems with too many edges per cell or vice versa
	  int iC0 = aFSEdges[i].iLeftCell();
	  int iC1 = aFSEdges[i + 1].iLeftCell();
	  Cell *pC0, *pC1;
	  if (iC0 >= 0)
	    pC0 = pCCell (iC0);
	  else {
	    BFace *pBF = pBFBFace (iBFace++);
	    pC0 = pBF;
	  }
	  if (iC1 >= 0)
	    pC1 = pCCell (iC1);
	  else {
	    BFace *pBF = pBFBFace (iBFace++);
	    pC1 = pBF;
	  }

	  Face *pF = pFNewFace();

	  pF->vAssign (pC0, pC1, pV0, pV1);
	  pC0->vAddFace (pF);
	  pC1->vAddFace (pF);
	}
	break;
      }
    default:
      // Handle non-manifold cases.
      vMessage(2, "Non-manifold surface along edge\n");
      vMessage(2, "         between vert %5d (%10g, %10g, %10g)\n",
	       iVertIndex(pV0), pV0->dX(), pV0->dY(), pV0->dZ());
      vMessage(2, "         and          %5d (%10g, %10g, %10g)\n",
	       iVertIndex(pV1), pV1->dX(), pV1->dY(), pV1->dZ());
      // Make an array with all the cells incident on this face.
      Cell **apCInc = new Cell*[iNumInc];
      MultiEdge *pME = pMENewMultiEdge();
      for (int ii = 0; ii < iNumInc; ii++) {
	iCell = aFSEdges[i + ii].iLeftCell();
	if (iCell < 0)
	  vFatalError("Illegal boundary connectivity (hanging edge)\n",
		      "SurfMesh::vConvertFromCellVert(...)");
	apCInc[ii] = pCCell(iCell);
	apCInc[ii]->vAddFace(pME);
      }
      pME->vAssign(pV0, pV1, apCInc, iNumInc);
      delete [] apCInc;
      break;
    }
    i += iNumInc;
  }
  // Done with this, so nuke it.
  delete [] aFSEdges;

  // Check whether all faces have the same orientation that they're
  // supposed to from the CellCV's.  Those that don't will have this
  // problem because of MultiEdge's.
  for (iCell = 0; iCell < iNCell; iCell++) {
    Vert *pV0CV = VATriCV[iCell].pVVert(0);
    Vert *pV1CV = VATriCV[iCell].pVVert(1);
    Vert *pV2CV = VATriCV[iCell].pVVert(2);

    Cell *pC = pCCell(iCell);
    Vert *pV0 = pC->pVVert(0);
    Vert *pV1 = pC->pVVert(1);
    Vert *pV2 = pC->pVVert(2);

    if (!( (pV0CV == pV0 && pV1CV == pV1 && pV2CV == pV2) ||
	   (pV0CV == pV1 && pV1CV == pV2 && pV2CV == pV0) ||
	   (pV0CV == pV2 && pV1CV == pV0 && pV2CV == pV1) )) {
      Face *pF0 = pC->pFFace(0);
      Face *pF1 = pC->pFFace(1);
      Face *pF2 = pC->pFFace(2);
      assert(pF0->eType() == Face::eMultiEdge);
      if (pF1->eType() == Face::eEdgeFace) {
	pC->vAssign(pF1, pF2, pF0);
      }
      else if (pF2->eType() == Face::eEdgeFace) {
	pC->vAssign(pF2, pF0, pF1);
      }
      else {
	vMessage(2, "%s %s\n",
		 "You have a face bounded by three multi-edges;",
		 "you're probably in trouble here.");
	vMessage(2, "  CellCV says verts are: %d %d %d\n",
		 iVertIndex(pV0CV), iVertIndex(pV1CV), iVertIndex(pV2CV));
	vMessage(2, "  SurfMesh says verts are: %d %d %d\n",
		 iVertIndex(pV0), iVertIndex(pV1), iVertIndex(pV2));
      }
    }
  }

#ifdef NDEBUG
  for (iCell = 0; iCell < iNCell; iCell++) {
    Vert *pV0CV = VATriCV[iCell].pVVert(0);
    Vert *pV1CV = VATriCV[iCell].pVVert(1);
    Vert *pV2CV = VATriCV[iCell].pVVert(2);

    Cell *pC = pCCell(iCell);
    Vert *pV0 = pC->pVVert(0);
    Vert *pV1 = pC->pVVert(1);
    Vert *pV2 = pC->pVVert(2);

    if (!( (pV0CV == pV0 && pV1CV == pV1 && pV2CV == pV2) ||
	   (pV0CV == pV1 && pV1CV == pV2 && pV2CV == pV0) ||
	   (pV0CV == pV2 && pV1CV == pV0 && pV2CV == pV1) )) {
      vMessage(2, "Possible problem with triangle orientation: \n");
      vMessage(2, "  CellCV says verts are: %d %d %d\n",
	       iVertIndex(pV0CV), iVertIndex(pV1CV), iVertIndex(pV2CV));
      vMessage(2, "  SurfMesh says verts are: %d %d %d\n",
	       iVertIndex(pV0), iVertIndex(pV1), iVertIndex(pV2));
    }
  }
#endif

  if (!qOrientationOK) {
    int aiErrorFreq[] = {0,0,0,0};
    for (i = 0; i < iNCell; i++) {
      aiErrorFreq[aiOrientErrors[i]]++;
    }
    vMessage(2, "Out of %d cells, %d had no orientation errors.\n",
	     iNCell, aiErrorFreq[0]);
    vMessage(2, "%d had one bad edge, %d had two, %d had three.\n",
	     aiErrorFreq[1], aiErrorFreq[2], aiErrorFreq[3]);
//      if (aiErrorFreq[2] + aiErrorFreq[3] > 0 &&
//	aiErrorFreq[2] != iNCell) {
//        // Re-orient anything with 2 or 3 bad edges.
//        vMessage
//	(1, "Trying heuristics to fix surface triangle orientation...\n");
//        vMessage(2, "Flipping all cells with 2 or 3 bad edges.\n");
//        int iNTriCV = VATriCV.iSize();
//        for (i = 0; i < iNCell; i++) {
//	if (aiOrientErrors[i] >= 2) {
//	  if (i < iNTriCV) {
//	    VATriCV.pTEntry(i)->vSwitchOrientation();
//	  }
//	  else {
//	    VAQuadCV.pTEntry(i-iNTriCV)->vSwitchOrientation();
//	  }
//	}
//        }
//        vMessage(2, "Trying conversion again.\n");
//        vConvertFromCellVert(VATriCV, VAQuadCV, VABFCV);
//      }
//      else {
      vFatalError ("Polygons do not all have the same orientation",
		   "SurfMesh::vConvertFromCellVert(...)");
//      }
  }
  delete [] aiOrientErrors;
  assert (iBFace == iNBEdge);
}

void SurfMesh::
vPurgeVerts ()
  // Renumber vertices, removing those that have been marked as deleted
{
  bool qLocalLog = false;
#ifdef SUMAA_LOG
  if (SUMAA_LOG_current != PURGE) {
    qLocalLog = true;
    SUMAA_LOG_EVENT_BEGIN(PURGE);
  }
#endif
  int *aiVertInd = new int[iNumVerts ()];
  int iV, iVActive;

  vSetAllHintFaces ();
  iVActive = 0;
  // Copy the vertex data down within the list and establish a lookup so
  // that the faces can figure out where their verts are now stored.
  for (iV = 0; iV < iNumVerts (); iV++) {
    Vert *pV = pVVert (iV);
    if (!pV->qDeleted ()) {
      aiVertInd[iV] = iVActive;
      Vert *pVA = pVVert (iVActive);
      // Copy the data from its old to its new location
      (*pVA) = (*pV);
      iVActive++;
    }
    else
      aiVertInd[iV] = -1;
  }

  // Note: Using vReplaceVert in the following loop fails because a face
  // can have two pointers to the same place, one before and one after
  // vertex list compression.  This would end up flipping the
  // orientation of the face, which is of course a no-no.
  for (int iF = iNumFaces () - 1; iF >= 0; iF--) {
    Face *pF = pFFace (iF);
    if (!pF->iFullCheck ())
      continue;
    Vert *pVNew0 = pVVert (aiVertInd[iVertIndex (pF->pVVert (0))]);
    Vert *pVNew1 = pVVert (aiVertInd[iVertIndex (pF->pVVert (1))]);

    // Passing the extra arguments does no harm, as they are ignored.
    (dynamic_cast<EdgeFace*>(pF))->vAssignVerts (pVNew0, pVNew1);
  }

  for (int iME = iNumMultiEdges () - 1; iME >= 0; iME--) {
    MultiEdge *pME = pMEMultiEdge (iME);
    if (!pME->iFullCheck ())
      continue;

    Vert *pVNew0 = pVVert (aiVertInd[iVertIndex (pME->pVVert (0))]);
    Vert *pVNew1 = pVVert (aiVertInd[iVertIndex (pME->pVVert (1))]);

    // Passing the extra arguments does no harm, as they are ignored.
    (dynamic_cast<EdgeFace*>(pME))->vAssignVerts (pVNew0, pVNew1);
  }

  delete[]aiVertInd;

  vSetupVerts (iVActive);
  vSetAllHintFaces ();
  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

SurfMesh::SurfMesh(const char strBaseFileName[], const int iQualMeas)
  : Mesh2D(iQualMeas), VAME(), iOnMultiEdge(0), pB3D(NULL), apBPatch(NULL),
    iBPatchSize(0)
{
  SUMAA_LOG_EVENT_BEGIN (INIT_TRI);
  SUMAA_LOG_EVENT_BEGIN (READ_MESH_2D);
  qLengthScaleFromCellSizes = true;

  //@@ Read data from a user-defined format input file and convert to Mesh2D
  int iNVerts, iNFaces, iNCells, iNIntBFaces;
  bool qFaceVert, qCellVert, qFaceCell, qCellFace, qCellRegion;
  bool qIntBFaceFace, qIntBFaceVert, qIntBFaceBC;
  int (*a2iFaceVert)[2], (*a2iFaceCell)[2];
  int (*a2iCellVert)[3], (*a2iCellFace)[3];
  int *aiCellRegion;
  int (*a2iIntBFaceFace)[2], *aiIntBFaceBC, (*a2iIntBFaceVert)[2];

  vReadFile_Surface(strBaseFileName, iNVerts, iNFaces, iNCells,
		    iNIntBFaces, qFaceVert, qCellVert,
		    qFaceCell, qCellFace, qCellRegion, qIntBFaceFace,
		    qIntBFaceVert, qIntBFaceBC, VAVerts,
		    a2iFaceVert, a2iFaceCell, a2iCellVert,
		    a2iCellFace, aiCellRegion, a2iIntBFaceFace,
		    aiIntBFaceBC, a2iIntBFaceVert);

  if (VAVerts.iSize () == 0)
    vFatalError ("No vertex coordinate data read", "surface mesh input");

  bool qConnectOK = false;
  if (qFaceVert && qFaceCell)
    qConnectOK = true;	// face-based structure

  else if (qCellVert)
    qConnectOK = true;	// FE structure

  if (!qConnectOK)
    vFatalError ("Inadequate connectivity data read; consult user's manual.",
		 "surface mesh input");

  //@@@ Convert face-based data to Mesh2D
  if (qFaceVert && qFaceCell) {
    // If the number of cells, faces, and boundary faces have been
    // provided, we can check for mesh validity (and the presence of
    // quads in the mesh; quads are not yet allowed).
    if (iNFaces == 0)
      vFatalError ("Mesh size unknown.", "surface mesh input");
    //@@@@ Verify mesh size and that numbering starts from 0
    {
      int iTmpNCells = 0, iTmpNBFaces = 0, iTmpNVerts = 0;
      for (int iF = 0; iF < iNFaces; iF++) {
	int iC0 = a2iFaceCell[iF][0];
	int iC1 = a2iFaceCell[iF][1];
#ifndef MAX3
#define MAX3(a,b,c) MAX((a), MAX((b),(c)))
#endif
	iTmpNCells = MAX3 (iTmpNCells, iC0, iC1);
	if (iC0 < 0 || iC1 < 0)
	  iTmpNBFaces++;
	int iV0 = a2iFaceVert[iF][0];
	int iV1 = a2iFaceVert[iF][1];
	iTmpNVerts = MAX3 (iTmpNVerts, iV0, iV1);
#undef MAX3
      }
      if (iNCells != 0 && iNCells != iTmpNCells + 1)
	vFatalError ("Number of cells in mesh doesn't match number given.",
		     "surface mesh input");
      if ((iNVerts != 0 && iNVerts != iTmpNVerts + 1)
	  || iTmpNVerts + 1 != VAVerts.iSize())
	vFatalError ("Number of verts in mesh doesn't match number given.",
		     "surface mesh input");
      iNCells = iTmpNCells + 1;
      iNVerts = iTmpNVerts + 1;
    }

    if (3 * iNCells != 2 * iNFaces)
      vFatalError ("Mesh is invalid or contains quadrilaterals.",
		   "surface mesh input");

    vSetupFaces (iNFaces);	// Ensure the existence of enough
    vSetupTriCells (iNCells);	// blank faces, bdry faces, and cells.

      //@@@ Assign cells to regions
      // The default region, set when the cell was created, is 1.
    if (qCellRegion) {
      assert(aiCellRegion != NULL);
      for (int iCell = 0; iCell < iNumCells (); iCell++)
	pCCell (iCell)->vSetRegion (aiCellRegion[iCell]);
      // Done with aiCellRegion, so get rid of it.  Not allocated unless
      // the cell region was actually read.
      free (aiCellRegion);
    }

    // Need to map vertex indices to point indices in the Bdry2D
    // struct, since Bdry2D no longer uses vertices directly.  This
    // mapping is stored in the following array.
    int *aiVertMap = new int[iNVerts];
    int iV;
    for (iV = 0; iV < iNVerts; aiVertMap[iV++] = -1) {}

    for (int i = 0; i < iNFaces; i++) {
      int iV0 = a2iFaceVert[i][0];
      int iV1 = a2iFaceVert[i][1];

      // Add entry to cell-face table or bdry face list on the fly
      int iC0 = a2iFaceCell[i][0];
      int iC1 = a2iFaceCell[i][1];

      Cell *pC0 = pCCell (iC0);
      Cell *pC1 = pCCell (iC1);

      // This needs to be done whether pC? is a real cell or a bdry face
      Face *pF = pFFace(i);
      pC0->vAddFace (pF);
      pC1->vAddFace (pF);

      // Tell the face everything it needs to know.
      pF->vAssign (pC0, pC1, pVVert (iV0), pVVert (iV1));
    } // Done assigning for this face.
    delete [] aiVertMap;
    vSetAllHintFaces ();
    if (!(qValid ()))
      vFoundBug ("surface mesh input from user-format face->cell-style file");
  }

  else if (qCellVert) {
    if (iNCells == 0)
      vFatalError ("Mesh size unknown.", "surface mesh input");
    //@@@@ Verify that number of verts is correct
    {
      int iTmpNVerts = 0;
      for (int iC = 0; iC < iNCells; iC++) {
	int iV0 = a2iCellVert[iC][0];
	int iV1 = a2iCellVert[iC][1];
	int iV2 = a2iCellVert[iC][2];
	iTmpNVerts = max(max(iTmpNVerts, iV0), max(iV1, iV2));
      }
      if (iNVerts != 0 && iNVerts != iTmpNVerts + 1)
	vFatalError ("Number of verts in mesh doesn't match number given.",
		     "surface mesh input");
    }
    // Store triangle data here
    VATriCellCV VACV;
    // If there were quad data, store it here
    VAQuadCellCV VAQuadCV;
    // Store boundary edge data here
    //      VABdryEdgeCV VABdryCV;
    VAEdgeBFaceCV VABdryCV;
    // Set up data for all the interior triangles
    for (int iC = 0; iC < iNCells; iC++) {
      Vert *pV0 = pVVert (a2iCellVert[iC][0]);
      Vert *pV1 = pVVert (a2iCellVert[iC][1]);
      Vert *pV2 = pVVert (a2iCellVert[iC][2]);
      TriCellCV *pC = VACV.pTNewEntry ();
      pC->vAssign (pV0, pV1, pV2);
    }
    vMessage (1, "Converting triangle list to edge-cell format.\n");
    //@@@@ Convert the triangles to edge-cell data
    vConvertFromCellVert (VACV, VAQuadCV, VABdryCV);
    vSetAllHintFaces ();
    if (!(qValid ()))
      vFoundBug ("surface mesh input from user-format FE-style file");
  } // Done translating cell-vert/bface-vert/bface-BC data

  //@@@ Tag faces and verts
  for (int iVert = 0; iVert < iNumVerts (); iVert++)
    pVVert (iVert)->vSetType (Vert::eInterior);
  for (int iFace = 0; iFace < iNumFaces (); iFace++) {
    Face *pF = pFFace (iFace);
    Cell *pCL = pF->pCCellLeft ();
    Cell *pCR = pF->pCCellRight ();
    if (pCL->eType () == Cell::eBdryEdge ||
	pCR->eType () == Cell::eBdryEdge)
      continue; // Do boundaries separately.
    // Mark internal boundaries

    if (pCL->iRegion () != pCR->iRegion ()) {
      pF->vSetFaceLoc (Face::eBdryTwoSide);
      pF->pVVert (0)->vSetType (Vert::eBdryTwoSide);
      pF->pVVert (1)->vSetType (Vert::eBdryTwoSide);
    }
    else // Not an internal boundary
      pF->vSetFaceLoc (Face::eInterior);
  }

  // Tag ordinary boundary verts, and tag verts that lie at the
  // intersection of internal and regular boundary faces as boundary
  // apexes, so that they can never be smoothed or removed.
  for (int iBFace = 0; iBFace < iNumBdryFaces (); iBFace++) {
    BFace *pBF = pBFBFace (iBFace);
    pBF->pFFace ()->vSetFaceLoc (Face::eBdryFace);
    for (int iV = 0; iV < pBF->iNumVerts (); iV++) {
      Vert *pV = pBF->pVVert (iV);
      switch (pV->iVertType ()) {
      case Vert::eBdryTwoSide:
	pV->vSetType (Vert::eBdryApex);
	break;
      case Vert::eInterior:
	pV->vSetType (Vert::eBdry);
	break;
      default:
	// Nothing needs to be done
	break;
      }
    }
  }

  qLengthScaleFromCellSizes = true;

  if (qFaceVert)
    free (a2iFaceVert);
  if (qCellVert)
    free (a2iCellVert);
  if (qFaceCell)
    free (a2iFaceCell);
  if (qCellFace)
    free (a2iCellFace);
  SUMAA_LOG_EVENT_END (READ_MESH_2D);
  SUMAA_LOG_EVENT_END (INIT_TRI);
}

Vert* SurfMesh::createVert(const double dX, const double dY, const double dZ)
{
  Vert *pV = pVNewVert();
  pV->vSetDefaultFlags();
  double adData[] = {dX, dY, dZ};
  pV->vSetCoords(3, adData);
  return pV;
}

