#include <stdlib.h>
#include "GR_assert.h"
#include "GR_misc.h"
#include "GR_events.h"
#include "GR_Mesh.h"
#include "GR_Vec.h"
#include "GR_VertConnect.h"

#define STRATEGY 2
#define REV_LEX 1
#undef ANISO_REINSERT

#include <set>
#include <algorithm>
using std::set;
using std::set_difference;

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

#ifdef ANISO_REINSERT
#include <vector>
using std::vector;

typedef struct {
  double adBdryLoc[3];
  vector<double> vecdXCoords, vecdYCoords, vecdZCoords;
  double dX(const unsigned int i) {
    assert(i >=0 && i < vecdXCoords.size());
    return vecdXCoords[i];
  }
  double dY(const unsigned int i) {
    assert(i >=0 && i < vecdYCoords.size());
    return vecdYCoords[i];
  }
  double dZ(const unsigned int i) {
    assert(i >=0 && i < vecdZCoords.size());
    return vecdZCoords[i];
  }
  void vAppend(const double adCoords[]) {
    vecdXCoords.push_back(adCoords[0]);
    vecdYCoords.push_back(adCoords[1]);
    vecdZCoords.push_back(adCoords[2]);
  }
  int iLength() {
    assert(vecdXCoords.size() == vecdYCoords.size());
    assert(vecdXCoords.size() == vecdZCoords.size());
    return vecdXCoords.size();
  }
} Chain;

static vector<Chain> vecCReinsert;
static Chain ChainNull;
#endif

int Mesh::iMarkChain(Vert * const pVSurf, Vert * const pVFirst,
		     const VertConnect aVC[],
		     const double dMaxDist, const int iSkipNum)
{
  pVSurf->vMarkToKeep(); // Be sure -not- to delete this vertex.
  bool qChainOK = true;
  int iRetVal = 0;
  int iCurr = iVertIndex(pVFirst);
  Vert *pVLastKept = pVSurf;
  Vert *pVCurr = pVFirst;
  Vert *pVPrev = pVSurf;
  Vert *pVNext = pVInvalidVert;
  Vert *apVSkipped[8]; // Store the intermediate skipped verts, so that
		       // those after the last kept vertex can be restored.

#ifdef ANISO_REINSERT
  // Set up a new re-insertion chain and initialize it with its starting
  // point.
  vecCReinsert.push_back(ChainNull);
  vecCReinsert.back().adBdryLoc[0] = pVSurf->adCoords()[0];
  vecCReinsert.back().adBdryLoc[1] = pVSurf->adCoords()[1];
  vecCReinsert.back().adBdryLoc[2] = pVSurf->adCoords()[2];
#endif /* ANISO_REINSERT */

  if (pVSurf->iSpaceDimen() == 2)
    vMessage(4, "Chain starting from location (%10g, %10g), type %d\n",
	     pVSurf->dX(), pVSurf->dY(), pVSurf->iVertType());
  else
    vMessage(4, "Chain starting from location (%10g, %10g, %10g), type %d\n",
	     pVSurf->dX(), pVSurf->dY(), pVSurf->dZ(), pVSurf->iVertType());

  while (1) { // Exit from this loop is via the return, below.
    int iSkipped = 0;
    while (iSkipped < iSkipNum && qChainOK) {
      assert(!pVPrev->qActive());
      assert(pVCurr->qActive());
      qChainOK = aVC[iCurr].qNormalNeighborExists(pVCurr, pVPrev, pVNext);
      pVCurr->vMarkStructured();
      // Make sure nothing happens to this vert later.
      pVCurr->vMarkInactive();
      apVSkipped[iSkipped] = pVCurr;
      if (pVCurr->iSpaceDimen() == 2)
	vMessage(4, "Chain skipping location (%10g, %10g), type %d\n",
		 pVCurr->dX(), pVCurr->dY(), pVCurr->iVertType());
      else
	vMessage(4, "Chain skipping location (%10g, %10g, %10g), type %d\n",
		 pVCurr->dX(), pVCurr->dY(), pVCurr->dZ(),
		 pVCurr->iVertType());

      if (!pVNext->qValid() || !pVNext->qActive()) qChainOK = false;
      if (qChainOK) {
	pVPrev = pVCurr;
	pVCurr = pVNext;
	iCurr = iVertIndex(pVCurr);
	iSkipped++;
      }
    }
    if (!qChainOK ||
	(dDistanceBetween(pVLastKept, pVCurr) > dMaxDist)) {
      vMessage(4, "Chain done\n");
      // Set length scale for last vertex and tag all its neighbors to
      // have their LS checked.  This prevents undesirable clustering
      // near the ends of chains.
      pVCurr->vSetLS(dMaxDist);
      std::set<Vert*> spVNeigh;
      std::set<Cell*> spCTmp;
      vNeighborhood(pVCurr, spCTmp, spVNeigh);
      spVUpdateLS.insert(spVNeigh.begin(), spVNeigh.end());

      // Restore the last batch of skipped verts to a state where they
      // can be included in the MIS for this dimension.
      for (int ii = iSkipped - 1; ii > (iSkipNum+1) / 2; ii++) {
	apVSkipped[ii]->vMarkActive();
      }
      return iRetVal;
    }
    // For chains lieing -on- the boundary, vertices are kept.  Chains
    // in the interior are deleted and re-inserted afterwards.
#ifdef ANISO_REINSERT
    if (pVCurr->qIsBdryVert())
      pVCurr->vMarkToKeep(); // Keep this vertex
    else {
      pVCurr->vMarkToDelete(); // Redundant
      vecCReinsert.back().vAppend(pVCurr->adCoords());
    }
#else
    pVCurr->vMarkToKeep(); // Keep this vertex
#endif

    pVLastKept = pVCurr;
    iRetVal ++;
    if (pVCurr->iSpaceDimen() == 2)
      vMessage(4, "Chain keeping location (%10g, %10g), type %d\n",
	       pVCurr->dX(), pVCurr->dY(), pVCurr->iVertType());
    else
      vMessage(4, "Chain keeping location (%10g, %10g, %10g), type %d\n",
	       pVCurr->dX(), pVCurr->dY(), pVCurr->dZ(),
	       pVCurr->iVertType());

    pVCurr->vMarkStructured();
    // Make sure nothing happens to this vert later.
    pVCurr->vMarkInactive();

    // Need to be pointing at the first vert to skip before returning to
    // the top of the loop.
    qChainOK = aVC[iCurr].qNormalNeighborExists(pVCurr, pVPrev, pVNext);
    if (!pVNext->qValid() || !pVNext->qActive()) qChainOK = false;
    if (qChainOK) {
      pVPrev = pVCurr;
      pVCurr = pVNext;
      iCurr = iVertIndex(pVCurr);
    }
    else {
      vMessage(4, "Chain done\n");
      // Set length scale for last vertex and tag all its neighbors to
      // have their LS checked.
      pVCurr->vSetLS(dMaxDist);
      std::set<Vert*> spVNeigh;
      std::set<Cell*> spCTmp;
      vNeighborhood(pVCurr, spCTmp, spVNeigh);
      spVUpdateLS.insert(spVNeigh.begin(), spVNeigh.end());

      // Restore the last batch of skipped verts to a state where they
      // can be included in the MIS for this dimension.
      for (int ii = iSkipped - 1; ii > (iSkipNum+1) / 2; ii++) {
	apVSkipped[ii]->vMarkActive();
      }
      return iRetVal;
    }
  }
}

void Mesh::vDeactivateChain(Vert * const pVSurf, Vert * const pVFirst,
			    const VertConnect aVC[],
			    const double dMaxDist) const
{
  bool qChainOK = true;
  int iCurr = iVertIndex(pVFirst);
  Vert *pVCurr = pVFirst;
  Vert *pVPrev = pVSurf;
  Vert *pVNext;

  do {
    assert(!pVPrev->qActive());
    pVCurr->vMarkStructured();
    pVCurr->vMarkInactive();
    qChainOK = aVC[iCurr].qNormalNeighborExists(pVCurr, pVPrev, pVNext);
    if (qChainOK) {
      pVPrev = pVCurr;
      pVCurr = pVNext;
      iCurr = iVertIndex(pVCurr);
    }
  } while (pVCurr->qActive() && qChainOK &&
	   dDistanceBetween(pVCurr, pVPrev) < dMaxDist);
}

void Mesh::vAddNormalVerts(const VertConnect aVC[],
			   const Vert::VertType VTVertTypeToStartFrom)
{
  int iV;
  int iAdded = 0;
  static const double dDistRatioLimit = 0.1;

  // Coarsen by choosing every iNSkip+1 vertex along valid lines.
  for (iV = 0; iV < iNumVerts(); ++ iV) {
    Vert *pV = pVVert(iV);
    // This vertex can be marched from only if (a) it is a sufficiently
    // important vertex (part of a sufficiently low dimension feature of
    // the mesh) and (b) the part of the mesh being marched into really
    // is anisotropic; this rule is actually a bit more complicated for
    // verts on BdryCurves.

    // Is this vertex part of a low-dimensional enough feature?
    if (pV->iVertType() <= VTVertTypeToStartFrom) {
      // Yes, it is.  Now find the distance to the nearest vertex that
      // is also low-D enough.  This is the long axis of the anisotropic
      // mesh, locally.
      double dLongDist =
	aVC[iV].dClosestLegitSurfPoint(pV, VTVertTypeToStartFrom);

      int iNSkip = iInteriorSkip();
      bool qInFill = qBoundaryFill();

      if (VTVertTypeToStartFrom <= Vert::eBdryCurve) {
	// Find the distance to the two closest surface points and the
	// closest interior point.
	double dDistS1 = LARGE_DBL, dDistS2 = LARGE_DBL, dDistInt = LARGE_DBL;

	if (dDistInt > dDistS1 && dDistS1 < dDistRatioLimit * dLongDist) {
	  // This is a case where an anisotropic mesh is present along
	  // one or both surfaces that meet at a fold.  March along
	  // lines on the surface keeping every fourth vertex to match
	  // what will happen in the interior of the mesh.  Only march
	  // from curve verts that are marked to be kept.
	  iNSkip = 3;
	  qInFill = false;
	}
	// For all other cases where the nearest surface vertex is
	// close, we're going to keep every second vertex on the
	// surface.
	iNSkip = 1;
	// If the two surfaces meeting at this curve are both
	// anisotropically meshed, we'll fill in along the curve.
	// Otherwise, we won't.
	qInFill = (dDistS2 < dDistRatioLimit * dLongDist);
      }

      // When qInFill is set, keep all verts that have a chain coming
      // out from them.
      bool qInitialDeleteRequest = pV->qDeletionRequested();
      int iNewlyAdded = 0;
      if (qInFill && qInitialDeleteRequest) {
	pV->vMarkToKeep();
	iNewlyAdded++;
      }

      // Only check neighbors as starting points of chains.
      std::set<Vert*> spVNeigh;
      {
	std::set<Cell*> spCTmp;
	vNeighborhood(pV, spCTmp, spVNeigh);
      }

      std::set<Vert*>::iterator iter;
      for (iter = spVNeigh.begin(); iter != spVNeigh.end(); iter++) {
	Vert *pVCand = *iter;
	if (!pVCand->qActive()) continue;
	int iAddedThisChain = 0;
	double dDistRatio = dDIST3D(pVCand->adCoords(), pV->adCoords())
	  / dLongDist;
	if (dDistRatio < dDistRatioLimit) {
	  pV->vMarkStructured();
	  if (!pV->qDeletionRequested()) {
	    iAddedThisChain = iMarkChain(pV, pVCand, aVC, dLongDist,
					 iNSkip);
	    iNewlyAdded += iAddedThisChain;
	  }
	} // Done running a chain out in one direction
	// The following clause obliterates chains that didn't manage
	// to build anything (including those where spacing was too
	// large initially), as well as chains that weren't meant to
	// build anything in the first place.
	if (iAddedThisChain == 0 && dDistRatio < 0.5) {
	  vDeactivateChain(pV, pVCand, aVC, dLongDist);
	}
      } // Done checking all active neighbors of this source vert
      // If the vertex was tagged for inclusion in the fine mesh only
      // because of a (non-existant) chain that could be built from it,
      // then re-mark the vertex for deletion.

      if (iNewlyAdded == 1 && qInitialDeleteRequest) {
	pV->vMarkToDelete ();
	iNewlyAdded--;
      }
      vMessage(4, "From vertex %d, added %d verts in chain.\n",
	       iV, iNewlyAdded);
      iAdded += iNewlyAdded;
    } // Done checking verts that can have a PSAM line built from them
  } // Done checking all verts
  vMessage(2, "Added %d verts in pseudo-structured anisotropic mesh\n",
	   iAdded);
}

void Mesh::vAddToCoarseVerts(const VertConnect aVC[]) const
{
  int iV, iChanges = 0;

  // The initial MIS is set up using a frontal advance through the mesh,
  // beginning from all kept verts, and working outward a layer at a
  // time.  For BdryCurves in 3D and Bdrys in 2D, this will produce
  // near-maximal size on the first pass, which lexicographic ordering
  // of verts does not accomplish and which subsequent improvement can
  // not achieve because there isn't enough wiggle room to move vertices
  // around.

  {
    set<Vert*> spVChecked;
    for (iV = 0; iV < iNumVerts(); ++ iV) {
      Vert *pV = pVVert(iV);
      if (!pV->qDeletionRequested()) {
	spVChecked.insert(pV);
      }
    }

    // Now do the frontal advance...
    //    int iOldLength = 0;
    set<Vert*> spVCurrentLayer;
    // Make an array copy of the latest layer of verts checked; this is
    // convenient to have, especially as the identity of new layers is
    // easiest to produce as an array.
    Vert** apVPrevLayer = new Vert*[spVChecked.size()];
    Vert** ppVEnd = std::copy(spVChecked.begin(), spVChecked.end(),
			      apVPrevLayer);
    assert(static_cast<unsigned int>(ppVEnd - apVPrevLayer) == spVChecked.size());
    int iPrevSize = spVChecked.size();
    do {
//        int iStop = iOldLength;
//        iOldLength = LpVChecked.iLength();
//        for (iV = LpVChecked.iLength() - 1; iV >= iStop; iV--) {
//	Vert *pV = LpVChecked[iV];
      spVCurrentLayer.clear();
      for (int ii = iPrevSize - 1; ii >= 0; ii--) {
	Vert *pV = apVPrevLayer[ii];
	int iVert = iVertIndex(pV);
	// Add all active neighbors of the current vert to the list, and
	// mark each of these verts for inclusion in the coarse mesh if
	// there isn't a conflicting vert already included.
	for (int iii = 0; iii < aVC[iVert].iSize(); iii++) {
	  Vert *pVNeigh = aVC[iVert].pVNeigh(iii);
	  if (pVNeigh->qActive() && pVNeigh->qDeletionRequested()) {
//	    LpVChecked.vAddItem(pVNeigh);
	    // Add to the list of all checked vertices and to the list
	    // of vertices checked on this pass.
	    spVCurrentLayer.insert(pVNeigh);
	    int iNeigh = iVertIndex(pVNeigh);
	    if (aVC[iNeigh].iCountKeptNeigh() == 0) {
	      pVNeigh->vMarkToKeep();
	      iChanges++;
	    }
	  }
	} // Done looping over one listed verts neighbors
      } // Done looping over all verts inserted in the last layer of the
	// advance

      // Put a list of the newest layer of verts into apVPrevLayer.
      delete [] apVPrevLayer;
      apVPrevLayer = new Vert*[spVCurrentLayer.size()];
      ppVEnd = set_difference(spVCurrentLayer.begin(),
			      spVCurrentLayer.end(),
			      spVChecked.begin(),
			      spVChecked.end(),
			      apVPrevLayer);
      iPrevSize = ppVEnd - apVPrevLayer;

      // Now update the list of all verts that have been checked.
      spVChecked.insert(apVPrevLayer, ppVEnd);
    } while (iPrevSize > 0);
  } // Done with first insertion pass over vertices.
  vMessage(2, "Initial site insertion pass for this phase added %7d verts\n",
	   iChanges);

  int iIter = 0;
  do {
    iIter++;
    iChanges = 0;
    for ( iV = 0; iV < iNumVerts() ; ++ iV ) {
      Vert *pV = pVVert(iV);
      if (pV->qActive() && pV->qDeletionRequested()) {
	int iNumKeptNeigh = aVC[iV].iCountKeptNeigh();

	if ( iNumKeptNeigh == 0 ) {
	  pV->vMarkToKeep();
	  ++ iChanges;
	}
	else if ( iNumKeptNeigh == 1) {
	  int iKept2nd = iCountKept2nd(iV, aVC, this);
	  for (int iNeigh = 0; iNeigh < aVC[iV].iSize(); iNeigh++) {
	    Vert *pVNeigh = aVC[iV].pVNeigh(iNeigh);
	    if (pVNeigh->qActive() && !pVNeigh->qDeletionRequested()) {
	      int iDiff = iKept2nd -
		iCountKept2nd(iVertIndex(pVNeigh), aVC, this);
	      double dThresh = atan( ( double(iDiff) -0.5 ) * 5. )
		/ M_PI + 0.5;
#ifndef WIN_VC32
	      if ( drand48() < dThresh ) {
#else
	      if ( rand()/double(RAND_MAX) < dThresh ) {
#endif
		assert(pV->qActive());
		assert(pVNeigh->qActive());
		pVNeigh->vMarkToDelete();
		pV->vMarkToKeep();
	      }
	      break;
	    } /* neighbor marked for keeping in the coarse mesh */
	  }   /* loop over first neighbors */
	}     /* exactly one kept neighbor */
      }       /* current vertex to be deleted */
    }         /* loop over all vertices */
    vMessage(2, "Site insertion pass%3d added%7d verts\n",
	     iIter, iChanges);
  } while (iIter < 20 && iChanges); /* iteration loop */

  // End with a pass over all vertices to ensure that the marked
  // vertices are a maximal independent set.  There is a small chance
  // that movement will result in a non-maximal set otherwise.
  for (iV = 0; iV < iNumVerts(); ++ iV) {
    Vert *pV = pVVert(iV);
    if (pV->qActive() && pV->qDeletionRequested()) {
      int iNumKeptNeigh = aVC[iV].iCountKeptNeigh();
      if (iNumKeptNeigh == 0) {
	pV->vMarkToKeep();
	++ iChanges;
      }
    }
  } // Done with final insertion pass over vertices.
  vMessage(2, "Final site insertion pass for this phase added %7d verts\n", iChanges);
}

void Mesh::vDeactivateAll() const
{
  for (int iV = 0; iV < iNumVerts(); iV++) {
    pVVert(iV)->vMarkInactive();
  }
}

int Mesh::iActivate(const Vert::VertType VT) const
{
  int iRetVal = 0;
  for (int iV = 0; iV < iNumVerts(); iV++) {
    Vert *pV = pVVert(iV);
    if (pV->iVertType() == VT) {
      pVVert(iV)->vMarkActive();
      iRetVal++;
    }
  }
  return iRetVal;
}

void Mesh::vSelectCoarseVerts(VertConnect aVC[])
  // This routine is designed to choose which vertices in a mesh should
  // remain when the mesh is coarsened.  These vertices are tagged using
  // Vert::vMarkToKeep; those that will be deleted are tagged using
  // Vert::vMarkToDelete.
  //
  // Selection is heirarchical, working first on bdry curves, then on
  // the entire boundary, then in the interior of the domain.  At each
  // level, first a check is made for stretched (pseudo-structured) mesh
  // regions, by building chains from verts already tagged for retention
  // into the next higher-dimensional entity.
{
  vMessage(1, "Deciding which verts to keep for the coarse mesh...\n");
  // Tag every vertex as not ACTIVE and DELETE_REQUESTED (except for
  // bdry apexes, which must always be keep for geometric integrity).
  int iV;
  for (iV = 0; iV < iNumVerts(); iV++) {
    Vert *pV = pVVert(iV);
    pV->vMarkInactive();
    if (pV->iVertType() == Vert::eBdryApex) {
      pV->vMarkToKeep();
      vMessage(4, "Vertex %d is an apex; kept.\n", iV);
    }
  }

  // First, find verts to add along bdry curves.  This will apply only
  // for 3D (and possible surface; that isn't clear yet), and checking
  // the number of Verts found on BdryCurves will short-circuit the real
  // work here when appropriate.
  int iActive = iActivate(Vert::eBdryCurve);
  if (iActive > 0) {
    vMessage(2, " Deciding which verts to keep along bdry curves...\n");
    // Add a MIS of maximal size within the active verts.
    vAddToCoarseVerts(aVC);

    // Lock all verts.
    vDeactivateAll();
  }

  // Now work with the rest of the boundary, including internal
  // boundaries and pseudo-surfaces.
  iActive = (iActivate(Vert::eBdryTwoSide) +
	     iActivate(Vert::eBdry) +
	     iActivate(Vert::ePseudoSurface));
  if (iActive > 0) {
    vMessage(2, " Deciding which verts to keep on boundaries...\n");
    vAddNormalVerts(aVC, Vert::eBdryCurve);
    vUpdateLengthScale();
    vUpdateConflictGraph(aVC);
    vAddToCoarseVerts(aVC);
    vDeactivateAll();
  }

  // In regions where the volume mesh is structured (STRUCTURED),
  // choose every fourth vertex along structured lines.
  iActive = iActivate(Vert::eInterior);
  if (iActive > 0) {
    vMessage(2, " Deciding which verts to keep in the interior...\n");
    // March out from surface-like stuff, skipping three vertices between
    // kept verts (keeping 1 out of 4) to equalize aspect ratios sooner.
    vAddNormalVerts(aVC, Vert::ePseudoSurface);
    vUpdateLengthScale();
    //    vUpdateConflictGraph(aVC);
    vAddToCoarseVerts(aVC);
    vDeactivateAll();
  }

//    vMessage(2, " Cleaning up along boundaries...\n");
//    // To improve coarse mesh quality near boundaries, retain any deleted
//    // non-STRUCTURED boundary vertex that is in conflict with a kept,
//    // non-STRUCTURED interior vertex.
//    for (iV = 0; iV < iNumVerts(); iV++) {
//      Vert *pV = pVVert(iV);
//      if (pV->qIsBdryVert() && pV->qDeletionRequested() &&
//	!pV->qIsStructured()) {
//        for (int ii = 0; ii < aVC[iV].iSize(); ii++) {
//	Vert *pVNeigh = aVC[iV].pVNeigh(ii);
//	if (!pVNeigh->qIsBdryVert() && !pVNeigh->qDeletionRequested() &&
//	    !pVNeigh->qIsStructured()) {
//	  pV->vMarkToKeep();
//	}
//        }
//      }
//    }
  vMessage(1, "Done deciding which verts to keep for the coarse mesh...\n");
}

void Mesh::vRemoveTaggedVerts()
{
  // Remove vertices from the mesh, with intermittent swapping to ease
  // the process and get good mesh quality.
  int iDeletedThisPass, iDeletedSincePurge, iTotalDeleted;
  int iPasses, iNumSwaps, iNMarkedForDeletion, iV;

  iDeletedSincePurge = iTotalDeleted = iPasses = iNumSwaps =
    iNMarkedForDeletion = 0;

//   If one were doing Delaunay vertex removal (which is non-robust),
//   then Delaunay swapping would be the right thing to do.  Otherwise,
//   it just wrecks the connectivity so that you can't remove verts as
//   well as you otherwise would.
//   vSetSwapType(eDelaunay);
  if (eSwapType() != eMinSine) {
    vSetSwapType(eMinSine);
    iSwap(1);
  }

  for (iV = 0; iV < iNumVerts(); iV++) {
    vMessage(4, "Vertex %d at (%f, %f) is%s tagged for deletion.\n",
	     iV, pVVert(iV)->dX(), pVVert(iV)->dY(),
	     pVVert(iV)->qDeletionRequested() ? "" : " not");
    if (pVVert(iV)->qDeletionRequested())
      iNMarkedForDeletion++;

  }

  vMessage(1, "Removing vertices from the mesh using strategy %d.\n",
	   STRATEGY);
  bool qBdryOnly = false;
  do {
    iPasses++;
    if (iPasses >= 3) qBdryOnly = false;
    vMessage(1, "  Pass %d:  ", iPasses);
    iDeletedThisPass = 0;
    int iPrintInt = max(25, iNumVerts() / 10 + 1);

#if (REV_LEX == 1) // Reverse lexicographic order
    for (iV = iNumVerts() - 1; iV >= 0; iV--) {
#else // Forward lexicographic order
    for (iV = 0; iV < iNumVerts(); iV++) {
#endif
      vMessage(4, "Checking vertex %d for removal\n", iV);
      Vert *pV = pVVert(iV);
#if (STRATEGY == 3)
      if (pV->qDeleted() ||
	  (!pV->qDeletionRequested()) ||
	  (!pV->qIsBdryVert() && qBdryOnly)) continue;
      // Contract only onto the nearest kept vertex.
      std::set<Cell*> spCDummy;
      std::set<Vert*> spVNeigh;
      vNeighborhood(pV, spCDummy, spVNeigh);

      // Find the nearest kept vertex.
      double dMinDist = LARGE_DBL;
      Vert *pVClosest = pVInvalidVert;
      for (std::set<Vert*>::iterator iterV = spVNeigh.begin();
	   iterV != spVNeigh.end(); iterV++) {
	Vert *pVCand = *iterV;
	if (!pVCand->qDeletionRequested()) {
	  double dDist = dDIST3D(pVCand->adCoords(), pV->adCoords());
	  if (dDist < dMinDist) {
	    dMinDist = dDist;
	    pVClosest = pVCand;
	  }
	}
      }
      // If there -is- a kept neighbor, try the deletion.
      int iNewSwaps;
      if (pVClosest->qValid() &&
	  qRemoveVertByContraction(pV, iNewSwaps, &pVClosest, 1)) {
	iNumSwaps += iNewSwaps;
	iTotalDeleted++;
	iDeletedThisPass++;
	iDeletedSincePurge++;
      } // Deleted a vert.
#elif (STRATEGY == 2)
      if (pV->qDeleted() ||
	  (!pV->qIsBdryVert() && qBdryOnly)) continue;
      int iNewSwaps;
      // Boundary verts are deleted using a different call than interior
      // verts, because the interior call panics when handed a bdry vert!
      if (pV->qDeletionRequested() &&
	  qRemoveVertByContraction(pV, iNewSwaps, NULL, 0)) {
//	    qRemoveVert(pV, iNewSwaps)) {
//	  ((pV->qIsBdryVert() &&
//	    qRemoveVertByContraction(pV, iNewSwaps, NULL, 0)) ||
//	   (!pV->qIsBdryVert() &&
//	    qRemoveVert(pV, iNewSwaps)))) {
	iNumSwaps += iNewSwaps;
	iTotalDeleted++;
	iDeletedThisPass++;
	iDeletedSincePurge++;
      } // Deleted a vert.
#elif (STRATEGY == 1)
    // A new approach: for each kept vertex, try to delete all deletable
    // neighbors.
      // Under odd circumstances, a vertex whose deletion was not
      // requested may actually be removed from the mesh, so we must
      // check for this...
      if (!pV->qDeletionRequested() && !pV->qDeleted()) {
	// Find all neighboring verts, then throw out the ones that
	// don't need to be deleted.
	std::set<Cell*> spCDummy;
	std::set<Vert*> spVNeigh;
	vNeighborhood(pV, spCDummy, spVNeigh);

	std::set<Vert*>::iterator iterV;
	for (iterV = spVNeigh.begin(); iterV != spVNeigh.end(); ) {
	  Vert *pVCand = *iterV;
	  assert(!pVCand->qDeleted());
	  if (!pVCand->qDeletionRequested()) {
	    std::set<Vert*>::iterator iterTmp = iterV++;
	    spVNeigh.erase(iterTmp);
	  }
	  else iterV++;
	}

	// Go through the list repeatedly until no more deletions are
	// made onto this vertex.
	bool qChange;
	do {
	  qChange = false;
	  for (iterV = spVNeigh.begin(); iterV != spVNeigh.end(); ) {
	    Vert *pVCand = *iterV;
	    std::set<Vert*>::iterator iterTmp = iterV++;
	    int iNewSwaps = 0;
	    if (pVCand->qDeleted() ||
		qRemoveVertByContraction(pVCand, iNewSwaps, &pV, 1)) {
	      qChange = true;
	      spVNeigh.erase(iterTmp);
	      iNumSwaps += iNewSwaps;
	      iTotalDeleted++;
	      iDeletedThisPass++;
	      iDeletedSincePurge++;
	    } // Deleted a vert.
	  } // Loop over neighboring verts.
	} while (qChange && !spVNeigh.empty());
      } // Is this vert being kept?
#endif
      if (iV % iPrintInt == 0 && iDeletedThisPass != 0) {
	vMessage(2, "At vertex %d out of %d; %d verts removed, %d swaps\n",
		 iV, iNumVerts(), iDeletedThisPass, iNumSwaps);
      }
    } // Loop over verts.
    assert(qValid());

    if (10 * iDeletedSincePurge > iNumVerts()) {
      vMessage(2, "Purging deleted mesh entities...\n");
      vPurge();
      iDeletedSincePurge = 0;
    }
    vMessage(2, "%6d vertices removed, %6d total.  %6d still to remove\n",
	     iDeletedThisPass, iTotalDeleted,
	     iNMarkedForDeletion - iTotalDeleted);
    // If nothing was deleted, this would be a waste of time.
    if (iDeletedThisPass != 0) {
      vMessage(2, "  Smoothing to improve cell shapes...\n");
      iSmooth(1);
      vMessage(2, "  Swapping to improve connectivity...\n");
      iNumSwaps += iSwap(1);
      // Using 3D bad cell repair techniques does no good here.
    }
  } while (iTotalDeleted < iNMarkedForDeletion
	   && iDeletedThisPass > 0 );
  vMessage(2, "Deleted %d vertices out of %d marked for deletion.\n",
	   iTotalDeleted, iNMarkedForDeletion);
  vMessage(2, "A total of %d swaps were performed %s\n",
	   iNumSwaps, "to allow deletion and improve mesh quality.");

  if (iDeletedSincePurge != 0) vPurge();

  // Info on undeleted verts.
  vMessage(2, "The following verts could not be deleted:\n");
  if (iTotalDeleted < iNMarkedForDeletion) {
    int iBdryCurve = 0, iBdryStruct = 0, iBdry = 0, iBdryTotal = 0;
    int iIntStruct = 0, iInt = 0, iIntTotal = 0, iGrandTotal = 0;
    for (iV = 0; iV < iNumVerts(); iV++) {
      Vert *pV = pVVert(iV);
      assert(!pV->qDeleted());
      if (pV->qDeletionRequested()) {
	switch (pV->iVertType()) {
	case Vert::eBdryCurve:
	  iBdryCurve++;
	  iBdryTotal++;
	  iGrandTotal++;
	  break;
	case Vert::eBdry:
	  if (pV->qIsStructured())
	    iBdryStruct++;
	  else
	    iBdry++;
	  iBdryTotal++;
	  iGrandTotal++;
	  break;
	case Vert::eInterior:
	  if (pV->qIsStructured())
	    iIntStruct++;
	  else
	    iInt++;
	  iIntTotal++;
	  iGrandTotal++;
	  break;
	default:
	  assert(0);
	  break;
	}
	if (pV->iSpaceDimen() == 2) {
	  vMessage(2, "Vert index %6d at (%13.6g, %13.6g). Type: %d\n",
		   iV, pV->dX(), pV->dY(), pV->iVertType());
	}
	else {
	  vMessage(2, "Vert index %6d at (%13.6g, %13.6g, %13.6g). Type: %d\n",
		   iV, pV->dX(), pV->dY(), pV->dZ(), pV->iVertType());
	}
      }
    }
    assert(iGrandTotal == iNMarkedForDeletion - iTotalDeleted);
    assert(iGrandTotal == iIntTotal + iBdryTotal);
    assert(iBdryTotal == iBdry + iBdryStruct + iBdryCurve);
    assert(iIntTotal == iInt + iIntStruct);
    vMessage(2, "Total verts not removed: %d\n", iGrandTotal);
    vMessage(2, " On the bdry: %3d on curves, %3d structured, %3d others. Total: %d\n",
	     iBdryCurve, iBdryStruct, iBdry, iBdryTotal);
    vMessage(2, " In interior:                %3d structured, %3d others. Total: %d\n",
	     iIntStruct, iInt, iIntTotal);
  }

#ifdef ANISO_REINSERT
//    {
//      //   An incredible hack to allow me to look at the surface mesh after
//      // all the interior verts have been pulled out.
//      VolMesh *pVM = static_cast<VolMesh*>(this);
//      vWriteFile_VolMesh(*pVM, "after-strip");
//    }

  // Re-insert interior structured verts, if needed.
  int iNChains = vecCReinsert.size(), iNVertsToReinsert = 0;
  int iC, ii;
  for (iC = 0; iC < iNChains; iC++) {
    iNVertsToReinsert += vecCReinsert[iC].iLength();
  }
  if (iNVertsToReinsert == 0) {
    vMessage(3, "No interior structured verts to re-insert.\n");
  }
  else {
    vMessage(2, "After stripping all interior structured verts, need to\n");
    vMessage(2, "reinsert %d, located along %d chains.\n",
	     iNVertsToReinsert, iNChains);

    SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
    int iTotalSwaps = 0;
    for (iC = 0; iC < iNChains; iC++) {
      Chain& C = vecCReinsert[iC];
      int iLen = C.iLength();

      // The following is a horrible hack, in that its running time is
      // basically O(NVerts * NBdryVerts).  It is, however, robust and
      // easy to code.
      Vert *pVLast;
      {
	for (iV = 0; iV < iNumVerts(); iV++) {
	  Vert *pVTmp = pVVert(iV);
	  if (pVTmp->qIsBdryVert() &&
	      pVTmp->dX() == C.adBdryLoc[0] &&
	      pVTmp->dY() == C.adBdryLoc[1] &&
	      (pVTmp->iSpaceDimen() == 2 ||
	       pVTmp->dZ() == C.adBdryLoc[2])) {
	    pVLast = pVTmp;
	    break;
	  }
	}
      }

      for (ii = 0; ii < iLen; ii++) {
	// Get location of new vert
	double adCoords[] = {C.dX(ii), C.dY(ii), C.dZ(ii)};
	// Create new vert
	Vert *pVNew = pVNewVert();
	pVNew->vSetCoords(pVLast->iSpaceDimen(), adCoords);
	// Find starting guess of cell (using BeginPipe, which should be
	// made a mesh member function and replicated for 2D (and
	// eventually for surfaces?))
	Cell *pCGuess = pCBeginPipe(pVLast, pVNew);
	// Insert (Lawson, probably)
	int iSwaps = 0;
	if (qInsertPoint(adCoords, pCGuess, &iSwaps, true, true,
			 pVNew)) {
	  pVNew->vMarkStructured();
	  iTotalSwaps += iSwaps;
	  pVLast = pVNew;
	}
      } // Done with inserting this chain
    } // Done inserting all chains
    SUMAA_LOG_EVENT_END(INSERT_VERTEX);

    // Try again to eliminate verts that weren't successfully removed
    // before?
    iSwap();
  } // Done with case where there were verts to re-insert
#endif

  // Finish up.
  if (!qValid())
    vFoundBug("removing verts to create a coarse mesh");

}

void Mesh::vUpdateConflictGraph(VertConnect aVC[]) const
  // The goal is to create a list, for each vertex i, of vertices j for
  // which the distance between V_i and V_j is smaller than the sum of
  // their length scales. An empirical scaling factor is used here to
  // make sure approximately the right number of verts remain after
  // coarsening.
  //
  // To accomplish this, a search marches outward from each vertex,
  // seeking all points that are topologically near a given vertex and
  // also geometrically nearby.  A strictly geometric search pierces
  // thin boundaries, which is bad, so that approach has been
  // abandoned.
  //
  // Only vertices that are current marked active are considered; others
  // are either already taken care of or not yet relevant.
{
  SUMAA_LOG_EVENT_BEGIN(PROBE1);
  vMessage(1, "Updating vertex-to-vertex conflict graph using length scale...\n");
  int iV, iInfoInterval = max(25, iNumVerts() / 10);
  double dBubbleSizeConst = dCoarsenConstant();
  int iUpdated = 0;
  unsigned int uiMaxCand = 100;
  Vert **apVCand = new Vert*[uiMaxCand];

  for (iV = 0; iV < iNumVerts(); iV++) {
    if (iV % iInfoInterval == 0) {
      vMessage(2, "  Updating conflicts for vertex %d...\n", iV);
    }
    Vert *pV = pVVert(iV);
    if (!pV->qActive()) continue;

    iUpdated++;


    // Initially, create a conflict set containing only pV, and a
    // candidate array containing the current conflict list for pV.
    // Also, initialize the set of tested verts to pV.

    set<Vert*> spVConflict, spVTested;
    spVConflict.insert(pV);
    spVTested.insert(pV);

    int iNCand = aVC[iV].iSize();
    if (uiMaxCand < static_cast<unsigned int>(abs(iNCand))) {
      uiMaxCand = iNCand;
      delete [] apVCand;
      apVCand = new Vert*[uiMaxCand];
    }
    Vert **ppVFirst = aVC[iV].ppVVerts();
    Vert **ppVLast = std::copy(ppVFirst, ppVFirst + iNCand,
			       apVCand);
    assert(ppVLast - apVCand == iNCand);

    int iPass = 0;
    set<Vert*> spVNewCand;
    do {
      spVNewCand.clear();
      spVTested.insert(apVCand, apVCand + iNCand);
      iPass ++;

      for (int i = 0; i < iNCand; i++) {
	Vert *pVCand = apVCand[i];
	if (!pVCand->qActive()) continue;

	double dDist = dDistanceBetween(pV, pVCand);

	// The definition of dMinLegalDist has been tweaked so that, in
	// 2D, the mesh is about one-fourth as dense as before
	// coarsening and the point distribution for coarse meshes is
	// esthetically reasonable.
	double dMinLegalDist = dBubbleSizeConst
	  * (pV->dLS() + pVCand->dLS());

	// If this is a bdry vert, then its neighbors on the bdry are
	// always considered to be in conflict.  This makes things much
	// nicer for 2D meshes, at least.
	if (dDist <= dMinLegalDist
	    ||
	    (iPass == 1 && pV->qIsBdryVert() && pVCand->qIsBdryVert())
	    ) {
	  // Conflict!  Add to conflict list and be sure to check
	  // anything in conflict with it.
	  vMessage(4, "Verts %d and %d. Dist = %12.6f, min legal = %12.6f, ratio = %12.6f\n",
		   iVertIndex(pV), iVertIndex(pVCand), dDist, dMinLegalDist,
		   dDist / dMinLegalDist);
	  spVConflict.insert(pVCand);
	  int iCand = iVertIndex(pVCand);
	  Vert **ppVBegin = aVC[iCand].ppVVerts();
	  spVNewCand.insert(ppVBegin, ppVBegin + aVC[iCand].iSize());
	} // Done with conflicting candidate
      } // Done checking all candidates
      if (spVNewCand.size() > uiMaxCand) {
	delete [] apVCand;
	uiMaxCand = spVNewCand.size()*2;
	apVCand = new Vert*[uiMaxCand];
      }
      ppVLast = std::set_difference(spVNewCand.begin(),
				    spVNewCand.end(),
				    spVTested.begin(),
				    spVTested.end(),
				    apVCand);
      iNCand = ppVLast - apVCand;
    } while (iNCand > 0);// No more candidates to check!
    spVConflict.erase(pV);

    // Construct a VertConnect entry for vertex iV
    aVC[iV].vSetSize(spVConflict.size());
    std::set<Vert*>::iterator it;
    int ii = 0;
    for (it = spVConflict.begin(); it != spVConflict.end(); it++, ii++) {
      aVC[iV].vSetNeigh(ii, *it);
    }
    assert(ii == aVC[iV].iSize());
  } // End of loop over verts
  vMessage(1, "Updated conflict graph for %d out of %d verts.\n",
	   iUpdated, iNumVerts());
  vMessage(2, "Maximum number of candidates considered at once was %d\n",
	   uiMaxCand);
  delete [] apVCand;
  SUMAA_LOG_EVENT_END(PROBE1);
}

void Mesh::vInitConflictGraph(VertConnect aVC[]) const
  // The full conflict graph is overkill for coarsening boundary curves
  // and for anisotropic coarsening.  In those cases, neighbor info is
  // all that's needed.  So compute neighbor info only, and later on
  // compute the true conflict graph incrementally using
  // vUpdateConflictGraph.
{
  vMessage(1, "Creating vertex-to-vertex conflict graph using neighbor data only...\n");
  int iV, iInfoInterval = max(25, iNumVerts() / 10);

  for (iV = 0; iV < iNumVerts(); iV++) {
    if (iV % iInfoInterval == 0) {
      vMessage(2, "  Checking conflicts for vertex %d...\n", iV);
    }
    Vert *pV = pVVert(iV);
    std::set<Vert*> spVConflict;
    {
      std::set<Cell*> spCTmp;
      vNeighborhood(pV, spCTmp, spVConflict);
    }
    spVConflict.erase(pV);

    // Construct a VertConnect entry for vertex iV
    aVC[iV].vSetSize(spVConflict.size());
    std::set<Vert*>::iterator iter;
    int ii = 0;
    for (iter = spVConflict.begin(); iter != spVConflict.end();
	 iter++, ii++) {
      aVC[iV].vSetNeigh(ii, *iter);
    }
  } // End of loop over verts
  vMessage(1, "Done initializing conflict graph.\n");
}

// Create a coarse unstructured mesh.
void Mesh::vCoarsen(const double dLengthRatio)
{
  assert(qSimplicial());
  // To improve vertex removal efficiency, pre-swapping has been tried.
  // These efforts were unsuccessful and have been scrapped.
  vMessage(2, "Increasing length scale at each vert by a factor of %.1f.\n",
	   dLengthRatio);
  int iV;
  for (iV = 0; iV < iNumVerts(); iV++) {
    pVVert(iV)->vSetLS(dLengthRatio * pVVert(iV)->dLS());
    spVUpdateLS.insert(pVVert(iV));
  }
  vMessage(2, "Done increasing length scale.\n");
  vCoarsenToLengthScale();
}

void Mesh::vCoarsenToLengthScale(const bool qHeal)
{
  assert(qSimplicial());
  bool qRebuildPatchesAfter = (iNumBdryFaces() == iNumBdryPatches());
  vMessage(1, "Coarsening mesh...\n");
  // Now make sure that grading, etc, is okay...
  vMessage(2, "Updating length scale to take grading into account...\n");
  int iV;
  for (iV = 0; iV < iNumVerts(); iV++) {
    spVUpdateLS.insert(pVVert(iV));
  }
  vUpdateLengthScale();

  vMessage(2, "Done with grading check.\n");

  // Classify verts on the basis of geometry only; this should prevent
  // degeneration of holes with multiple coarsening.

  // The following code is a little bit dubious, in that it could cause
  // problems with points that lie on patch boundaries, especially in
  // the presence of re-reading geometry for periodic cases.  So for now
  // just make sure that anything definitively ID'd as an apex from the
  // bdry patches isn't reclassified. 
  // FIX ME 0.5.0
  vIdentifyVertexTypesGeometrically();
  for (iV = 0; iV < iNumVerts(); iV++) {
    vMessage(4, "Vert: %5d  Type: %d\n", iV, pVVert(iV)->iVertType());
  }

  // Always keep bdry apexes.  Also, for the mesh healing case, keep all
  // boundary verts.  This should be okay unless the boundary itself has
  // changed shape rapidly, but it is a bit of a hack. CFO-G 10/17/02
  for (iV = 0; iV < iNumVerts(); iV++) {
    Vert *pV = pVVert(iV);
    if (pV->iVertType() == Vert::eBdryApex) {
      pV->vMarkToKeep();
    }
    else if (pV->qIsBdryVert() && qHeal) {
      pV->vMarkToKeep();
    }
    else {
      pV->vMarkToDelete();
    }
  }

  VertConnect *aVC = new VertConnect[iNumVerts()];
  vInitConflictGraph(aVC);

  vSelectCoarseVerts(aVC);

  // Remove connectivity information
  delete [] aVC;

  vRemoveTaggedVerts();
  if (qRebuildPatchesAfter)
    vBuildBdryPatches();

  vMessage(1, "Done coarsening mesh\n");
}
