#include <stdlib.h>
#include "GR_Geometry.h"
#include "GR_Mesh2D.h"
#include "SMsmooth.h"
#include "GR_events.h"
#include "GR_BFace.h"
#include "GR_Vec.h"

#define iMaxVertNeigh MAX_NUM_PTS

bool Mesh2D::
qRemoveVertCheaply (Vert * const pV, int& iSwaps)
{
  // Simple:  remove edges until only three remain.  Then the vertex to
  // be removed is surrounded by three triangles which can be combined
  // to form one.  Afterwards, try to swap each edge that once shared a
  // triangle with the vertex.
  //
  // For boundary edges, when three edges remain, two are boundary
  // edges.  Those two edges are merged, the third disappears.

  if (pV->iVertType () == Vert::eBdryApex ||
      pV->iVertType () == Vert::eInteriorFixed)
    return false;

  // Can't remove bdry verts when the boundary is precious.
  if (!qAllowBdryChanges && (pV->iVertType () == Vert::eBdryApex ||
			     pV->iVertType () == Vert::eBdryCurve ||
			     pV->iVertType () == Vert::eBdryTwoSide ||
			     pV->iVertType () == Vert::eBdry))
    return false;

  bool qRetVal = true;
  iSwaps = 0;

  std::set<Face*> spFSwapEdges;
  {
    std::set<Vert*> spVTmp;
    std::set<Cell*> spCTmp;
    vNeighborhood(pV, spCTmp, spVTmp, NULL, NULL, &spFSwapEdges);

    std::set<Cell*>::iterator iter, iterEnd = spCTmp.end();
    for (iter = spCTmp.begin(); iter != iterEnd; iter++) {
      if ((*iter)->eType() == Cell::eQuadCell)
	return false;
    }
  }
  std::set<Face*> spFIncEdges;
  std::set<Face*>::iterator iterF;
  for (iterF = spFSwapEdges.begin(); iterF != spFSwapEdges.end(); iterF++) {
    Face *pF = *iterF;
    if (pF->qHasVert(pV)) spFIncEdges.insert(pF);
  }

  // Now have a list of every edge incident on the vertex.  Cycle
  // through the list, removing edges, until only three remain.
  int iLastLength = -1;
  while (spFIncEdges.size () > 3
	 && int(spFIncEdges.size ()) != iLastLength) {
    iLastLength = spFIncEdges.size ();

    // The List<> version of this loop was doubled (working out from the
    // middle ) to try to fix a prob lem with vertex deletion for
    // structured meshes.
    for (iterF = spFIncEdges.begin(); iterF != spFIncEdges.end(); ) {
      Face *pF = *iterF;
      std::set<Face*>::iterator iterTmp = iterF++;
      if (qRemoveEdge (pF)) {
	spFIncEdges.erase(iterTmp);
      }
    }
  }
  int iNEdges = spFIncEdges.size();
  if (iNEdges < 2 || iNEdges > 4) return false;


  std::vector<Face*> vecpFIncEdges(iNEdges);
  std::copy(spFIncEdges.begin(), spFIncEdges.end(), vecpFIncEdges.begin());
  
  // List the cells incident on these edges
  std::set<Cell*> spCIncCells;
  for (int i = 0; i < iNEdges; i++) {
    spCIncCells.insert (vecpFIncEdges[i]->pCCellLeft ());
    spCIncCells.insert (vecpFIncEdges[i]->pCCellRight ());
  }

  if (pV->iVertType () == Vert::eInterior) {
    if (iNEdges == 4) return false;
    assert(iNEdges == 3 && spCIncCells.size() == 3);
    // Identify all the necessary parts

    Face *pFA1 = vecpFIncEdges[0];
    Face *pFB1 = vecpFIncEdges[1];
    Face *pFC1 = vecpFIncEdges[2];

    Cell *pCA = pCCommonCell (pFB1, pFC1);
    Cell *pCB = pCCommonCell (pFC1, pFA1);
    Cell *pCC = pCCommonCell (pFA1, pFB1);
    assert (spCIncCells.count(pCA) == 1);
    assert (spCIncCells.count(pCB) == 1);
    assert (spCIncCells.count(pCC) == 1);

    Vert *pVA = dynamic_cast<TriCell*>( pCB)->pVVertOpposite (pFC1);
    Vert *pVB = dynamic_cast<TriCell*>( pCC)->pVVertOpposite (pFA1);
    Vert *pVC = dynamic_cast<TriCell*>( pCA)->pVVertOpposite (pFB1);
    spVUpdateLS.insert(pVA);
    spVUpdateLS.insert(pVB);
    spVUpdateLS.insert(pVC);

    Face *pFB2 = dynamic_cast<TriCell*>( pCB)->pFFaceOpposite (pV);
    Face *pFC2 = dynamic_cast<TriCell*>( pCC)->pFFaceOpposite (pV);

    pCA->vReplaceFace (pFB1, pFC2);	// sic
    pCA->vReplaceFace (pFC1, pFB2);	// sic

    pFB2->vReplaceCell (pCB, pCA);
    pFC2->vReplaceCell (pCC, pCA);

    pFA1->vMarkDeleted ();
    pFB1->vMarkDeleted ();
    pFC1->vMarkDeleted ();

    pCB->vMarkDeleted ();
    pCC->vMarkDeleted ();

    pV->vMarkDeleted ();

    pVA->vSetHintFace (pFB2);
    pVB->vSetHintFace (pFC2);
    pVC->vSetHintFace (pFB2);
    assert (pCA->iFullCheck ());
    assert (pFB2->iFullCheck ());
    assert (pFC2->iFullCheck ());
  }				// Done deleting interior vert

  else if (pV->iVertType () == Vert::eBdryTwoSide) {
    // Delete internal boundary vert
    // Version 0.1.8: Implement this
    return false;
  }				// Done deleting internal boundary vert

  else {			// Delete boundary vert
    if (iNEdges == 3) {
      // Identify all the parts
      Cell* apCIncCells[4];
      std::copy(spCIncCells.begin(), spCIncCells.end(), apCIncCells);
      Cell *pCA, *pCB;
      if (apCIncCells[0]->eType () == Cell::eTriCell) {
	pCA = apCIncCells[0];
	if (apCIncCells[1]->eType () == Cell::eTriCell)
	  pCB = apCIncCells[1];
	else if (apCIncCells[2]->eType () == Cell::eTriCell)
	  pCB = apCIncCells[2];
	else
	  pCB = apCIncCells[3];
      }
      else if (apCIncCells[1]->eType () == Cell::eTriCell) {
	pCA = apCIncCells[1];
	if (apCIncCells[2]->eType () == Cell::eTriCell)
	  pCB = apCIncCells[2];
	else
	  pCB = apCIncCells[3];
      }
      else {
	pCA = apCIncCells[2];
	pCB = apCIncCells[3];
      }
      assert (pCA->eType () == Cell::eTriCell);
      assert (pCB->eType () == Cell::eTriCell);
      assert (pCA != pCB);
      
      Face *pFD = pFInvalidFace;
      for (int ii = 0; ii < 3; ii++) {
	pFD = vecpFIncEdges[ii];
	if (pFD->qHasCell (pCA) && pFD->qHasCell (pCB))
	  break;
      }
      assert (pFD->qHasCell (pCA) && pFD->qHasCell (pCB));
      
      Vert *pVA = dynamic_cast<TriCell*>( pCB)->pVVertOpposite (pFD);
      Vert *pVC = (pFD->pVVert (0) == pV) ?
	pFD->pVVert (1) : pFD->pVVert (0);
      spVUpdateLS.insert(pVA);
      spVUpdateLS.insert(pVC);
      
      Face *pFB = dynamic_cast<TriCell*>( pCB)->pFFaceOpposite (pV);
      Face *pFC1 = dynamic_cast<TriCell*>( pCB)->pFFaceOpposite (pVC);
      Face *pFC2 = dynamic_cast<TriCell*>( pCA)->pFFaceOpposite (pVC);
      
      assert (pFC1->pCCellOpposite (pCB)->eType () == Cell::eBdryEdge);
      assert (pFC2->pCCellOpposite (pCA)->eType () == Cell::eBdryEdge);
      
      BdryEdge *pBE1 = dynamic_cast<BdryEdge *>(pFC1->pCCellOpposite (pCB));
      BdryEdge *pBE2 = dynamic_cast<BdryEdge *>(pFC2->pCCellOpposite (pCA));
      
      // Only delete the vertex if the boundary conditions are the same
      if (pBE1->iBdryCond () == pBE2->iBdryCond ()) {
	pCA->vReplaceFace (pFD, pFB);
	pFB->vReplaceCell (pCB, pCA);
	pFC2->vReplaceVert (pV, pVA);
	pFC1->vMarkDeleted ();
	pBE1->vMarkDeleted ();
	pFD->vMarkDeleted ();
	pCB->vMarkDeleted ();
	pV->vMarkDeleted ();
	pVA->vSetHintFace (pFB);
	pVC->vSetHintFace (pFB);
	assert (pCA->iFullCheck ());
      }
      else
	qRetVal = false;
    } // Done with three edges incident on a bdry vert
    else {
      // Two edges incident on a bdry vert
      // Identify all the parts
      assert(iNEdges == 2);
      assert(spCIncCells.size() == 3);
      Cell* apCIncCells[3];
      std::copy(spCIncCells.begin(), spCIncCells.end(), apCIncCells);
      Cell *pCA;
      if (apCIncCells[0]->eType () == Cell::eTriCell) {
	pCA = apCIncCells[0];
      }
      else if (apCIncCells[1]->eType () == Cell::eTriCell) {
	pCA = apCIncCells[1];
      }
      else {
	pCA = apCIncCells[2];
      }
      assert (pCA->eType () == Cell::eTriCell);
      
      Face *pF0 = vecpFIncEdges[0];
      Face *pF1 = vecpFIncEdges[1];
      
      Vert *pVA = (pF0->pVVert (0) == pV) ?
	pF0->pVVert (1) : pF0->pVVert (0);
      Vert *pVC = (pF1->pVVert (0) == pV) ?
	pF1->pVVert (1) : pF1->pVVert (0);
      spVUpdateLS.insert(pVA);
      spVUpdateLS.insert(pVC);
      
      Face *pFA = dynamic_cast<TriCell*>( pCA)->pFFaceOpposite (pV);

      assert (pF0->pCCellOpposite (pCA)->eType () == Cell::eBdryEdge);
      assert (pF1->pCCellOpposite (pCA)->eType () == Cell::eBdryEdge);
      
      BdryEdge *pBE0 = dynamic_cast<BdryEdge *>(pF0->pCCellOpposite (pCA));
      BdryEdge *pBE1 = dynamic_cast<BdryEdge *>(pF1->pCCellOpposite (pCA));
      
      // Only delete the vertex if the boundary conditions are the same
      if (pBE0->iBdryCond () == pBE1->iBdryCond ()) {
 	pFA->vReplaceCell (pCA, pBE0);
	pBE0->vReplaceFace(pF0, pFA);
	pF0->vMarkDeleted ();
	pF1->vMarkDeleted ();
	pBE1->vMarkDeleted ();
	pCA->vMarkDeleted ();
	pV->vMarkDeleted ();
	pVA->vSetHintFace (pFA);
	pVC->vSetHintFace (pFA);
      }
      else
	qRetVal = false;
    } // Done with two edges incident on a bdry vert
  }				// Done deleting boundary vert

  SUMAA_LOG_EVENT_BEGIN (NET_SWAPPING);
  for (iterF = spFSwapEdges.begin();
       iterF != spFSwapEdges.end(); iterF++) {
    EdgeFace *pEF = dynamic_cast<EdgeFace*>(*iterF);
    if (pEF->qValid () && !pEF->qDeleted ()) {
      pEF->pVVert (0)->vMarkIllShaped ();
      pEF->pVVert (1)->vMarkIllShaped ();
      iSwaps = iSwaps + iFaceSwap (pEF);
    }
  }

  SUMAA_LOG_EVENT_END (NET_SWAPPING);
  vUpdateLengthScale();
  return qRetVal;
}

bool Mesh2D::
qBdrySmooth (Vert * const pVBdry, Face * const pFBdry,
	     BdryEdge * const pBE_Init)
{
  assert (qSimplicial ());
  assert (pBE_Init->qHasFace (pFBdry));
  assert (pFBdry->qHasVert (pVBdry));

  if (!qAllowBdryChanges || pVBdry->iVertType() == Vert::eBdryApex)
    return false;

  // This code doesn't behave properly for internal bdrys, because it
  // doesn't wrap around to include both regions in the data passed to
  // OptMS.
  if (pBE_Init->eType() == Cell::eIntBdryEdge)
    return false;

  // Also misbehaves for periodic boundaries, for the same reasons as
  // for internal boundaries.
  if (pBE_Init->pPatchPointer()->qIsPeriodic()) return false;

  bool qRetVal = false;
  double **a2dNeighPts = new double *[iMaxVertNeigh];
  int **a2iFaceConn = new int *[iMaxVertNeigh];
  int i;
  for (i = 0; i < iMaxVertNeigh; i++) {
    a2dNeighPts[i] = new double[2];
    a2iFaceConn[i] = new int[2];
  }
  Face **apFIncFaces = new Face *[iMaxVertNeigh];

  int iNeigh = 0;

  Cell *pC = pBE_Init;
  Face *pF = pFBdry;
  do {
    if (iNeigh < iMaxVertNeigh) {
      Vert *pVOpp;
      if (pF->pVVert (0) == pVBdry)
	pVOpp = pF->pVVert (1);
      else
	pVOpp = pF->pVVert (0);
      a2dNeighPts[iNeigh][0] = pVOpp->dX ();
      a2dNeighPts[iNeigh][1] = pVOpp->dY ();
      apFIncFaces[iNeigh] = pF;
      iNeigh++;
    }

    pC = pF->pCCellOpposite (pC);
    if (pC->eType () == Cell::eTriCell) {
      for (int iF = 0; iF < 3; iF++) {
	Face *pFEdge = pC->pFFace (iF);
	assert (pFEdge->qValid ());
	if (pF == pFEdge)
	  continue;
	if (pFEdge->pVVert (0) == pVBdry ||
	    pFEdge->pVVert (1) == pVBdry) {
	  pF = pFEdge;
	  break;
	}
      }
    }
  }
  while (pC->eType () == Cell::eTriCell);

  // Modified 21/09/99 by Charles Boivin
  // Added the possibility of an internal boundary
  assert ((pC->eType () == Cell:: eBdryEdge) || (pC->eType () == Cell::eIntBdryEdge));
  BdryEdgeBase *pBE_Final = dynamic_cast<BdryEdgeBase *>(pC);

  // Check whether it's okay to smooth this vertex.
  //   Are the boundary conditions the same?
  if (pBE_Init->iBdryCond () != pBE_Final->iBdryCond ())
    goto Return;
  //   Are the faces co-linear?  If so, the unit normals will have a
  //   dot product which has magnitude of almost precisely 1.
  {
    double adNorm_Init[3], adNorm_Final[3];
    pFBdry->vUnitNormal (adNorm_Init);
    pF->vUnitNormal (adNorm_Final);
    double dDot = fabs (dDOT2D (adNorm_Init, adNorm_Final));
    if (iFuzzyComp (dDot, 1.) != 0)
      goto Return;
    // Co-linear!  Set the surface normal accordingly.  Orientation of
    // the normal (inward vs. outward) doesn't matter here.
    SMsetNormal (adNorm_Init);
  }

  // Now that we know we're going to smooth, set up the connectivity,
  // etc and get on with it.

  // List all opposite faces (w/ correct orientation) using indices
  // into the location array. There are iNeigh-1 such faces.
  if ((pFBdry->pCCellRight () == static_cast<Cell *>(pBE_Init) &&
       pFBdry->pVVert (1) == pVBdry)
      ||
      (pFBdry->pCCellLeft () == static_cast<Cell *>(pBE_Init) &&
       pFBdry->pVVert (0) == pVBdry)) {
    for (int ii = 0; ii < iNeigh - 1; ii++) {
      a2iFaceConn[ii][0] = ii + 1;
      a2iFaceConn[ii][1] = ii;
    }
  }
  else {
    for (int ii = 0; ii < iNeigh - 1; ii++) {
      a2iFaceConn[ii][0] = ii;
      a2iFaceConn[ii][1] = ii + 1;
    }
  }

  if (iNeigh <= iMaxVertNeigh) {
    double adCurrPt[] =
    {pVBdry->dX (), pVBdry->dY ()};
    double adLoc[] =
    {adCurrPt[0], adCurrPt[1]};
    SMsmooth (iNeigh, iNeigh - 1, adCurrPt, a2dNeighPts, a2iFaceConn,
	      pvSmoothData, 1);
    if (dDIST2D (adLoc, adCurrPt) > 1.e-8) {
      pVBdry->vSetCoords (2, adCurrPt);
    }

    // Only mark this vertex and its neighborhood for future
    // improvement if the worst angle in the submesh is less than 40
    // degrees or more than 100 degrees.
    double dQualValue =
      static_cast<SMsmooth_data *>(pvSmoothData)->local_mesh->current_active_value;
    SMconvertToDegrees (static_cast<SMsmooth_data *>(pvSmoothData)->smooth_param
			->function_id, &dQualValue);
    if (dQualValue < 40 || dQualValue > 100) {
      pVBdry->vMarkIllShaped ();
      // Also, mark all incident edges as dirty.
      for (int ii = 0; ii < iNeigh; ii++)
	apFIncFaces[ii]->vMarkForSwap ();
    }
    else {
      pVBdry->vMarkWellShaped ();
    }
  }
 Return:
  for (i = 0; i < iMaxVertNeigh; i++) {
    delete[]a2dNeighPts[i];
    delete[]a2iFaceConn[i];
  }
  delete[]a2dNeighPts;
  delete[]a2iFaceConn;
  delete[]apFIncFaces;
  return (qRetVal);
}

//@ Optimization driver for 2D mesh smoothing
int Mesh2D::
iSmooth (const int iMaxPasses)
{
  assert (qSimplicial ());
  int iSmooths = 1;
  int iTotal = 0;
  int iPasses = 0;

  // Set up smoothing stuff that needs doing only once
  Face **apFIncFaces = new Face *[iMaxVertNeigh];
  double **a2dNeighPts = new double *[iMaxVertNeigh];
  int **a2iFaceConn = new int *[iMaxVertNeigh];
  int i;
  for (i = 0; i < iMaxVertNeigh; i++) {
    a2dNeighPts[i] = new double[2];
    a2iFaceConn[i] = new int[2];
  }

  while (iSmooths && iPasses < iMaxPasses) {
    double adCurrPt[2];

    SMinitSmoothStats (pvSmoothData);
    ++iPasses;
    iSmooths = 0;
    for (i = 0; i < iNumVerts (); i++) {
      // Set up coordinates of the current vertex
      Vert *pV = pVVert (i);
      // If this vertex is happy where it is, don't bother trying to
      // smooth it.  Also, if the vert is in a structured part of the
      // mesh, don't touch it.
      if (pV->qWellShaped () || pV->qIsStructured ())
	continue;

      adCurrPt[0] = pV->dX ();
      adCurrPt[1] = pV->dY ();

      bool qBdryVert = false;
      int iNeigh = 0;
      // Don't try to smooth the vertex if it's been removed from the
      // mesh, and don't bother if the neighborhood is known to be
      // well-shaped already.
      if (pV->qDeleted () || pV->qWellShaped ())
	continue;
      Face *pF = pV->pFHintFace ();
      Cell *pC;
      if (pF->pVVert (0) == pV)
	pC = pF->pCCellLeft ();
      else
	pC = pF->pCCellRight ();

      Face *pFStart = pF;
      do {
	if (pC->eType () == Cell::eTriCell && iNeigh < iMaxVertNeigh) {
	  Vert *pVOpp;
	  if (pF->pVVert (0) == pV)
	    pVOpp = pF->pVVert (1);
	  else
	    pVOpp = pF->pVVert (0);
	  a2dNeighPts[iNeigh][0] = pVOpp->dX ();
	  a2dNeighPts[iNeigh][1] = pVOpp->dY ();
	  apFIncFaces[iNeigh] = pF;
	  iNeigh++;
	}
	else {
	  qBdryVert = true;
	  break;
	}

	pC = pF->pCCellOpposite (pC);
	if (pC->eType () == Cell::eTriCell) {
	  for (int iF = 0; iF < 3; iF++) {
	    Face *pFEdge = pC->pFFace (iF);
	    assert (pFEdge->qValid ());
	    if (pF == pFEdge)
	      continue;
	    if (pFEdge->pVVert (0) == pV ||
		pFEdge->pVVert (1) == pV) {
	      pF = pFEdge;
	      break;
	    }
	  }
	}
	else {
	  qBdryVert = true;
	}
      }
      while (pF != pFStart && !qBdryVert);

      // Boundary vertex smoothing is handled in a separate routine.
      if (pV->iVertType () == Vert::eBdryTwoSide) {
	// Version 0.1.8: Implement this
	continue;
      }
      else if (pV->iVertType () == Vert::eBdryApex ||
	       pV->iVertType () == Vert::eInteriorFixed) {
	// Can't smooth these
	continue;
      }
      else if (qBdryVert) {
	// Note by Charles Boivin - 21/09/99
	// This seems to be reached for vertices on internal boundaries too
	// Made sure the assertion was not triggered for no reason
	assert ((pC->eType () == Cell:: eBdryEdge) ||
		(pC->eType () == Cell::eIntBdryEdge));
	if (pC->eType() == Cell::eBdryEdge) {
	  BdryEdge *pBE = dynamic_cast<BdryEdge *>(pC);
	  if (qBdrySmooth (pV, pF, pBE))
	    iSmooths++;
	  continue;
	}
	else {
	  // Can't currently smooth points on internal bdry edges.
	}
      }

      // List all opposite faces (w/ correct orientation) using indices
      // into the location array
      a2iFaceConn[0][0] = 0;
      a2iFaceConn[0][1] = iNeigh - 1;
      int ii;
      for (ii = 1; ii < iNeigh; ii++) {
	a2iFaceConn[ii][0] = ii;
	a2iFaceConn[ii][1] = ii - 1;
      }

      if (iNeigh <= iMaxVertNeigh) {
	double adLoc[] =
	{adCurrPt[0], adCurrPt[1]};
	SMsmooth (iNeigh, iNeigh, adCurrPt, a2dNeighPts, a2iFaceConn,
		  pvSmoothData, 0);
	if (dDIST2D (adLoc, adCurrPt) > 1.e-8) {
	  pV->vSetCoords (2, adCurrPt);
	  iSmooths++;
	}

	// Only mark this vertex and its neighborhood for future
	// improvement if the worst angle in the submesh is less than 30
	// degrees or more than 150 degrees.
	double dQualValue =
	  static_cast<SMsmooth_data *>(pvSmoothData)->local_mesh->current_active_value;
	SMconvertToDegrees (static_cast<SMsmooth_data *>(pvSmoothData)
			    ->smooth_param->function_id, &dQualValue);
	if (dQualValue < 40 || dQualValue > 100) {
	  pV->vMarkIllShaped ();
	  // Also, mark all incident edges as dirty.
	  for (ii = 0; ii < iNeigh; ii++)
	    apFIncFaces[ii]->vMarkForSwap ();
	}
	else {
	  pV->vMarkWellShaped ();
	}
      }
    }
    SMprintSmoothStats (pvSmoothData);
    vResetThreshold();

    iTotal += iSmooths;
    vMessage (1, "Pass %d:  %5d vertices moved by smoothing, %6d total\n",
	      iPasses, iSmooths, iTotal);

  }
  for (i = 0; i < iMaxVertNeigh; i++) {
    delete[]a2dNeighPts[i];
    delete[]a2iFaceConn[i];
  }
  delete[]a2dNeighPts;
  delete[]a2iFaceConn;
  delete[]apFIncFaces;
  vMessage (1, "Vertex smoothing complete.\n");
  return iTotal;
}

//@ Optimization driver for face swapping.  Swap measure is set elsewhere.
int Mesh2D::
iSwap (const int iMaxPasses, const bool qAlreadyMarked)
{
  SUMAA_LOG_EVENT_BEGIN (SWAP_PASSES);
  SUMAA_LOG_EVENT_BEGIN (NET_SWAPPING);
  int iSwaps = 1;
  assert (eSwapMeasure != eNone);

  if (!qAlreadyMarked) {
    for (int iFace = 0; iFace < iNumFaces (); iFace++)
      pFFace (iFace)->vMarkForSwap ();
  }

  int iTotal = 0;
#ifndef NDEBUG
  int iPasses = 0;

  while (iSwaps && iPasses < iMaxPasses) {
    ++iPasses;
    if (iMaxPasses > 1)
      vMessage (2, "   Pass %d:", iPasses);
#endif
    iSwaps = 0;
    for (int i = 0; i < iNumFaces (); i++) {
      // Always swap recursively
      iSwaps += iFaceSwap (pFFace (i));
    }
    iTotal += iSwaps;
    vMessage (2, "%6d swaps", iSwaps);
#ifndef NDEBUG
    if (iMaxPasses > 1)
      vMessage (2, ",%6d total", iTotal);
#endif
    vMessage (2, "\n");
#ifndef NDEBUG
  }
  if (iMaxPasses > 1)
    vMessage (2, "Face swapping complete.\n");
#endif

  SUMAA_LOG_EVENT_END (NET_SWAPPING);
  // No need to purge faces and cells after 2D swapping
  vSetAllHintFaces ();
  SUMAA_LOG_EVENT_END (SWAP_PASSES);
  return iTotal;
}

bool Mesh2D::qDelaunayize()
{
  assert(qSimplicial());
  // Do incirclesphere face swapping; this gives a Delaunay mesh without
  // further ado.
  vAllowSwapRecursion();
  vSetSwapType(eDelaunay);
  iSwap();
  return true;
}
