#include <string.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <map>
using std::multimap;
using std::pair;

#include "GR_Bdry2D.h"
#include "GR_BdryPatch2D.h"
#include "GR_Mesh2D.h"
#include "GR_CellCV.h"
#include "GR_Geometry.h"
#include "GR_InsertionQueue.h"
#include "GR_Vec.h"
#include "GR_misc.h"
#include "GR_List.h"
#include "GR_events.h"
#include "GR_BFace.h"
#include "GR_BFaceCV.h"
#include "GR_Sort.h"

#include "OptMS.h"
#include "SMsmooth.h"

#ifndef WIN_VC32
#include <unistd.h>
#else
#pragma warning(disable:4786) // DEBUG symbol > 255 chars
#endif

// Someplace to stash the info we need to correctly disconnect the
// exterior part of the mesh and replace it with boundary conditions.
struct FaceInfo {
  Face *pF;
  BFace *pBF;
  Cell *pC;
  FaceInfo () :
    pF(pFInvalidFace),
    pBF(pBFInvalidBFace),
    pC(pCInvalidCell) {}
};

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

//@ Default constructor
Mesh2D::Mesh2D(const int iQualMeas)
  : Mesh (), VAEdgeF(), VAEdgeBF(), VATri(), VAQuad(),
    VAIntEdgeBF(), iOnEdge(0), iOnBdryEdge(0), iOnBdryEdge2(0),
    iInterior(0), iInternalBdry(0), iInternalBdry2(0), B2D()
{
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(2, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);
}

static void
vMarkCleanNeighbors (Cell * pC)
{
  int iReg = pC->iRegion ();
  assert (iReg != iInvalidRegion);
  // Cells in the exterior do not get marked with a region.  Instead,
  // they retain their invalid region tag and no traversal is done.
  for (int iF = 0; iF < pC->iNumFaces (); iF++) {
    Face *pF = pC->pFFace (iF);
    if (pF->qValid ()) {
      Cell *pCCand = pF->pCCellOpposite (pC);
      assert (pCCand->qValid ());
      if (pCCand->eType () == Cell::eBdryEdge)
	continue;
      if (pCCand->iRegion () == iInvalidRegion) {
	// If the opposite cell hasn't been marked, mark it the same as
	// this one.
	pCCand->vSetRegion (iReg);
	vMarkCleanNeighbors (pCCand);
      }				// Done tagging candidate cell

    }				// Face is OK

  }				// Loop over faces

}

static Face *
pFFaceBetween (Vert * const pV0, Vert * const pV1)
{
  // Find all the faces in cells incident on pV0
  std::set<Cell*> spCInc;
  std::set<Vert*> spVNeigh;
  std::set<BFace*> spBFInc;
  std::set<Face*> spFNearby;
  bool qIsBdryVert;

  vNeighborhood(pV0, spCInc, spVNeigh, &spBFInc, &qIsBdryVert, &spFNearby);

  // Find the one face that is incident on both pV0 and pV1
  std::set<Face*>::iterator iterF;
  for (iterF = spFNearby.begin(); iterF != spFNearby.end(); iterF++) {
    Face *pFCand = *iterF;
    if (pFCand->qHasVert (pV0) && pFCand->qHasVert (pV1))
      return pFCand;
  }
  // If we got to here, then the two verts are not in fact adjacent in
  // the mesh.
  return pFInvalidFace;
}


//@ Delaunay triangulation of a point set, with optional recovery of perimeter
Mesh2D::Mesh2D (const VAVertex & VA, const int iQualMeas,
		const int iNumBdryEdges, const int (*const a2iConnectList)[2])
  : Mesh(VA), VAEdgeF(), VAEdgeBF(), VATri(),
    VAQuad(), VAIntEdgeBF(), iOnEdge(0), iOnBdryEdge(0), iOnBdryEdge2(0),
    iInterior(0), iInternalBdry(0), iInternalBdry2(0), B2D()
{
  SUMAA_LOG_EVENT_BEGIN (INIT_TRI);
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(2, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);

  // Figure out which vertices are really used in the 2D mesh
  bool *aqVertUsed = new bool[iNumVerts ()];
  int iVert, iBdryEdge;
  if (iNumBdryEdges > 0) {
    for (iVert = 0; iVert < iNumVerts (); iVert++)
      aqVertUsed[iVert] = false;
    for (iBdryEdge = 0; iBdryEdge < iNumBdryEdges; iBdryEdge++) {
      aqVertUsed[a2iConnectList[iBdryEdge][0]] = true;
      aqVertUsed[a2iConnectList[iBdryEdge][1]] = true;
    }
  }
  else {
    // In this case, triangulate them all.
    for (iVert = 0; iVert < iNumVerts (); iVert++)
      aqVertUsed[iVert] = true;
  }
  
  // Set up an outer triangulation of a big square
  double dXMin, dXMax, dYMin, dYMax;
  dXMin = dYMin = 1.e100;
  dXMax = dYMax = -1.e100;

  for (iVert = 0; iVert < iNumVerts (); iVert++) {
    if (aqVertUsed[iVert]) {
      double dX = pVVert (iVert)->dX ();
      double dY = pVVert (iVert)->dY ();
      dXMin = min (dXMin, dX);
      dXMax = max (dXMax, dX);
      dYMin = min (dYMin, dY);
      dYMax = max (dYMax, dY);
    }
  }

  double dXLo = dXMin - (dXMax - dXMin);
  double dXHi = dXMax + (dXMax - dXMin);
  double dYLo = dYMin - (dYMax - dYMin);
  double dYHi = dYMax + (dYMax - dYMin);

  int iNumOrigVerts = iNumVerts ();

  Vert *pVLoLeft = pVNewVert ();
  Vert *pVLoRight = pVNewVert ();
  Vert *pVHiLeft = pVNewVert ();
  Vert *pVHiRight = pVNewVert ();

  double adLoc[2];
  adLoc[0] = dXLo;
  adLoc[1] = dYLo;
  pVLoLeft->vSetCoords (2, adLoc);
  adLoc[1] = dYHi;
  pVHiLeft->vSetCoords (2, adLoc);
  adLoc[0] = dXHi;
  pVHiRight->vSetCoords (2, adLoc);
  adLoc[1] = dYLo;
  pVLoRight->vSetCoords (2, adLoc);

  // Set up initial mesh
  //
  //    HL-------0--------HR
  //     |              _/|
  //     |     0     __/  |
  //     |         _/     |
  //     1      _4/       3
  //     |    _/     1    |
  //     | __/            |
  //     |/               |
  //    LL-------2--------LR
  //

  vSetupFaces (5);
  vSetupBFaces (4);
  vSetupTriCells (2);
  pFFace (0)->vAssign (pCCell (0), pBFBFace (0), pVHiRight, pVHiLeft);
  pBFBFace (0)->vAssign (pFFace (0));
  pFFace(0)->vSetFaceLoc(Face::eBdryFace);

  pFFace (1)->vAssign (pCCell (0), pBFBFace (1), pVHiLeft, pVLoLeft);
  pBFBFace (1)->vAssign (pFFace (1));
  pFFace(1)->vSetFaceLoc(Face::eBdryFace);

  pFFace (2)->vAssign (pCCell (1), pBFBFace (2), pVLoLeft, pVLoRight);
  pBFBFace (2)->vAssign (pFFace (2));
  pFFace(2)->vSetFaceLoc(Face::eBdryFace);

  pFFace (3)->vAssign (pCCell (1), pBFBFace (3), pVLoRight, pVHiRight);
  pBFBFace (3)->vAssign (pFFace (3));
  pFFace(3)->vSetFaceLoc(Face::eBdryFace);

  pFFace (4)->vAssign (pCCell (0), pCCell (1), pVLoLeft, pVHiRight);

  pCCell (0)->vAssign (pFFace (0), pFFace (1), pFFace (4));
  pCCell (1)->vAssign (pFFace (2), pFFace (3), pFFace (4));

  vSetAllHintFaces ();
  assert (qValid ());

  // Do Delaunay insertion of all the existing verts in the mesh
  vSetSwapType (eDelaunay);
  int iNumSwaps = 0;
  for (iVert = 0; iVert < iNumOrigVerts; iVert++) {
    Vert *pV = pVVert (iVert);
    if (aqVertUsed[iVert]) {
      double adCoords[2];
      adCoords[0] = pV->dX ();
      adCoords[1] = pV->dY ();
      int iNewSwaps;
      // Insert, then swap.  Force insertion at the specified point and
      // use the existing vertex location.
      SUMAA_LOG_EVENT_BEGIN (INSERT_VERTEX);
#ifndef NDEBUG
      bool qSuccess = qInsertPoint (adCoords, pCCell (0), &iNewSwaps, true,
				    true, pV);
      assert (qSuccess);
#else
      (void) qInsertPoint (adCoords, pCCell (0), &iNewSwaps, true, true, pV);
#endif
      SUMAA_LOG_EVENT_END (INSERT_VERTEX);
      iNumSwaps += iNewSwaps;
    }
  }
  delete[]aqVertUsed;
  vSetAllHintFaces ();
  assert (qValid ());

  // Force swapping to recover the requested perimeter
  assert ((iNumBdryEdges == 0 && a2iConnectList == NULL) ||
	  (iNumBdryEdges != 0 && a2iConnectList != NULL));

  if (iNumBdryEdges == 0) {
    // No perimeter specified
    // Delete those four outer points and anything incident on them
    pVLoLeft->vMarkDeleted ();
    pVLoRight->vMarkDeleted ();
    pVHiLeft->vMarkDeleted ();
    pVHiRight->vMarkDeleted ();

    pBFBFace (0)->vMarkDeleted ();
    pBFBFace (1)->vMarkDeleted ();
    pBFBFace (2)->vMarkDeleted ();
    pBFBFace (3)->vMarkDeleted ();
    // Do this the easy way, in which every cell is checked for incidence
    // on one of the deleted verts, then local marking is done.
    int iCell, iFace;
    for (iCell = 0; iCell < iNumCells (); iCell++) {
      Cell *pC = pCCell (iCell);
      if (pC->pVVert (0)->qDeleted () ||
	  pC->pVVert (1)->qDeleted () ||
	  pC->pVVert (2)->qDeleted ()) {
	pC->vMarkDeleted ();
	// Find the face that is -not- incident on a deleted vert.  Set it
	// up as a new boundary vert.
	for (iFace = 0; iFace < 3; iFace++) {
	  Face *pF = pC->pFFace (iFace);
	  if (!pF->pVVert (0)->qDeleted () &&
	      !pF->pVVert (1)->qDeleted ()) {
	    BFace *pBF = pBFNewBFace ();
	    pBF->vAssign (pF);
	    pF->vReplaceCell (pC, pBF);
	  }
	}
      }
    }

    for (iFace = 0; iFace < iNumFaces (); iFace++) {
      Face *pF = pFFace (iFace);
      if (pF->pVVert (0)->qDeleted () ||
	  pF->pVVert (1)->qDeleted ())
	pF->vMarkDeleted ();
    }
  }				// Done with the case of no perimeter

  else {
    int iEdge;
    // Recover all edges; this is trivial
    for (iEdge = 0; iEdge < iNumBdryEdges; iEdge++) {
      int iV0 = a2iConnectList[iEdge][0];
      int iV1 = a2iConnectList[iEdge][1];
#ifndef NDEBUG
      bool qOK = qRecoverEdge (iV0, iV1);
      assert (qOK);
#else
      (void) qRecoverEdge (iV0, iV1);
#endif
    }

    // Mark all cells as having an invalid region, since we don't yet know
    // where any of them are.
    int iCell;
    for (iCell = 0; iCell < iNumCells (); iCell++)
      pCCell (iCell)->vSetRegion (iInvalidRegion);

    // Someplace to stash the info we need to correctly disconnect the
    // exterior part of the mesh and replace it with boundary conditions.
    struct FaceInfo *aFI = new struct FaceInfo[iNumBdryEdges];

    // For each bdry edge:
    //   1. Mark cells that are in the interior of the domain to be meshed
    //   with their region.
    //   2. Mark cells that are in the exterior for deletion.
    for (iEdge = 0; iEdge < iNumBdryEdges; iEdge++) {
      int iV0 = a2iConnectList[iEdge][0];
      int iV1 = a2iConnectList[iEdge][1];
      Vert *pV0 = pVVert (iV0);
      Vert *pV1 = pVVert (iV1);

      // Find the edge that connects these two vertices.
      Face *pF = pFFaceBetween (pV0, pV1);
      assert (pF->qValid ());
      assert (pF->pVVert (0) == pV0 || pF->pVVert (0) == pV1);
      assert (pF->pVVert (1) == pV0 || pF->pVVert (1) == pV1);

      // Mark the interior cells
      bool qSameOrient = (pF->pVVert (0) == pV0);

      Cell *pC;
      if (qSameOrient) {
	pF->pCCellLeft ()->vSetRegion (1);
	pC = pF->pCCellRight ();
      }
      else {
	pF->pCCellRight ()->vSetRegion (1);
	pC = pF->pCCellLeft ();
      }

      // Set up a new bdry face
      BFace *pBF = pBFNewBFace ();
      pBF->vAssign (pF);

      pF->vSetFaceLoc (Face::eBdryFace);

      aFI[iEdge].pC = pC;
      aFI[iEdge].pF = pF;
      aFI[iEdge].pBF = pBF;
    }

    // For every edge separating interior and exterior, replace the
    // exterior cell with a boundary face (which has a pointer to the
    // correct boundary patch so that geometric info is not lost).
    for (iEdge = 0; iEdge < iNumBdryEdges; iEdge++) {
      // Only do this stuff if this is -really- a boundary face.
      if (aFI[iEdge].pF->qValid ()) {
	Cell *pC = aFI[iEdge].pC;
	Face *pF = aFI[iEdge].pF;
	BFace *pBF = aFI[iEdge].pBF;

	// Disconnect the old cell.
	pC->vRemoveFace (pF);

	// Replace the old cell with the new bdry face in the face
	// connectivity.
	pF->vReplaceCell (pC, pBF);

	// And that's it.  The rest has already been taken care of.
      }
    }
    delete[]aFI;

    // March out from marked cells, filling all unmarked neighbors with
    // correct region tag (interior) or marking for deletion (exterior) in
    // a second pass.
    for (iCell = 0; iCell < iNumCells (); iCell++) {
      Cell *pC = pCCell (iCell);
      // Can't march out from cells that don't know where they are, and
      // don't need to march out from cells in the exterior.
      if (pC->iRegion () == iInvalidRegion)
	continue;
      vMarkCleanNeighbors (pC);
    }

    // Mark exterior cells for deletion.
    for (iCell = 0; iCell < iNumCells (); iCell++) {
      Cell *pC = pCCell (iCell);
      if (pC->iRegion () == iInvalidRegion)
	pC->vMarkDeleted ();
    }

    // Mark four farfield vertices for deletion.
    for (iVert = iNumOrigVerts; iVert < iNumVerts (); iVert++)
      pVVert (iVert)->vMarkDeleted ();

    // Mark farfield boundary faces for deletion.  These are always the
    // first four BFaces.
    int iBFace;
    for (iBFace = 0; iBFace < 4; iBFace++)
      pBFBFace (iBFace)->vMarkDeleted ();

    // Mark all faces bounded by two deleted cells for deletion. This must
    // come -after- bdry face deletion so that the farfield faces are
    // deleted properly.
    int iFace;
    for (iFace = 0; iFace < iNumFaces (); iFace++) {
      Face *pF = pFFace (iFace);
      if (pF->pCCellLeft ()->qDeleted () &&
	  pF->pCCellRight ()->qDeleted ())
	pF->vMarkDeleted ();
    }
  }

  // Clean up and exit
  vPurge ();
  assert (qValid ());
  SUMAA_LOG_EVENT_END (INIT_TRI);
}

void Mesh2D::
vCheckPatches ()
{
  // Used to check if the patches are ok...
  BdryPatch2D *pPatch;
  double adCoords1[2], adCoords2[2];

  fprintf (stdout, "Now checking the patches...\n-------------------------------------\n");
  fprintf (stdout, "Number of patches: %d\n", B2D.iNumPatches ());
  fprintf (stdout, "Number of verts: %d\n", B2D.iNumPoints ());
  /*
    for (i = 0; i < B2D.iNumBdryPatches (); i++) {
    fprintf (stdout, "\nNow accessing patch # %d\n", i);
    pV1 = B2D.pBdry (i)->pVBeginVert ();
    pV2 = B2D.pBdry (i)->pVEndVert ();
    fprintf (stdout, "   Begin Vertex coords are: (%f, %f)\n", pV1->dX (), pV1->dY ());
    fprintf (stdout, "   This is the index of this vertex: %d\n", B2D.iBeginVert (i));
    fprintf (stdout, "   End Vertex coords are: (%f, %f)\n", pV2->dX (), pV2->dY ());
    fprintf (stdout, "   This is the index of this vertex: %d\n", B2D.iEndVert (i));
    fprintf (stdout, "   ArcLength is %f\n", B2D.pBdry (i)->dArcLength ());
    fprintf (stdout, "   The BC is on side: %d\n", B2D.pBdry (i)->iBCWhichSide ());
    fprintf (stdout, "   And it has this value: %d\n", B2D.pBdry (i)->iBdryCond ());
    fprintf (stdout, "   The region is on side: %d\n", B2D.pBdry (i)->iRegionWhichSide ());
    B2D.pBdry (i)->vPointAtParameterRatio (0.5, adCoords);
    fprintf (stdout, "   This would be the midpoint: (%f, %f)\n\n", adCoords[0], adCoords[1]);
    }
  */
  for (int iBF = iNumBdryFaces () - 1; iBF >= 0; iBF--) {
    BFace *pBF = pBFBFace (iBF);
    Vert *pV1 = pBF->pVVert (0);
    Vert *pV2 = pBF->pVVert (1);


    pPatch = B2D.pBdry (dynamic_cast<BdryEdgeBase*>(pBF)->iPatchIndex ());

    fprintf (stdout, "Now checking BFace # %d\n", iBF);
    //        fprintf(stdout, "This is the index to the bpatch: %d\n", ((BdryEdge*)pBF)->iPatchIndex());
    adCoords1[0] = pV1->dX ();
    adCoords1[1] = pV1->dY ();
    //        fprintf(stdout, "This is vertex coords 1: (%f, %f)\n", adCoords1[0], adCoords1[1]);
    adCoords2[0] = pV2->dX ();
    adCoords2[1] = pV2->dY ();
    //        fprintf(stdout, "This is vertex coords 2: (%f, %f)\n", adCoords2[0], adCoords2[1]);
    //        fprintf(stdout, "Patch coords 1: (%f, %f)\n",
    //                        B2D.pBdry(((BdryEdge*)pBF)->iPatchIndex())->pVBeginVert()->dX(),
    //                        B2D.pBdry(((BdryEdge*)pBF)->iPatchIndex())->pVBeginVert()->dY() );
    //        fprintf(stdout, "Patch coords 2: (%f, %f)\n",
    //                        B2D.pBdry(((BdryEdge*)pBF)->iPatchIndex())->pVEndVert()->dX(),
    //                        B2D.pBdry(((BdryEdge*)pBF)->iPatchIndex())->pVEndVert()->dY() );
    fprintf (stdout, "BFace coords: (%f, %f) - (%f, %f)\n",
	     adCoords1[0], adCoords1[1],
	     adCoords2[0], adCoords2[1]);
    fprintf (stdout, "Arclengths 1: %f, 2: %f, Total: %f\n",
	     pPatch->dParameterAtPoint (adCoords1),
	     pPatch->dParameterAtPoint (adCoords2),
	     pPatch->dMaxParameter()
	     );
  }


  fprintf (stdout, "\nDone checking the patches\n----------------------------------------\n\n");
  //    scanf("%c\n", &cDummy);
}


void Mesh2D::
vGetBdry2DInfo (const Bdry2D & B2DInfo)
{
  // Copy the Bdry2D info in the internal structure...
  int i;
  int iNumSegs = 0;
  BdryPatch2D *pPatch;
  structVertConnect sVConnect;
  double adC[2];

  vMessage (4, "Copying the Bdry2D information\n");

  // Initialize size..
  //    B2D.vSetSize(B2DInfo.iNumBdryPatches(), B2DInfo.iNumVerts());
  // Should set the size of the num of patches to the num of patches used!
  //  vMessage (3, "Initializing the new Bbdry2D class with %d patches and %d verts\n",
  //	    B2DInfo.iNumBdryPatchesUsed (), B2DInfo.iNumPoints ());
  int iNumPatch = B2DInfo.iNumPatchesUsed();
  int iNumPts = B2DInfo.iNumPoints();
  int iNumPer = B2DInfo.iNumPeriodic();
  B2D.vSetSize (iNumPatch, iNumPts, iNumPer);

  // Add each vertex...
  for (i = 0; i < iNumPts; i++) {
    //    vMessage(2, "   Adding vertex # %d ...", i);
    B2DInfo.vGetPoint(i, adC);
    B2D.vAddPoint(adC);
    B2D.vSetBdryPoint (i, B2DInfo.qIsBdryPoint (i));
    B2D.vSetNumPatchesConnectedToPoint(i, B2DInfo.iNumPatchesConnectedToPoint(i));
    //    vMessage(2, " (%f, %f) - Used: %d\n", adC[0], adC[1], B2D.qIsBdryPoint (i));
  }

  // Copy vert connectivity stuff too
  for (i = 0; i < B2DInfo.iConnectivityListLength(); i++) {
    sVConnect = B2DInfo.sVCConnectivityListContent(i);
    structVertConnect sVC2;
    //		sVC2 = new structVertConnect;
    sVC2.iVert  = sVConnect.iVert;
    sVC2.iPatch = sVConnect.iPatch;
    sVC2.dAngle = sVConnect.dAngle;
    B2D.vAddToConnectivityList(sVC2);
  }

  // Copy small angles info
  for (i = 0; i < iNumPts; i++) {
    B2D.vSetSmallAngle(i, B2DInfo.qVertSmallAngle(i));
    int iStart = B2DInfo.iStartVert(i);
    B2D.vSetStartVert(i, iStart);
  }

  // Then add each patch...
  for (i = 0; i < B2DInfo.iNumPatches (); i++) {
    // Must make sure we only copy USED patches...
    //    if (B2DInfo.qIsPatchUsed (i)) {
    pPatch = B2DInfo.pBdry (i);
    // Should really check the type of patch first..!!!!!
    switch (pPatch->eGeom2DType) {
    case BdryPatch2D::ePolyLine:
      {
	vMessage(4, "   Creating BdrySeg #%d (was %d before) ...", iNumSegs, i);

	BdryPatch2D *pBChk = new BdrySeg (&B2D, pPatch->iLeftBdryCond (),
					  pPatch->iRightBdryCond (),
					  pPatch->iLeftRegion (),
					  pPatch->iRightRegion (),
					  pPatch->iBeginPointIndex (),
					  pPatch->iEndPointIndex ()
					  );
	B2D.vAddPatchToList (pBChk);
	iNumSegs++;
	/*
	  B2D.vAddPatchToList(BdryPatch2D::ePolyLine,
	  B2DInfo.pBdry(i)->iLeftBdryCond(),
	  B2DInfo.pBdry(i)->iRightBdryCond(),
	  B2DInfo.pBdry(i)->iLeftRegion(),
	  B2DInfo.pBdry(i)->iRightRegion(),
	  B2DInfo.iBeginVert(i),
	  B2DInfo.iEndVert(i)
	  );
	*/
	vMessage(4, " and added\n");
      }
      break;
    case BdryPatch2D::eCircle:
      {
				// For circles...
	vMessage(4, "   Creating circle BdryArc #%d (was %d before) ...", iNumSegs, i);

	dynamic_cast<BdryArc*>(pPatch)->vGetCenter(adC);

	BdryPatch2D *pBChk = new BdryArc (&B2D,
					  pPatch->iLeftBdryCond (),
					  pPatch->iRightBdryCond (),
					  pPatch->iLeftRegion (),
					  pPatch->iRightRegion (),
					  dynamic_cast<BdryArc*>(pPatch)->dGetRadius(),
					  dynamic_cast<BdryArc*>(pPatch)->iCenterPointIndex ()
					  );
	B2D.vAddPatchToList (pBChk);
	iNumSegs++;

	vMessage(4, " and added\n");
      }
      break;
    case BdryPatch2D::eCircularArc:
      {
				// For arcs....
	vMessage(4, "   Creating arc BdryArc #%d (was %d before) ...", iNumSegs, i);

	dynamic_cast<BdryArc*>(pPatch)->vGetCenter(adC);

	BdryPatch2D *pBChk = new BdryArc ( &B2D,
					   pPatch->iLeftBdryCond (),
					   pPatch->iRightBdryCond (),
					   pPatch->iLeftRegion (),
					   pPatch->iRightRegion (),
					   pPatch->iBeginPointIndex (),
					   pPatch->iEndPointIndex (),
					   dynamic_cast<BdryArc*>(pPatch)->dGetRadius()
					   );
	B2D.vAddPatchToList (pBChk);
	iNumSegs++;

	vMessage(4, " and added\n");
      }
      break;
    case BdryPatch2D::eLongArc:
      {
				// For arcs....
	vMessage(4, "   Creating long arc BdryArc #%d (was %d before) ...", iNumSegs, i);

	dynamic_cast<BdryArc*>(pPatch)->vGetCenter(adC);

	BdryPatch2D *pBChk = new BdryArc ( &B2D,
					   pPatch->iLeftBdryCond (),
					   pPatch->iRightBdryCond (),
					   pPatch->iLeftRegion (),
					   pPatch->iRightRegion (),
					   pPatch->iBeginPointIndex (),
					   pPatch->iEndPointIndex (),
					   dynamic_cast<BdryArc*>(pPatch)->dGetRadius(),
					   true
					   );
	B2D.vAddPatchToList (pBChk);
	iNumSegs++;

	vMessage(4, " and added\n");
      }
      break;
    case BdryPatch2D::eCubicParam:
      {
				// Cubic parametric curves
	vMessage(4, "   Creating cubic param curve #%d (was %d before) ...", iNumSegs, i);

	double dax, dbx, dcx, ddx, day, dby, dcy, ddy;
				// Get the parameters
	dynamic_cast<BdryCubicParam2D*>(pPatch)->vGetAllCoefficients(&dax, &dbx, &dcx, &ddx,
							 &day, &dby, &dcy, &ddy);
				// Create the new patch
	BdryPatch2D *pBChk = new BdryCubicParam2D(
						  &B2D,
						  pPatch->iLeftBdryCond (),
						  pPatch->iRightBdryCond (),
						  pPatch->iLeftRegion (),
						  pPatch->iRightRegion (),
						  pPatch->iBeginPointIndex (),
						  pPatch->iEndPointIndex (),
						  dax, dbx, dcx, ddx, day, dby, dcy, ddy
						  );

	B2D.vAddPatchToList(pBChk);
	iNumSegs++;


	vMessage(4, " and added\n");
      }
      break;

    case BdryPatch2D::eBezier:
      {
	vMessage(4, "   Creating bezier curve #%d (was %d before) ...", iNumSegs, i);

				// Create the new patch
	BdryPatch2D *pBChk = new BdryBezier2D(
					      &B2D,
					      pPatch->iLeftBdryCond (),
					      pPatch->iRightBdryCond (),
					      pPatch->iLeftRegion (),
					      pPatch->iRightRegion (),
					      pPatch->iBeginPointIndex (),
					      pPatch->iEndPointIndex (),
					      dynamic_cast<BdryBezier2D*>(pPatch)->iControlPointIndex(1),
					      dynamic_cast<BdryBezier2D*>(pPatch)->iControlPointIndex(2)
					      );

	B2D.vAddPatchToList(pBChk);
	iNumSegs++;


	vMessage(4, " and added\n");

      }

      break;
    case BdryPatch2D::eInterpolated:
      {
	vMessage(4, "   Creating interpolated spline #%d (was %d before) ...", iNumSegs, i);

				// Get all the points in the list
	List<int> LdPts;
	for (int ii = 0;
	     ii < dynamic_cast<BdryInterpolated2D*>(pPatch)->iNumPoints(); ii++) {
	  int iDum = dynamic_cast<BdryInterpolated2D*>(pPatch)->iGetPoint(ii);
	  LdPts.vAppendItem( iDum );
	}

				// Create the new patch
	BdryPatch2D *pBChk = new BdryInterpolated2D(
						    &B2D,
						    pPatch->iLeftBdryCond (),
						    pPatch->iRightBdryCond (),
						    pPatch->iLeftRegion (),
						    pPatch->iRightRegion (),
						    LdPts
						    );

	B2D.vAddPatchToList(pBChk);
	iNumSegs++;

				// Now must set the coefficients manually in case they were changed
				// because of a split or something
	double d1, d2, d3, d4, d5, d6, d7, d8;

	assert( dynamic_cast<BdryInterpolated2D*>(pPatch)->iNumSplines() ==
		dynamic_cast<BdryInterpolated2D*>(pBChk)->iNumSplines());

	for (int iCP = 0; iCP < dynamic_cast<BdryInterpolated2D*>(pBChk)->iNumSplines(); iCP++) {
	  dynamic_cast<BdryInterpolated2D*>(pPatch)->vGetSplineCoefficients(iCP, &d1, &d2, &d3, &d4, &d5, &d6, &d7, &d8);
	  dynamic_cast<BdryInterpolated2D*>(pBChk)->vSetSplineCoefficients(iCP, d1, d2, d3, d4, d5, d6, d7, d8);
	}

	dynamic_cast<BdryInterpolated2D*>(pBChk)->vRefresh();

	vMessage(4, " and added\n");

      }

      break;
    case BdryPatch2D::eUnknown:
      assert(0);
      break;
    } // end switch..
    //    }
  }

    // Now copy the periodic data stuff.
  for (int iP = 0; iP < iNumPer; iP++) {
    const Periodic *pP = B2DInfo.pPeriodicFace(iP);
    BdryPatch2D *pBP0 =
      B2D.pBdry(B2DInfo.iPatchIndex(dynamic_cast<BdryPatch2D*>(pP->pBPPatch(0))));
    BdryPatch2D *pBP1 =
      B2D.pBdry(B2DInfo.iPatchIndex(dynamic_cast<BdryPatch2D*>(pP->pBPPatch(1))));
    
    Periodic* pPNew = B2D.pPeriodicFace(iP);
    pPNew->vSetPatches(pBP0, pBP1);
  }
    
  vMessage(4, "Finished copying the Bdry2D info\n");
}

//////////////////////////////////////////////////////////////
//
//@ Triangulation of a generalized boundary geometry
//
//////////////////////////////////////////////////////////////
Mesh2D::Mesh2D (const Bdry2D & B2DInfo, const int iQualMeas,
		bool qInsertExtraVerts):
  Mesh (), VAEdgeF (), VAEdgeBF (), VATri (),
  VAQuad (), VAIntEdgeBF(), iOnEdge(0), iOnBdryEdge(0), iOnBdryEdge2(0),
  iInterior(0), iInternalBdry(0), iInternalBdry2(0), B2D()
{
  SUMAA_LOG_EVENT_BEGIN (INIT_TRI);
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(2, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);

  // List of edges to recover
  List <structTrackEdges*> LsTrackEdges;
  structTrackEdges *pEdge;
  int EdgeiV0 = -1, EdgeiV1 = -1, EdgeiPatch;

  // Concurrent list to keep track of what patch a vert is associated
  // with.
  // LsVertsPatch[0] refers to vert 0, etc...
  List <int> LsVertsPatch;


  dLipschitzAlpha = 1.0;
  dMeshResolution = 1.0;

  // First make a copy of the Bdry Info...
  vGetBdry2DInfo (B2DInfo);

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

  // Find max/min of segments
  int iPatch, iNBC = B2DInfo.iNumPatches();
  for (iPatch = 0; iPatch < iNBC; iPatch++) {
    double dSegXMin = B2D[iPatch]->dXMin();
    double dSegXMax = B2D[iPatch]->dXMax();
    double dSegYMin = B2D[iPatch]->dYMin();
    double dSegYMax = B2D[iPatch]->dYMax();

    dXMin = min (dXMin, dSegXMin);
    dXMax = max (dXMax, dSegXMax);
    dYMin = min (dYMin, dSegYMin);
    dYMax = max (dYMax, dSegYMax);
  }

  double dXLo = dXMin - (dXMax - dXMin);
  double dXHi = dXMax + (dXMax - dXMin);
  double dYLo = dYMin - (dYMax - dYMin);
  double dYHi = dYMax + (dYMax - dYMin);

  vMessage(3, "dXLo: %f dXHi: %f dYLo: %f dYHi: %f\n",
	   dXLo, dXHi, dYLo, dYHi);

  int iNumOrigVerts = iNumVerts ();

  Vert *pVLoLeft = pVNewVert ();
  Vert *pVLoRight = pVNewVert ();
  Vert *pVHiLeft = pVNewVert ();
  Vert *pVHiRight = pVNewVert ();
  // add dummy entries..
  for (int iDum = 0; iDum < 4; iDum++) {
    LsVertsPatch.vAppendItem(-100);
  }


  double adLoc[2];
  adLoc[0] = dXLo;
  adLoc[1] = dYLo;
  pVLoLeft->vSetCoords (2, adLoc);
  adLoc[1] = dYHi;
  pVHiLeft->vSetCoords (2, adLoc);
  adLoc[0] = dXHi;
  pVHiRight->vSetCoords (2, adLoc);
  adLoc[1] = dYLo;
  pVLoRight->vSetCoords (2, adLoc);

  // Set up initial mesh
  //
  //    HL-------0--------HR
  //     |              _/|
  //     |     0     __/  |
  //     |         _/     |
  //     1      _4/       3
  //     |    _/     1    |
  //     | __/            |
  //     |/               |
  //    LL-------2--------LR
  //

  // Setup sizes for initial connectivity arrays
  vSetupFaces (5);
  vSetupBFaces (4);
  vSetupTriCells (2);

  pFFace (0)->vAssign (pCCell (0), pBFBFace (0), pVHiRight, pVHiLeft);
  pBFBFace (0)->vAssign (pFFace (0));

  pFFace (1)->vAssign (pCCell (0), pBFBFace (1), pVHiLeft, pVLoLeft);
  pBFBFace (1)->vAssign (pFFace (1));

  pFFace (2)->vAssign (pCCell (1), pBFBFace (2), pVLoLeft, pVLoRight);
  pBFBFace (2)->vAssign (pFFace (2));

  pFFace (3)->vAssign (pCCell (1), pBFBFace (3), pVLoRight, pVHiRight);
  pBFBFace (3)->vAssign (pFFace (3));

  // Mark these as bdry faces
  for (int iF = 0; iF <= 3; iF++)
    pFFace(iF)->vSetFaceLoc(Face::eBdryFace);

  pFFace (4)->vAssign (pCCell (0), pCCell (1), pVLoLeft, pVHiRight);

  pCCell (0)->vAssign (pFFace (0), pFFace (1), pFFace (4));
  pCCell (1)->vAssign (pFFace (2), pFFace (3), pFFace (4));

  vSetAllHintFaces ();
  assert (qValid ());

  // Do Delaunay insertion of all the existing verts in the mesh
  vSetSwapType (eDelaunay);

  //	vMessage(2, "Number of verts now: %d\n", iNumVerts());
  //	for (int kk = 0; kk < iNumVerts(); kk++) {
  //		vMessage(2, "  %d -- (%f, %f)\n", kk, pVVert(kk)->dX(), pVVert(kk)->dY());
  //	}

  vMessage(2, "Going through patches to determine original discretization... ");

  int iOriginalV = -1;
  BdryPatch2D *pPatch;
  List <int> LiInsertedVerts;
  List <int> LiMeshToBdryVerts;
  LiInsertedVerts.vClear();
  LiInsertedVerts.vDisallowSorting();
  LiMeshToBdryVerts.vClear();
  LiMeshToBdryVerts.vDisallowSorting();

  for (iPatch = 0; iPatch < B2DInfo.iNumPatches(); iPatch++) {
    // Get the patch
    pPatch = B2D[iPatch];
    // Call the original discretization info
    List <double> La2dOtherV;
    bool qInsertPointers;
    bool qLoop;
    double adCoords[2];
    int iNewSwaps;

    La2dOtherV.vClear();

    La2dOtherV = pPatch->LdOriginalDiscretization(&qInsertPointers, &qLoop);

    // Insert the necessary verts, and build edge tracking list
    // at same time...

    bool qOriginalV = true;
    int iVert;
    int iMeshVert;
    Vert *pV;


    if (qInsertPointers) {
      iVert = pPatch->iBeginPointIndex();
      iMeshVert = iNumVerts();
      //      pV = pVVert(iVert);


      // We know iVert is associated to iPatch
      LsVertsPatch.vAppendItem(iPatch);
      pPatch->vGetBeginPoint(adCoords);

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

      if (qOriginalV) {
	iOriginalV = iMeshVert;
	qOriginalV = false;
      }

      // Insert pV0 for this patch
      // pV->vSetType(Vert::eBdryApex);

      //			vMessage(2, "iVert is: %d\n", iVert);
      if (!LiInsertedVerts.qInList(iVert)) {
	vMessage(3, "Inserting vert (%f, %f)\n", adCoords[0], adCoords[1]);
				// Must create the vertex
	pV = pVNewVert();
	if (B2D.qVertSmallAngle(iVert)) {
	  pV->vSetSmallAngleCorner(true);
	}
	pV->vSetCoords(2, adCoords);

	SUMAA_LOG_EVENT_BEGIN (INSERT_VERTEX);
#ifndef NDEBUG
	bool qSuccess = qInsertPoint (adCoords, pCCell (0), &iNewSwaps,
				      true, true, pV);
	assert (qSuccess);
#else
	(void) qInsertPoint (adCoords, pCCell (0), &iNewSwaps, true, true, pV);
#endif
	SUMAA_LOG_EVENT_END (INSERT_VERTEX);
	pV->vSetType(Vert::eBdryApex);


	LiInsertedVerts.vAppendItem(iVert);
	LiMeshToBdryVerts.vAppendItem(iMeshVert);

	//	vMessage(2, "Number of verts now: %d\n", iNumVerts());
	//	for (int kk = 0; kk < iNumVerts(); kk++) {
	//	  vMessage(2, "  %d -- (%f, %f)\n", kk, pVVert(kk)->dX(), pVVert(kk)->dY());
	//	}
	EdgeiV0 = iMeshVert;

      }
      else {
	// iV0 was already inserted. Find out which one it is in the mesh
	int iTarget = LiInsertedVerts.iSeek(iVert);
	EdgeiV0 = LiMeshToBdryVerts[iTarget];
      }

    } // if pointers are supposed to be inserted

    // now insert the middle "extra" verts specificied in La2dOtherV

    if (La2dOtherV.iLength() > 0) {
      for (int i = 0; i < La2dOtherV.iLength()/2; i++) {

				// insert point
	adCoords[0] = La2dOtherV[2*i];
	adCoords[1] = La2dOtherV[(2*i) + 1];

	iMeshVert = iNumVerts();

				// Create new vert
	pV = pVNewVert();
				// we know the patch for this vert...
	LsVertsPatch.vAppendItem(iPatch);

				// Set vert info
	pV->vSetCoords(2, adCoords);

	if (!qOriginalV) {
	  EdgeiV1 = iNumVerts()-1;
	  EdgeiPatch = iPatch;

	  pEdge = new structTrackEdges;
	  pEdge->iV0 = EdgeiV0;
	  pEdge->iV1 = EdgeiV1;
	  pEdge->iPatch = EdgeiPatch;
	  pEdge->pF = pFInvalidFace;
	  pEdge->qTriedRecover = false;
	  LsTrackEdges.vAppendItem(pEdge);

				// prepare next one..
	  EdgeiV0 = EdgeiV1;

	} else {
	  qOriginalV = false;
	  EdgeiV0 = iNumVerts()-1;
	  iOriginalV = EdgeiV0;
	}

	vMessage(3, "Inserting vert (%f, %f)\n", adCoords[0], adCoords[1]);

#ifndef NDEBUG
	bool qSuccess = qInsertPoint (adCoords, pCCell (0), &iNewSwaps,
				      true, true, pV);
	assert (qSuccess);
#else
	(void) qInsertPoint (adCoords, pCCell (0), &iNewSwaps, true, true, pV);
#endif
	pV->vSetType(Vert::eBdryApex);
				// This vertex is not in the Bdry2D list of points..
	int iDum = -1;
	LiInsertedVerts.vAppendItem(iDum);
	LiMeshToBdryVerts.vAppendItem(iMeshVert);

      } // loop over "middle" verts
    } // if there were any middle verts

    if (qInsertPointers) {
      // insert the "end" vert now
      iVert = pPatch->iEndPointIndex();
      iMeshVert = iNumVerts();
      //      pV = pVVert(iVert);

      LsVertsPatch.vAppendItem(iPatch);

      pPatch->vGetEndPoint(adCoords);
      //      adCoords[0] = pV->dX();
      //      adCoords[1] = pV->dY();



      //      vMessage(2, "Adding Patch: %d, iV0: %d, iV1: %d to list\n",
      //				pEdge->iPatch, pEdge->iV0, pEdge->iV1);
      //			vMessage(2, "(%f, %f) -- (%f, %f)\n",
      //				pVVert(pEdge->iV0)->dX(), pVVert(pEdge->iV0)->dY(),
      //				pVVert(pEdge->iV1)->dX(), pVVert(pEdge->iV1)->dY());

      // Insert pV1 for this patch
      //      pV->vSetType(Vert::eBdryApex);

      if (!LiInsertedVerts.qInList(iVert)) {
	vMessage(3, "Inserting vert (%f, %f)\n", adCoords[0], adCoords[1]);

				// Must create the vertex
	pV = pVNewVert();
	if (B2D.qVertSmallAngle(iVert)) {
	  pV->vSetSmallAngleCorner(true);
	}
	pV->vSetCoords(2, adCoords);


#ifndef NDEBUG
	bool qSuccess = qInsertPoint (adCoords, pCCell (0), &iNewSwaps,
				      true, true, pV);
	assert (qSuccess);
#else
	(void) qInsertPoint (adCoords, pCCell (0), &iNewSwaps, true, true, pV);
#endif
	pV->vSetType(Vert::eBdryApex);

	LiInsertedVerts.vAppendItem(iVert);
	LiMeshToBdryVerts.vAppendItem(iMeshVert);

				//
	EdgeiV1 = iMeshVert;

      }
      else {
				// This vertex was already inserted. Use the
				// mesh index stored in LiMeshToBdryVerts
	int iTarget = LiInsertedVerts.iSeek(iVert);
	EdgeiV1 = LiMeshToBdryVerts[iTarget];

      }

      EdgeiPatch = iPatch;

      pEdge = new structTrackEdges;
      pEdge->iV0 = EdgeiV0;
      pEdge->iV1 = EdgeiV1;
      pEdge->iPatch = EdgeiPatch;
      pEdge->pF = pFInvalidFace;
      pEdge->qTriedRecover = false;
      LsTrackEdges.vAppendItem(pEdge);


    }		// if end pointer was to be inserted...

    // Now check if we must loop
    if (qLoop) {
      EdgeiV0 = EdgeiV1;
      EdgeiV1 = iOriginalV;
      EdgeiPatch = iPatch;

      pEdge = new structTrackEdges;
      pEdge->iV0 = EdgeiV0;
      pEdge->iV1 = EdgeiV1;
      pEdge->iPatch = EdgeiPatch;

      pEdge->pF = pFInvalidFace;
      pEdge->qTriedRecover = false;
      LsTrackEdges.vAppendItem(pEdge);

    }

  } // loop over all patches

  iSwap(100);

  vMessage(2, "done.\nNow trying to recover all of the edges... ");

  vMessage(3, "List of edges to be recovered contains %d items\n", LsTrackEdges.iLength());

  for (int j = 0; j < LsTrackEdges.iLength(); j++) {
    pEdge = LsTrackEdges[j];
    vMessage(3, "  %d) iV0: %d, iV1: %d, iPatch: %d\n", j, pEdge->iV0, pEdge->iV1, pEdge->iPatch);
  }

  vSetAllHintFaces ();
  assert (qValid ());

  // Another step must be performed here.
  // We will make sure that all the boundary edges can be recovered
  // by using something similar to stitching.
  //
  // If an edge that is in the list of edges does NOT exist in the
  // triangulation then insert a vertex at its midpoint (making sure the
  // list gets updated properly).
  //
  // This fixes problems that arise when two boundary edges (wrongly) cross
  // each other.
  //
  // Now, everytime an edge that IS present in the mesh is found,
  // check for the following:
  //  - check for encroachment (2 opposite verts)
  //  - If it is encroached on, check distance to that bdry patch
  //  - If distance is positive, split.
  //
  // When this is over, the encroachment check step can be skipped
  // altogether..!

  bool qAllEdgesPresent = false;
  int iEdge = 0;
  int iNumEdgesPresent = 0;

  while (!qAllEdgesPresent) {
    // Look at iEdge
    // if it is OK, iNumEdgesPresent++
    // iEdge++ <-- Make sure it does not go too far

    // Figure out if the edge is present in mesh
    bool qEdgePresent = false;
    // Get the edge info
    pEdge = LsTrackEdges[iEdge];
    int iV0, iV1;
    Vert *pV0, *pV1;
    // Indices
    iV0 = pEdge->iV0;
    iV1 = pEdge->iV1;
    iPatch = pEdge->iPatch;
    // Get pointers
    pV0 = pVVert(iV0);
    pV1 = pVVert(iV1);
    vMessage(2, "Looking at edge %d -- (%f, %f)->(%f, %f) -- Patch #%d\n",
	     iEdge, pV0->dX(), pV0->dY(), pV1->dX(), pV1->dY(), iPatch);


    // Can save some time if we already have a face pointer
    if (pEdge->pF != pFInvalidFace) {
				// Edge was at some point present. Make sure it
				// still is.
      if (!pEdge->pF->qDeleted()) {
	Vert *pFV0 = pEdge->pF->pVVert(0);
	Vert *pFV1 = pEdge->pF->pVVert(1);

				// Make sure vertices are the same
	qEdgePresent =
	  ( ((pFV0 == pV0) && (pFV1 == pV1)) ||
	    ((pFV0 == pV1) && (pFV1 == pV0)) );
      }
      else {
				// Edge was deleted..
	qEdgePresent = false;
      }

    }
    else {
				// No pointer to the face -- Have to determine if
				// the edge is present in the mesh

				// Neighbor lists
      List<Vert*> LpVNear;
      List<Face*> LpFNear;
				// Get the neighbors
      vGetNeighborVertsAndFaces(pV0, LpVNear, LpFNear);
				// Go through vertex list -- is pV1 in there?
      int iV = 0;
      while ( (iV < LpVNear.iLength()) && (qEdgePresent == false) ) {
	if (LpVNear[iV] == pV1) qEdgePresent = true;
	iV++;
      }

      if (qEdgePresent) {
				// Set face pointer -- will save time later
	pEdge->pF = LpFNear[iV-1];
      }
    } // whether or not a pointer to the face existed already

    // Check for encroachment if edge exists
    bool qEncroach = false;
    if (qEdgePresent) {
      // Make sure it is not encroached on if it's a curved boundary edge
      // if (B2D[iPatch]->eGeom2DType != BdryPatch2D::ePolyLine) {
      // Check for encroachment. Find 2 opposite verts.
      Cell *pCL, *pCR;
      Vert *pVL, *pVR;

      double adC0[2], adC1[2];

      adC0[0] = pV0->dX();
      adC0[1] = pV0->dY();
      adC1[0] = pV1->dX();
      adC1[1] = pV1->dY();

      // Get the midpoint of the edge...
      double adMid[2];
      B2D[iPatch]->vRatioBetweenTwoPoints(0.5, adC0, adC1, adMid);


      // Is the left cell *really* to the left
      // of the patch..?

      pCL = pEdge->pF->pCCellLeft();
      pCR = pEdge->pF->pCCellRight();
      pVL = pCL->pVVertOpposite(pEdge->pF);
      pVR = pCR->pVVertOpposite(pEdge->pF);

      // Get orientations...
      int iOldL, iOldR, iNewL, iNewR;
      iOldL = iOrient2D(pV0, pV1, pVL);
      iOldR = iOrient2D(pV0, pV1, pVR);
      iNewL = iOrient2D(adC0, adMid, pVL->adCoords());
      iNewR = iOrient2D(adC0, adMid, pVR->adCoords());

      // Check difference in orientations
      if ( (iOldL != iNewL) || (iOldR != iNewR) ) {
	// Hmm.. Something wrong, I think....
	// vMessage(2, "NEW METHOD: Wrong orientations!\n");
	qEncroach = true;
      }
      /*
	iPL = LsVertsPatch[iVertIndex(pVL)];
	iPR = LsVertsPatch[iVertIndex(pVR)];
	
	B2D[iPatch]->vDistWhichSide(pVL->adCoords(), &dDL, &iSideL);
	B2D[iPatch]->vDistWhichSide(pVR->adCoords(), &dDR, &iSideR);
	
	vMessage(2, "Left vert: (%f, %f) -- Dist: %f -- Side: %d\n",
	pVL->dX(), pVL->dY(), dDL, iSideL);
	vMessage(2, "Right vert: (%f, %f) -- Dist: %f -- Side: %d\n",
	pVR->dX(), pVR->dY(), dDR, iSideR);
	
	
	// Make sure other vert is not part of the same patch..
	if ( (iFuzzyComp(dDR, 0.0) > 0) && (iFuzzyComp(dDL, 0.0) > 0) ) {
	// If the side is the same... there is something wrong!
	if (iSideL == iSideR) {
	qEncroach = true;
	vMessage(2, "OLD METHOD: Wrong orientations!\n");
	}
	}
      */
    }
    
    // if these two are OK, then increment, otherwise, split!
    if ( (qEdgePresent) &&  (!qEncroach) ) {
      vMessage(2, "Edge %d was present.\n", iEdge);
      iNumEdgesPresent++;
    }
    else {
				// Reset counter
      iNumEdgesPresent = 0;
				// Make sure the face pointer is invalid
      pEdge->pF = pFInvalidFace;
				// Get the patch
      pPatch = B2D[iPatch];

      if ( ( pPatch->eGeom2DType == BdryPatch2D::ePolyLine ) ||
	   (    (!pEdge->qTriedRecover) && (!qEncroach)    )
	   ) {
				// Try to recover the old way...
	vMessage(4, "Trying to recover that edge...\n");
	qRecoverEdge(pEdge->iV0, pEdge->iV1);
	pEdge->qTriedRecover = true;
      } else {
	// Already tried to recover the "old" way, so do something else...
	// If it's not a line,then insert a point...
	double adLoc1[2], adLoc2[2], adNewPoint[2];
	adLoc1[0] = pV0->dX();
	adLoc1[1] = pV0->dY();
	adLoc2[0] = pV1->dX();
	adLoc2[1] = pV1->dY();
	// Get the point
	pPatch->vRatioBetweenTwoPoints(0.5, adLoc1, adLoc2, adNewPoint);
	// insert it!
	int iVNew = iNumVerts();
	Vert *pV = pVNewVert();
	// Associate patch with it
	LsVertsPatch.vAppendItem(iPatch);

	// Set type;
	pV->vSetCoords(2, adNewPoint);
	pV->vSetType(Vert::eBdryApex);
				//
	int iNewSwaps = 0;

	SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
#ifndef NDEBUG
	bool qSuccess = qInsertPoint (adNewPoint, pCCell (0), &iNewSwaps,
				      true, true, pV);
	assert (qSuccess);
#else
	(void) qInsertPoint (adNewPoint, pCCell (0), &iNewSwaps, true, true, pV);
#endif
	SUMAA_LOG_EVENT_END(INSERT_VERTEX);
	// Delete the old edge from the tracking list
	LsTrackEdges.vDelete(pEdge);

	// add the two new ones
	EdgeiV0 = iV0;
	EdgeiV1 = iVNew;
	EdgeiPatch = iPatch;

	pEdge = new structTrackEdges;
	pEdge->iV0 = EdgeiV0;
	pEdge->iV1 = EdgeiV1;
	pEdge->iPatch = EdgeiPatch;
	pEdge->pF = pFInvalidFace;
	pEdge->qTriedRecover = false;
	LsTrackEdges.vAppendItem(pEdge);

	EdgeiV0 = iVNew;
	EdgeiV1 = iV1;

	pEdge = new structTrackEdges;
	pEdge->iV0 = EdgeiV0;
	pEdge->iV1 = EdgeiV1;
	pEdge->iPatch = EdgeiPatch;
	pEdge->pF = pFInvalidFace;
	pEdge->qTriedRecover = false;
	LsTrackEdges.vAppendItem(pEdge);
      }
    }

    // Do the iEdge increment properly
    if (iEdge < LsTrackEdges.iLength()-1) {
				// Can increment
      iEdge++;
    }
    else {
				// Loop back
      vMessage(2, "Looping back; %d edges present\n", iNumEdgesPresent);
      iEdge = 0;
    }

    // Are we there yet?
    if (iNumEdgesPresent == LsTrackEdges.iLength()) qAllEdgesPresent = true;
  }

//   vMessage(4, "Now marking all verts at the ends of curved patches as apexes\n");
//   for (iPatch = 0; iPatch < iNumPatch; iPatch++) {
//     // Get the patch
//     pPatch = B2D[iPatch];

//     // Do nothing for polylines and patches of unknown type (the latter
//     // are impossible anyway).
//     if (pPatch->eGeom2DType == BdryPatch2D::eUnknown ||
//	pPatch->eGeom2DType == BdryPatch2D::ePolyLine) continue;

//     // For other patches, mark their endpoints as apexes to prevent
//     // inappropriate surface smoothing.  The "if" is there because
//     // circles have -1 as their begin and end indices.
//     int iBeginVert = pPatch->iBeginPointIndex();
//     if (iBeginVert >= 0) pVVert(iBeginVert)->vSetType(Vert::eBdryApex);
//     int iEndVert   = pPatch->iEndPointIndex();
//     if (iEndVert >= 0) pVVert(iEndVert)->vSetType(Vert::eBdryApex);
//   }

  // Mark all cells as having an invalid region, since we don't yet know
  // where any of them are.
  int iCell;
  for (iCell = 0; iCell < iNumCells (); iCell++)
    pCCell (iCell)->vSetRegion (iInvalidRegion);

  struct FaceInfo *aFI = new struct FaceInfo[LsTrackEdges.iLength()];

  // For each boundary patch:
  //   1. Mark cells that are in the interior of the domain to be meshed
  //   with their region.
  //   2. Mark cells that are in the exterior for deletion.

  int iNumOneBdryPatches = 0;
  vMessage(4, "Assign proper boundary patches data to neighbor cells\n");
  for (iEdge = 0; iEdge < LsTrackEdges.iLength(); iEdge++) {

    vMessage(4, "Edge %d\n", iEdge);
    pEdge = LsTrackEdges[iEdge];

    int iV0 = pEdge->iV0;
    iPatch = pEdge->iPatch;
    Vert *pV0 = pVVert (iV0);
    BFace *pBF;

    // Edge was already found before
    Face *pF = pEdge->pF;
    vMessage(4, "Face pointer is %p\n", pF);

    assert (pF->qValid ());
#ifndef NDEBUG
    int iV1 = pEdge->iV1;
    Vert *pV1 = pVVert (iV1);
    assert (pF->pVVert (0) == pV0 || pF->pVVert (0) == pV1);
    assert (pF->pVVert (1) == pV0 || pF->pVVert (1) == pV1);
#endif

    vMessage(4, "   pVVert(0): (%f, %f) -- pVVert(1) (%f, %f)\n",
	     pF->pVVert(0)->dX(), pF->pVVert(0)->dY(),
	     pF->pVVert(1)->dX(), pF->pVVert(1)->dY());

    // Determine whether the edge and the bdry patch are pointing the
    // same direction.
    bool qSameOrient = (pF->pVVert (0) == pV0);

    vMessage(4, "Marking regions..\n");

    BdryPatch2D *pBPatch = B2D[iPatch];
    // 1. Mark regions.
    if (pBPatch->iLeftRegion () != iInvalidRegion) {
      if (qSameOrient) {
	pF->pCCellLeft ()->vSetRegion (pBPatch->iLeftRegion ());
	vMessage(4, "Marking left cell with left region number\n");
      }
      else {
	pF->pCCellRight ()->vSetRegion (pBPatch->iLeftRegion ());
	vMessage(4, "Marking right cell with left region number\n");
      }
    }
    if (pBPatch->iRightRegion () != iInvalidRegion) {
      if (qSameOrient) {
	pF->pCCellRight ()->vSetRegion (pBPatch->iRightRegion ());
	vMessage(4, "Marking right cell with right region number\n");
      }
      else {
	pF->pCCellLeft ()->vSetRegion (pBPatch->iRightRegion ());
	vMessage(4, "Marking left cell with right region number\n");
      }
    }

    // 2. Mark boundaries.
    // Is one side of the edge a boundary?
    if (pBPatch->iBCWhichSide () != BC_SIDE_NONE) {
      // Make sure this is only done when it's an external boundary

      vMessage(4, "   This patch is an external boundary\n");
      Cell *pC;
      if (pBPatch->iBCWhichSide () == BC_SIDE_LEFT) {
				// BC is left side of the patch
	if (qSameOrient)
	  pC = pF->pCCellLeft ();
	else
	  pC = pF->pCCellRight ();
      }
      else {
				// BC is right side of the patch...
	if (qSameOrient)
	  pC = pF->pCCellRight ();
	else
	  pC = pF->pCCellLeft ();
      }

      pBF = pBFNewBFace ();

      // Make sure the face is set properly.
      pBF->vAssign (pF);
      pF->vSetFaceLoc (Face::eBdryFace);

      // This is used later on to delete the unused exterior cells
      aFI[iNumOneBdryPatches].pC = pC;
      aFI[iNumOneBdryPatches].pF = pF;
      aFI[iNumOneBdryPatches].pBF = pBF;
      iNumOneBdryPatches++;
    }
    else {
      // This is an internal boundary

      vMessage(4, "   This patch is an internal boundary\n");

      // Must create...
      //  1. An internal boundary face..
      //  2. Another face. The IBF is located "between" the two faces..

      vMessage(4, "   Creating the extra face\n");

      Face *pF2 = pFNewFace ();

      vMessage(4, "   Extra face is: %d (%p)\n", iFaceIndex(pF2), pF2);

      vMessage(4, "   Creating the new IntBFace..\n");

      pBF = pIBENewIntBdryEdge ();


      // The new face will be the "right" face of the IBF.
      // This means that the IBF is located on the left side
      // of the new face..

      vMessage(4, "   Assigning the cells and verts to the faces\n");


      Cell *pCOther;
      if (qSameOrient) {
	pCOther = pF->pCCellRight ();
	pF2->vAssign (pBF, pCOther, pF->pVVert (0), pF->pVVert (1));
      }
      else {
	pCOther = pF->pCCellLeft ();
	pF2->vAssign (pCOther, pBF, pF->pVVert (0), pF->pVVert (1));
      }

      pF->vReplaceCell (pCOther, pBF);
      pCOther->vReplaceFace (pF, pF2);

      vMessage(4, "   Setting faces types\n");

      // Make sure they know they are internal boundaries
      pF2->vSetFaceLoc (Face::eBdryTwoSide);
      pF->vSetFaceLoc (Face::eBdryTwoSide);


      vMessage(4, "   Assigning faces to the IBF\n");

      // Make sure IBF knows about the two faces...
      //
      // Is the order really important here..?
      //
      dynamic_cast<IntBdryEdge*>( pBF)->vAssign (pF, pF2);
    }

    // Set the patch index
    dynamic_cast<BdryEdgeBase *>(pBF)->vSetPatch (iPatch, &B2D);

  }				// for loop over all patches

  vMessage(2, "done.\nNow deleting unnecessary cells and faces... ");

  // For every edge separating interior and exterior, replace the
  // exterior cell with a boundary face (which has a pointer to the
  // correct boundary patch so that geometric info is not lost).

  vMessage(4, "Now deleting the exterior cells of Bdry Faces...\n");
  for (iPatch = 0; iPatch < iNumOneBdryPatches; iPatch++) {
    // Only do this stuff if this is -really- a boundary face.
    vMessage(4, "   Checking exterior patch %d\n", iPatch);
    if (aFI[iPatch].pF->qValid ()) {
      Cell *pC = aFI[iPatch].pC;
      Face *pF = aFI[iPatch].pF;
      BFace *pBF = aFI[iPatch].pBF;

      // Disconnect the old cell.
      pC->vRemoveFace (pF);

      // Replace the old cell with the new bdry face in the face
      // connectivity.
      pF->vReplaceCell (pC, pBF);

      pC->vMarkDeleted();
      // And that's it.  The rest has already been taken care of.
    }
  }
  delete[]aFI;


  vMessage(4, "Marking regions\n");

  // March out from marked cells, filling all unmarked neighbors with
  // correct region tag (interior) or marking for deletion (exterior) in
  // a second pass.
  for (iCell = 0; iCell < iNumCells (); iCell++) {
    Cell *pC = pCCell (iCell);
    // Can't march out from cells that don't know where they are, and
    // don't need to march out from cells in the exterior.
    if (pC->iRegion () == iInvalidRegion)
      continue;
    vMarkCleanNeighbors (pC);
  }

  vMessage (4, "Marking exterior cells for deletion\n");

  // Mark exterior cells for deletion.
  for (iCell = 0; iCell < iNumCells (); iCell++) {
    Cell *pC = pCCell (iCell);
    if (pC->iRegion () == iInvalidRegion)
      pC->vMarkDeleted ();
  }

  vMessage (4, "Marking exterior verts for deletion\n");

  // Mark four farfield vertices for deletion.
  int iVert;
  for (iVert = iNumOrigVerts; iVert < iNumOrigVerts+4; iVert++)
    pVVert (iVert)->vMarkDeleted ();

  // Doublecheck to be sure that no kept cells have a deleted vertex.
  // If so, this is an input error.  Further diagnosis requires throwing
  // an exception.
  // FIX ME: 0.2.2.
  for (iCell = 0; iCell < iNumCells(); iCell++) {
    Cell *pC = pCCell(iCell);
    if (!pC->qDeleted()) {
      for (iVert = 0; iVert < pC->iNumVerts(); iVert++) {
	if (pC->pVVert(iVert)->qDeleted()) {
	  vFatalError("Problem with polygon orientation in input file.",
		      "Mesh2D::Mesh2D(Bdry2D)");
	}
      }
    }
  }

  // Mark farfield boundary faces for deletion.  These are always the
  // first four BFaces.
  int iBFace;
  for (iBFace = 0; iBFace < 4; iBFace++)
    pBFBFace (iBFace)->vMarkDeleted ();

  // Mark all faces bounded by two deleted cells for deletion. This must
  // come -after- bdry face deletion so that the farfield faces are
  // deleted properly.
  int iFace;
  for (iFace = 0; iFace < iNumFaces (); iFace++) {
    Face *pF = pFFace (iFace);
    if (pF->pCCellLeft ()->qDeleted () &&
	pF->pCCellRight ()->qDeleted ())
      pF->vMarkDeleted ();
  }

  // Verify that markings are consistent in that all interior regions
  // are bounded away from each other by region boundaries.
  for (iFace = 0; iFace < iNumFaces (); iFace++) {
    Face *pF = pFFace (iFace);
    if (pF->qDeleted ())
      continue;			// Don't bother about deleted edges.

    if (pF->iFaceLoc () == Face::eBdryFace)
      continue;			// Ditto for boundaries.

    int iRegL = pF->pCCellLeft ()->iRegion ();
    int iRegR = pF->pCCellRight ()->iRegion ();
    if ((iRegL != iRegR) && (pF->iFaceLoc () != Face::eBdryTwoSide))
      vFatalError ("Orientation of an internal boundary is incorrect.",
		   "2D initial meshing from boundary geometry.");
  }

  // Verify that cells aren't screwed up; this can happen with bad
  // boundary orientation.
  for (iCell = 0; iCell < iNumCells (); iCell++) {
    Cell *pC = pCCell (iCell);
    if (pC->qDeleted ())
      continue;
    if (!(pC->pFFace (0)->qValid () &&
	  pC->pFFace (1)->qValid () &&
	  pC->pFFace (2)->qValid ()))
      vFatalError ("Orientation of an internal boundary is incorrect.",
		   "2D initial meshing from boundary geometry.");
  }


  vMessage(4, "About to purge mesh\n");
  // Do the actual purge to eliminate that stuff from the database.
  vPurge ();
  if (3 * iNumCells () + iNumBdryFaces () + 2 * iNumIntBdryFaces () !=
      2 * iNumFaces ()) {
    vMessage(1, "NumCells: %d, NumIntBdryFaces: %d, NumBdryFaces: %d, NumFaces: %d\n",
	     iNumCells (), iNumIntBdryFaces (), iNumBdryFaces (), iNumFaces ());
    vFatalError ("Incompatible number of cells, faces, and boundary faces.\nProbably a boundary is oriented incorrectly.",
		 "2D initial meshing from boundary geometry.");
  }

  // Mark vertices on internal boundaries.
  vMessage(4, "Now marking vertices on internal boundaries\n");
  for (iFace = 0; iFace < iNumFaces (); iFace++) {
    Face *pF = pFFace (iFace);
    if (pF->iFaceLoc () == Face::eBdryTwoSide) {
      vMessage(4, "  ... for Face %d\n", iFace);
      pF->pVVert (0)->vSetType (Vert::eBdryTwoSide);
      pF->pVVert (1)->vSetType (Vert::eBdryTwoSide);
    }
  }

  // Mark vertices on external boundaries.  Must be done after internal
  // boundaries.
  vMessage(4, "Now marking vertices on external boundaries...\n");
  for (iBFace = 0; iBFace < iNumBdryFaces (); iBFace++) {
    vMessage(4, "  ... for BFace %d\n", iBFace);
    Face *pF = pBFBFace (iBFace)->pFFace (0);
    //      pF->pVVert(0)->vSetType(Vert::eBdryApex);
    //      pF->pVVert(1)->vSetType(Vert::eBdryApex);
    switch (pF->pVVert (0)->iVertType ()) {
    case Vert::eBdryTwoSide:
      pF->pVVert (0)->vSetType (Vert::eBdryApex);
      break;
    case Vert::eBdryApex:
      break;
    default:
      pF->pVVert (0)->vSetType (Vert::eBdry);
      break;
    }
    switch (pF->pVVert (1)->iVertType ()) {
    case Vert::eBdryTwoSide:
      pF->pVVert (1)->vSetType (Vert::eBdryApex);
      break;
    case Vert::eBdryApex:
      break;
    default:
      pF->pVVert (1)->vSetType (Vert::eBdry);
      break;
    }
  }

  vMessage (4, "Getting rid of list of edges to be recovered\n");
  for (int kk = LsTrackEdges.iLength()-1; kk >= 0; kk--) {
    pEdge = LsTrackEdges[kk];
    delete pEdge;
  }
  LsTrackEdges.vClear();

  vMessage(2, "done.\n ");

  if (qInsertExtraVerts) {
    // Insert remaining points from Bdry2D object.
    int iTotalSwaps = 0, iAddedPoints = 0;
    for (int iP = 0; iP < B2DInfo.iNumPoints(); iP++) {
      if (B2DInfo.iNumPatchesConnectedToPoint(iP) == 0) {
	B2DInfo.vGetPoint(iP, adLoc);
	int iSwaps = 0;
	vMessage(2, "   Inserting a point at (%f, %f)\n",
		 adLoc[0], adLoc[1]);
	// Be sure to force insertion at the given point location...
	if (qInsertPoint(adLoc, pCCell(0), &iSwaps, true, true)) {
	  Vert *pVNew = pVVert(iNumVerts() - 1);
	  vMessage(2, "    Actually at (%f, %f)\n",
		   pVNew->dX(), pVNew->dY());
	  pVNew->vSetType(Vert::eInteriorFixed);
	  iAddedPoints++;
	  iTotalSwaps += iSwaps;
	}
      }
    }
    if (iAddedPoints > 0) {
      vMessage(2, "Added %d additional points from the input file; %d swaps.\n",
	       iAddedPoints, iTotalSwaps);
    }
  }

  // Don't know if the following is necessary...?
  vPurge ();
  vSetAllHintFaces ();

  assert (qValid ());
  SUMAA_LOG_EVENT_END (INIT_TRI);

  vMessage(2, "Swapping to make sure everything is Delaunay...\n");

  // try to swap
  int iNSwaps = iSwap();
  vMessage(2, "done; %d swaps.\n", iNSwaps);

  vResetPeriodicData();
}

void Mesh2D::vResetPeriodicData() 
{
  if (B2D.iNumPeriodic() != 0) {
    vMessage(2, "Resetting data for periodic faces...");
    B2D.vClearPeriodicFaceData();
    int iNP = 0;

    for (int iBF = 0; iBF < iNumBdryFaces(); iBF++) {
      BFace* pBF = pBFBFace(iBF);
      if (pBF->qDeleted()) continue;
      BdryPatch* pBP = pBF->pPatchPointer();
      if (pBP->qIsPeriodic()) {
	iNP++;
	Periodic* pP = pBP->pPGetPeriodicData();
	pP->vAddFace(pBP, pBF->pFFace());
      }
    }
    vMessage(2, " %d periodic faces.\n", iNP);
  }
}

//@ Copy constructor
Mesh2D::Mesh2D (const Mesh2D & M)
  : Mesh (M), VAEdgeF (M.iNumFaces()),
    VAEdgeBF(M.iNumBdryFaces()), VATri(M.iNumTriCells()),
    VAQuad(M.iNumQuadCells()), VAIntEdgeBF(M.iNumIntBdryFaces()),
    iOnEdge(0), iOnBdryEdge(0), iOnBdryEdge2(0), iInterior(0),
    iInternalBdry(0), iInternalBdry2(0), B2D()
{
  // Set up the pointers for the new edges.
  int iF;
  for (iF = iNumFaces () - 1; iF >= 0; iF--) {
    Face *pF_Old = M.pFFace (iF);
    int iV0 = M.iVertIndex (pF_Old->pVVert (0));
    int iV1 = M.iVertIndex (pF_Old->pVVert (1));

    Cell *pCL_Old = pF_Old->pCCellLeft ();
    Cell *pCL_New;
    if (!pCL_Old->qValid())
      pCL_New = pCInvalidCell;
    else if (pCL_Old->eType () == Cell::eBdryEdge)
      pCL_New = pBFBFace (M.iBFaceIndex (dynamic_cast<BdryEdgeBase *>(pCL_Old)));
    else
      pCL_New = pCCell (M.iCellIndex (pCL_Old));

    Cell *pCR_Old = pF_Old->pCCellRight ();
    Cell *pCR_New;
    if (!pCR_Old->qValid())
      pCR_New = pCInvalidCell;
    else if (pCR_Old->eType () == Cell::eBdryEdge)
      pCR_New = pBFBFace (M.iBFaceIndex (dynamic_cast<BdryEdgeBase *>(pCR_Old)));
    else
      pCR_New = pCCell (M.iCellIndex (pCR_Old));

    pFFace (iF)->vAssign (pCL_New, pCR_New,
			  pVVert (iV0), pVVert (iV1));
    pFFace (iF)->vCopyAllFlags (*pF_Old);
  }

  // Set up the pointers for the new cells, both quads and tris.
  for (int iC = iNumCells () - 1; iC >= 0; iC--) {
    Cell *pC_Old = M.pCCell (iC);
    if (!pC_Old->qValid()) continue;
    Face *apF[4];
    for (iF = 0; iF < pC_Old->iNumFaces (); iF++)
      apF[iF] = pFFace (M.iFaceIndex (pC_Old->pFFace (iF)));
    pCCell (iC)->vAssign (apF);
    pCCell (iC)->vCopyAllFlags (*pC_Old);
  }

  // Set up pointers for new bdry faces.
  // Should probably have the pointers to the patches too, huh?
  for (int iBF = iNumBdryFaces () - 1; iBF >= 0; iBF--) {
    BFace *pBF_Old = M.pBFBFace (iBF);
    BFace *pBF_New = pBFBFace (iBF);
    pBF_New->vAssign (pFFace (M.iFaceIndex (pBF_Old->pFFace (0))));

    dynamic_cast<BdryEdgeBase*>( pBF_New)->vSetPatch (dynamic_cast<BdryEdgeBase*>( pBF_Old)->iPatchIndex (),
				       &B2D);
    pBF_New->vCopyAllFlags (*pBF_Old);
  }

  // Charles added this.. Is it correct?
  // B2D = M.B2D;
  // End of addition

  // Carl's guess at what is correct.  Need a real operator= for Bdry2D
  // objects, though, and perhaps a copy constructor
  vGetBdry2DInfo (M.B2D);

  // Mesh doesn't have its own copy constructor, so this stuff has to be
  // copied here.
  qSwapRecur = M.qSwapRecur;
  qAllowBdryChanges = M.qAllowBdryChanges;
  eSwapMeasure = M.eSwapMeasure;
  dLengthFactor = M.dLengthFactor;
  dfLengthScale = M.dfLengthScale;
  dLipschitzAlpha = M.dLipschitzAlpha;
  dMeshResolution = M.dMeshResolution;

  // Mesh2D scalar data
  qSimplex = M.qSimplex;
  iOnEdge = M.iOnEdge;
  iOnBdryEdge = M.iOnBdryEdge;
  iOnBdryEdge2 = M.iOnBdryEdge2;
  iInterior = M.iInterior;
  iInternalBdry = M.iInternalBdry;
  iInternalBdry2 = M.iInternalBdry2;
  qLengthScaleFromCellSizes = M.qLengthScaleFromCellSizes;

  // Set hint faces.
//   vSetAllHintFaces ();
  for (int iV = 0; iV < iNumVerts(); iV++) {
    Vert *pVNew = pVVert(iV);
    Vert *pVOld = M.pVVert(iV);
    Face *pFOldHint = pVOld->pFHintFace();
    if (pFOldHint->qValid() && !pFOldHint->qDeleted() &&
	pFOldHint->qHasVert(pVOld)) {
      // It's a go.
      int iFace = M.iFaceIndex(pFOldHint);
      pVNew->vSetHintFace(pFFace(iFace));
    }
    else {
      pVNew->vClearHintFace();
    }
  }
  assert (qValid ());
}

Mesh2D& Mesh2D::
operator= (const Mesh2D & M)
{
  assert (M.qValid ());
  if (this != &M) {
    qSimplex = M.qSimplex;
    qBdryFill = M.qBdryFill;
    iIntSkip = M.iIntSkip;
    eEncType = M.eEncType;

    // Copy the essential smoothing data over existing data.
    {
      SMsmooth_data *pSMSmooth = static_cast<SMsmooth_data*>(M.pvSmoothData);
      int iFunc = pSMSmooth->smooth_param->function_id;
      int iTech = pSMSmooth->smooth_param->smooth_technique;
      double dAcceptValue = pSMSmooth->smooth_param->lap_accept_value;
      SMconvertToDegrees(static_cast<SMsmooth_data *>(pSMSmooth)->smooth_param
			 ->function_id, &dAcceptValue);

      SMinitSmoothing(2, iTech, iFunc, dAcceptValue, &pvSmoothData);
      SMinitGlobalMinValue(pvSmoothData);
    }

    // This one copies easily because there are no pointers (hint
    // face is cleared.
    VAVerts = M.VAVerts;

    // The remainder are trickier.  First, set up sizes.
    VAEdgeF.vSetup (M.VAEdgeF.iSize ());
    VAEdgeBF.vSetup (M.VAEdgeBF.iSize ());
    VATri.vSetup (M.VATri.iSize ());
    VAQuad.vSetup (M.VAQuad.iSize ());

    // Set up the pointers for the new edges.
    int iF;
    for (iF = iNumFaces () - 1; iF >= 0; iF--) {
      Face *pF_Old = M.pFFace (iF);
      int iV0 = M.iVertIndex (pF_Old->pVVert (0));
      int iV1 = M.iVertIndex (pF_Old->pVVert (1));

      Cell *pCL_Old = pF_Old->pCCellLeft ();
      Cell *pCL_New;
      if (pCL_Old->eType () == Cell::eBdryEdge)
      pCL_New = pBFBFace (M.iBFaceIndex (dynamic_cast<BdryEdgeBase*>( pCL_Old)));
      else
      pCL_New = pCCell (M.iCellIndex (pCL_Old));

      Cell *pCR_Old = pF_Old->pCCellRight ();
      Cell *pCR_New;
      if (pCR_Old->eType () == Cell::eBdryEdge)
      pCR_New = pBFBFace (M.iBFaceIndex (dynamic_cast<BdryEdgeBase*>( pCR_Old)));
      else
      pCR_New = pCCell (M.iCellIndex (pCR_Old));

      pFFace (iF)->vAssign (pCL_New, pCR_New,
			    pVVert (iV0), pVVert (iV1));
      pFFace (iF)->vCopyAllFlags (*pF_Old);
    }

    // Set up the pointers for the new cells, both quads and tris.
    for (int iC = iNumCells () - 1; iC >= 0; iC--) {
      Cell *pC_Old = M.pCCell (iC);
      Face *apF[4];
      for (iF = 0; iF < pC_Old->iNumFaces (); iF++)
      apF[iF] = pFFace (M.iFaceIndex (pC_Old->pFFace (iF)));
      pCCell (iC)->vAssign (apF);
      pCCell (iC)->vCopyAllFlags (*pC_Old);
    }

    // Set up pointers for new bdry faces.
    // Should probably have the pointers to the patches too, huh?
    for (int iBF = iNumBdryFaces () - 1; iBF >= 0; iBF--) {
//        int iPatchGoneCorrection = 0;
      BFace *pBF_Old = M.pBFBFace (iBF);
      BFace *pBF_New = pBFBFace (iBF);
      pBF_New->vAssign (pFFace (M.iFaceIndex (pBF_Old->pFFace (0))));

//        // Should set the pointer to pPatch instead...
//        for (int iChk = 0; iChk <= dynamic_cast<BdryEdgeBase*>( pBF_Old)->iPatchIndex ();
//	   iChk++) {
//	     // Check if the old patch is used...
//	     if (!M.B2D.qIsPatchUsed (iChk))
//	     iPatchGoneCorrection++;
//	   }
//        dynamic_cast<BdryEdgeBase*>( pBF_New)->vSetPatch (
//					 dynamic_cast<BdryEdgeBase*>( pBF_Old)->iPatchIndex () - iPatchGoneCorrection,
//					 &B2D);
      dynamic_cast<BdryEdgeBase*>(pBF_New)->
	vSetPatch(dynamic_cast<BdryEdgeBase*>(pBF_Old)->iPatchIndex(), &B2D);
      pBF_New->vCopyAllFlags (*pBF_Old);
    }

    // Added this part to copy the Bdry2D object..
    vGetBdry2DInfo (M.B2D);

    iOnEdge = M.iOnEdge;
    iOnBdryEdge = M.iOnBdryEdge;
    iOnBdryEdge2 = M.iOnBdryEdge2;
    iInterior = M.iInterior;
    iInternalBdry = M.iInternalBdry;
    iInternalBdry2 = M.iInternalBdry2;
    qLengthScaleFromCellSizes = M.qLengthScaleFromCellSizes;
    vSetAllHintFaces ();
  }
  assert (qValid ());
  return (*this);
}

static Face *
pFFindEdge (const Vert * const pV0, const Vert * const pV1,
	    const VAEdgeFace * const pVAEF)
{
  assert (pV0->qValid ());
  assert (pV1->qValid ());

  for (int i = 0; i < pVAEF->iSize (); i++) {
    Face *pF = pVAEF->pTEntry (i);
    if (pF->qHasVert (pV0) && pF->qHasVert (pV1) &&
	(pF->pCCellLeft () == pCInvalidCell ||
	 pF->pCCellRight () == pCInvalidCell))
      return pF;
  }
  return (pFInvalidFace);
}

//@ Merge together two connectivities which are topologically two-dimensional
void Mesh2D::
vMerge (const Mesh2D & MeshA)
  // The two connectivities being merged are assumed to be disjoint.
  // If they aren't, assertions will alert one to this fact.
  //
  // There are two choices for doing this:
  //
  // If one connectivity is much smaller than the other or both are
  // small, the merge will be done by explicit matching of faces.
  //
  // If the two are comparable in size and large, the whole thing
  // will be done from scratch, with edges paired up by sorting.
{
  // For now, just merge everything.
  //   // Verify that the respective vertices of the two meshes are identical.
  //   // This is critical to correct function of code later in the routine.
  // #ifndef NDEBUG
  //   assert(iNumVerts() == MeshA.iNumVerts());
  //   for (int iTest = 0; iTest < iNumVerts(); iTest++) {
  //     const double * const adC = pVVert(iTest)->adCoords();
  //     const double * const adCA = MeshA.pVVert(iTest)->adCoords();
  //     assert(dDIST2D(adC, adCA) < 1.e-8);
  //   }
  // #endif

  // All cells must be added to *this.  The new cells are initialized to
  // be clones of the old cells, and a note is made of where the new one
  // goes so that the connectivity can be adjusted later.
  Cell **apCNew = new Cell *[MeshA.iNumCells ()];
  int i, iFacesOnEntry = iNumFaces ();
  for (i = 0; i < MeshA.iNumCells (); apCNew[i++] = pCInvalidCell) {}
  for (i = 0; i < MeshA.iNumCells (); i++) {
    Cell *pC = MeshA.pCCell (i);
    apCNew[i] = pCNewCell ();
    // Only copy the face pointers over; don't muck with updating which
    // faces and cells point to each other, since that is done below.
    {
      Face *apF[4];
      for (int ii = pC->iNumFaces () - 1; ii >= 0; ii--)
	apF[ii] = pC->pFFace (ii);
      apCNew[i]->vAssign (apF);
    }
  }

  // For each face(edge) determine whether it coincides with an existing
  // one.  If not create the edge.  In any case, adjust the connectivity
  // in *this.
  for (i = 0; i < MeshA.iNumFaces (); i++) {
    Face *pF = MeshA.pFFace (i);
    Vert *pV0 = pVVert (MeshA.iVertIndex (pF->pVVert (0)));
    Vert *pV1 = pVVert (MeshA.iVertIndex (pF->pVVert (1)));
    Face *pFMatch = pFFindEdge (pV0, pV1, &VAEdgeF);
    if (pFMatch == pFInvalidFace) {
      // Create new face.
      Face *pFNew = pFNewFace ();

      // If both cells bordering the face are interior cells, reconnect
      // them both.  If not, make sure that boundary faces have a raw
      // cell left over by forcing the face to point to a null cell.
      Cell *pCLNew = pCInvalidCell, *pCRNew = pCInvalidCell;
      assert (pF->pCCellLeft ()->qValid () && pF->pCCellRight ()->qValid ());
      if (pF->pCCellLeft ()->eType () != Cell::eBdryEdge) {
	pCLNew = apCNew[MeshA.iCellIndex (pF->pCCellLeft ())];
	pCLNew->vReplaceFace (pF, pFNew);
      }
      if (pF->pCCellRight ()->eType () != Cell::eBdryEdge) {
	pCRNew = apCNew[MeshA.iCellIndex (pF->pCCellRight ())];
	pCRNew->vReplaceFace (pF, pFNew);
      }
      // Set up face connectivity.
      pFNew->vAssign (pCLNew, pCRNew, pV0, pV1);

      // Correct the connectivity of the (guaranteed to be) new cells on
      // either side of it.

    }				// Face in MeshA does not match face in *this

    else {
      // Ensure that the faces in *this and MeshA don't have neighboring
      // cells which conflict with each others' locations.
      Cell *pCLthis = pFMatch->pCCellLeft ();
      Cell *pCRthis = pFMatch->pCCellRight ();
      Cell *pCLOld = pF->pCCellLeft ();
      Cell *pCROld = pF->pCCellRight ();

      Cell *pCLNew = pCInvalidCell, *pCRNew = pCInvalidCell;

      if (pCLOld->eType () != Cell::eBdryEdge)
	pCLNew = apCNew[MeshA.iCellIndex (pCLOld)];

      if (pCROld->eType () != Cell::eBdryEdge)
	pCRNew = apCNew[MeshA.iCellIndex (pCROld)];

      if (pV0 == pFMatch->pVVert (0)) {
      }
      else {
	Cell *pCTemp = pCRNew;
	pCRNew = pCLNew;
	pCLNew = pCTemp;
      }
      // Exactly one of each of these pairs must be invalid, unless the
      // hull being added to consists of only a bare edge.
      if (iFacesOnEntry != 1) {
	assert (XOR (pCLthis->qValid (), pCRthis->qValid ()));
	assert (XOR (pCLNew->qValid (), pCRNew->qValid ()));
	assert (XOR (pCLthis->qValid (), pCLNew->qValid ()));
	assert (XOR (pCRthis->qValid (), pCRNew->qValid ()));
	// Update the connectivity of the existing face in *this.
	// Correct the connectivity of the (guaranteed to be) one new cell
	// for this face.
	if (pCLthis->qValid ()) {
	  pFMatch->vAddCell (pCRNew);
	  pCRNew->vReplaceFace (pF, pFMatch);
	}
	else {
	  pFMatch->vAddCell (pCLNew);
	  pCLNew->vReplaceFace (pF, pFMatch);
	}
      }
      else {
	// Exactly one of the new cells should be valid
	assert (!pCLthis->qValid () && !pCRthis->qValid ());
	assert (XOR (pCLNew->qValid (), pCRNew->qValid ()));
	// Modify the original bare edge to point to the new cells
	pFMatch->vAssign (pCLNew, pCRNew, pFMatch->pVVert (0),
			  pFMatch->pVVert (1));
	// Correct the connectivity of the (guaranteed to be) one new cell
	// for this face.
	if (pCLthis->qValid ())
	  pCRNew->vReplaceFace (pF, pFMatch);
	else
	  pCLNew->vReplaceFace (pF, pFMatch);
      }
    }				// Face in MeshA matches face in *this

  }				// Done adding all new faces

  delete[]apCNew;
}

//@ Retrieval by index
Face* Mesh2D::
pFFace (const int i) const
{
  assert (i >= 0 && i < iNumFaces ());
  return VAEdgeF.pTEntry (i);
}

BFace* Mesh2D::
pBFBFace (const int i) const
{
  assert (i >= 0 && i < iNumTotalBdryFaces ());
  if (i < iNumBdryFaces())
    return VAEdgeBF.pTEntry(i);
  else
    return VAIntEdgeBF.pTEntry(i - iNumBdryFaces());
}

IntBdryEdge* Mesh2D::
pIBEIntBdryEdge (const int i) const
{
  assert (i >= 0 && i < iNumIntBdryFaces ());
  return VAIntEdgeBF.pTEntry (i);
}

Cell* Mesh2D::
pCCell (const int i) const
{
  assert (i >= 0 && i < iNumCells ());
  if (i < VATri.iSize ())
    return VATri.pTEntry (i);
  else
    return VAQuad.pTEntry (i - VATri.iSize ());
}

int Mesh2D::
iFaceIndex (const Face * const pF) const
{
  assert (pF->eType () == Face::eEdgeFace);
  return (VAEdgeF.iIndex (dynamic_cast<const EdgeFace *>(pF)));
}

int Mesh2D::
iBFaceIndex (const BFace * const pBF) const
{
  switch (pBF->eType()) {
  case Cell::eBdryEdge:
    return VAEdgeBF.iIndex(dynamic_cast<const BdryEdge*>(pBF));
  case Cell::eIntBdryEdge:
    return iNumBdryFaces() + VAIntEdgeBF.iIndex(dynamic_cast<const IntBdryEdge*>(pBF));
  default:
    assert(0);
    return -1;
  }
}

int Mesh2D::
iCellIndex (const Cell * const pC) const
{
  assert ((pC->eType () == Cell::eQuadCell && !qSimplex) ||
	  pC->eType () == Cell::eTriCell);
  if (pC->eType () == Cell::eTriCell)
    return VATri.iIndex (dynamic_cast<const TriCell *>(pC));
  else
    {
      return VAQuad.iIndex (dynamic_cast<const QuadCell *>(pC)) + VATri.iSize ();
    }
}				//@ Creation of new entities

bool Mesh2D::
qIsValidEnt (void* pvEnt) const
  // Determines if the memory location pointed to by the pointer
  // given is valid (within the range) for any of the entity types
  // for this mesh
{
  return ( VAEdgeF.qValid( static_cast<EdgeFace*>(pvEnt) ) ||
	   VATri.qValid( static_cast<TriCell*>(pvEnt) )    ||
	   VAQuad.qValid( static_cast<QuadCell*>(pvEnt) )  ||
	   VAVerts.qValid( static_cast<Vert*>(pvEnt) )     );
}

Face* Mesh2D::
pFNewFace (const int iNV)
{
  assert (iNV == 2);
  Face *pF = VAEdgeF.pTNewEntry ();
  pF->vSetDefaultFlags ();
  return pF;
}

BFace* Mesh2D::
pBFNewBFace (const int iNBV)
{
  assert (iNBV == 2);
  BFace *pBF = VAEdgeBF.pTNewEntry ();
  pBF->vSetDefaultFlags ();
  return pBF;
}

IntBdryEdge* Mesh2D::
pIBENewIntBdryEdge (const int iNBV)
{
  assert (iNBV == 2);
  IntBdryEdge *pIBE = VAIntEdgeBF.pTNewEntry ();
  pIBE->vSetDefaultFlags ();
  return pIBE;
}

Cell* Mesh2D::
pCNewCell (const int iNF)
{
  assert (iNF == 3 || (iNF == 4 && !qSimplex));
  if (iNF == 3) {
    TriCell *pTriC = VATri.pTNewEntry ();
    pTriC->vSetDefaultFlags ();
    pTriC->vSetType (Cell::eTriCell);
    return pTriC;
  }
  else {
    QuadCell *pQuadC = VAQuad.pTNewEntry ();
    pQuadC->vSetDefaultFlags ();
    pQuadC->vSetType (Cell::eQuadCell);
    return pQuadC;
  }
}

//@ Purge unused entries
void Mesh2D::
vPurge ()
{
  SUMAA_LOG_EVENT_BEGIN(PURGE);
  vPurgeCells ();
  vPurgeBdryFaces ();
  vPurgeFaces ();
  vPurgeVerts ();
  SUMAA_LOG_EVENT_END(PURGE);
}

void Mesh2D::
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()) {// If vert is deleted, copy over it.
      aiVertInd[iV] = -1;
    }
    else {
      aiVertInd[iV] = iVActive;
      Vert *pVA = pVVert(iVActive);
      // Copy the data from its old to its new location
      (*pVA) = (*pV);
      iVActive++;
    }
  }
  for (int ii = iVActive ; ii < iNumVerts(); ii++) {
    // Mark all the rest as deleted, just to be on the safe side.
    Vert *pV = pVVert(ii);
    pV->vMarkDeleted();
  }

  // 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 *apVNew[4];
    for (int ii = pF->iNumVerts() - 1; ii >= 0; ii--)
      apVNew[ii] = pVVert(aiVertInd[iVertIndex(pF->pVVert(ii))]);
    // Passing the extra arguments does no harm, as they are ignored.
    pF->vAssign(pF->pCCellLeft(), pF->pCCellRight(), apVNew[0],
		apVNew[1], apVNew[2], apVNew[3]);
  }
  vSetAllHintFaces();

  // Now change the influencing verts, where they occur, to match the
  // new vertex numbering.

  for (iV = iVActive - 1; iV >= 0; iV--) {
    Vert *pV = pVVert(iV);
    assert(pV->qValid() && !pV->qDeleted());
    if (pV->dLS() < 0) {
      // Don't bother copying data if no length scale exists yet.
      pV->vClearLS();
      continue;
    }

    Vert *pVInf;
    BdryPatch *pBPJunk;
    pV->vGetClosest(&pVInf, &pBPJunk);
    if (pVInf->qValid()) {
      pVInf = pVVert(aiVertInd[iVertIndex(pVInf)]);
      pV->vSetClosest(pVInf);
    }

    pV->vGetSecond(&pVInf, &pBPJunk);
    if (pVInf->qValid()) {
      pVInf = pVVert(aiVertInd[iVertIndex(pVInf)]);
      pV->vSetSecond(pVInf);
    }
  }

  delete [] aiVertInd;

  vSetupVerts(iVActive);
//   vUpdateInfluence();
  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

void Mesh2D::
vPurgeFaces ()
  // Get rid of all unused faces in the connectivity table.
{
  bool qLocalLog = false;
#ifdef SUMAA_LOG
  if (SUMAA_LOG_current != PURGE) {
    qLocalLog = true;
    SUMAA_LOG_EVENT_BEGIN(PURGE);
  }
#endif
  EdgeFace *pFFor = VAEdgeF.pTFirst ();
  EdgeFace *pFBack = VAEdgeF.pTLast ();

  int iNewSize = 0;
  // The case of exactly one valid face
  if (pFFor == pFBack &&
      pFFor != pFInvalidFace &&
      !pFFor->qDeleted ())
    iNewSize++;

  while (pFFor != pFBack) {
    while (!pFFor->qDeleted () && pFFor != pFBack) {
      pFFor = VAEdgeF.pTNext (pFFor);
      iNewSize++;
    }

    while (pFBack->qDeleted () && pFFor != pFBack)
      pFBack = VAEdgeF.pTPrev (pFBack);

    if (pFFor != pFBack) {
      assert (pFFor->qDeleted ());
      iNewSize++;
      (*pFFor) = (*pFBack);
      assert (!pFFor->qDeleted ());
      pFBack->vMarkDeleted ();
      pFFor = VAEdgeF.pTNext (pFFor);
    }
    // Allow for the case where nothing has been deleted
    else if (!pFFor->qDeleted ()) {
      iNewSize++;
    }
  }				// Loop over all faces
  // This resets the number of entries in use to iNewSize.

  VAEdgeF.vSetup (iNewSize);
  vResetPeriodicData();
  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

void Mesh2D::
vPurgeCells ()
  // Get rid of all unused cells in the connectivity table.
{
  bool qLocalLog = false;
#ifdef SUMAA_LOG
  if (SUMAA_LOG_current != PURGE) {
    qLocalLog = true;
    SUMAA_LOG_EVENT_BEGIN(PURGE);
  }
#endif
  {
    TriCell *pCFor = VATri.pTFirst ();
    TriCell *pCBack = VATri.pTLast ();
    int iNewSize = 0;
    // The case of exactly one valid cell
    if (pCFor == pCBack &&
	pCFor != pCInvalidCell &&
	!pCFor->qDeleted ()) {
      iNewSize = 1;
    }

    while (pCFor != pCBack) {
      while (!pCFor->qDeleted () && pCFor != pCBack) {
	pCFor = VATri.pTNext (pCFor);
	iNewSize++;
      }

      while (pCBack->qDeleted () && pCFor != pCBack)
	pCBack = VATri.pTPrev (pCBack);

      if (pCFor != pCBack) {
	assert (pCFor->qDeleted ());
	iNewSize++;
	(*pCFor) = (*pCBack);
	pCBack->vMarkDeleted ();
#ifndef NDEBUG
	for (int i = 0; i < pCFor->iNumFaces (); i++)
	  assert (pCFor->pFFace (i)->iFullCheck ());
#endif
	pCFor = VATri.pTNext (pCFor);
      }
      // Allow for the case where nothing has been deleted
      else if (!pCFor->qDeleted ())
	iNewSize++;
    }				// Loop over triangles
    // Reset the number of tri cells in use

    VATri.vSetup (iNewSize);
  }				// Contex for loop over triangles

  {
    QuadCell *pCFor = VAQuad.pTFirst ();
    QuadCell *pCBack = VAQuad.pTLast ();
    int iNewSize = 0;
    // The case of exactly one valid cell
    if (pCFor == pCBack &&
	pCFor != pCInvalidCell &&
	!pCFor->qDeleted ())
      iNewSize++;
    while (pCFor != pCBack) {
      while (!pCFor->qDeleted () && pCFor != pCBack) {
	pCFor = VAQuad.pTNext (pCFor);
	iNewSize++;
      }

      while (pCBack->qDeleted () && pCFor != pCBack)
	pCBack = VAQuad.pTPrev (pCBack);

      if (pCFor != pCBack) {
	assert (pCFor->qDeleted ());
	iNewSize++;

	(*pCFor) = (*pCBack);
	pCBack->vMarkDeleted ();
#ifndef NDEBUG
	for (int i = 0; i < pCFor->iNumFaces (); i++)
	  assert (pCFor->pFFace (i)->iFullCheck ());
#endif
	pCFor = VAQuad.pTNext (pCFor);
      }
      // Allow for the case where nothing has been deleted
      else if (!pCFor->qDeleted ())
	iNewSize++;
    }				// Loop over quads
    // Reset the number of quad cells in use

    VAQuad.vSetup (iNewSize);
  }				// Context for loop over quads
  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

void Mesh2D::
vPurgeBdryFaces ()
  // Get rid of all unused bdry faces in the connectivity table.
{
  bool qLocalLog = false;
#ifdef SUMAA_LOG
  if (SUMAA_LOG_current != PURGE) {
    qLocalLog = true;
    SUMAA_LOG_EVENT_BEGIN(PURGE);
  }
#endif
  // First the "real" bdry edges
  {
    BdryEdge *pEBFFor = VAEdgeBF.pTFirst ();
    BdryEdge *pEBFBack = VAEdgeBF.pTLast ();

    int iNewSize = 0;

    while (pEBFFor != pEBFBack) {
      while (!pEBFFor->qDeleted () && pEBFFor != pEBFBack) {
	pEBFFor = VAEdgeBF.pTNext (pEBFFor);
	iNewSize++;
      }

      while (pEBFBack->qDeleted () && pEBFFor != pEBFBack)
	pEBFBack = VAEdgeBF.pTPrev (pEBFBack);

      if (pEBFFor != pEBFBack) {
	assert (pEBFFor->qDeleted ());
	iNewSize++;
	(*pEBFFor) = (*pEBFBack);
	pEBFBack->vMarkDeleted ();
#ifndef NDEBUG
	for (int i = 0; i < pEBFFor->iNumFaces (); i++)
	  assert (pEBFFor->pFFace (i)->iFullCheck ());
#endif
	pEBFFor = VAEdgeBF.pTNext (pEBFFor);
      }
      else if (!pEBFFor->qDeleted ())
	iNewSize++;
    }
    // Reset the number of bdry faces in use
    VAEdgeBF.vSetup (iNewSize);
  }

  // Now the internal bdry faces
  {
    int iNewSize = 0;
    IntBdryEdge *pIEBFFor = VAIntEdgeBF.pTFirst ();
    IntBdryEdge *pIEBFBack = VAIntEdgeBF.pTLast ();
    // The case of exactly one valid internal bdry edge
    if (pIEBFFor == pIEBFBack &&
	pIEBFFor != pBFInvalidBFace &&
	!pIEBFFor->qDeleted ()) {
      iNewSize = 1;
    }

    while (pIEBFFor != pIEBFBack) {
      while (!pIEBFFor->qDeleted () && pIEBFFor != pIEBFBack) {
	pIEBFFor = VAIntEdgeBF.pTNext (pIEBFFor);
	iNewSize++;
      }

      while (pIEBFBack->qDeleted () && pIEBFFor != pIEBFBack)
	pIEBFBack = VAIntEdgeBF.pTPrev (pIEBFBack);

      if (pIEBFFor != pIEBFBack) {
	assert (pIEBFFor->qDeleted ());
	iNewSize++;
	(*pIEBFFor) = (*pIEBFBack);
	pIEBFBack->vMarkDeleted ();
#ifndef NDEBUG
	for (int i = 0; i < pIEBFFor->iNumFaces (); i++)
	  assert (pIEBFFor->pFFace (i)->iFullCheck ());
#endif
	pIEBFFor = VAIntEdgeBF.pTNext (pIEBFFor);
      }
      else if (!pIEBFFor->qDeleted ())
	iNewSize++;
    }
    // Reset the number of bdry faces in use
    VAIntEdgeBF.vSetup (iNewSize);
  }

  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

//@ Validity checking of entire mesh
bool Mesh2D::
qValid () const
{
  bool qRetVal = true;
  int i;

  for (i = iNumVerts () - 1; i >= 0 && qRetVal; i--)
    {
      Vert *pV = pVVert (i);
      qRetVal = qRetVal && pV->qValid ();
      if (!pV->qDeleted ()) {		// If the vertex still belongs to the mesh...

	Face *pF = pV->pFHintFace ();
	if (pF->qValid ())		// and the vertex has a hint face defined...

	  qRetVal = qRetVal && pF->qHasVert (pV);	// that hint has to be correct.

      }
    }

  for (i = iNumFaces () - 1; i >= 0 && qRetVal; i--) {
    Face *pF = pFFace (i);
    if (!(pF->qValid ())) {
      qRetVal = false;
      break;
    }
    if (pF->qDeleted () || !pF->pCCellLeft ()->qValid () ||
	!pF->pCCellRight ()->qValid ())
      continue;
    int iLoc = pF->iFaceLoc ();
    switch (iLoc) {
    case (Face:: eInterior):
      if (pF->pCCellLeft ()->iRegion () !=
	  pF->pCCellRight ()->iRegion ())
	qRetVal = false;
      break;
    case (Face:: eBdryTwoSide):
      //      if (pF->pCCellLeft()->iRegion() == pF->pCCellRight()->iRegion())
      //      qRetVal = false;
      // Must make sure that at least one side is an IntBdryEdge...
      if ((pF->pCCellLeft ()->eType () != Cell::eIntBdryEdge) &&
	  (pF->pCCellRight ()->eType () != Cell::eIntBdryEdge)) {
	qRetVal = false;
      }
      break;
    default:
      break;
    }
    if (pF->qDeleted ()) continue;
    //      if (!pF->pCCellLeft ()->qValid () ||
    //	!pF->pCCellRight ()->qValid ())
    //        continue;
    if (pF->iFullCheck() == 0) qRetVal = false;
  }

  for (i = iNumBdryFaces () - 1; i >= 0 && qRetVal; i--)
    qRetVal = qRetVal && pBFBFace (i)->qValid ();

  for (i = iNumCells () - 1; i >= 0 && qRetVal; i--) {
    Cell *pC = pCCell (i);
    qRetVal = qRetVal && pC->qValid ();
    if (pC->qDeleted ())
      continue;
    bool qRegOK = (pC->iRegion () != iInvalidRegion);
    qRetVal = qRetVal && qRegOK;
    // Check that the cell is actually closed
    qRetVal = qRetVal && pC->qIsClosed ();
  }

  //   qRetVal = qRetVal && qWatertight();
  return qRetVal;
}

//@ Computation of perimeter, total surface area, and watertightness of surface
double Mesh2D::
dBoundarySize () const
{
  double dRetVal = 0;
  for (int i = iNumBdryFaces () - 1; i >= 0; i--)
    dRetVal += pBFBFace (i)->dSize ();
  return dRetVal;
}

double Mesh2D::
dInteriorSize () const
{
  double dRetVal = 0;
  for (int i = iNumCells () - 1; i >= 0; i--)
    dRetVal += pCCell (i)->dSize ();
  return dRetVal;
}

bool Mesh2D::
qWatertight () const
{
  double adSum[] =
  {0., 0., 0.};
  double dArea = 0;
  bool qAnyCells = false;
  for (int i = iNumCells () - 1; i >= 0; i--)
    {
      Cell *pC = pCCell (i);
      // Only compute the directed surface area for cells which have valid
      // connectivity.  This way, surface meshes with deleted cells (ie,
      // no valid faces) can be checked for watertightness.
      if (pC->iFullCheck () == 1) {
	qAnyCells = true;
	double adTemp[3];
	pC->vVecSize (adTemp);
	adSum[0] += adTemp[0];
	adSum[1] += adTemp[1];
	adSum[2] += adTemp[2];
	dArea += pC->dSize ();
      }
    } if (!qAnyCells)
      return true;
  double dInvArea = 1. / dArea;
  vSCALE3D (adSum, dInvArea);
  if (iFuzzyComp (adSum[0], 0) == 0 &&
      iFuzzyComp (adSum[1], 0) == 0 &&
      iFuzzyComp (adSum[2], 0) == 0)
    return (true);
  else {
    vMessage (2, "Directed surface areas are non-zero:");
    vMessage (2, "  x: %12.5g\n  y: %12.5g\n  z: %12.5g\n", adSum[0],
	      adSum[1], adSum[2]);
    return (false);
  }
}

//@ Convert from cell-vert data format to face-cell

//@@ This is a general conversion, for [mixed] 2D meshes
// Meshes with only triangles will simply have no quads.
void Mesh2D::
vConvertFromCellVert (const VATriCellCV & VATriCV,
		      const VAQuadCellCV & VAQuadCV,
		      const VAEdgeBFaceCV & VABFCV)
{
  int iNTri = VATriCV.iSize ();
  int iNQuad = VAQuadCV.iSize ();
#ifndef NDEBUG
  int iNCell = iNTri + iNQuad;
#endif
  int iNBEdge = VABFCV.iSize();
  // 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
  int iNEdge = 3 * iNTri + 4 *iNQuad + iNBEdge;
  FaceSort *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);
#ifndef NDEBUG
  assert (iCell == iNCell);
#endif
  std::sort(aFSEdges, aFSEdges + iEdge);

  // If not bdry info was given, count the bdry edges.
  if (iNBEdge == 0) {
    for (i = 0; i < iNEdge; ) {
      // How many faces are incident on this edge?
      int iNumInc = 1;
      while (i + iNumInc < iNEdge &&
	     qIdenticalVerts(aFSEdges[i + iNumInc], aFSEdges[i]))
	iNumInc++;

      if (iNumInc == 1) iNBEdge++;
      i += iNumInc;
    }
  }

  // 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);
  VAEdgeF.vSetup ((3*iNTri + iNBEdge + 4*iNQuad) / 2);
  VAEdgeBF.vSetup (iNBEdge);

  // Read the data off
  int iBFace = 0;
  int iFace = 0;
  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 + iNumInc], aFSEdges[i]))
      iNumInc++;
    assert(i + iNumInc <= iNEdge);

    if (iNumInc == 1) {
      // Create a bdry face out of thin air.
      FaceSort FS = aFSEdges[i];
      assert(FS.iLeftCell() >= 0);
      Face *pF = pFFace(iFace);
      Cell *pC = pCCell(FS.iLeftCell());
      pC->vAddFace(pF);
      BFace *pBF = pBFBFace(iBFace);
      pBF->vAssign(pF);
      pF->vAssign(pBF, pC, pV0, pV1);
      iFace++; iBFace++;
    }
    else {
      if ((pV0 != aFSEdges[i + 1].pVVert(1)) ||
	  (pV1 != aFSEdges[i + 1].pVVert(0))) {
	qOrientationOK = false;
	assert(pV0->iSpaceDimen() == 2);
	vMessage(0, "Warning: Inconsistent orientation of input polygons along edge\n");
	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());
      }

      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 = pFFace(iFace++);

	pF->vAssign (pC0, pC1, pV0, pV1);
	pC0->vAddFace (pF);
	pC1->vAddFace (pF);
      }
    } // Done with matching edges.
    i += iNumInc;
  }
  if (!qOrientationOK) {
    vFatalError ("Polygons do not all have the same orientation",
		 "Mesh2D::vConvertFromCellVert(...)");
  }
  assert (2*iFace == iNBEdge + 3*iNTri + 4*iNQuad);
  assert (iBFace == iNBEdge);
}

typedef struct _vert_vert_ {
  int iV0, iV1;
  _vert_vert_ (const int i0 = -1, const int i1 = -1) :
    iV0(i0), iV1(i1) {}
}
VertVert;

extern "C"
{
  static int qVVComp (const void *pv0, const void *pv1) {
    const VertVert *const pVV0 = (const VertVert *) pv0;
    const VertVert *const pVV1 = (const VertVert *) pv1;
    return (pVV0->iV0 < pVV1->iV0 ||
	    (pVV0->iV0 == pVV1->iV0 && pVV0->iV1 < pVV1->iV1)) ? -1 : 1;
  }
}

typedef struct _renum_ {
  int iV, iDeg;
  _renum_ (const int iVIn = -1, const int iDegIn = 0) :
    iV(iVIn), iDeg(iDegIn) {}
} Renum;

extern "C"
{
  static int qRenumComp2D (const void *pv0, const void *pv1) {
    const Renum *const pR0 = (const Renum *) pv0;
    const Renum *const pR1 = (const Renum *) pv1;
    return ((pR0->iDeg < pR1->iDeg) ? -1 :
	    ((pR0->iDeg == pR1->iDeg) ? 0 : 1));
  }
}


//@ Reorder vertices in a 2D mesh using the reverse Cuthill-McKee algorithm
void Mesh2D::
vReorderVerts_RCM (void)
{
  vMessage(2, "Renumbering vertices...");
  int iVert, iFace, iNVerts = iNumVerts (), iNFaces = iNumFaces ();
  //@@ Construct the vert-vert connectivity
  // Make a list containing each edge (once in each direction) and sort
  VertVert *aVV = new VertVert[2 * iNFaces];
  for (iFace = iNFaces - 1; iFace >= 0; iFace--) {
    Face *pF = pFFace (iFace);
    int iV0 = iVertIndex (pF->pVVert (0));
    int iV1 = iVertIndex (pF->pVVert (1));
    aVV[2 * iFace + 1] = VertVert (iV0, iV1);
    aVV[2 * iFace] = VertVert (iV1, iV0);
  }
  qsort (aVV, 2 * iNFaces, sizeof (VertVert), qVVComp);
  // Tabulate starting locations and degree for each vertex.
  int *aiStart = new int[iNVerts + 1];
  int *aiDegree = new int[iNVerts];
  bool *aqUsed = new bool[iNVerts];
  int iStart = 0, iNAttachedVerts = 0;
  for (iVert = 0; iVert < iNVerts; iVert++) {
    aiStart[iVert] = iStart;
    for (; (iStart < 2*iNFaces) && (aVV[iStart].iV0 == iVert); iStart++) {}
    if (iStart != aiStart[iVert])
      iNAttachedVerts++;
    aiDegree[iVert] = iStart - aiStart[iVert];
    aqUsed[iVert] = false;
  }
  aiStart[iNVerts] = 2 * iNFaces;

  // Set up space to compute re-ordering
  Renum *aRReorder = new Renum[iNVerts];

  int iNumRenumbered = 0;
  while (iNumRenumbered < iNAttachedVerts) {
    //@@ Initialize the front to be a single unrenumbered vertex with lowest degree
    // Find lowest degree vertex
    int iMinDeg = 100, iMinVert = -1;
    for (iVert = 0; iVert < iNAttachedVerts; iVert++) {
      int iDeg = aiStart[iVert + 1] - aiStart[iVert];
      if (iDeg < iMinDeg && iDeg > 0 && !aqUsed[iVert]) {
	iMinDeg = iDeg;
	iMinVert = iVert;
	if (iMinDeg == 2)
	  break;	// Can never be lower than this, so why go on?
      }
    }
    assert (iMinVert != -1);

    // Initialize
    aqUsed[iMinVert] = true;
    aRReorder[iNumRenumbered] = Renum (iMinVert, iMinDeg);
    iNumRenumbered++;
    int iEnd = iNumRenumbered, iBegin = iNumRenumbered-1;

    //@@ Propagate across the mesh
    do {
      //@@@ Add all the next layer of vertices
      int iEnt;
      for (iEnt = iBegin; iEnt < iEnd; iEnt++) {
	int iV = aRReorder[iEnt].iV;
	for (int iNeigh = aiStart[iV]; iNeigh < aiStart[iV + 1]; iNeigh++) {
	  // Grab a candidate
	  assert (aVV[iNeigh].iV0 == iV);
	  int iCand = aVV[iNeigh].iV1;
	  if (!aqUsed[iCand]) {
	    // Add it if it isn't already in the list
	    aRReorder[iNumRenumbered].iV = iCand;
	    aRReorder[iNumRenumbered].iDeg = aiDegree[iCand];
	    aqUsed[iCand] = true;
	    iNumRenumbered++;
	  }
	}			// Loop over neighbors of a given vertex
      }				// Loop over all verts in the most recent layer

      //@@@ Sort the new layer in order of increasing degree
      iBegin = iEnd;
      iEnd = iNumRenumbered;
      qsort (&(aRReorder[iBegin]), static_cast<size_t>(iEnd - iBegin),
	     sizeof (Renum), qRenumComp2D);
    } while (iEnd != iBegin);
  }
  delete [] aVV;
  delete [] aiStart;
  delete [] aiDegree;
  delete [] aqUsed;

  for (iVert = 0; iNumRenumbered != iNVerts; iVert++) {
    if (!aqUsed[iVert]) {
      aRReorder[iNumRenumbered].iV = iVert;
      aqUsed[iVert] = true;
      iNumRenumbered++;
    }
  }

  //@@ Reverse the ordering
  for (iVert = 0; iVert < iNVerts / 2; iVert++) {
    int iTemp = aRReorder[iVert].iV;
    aRReorder[iVert].iV = aRReorder[iNVerts - iVert - 1].iV;
    aRReorder[iNVerts - iVert - 1].iV = iTemp;
  }

  //@@ Move vertex data
  // Create a new copy
  VAVertex VATemp;
  VATemp.vSetup (iNVerts);
  for (iVert = 0; iVert < iNVerts; iVert++) {
    // Copy both coordinates and vertex tags
    Vert *pVTmp = pVVert (iVert);
    VATemp[iVert].vSetCoords (pVTmp->iSpaceDimen (), pVTmp->adCoords ());
    VATemp[iVert].vCopyAllFlags (pVTmp);
    VATemp[iVert].vCopyInfluence(pVTmp);
    VATemp[iVert].vSetLS(pVTmp->dLS());
  }
  for (iVert = 0; iVert < iNVerts; iVert++) {
    // Copy back into re-ordered locations
    Vert *pVSource = &(VATemp[aRReorder[iVert].iV]);
    pVVert (iVert)->vSetCoords (pVSource->iSpaceDimen (), pVSource->adCoords ());
    pVVert (iVert)->vCopyAllFlags (pVSource);
    pVVert (iVert)->vCopyInfluence(pVSource);
    pVVert (iVert)->vSetLS(pVSource->dLS());
  }

  //@@ Re-label all the face data
  int *aiOld2New = new int[iNVerts];
  for (iVert = 0; iVert < iNVerts; iVert++)
    aiOld2New[aRReorder[iVert].iV] = iVert;
  delete[]aRReorder;

  for (iFace = 0; iFace < iNFaces; iFace++) {
    Face *pF = pFFace (iFace);
    Vert *apVNew[2];		// Only two verts to an edge

    apVNew[0] = pVVert (aiOld2New[iVertIndex (pF->pVVert (0))]);
    apVNew[1] = pVVert (aiOld2New[iVertIndex (pF->pVVert (1))]);
    pF->vAssign (pF->pCCellLeft (), pF->pCCellRight (), apVNew[0],
		 apVNew[1]);
  }
  delete[]aiOld2New;

  //@@ Reset vertex hints
  vSetAllHintFaces ();
  if (!(qValid ())) {
    vMessage(2, "aborted!\n");
    vFoundBug ("Vertex re-ordering for 2D meshes");
  }
  vMessage(2, "done.\n");
}

void Mesh2D::vReorderCells()
  // Reorder cells according to vertex indices, so that cells with the
  // lowest sum of vertex indices are first in order.  The hardest part
  // is actually to get the faces and cells reconnected afterwards.
{
  vMessage(2, "Renumbering cells...");
  assert(qSimplicial());
  multimap<int, int> mOrderedCells;
  int iC;
  for (iC = 0; iC < iNumCells(); iC++) {
    Cell *pC = pCCell(iC);
    int iSum = 0, iNV = pC->iNumVerts();
    for (int iV = 0; iV < iNV; iV++) {
      iSum += iVertIndex(pC->pVVert(iV));
    }
    int iMean = iSum / iNV;
    mOrderedCells.insert(pair<int, int>(iMean, iC));
  }
  // Now we have a list of all the cells, sorted in order by smallest
  // average of vertex indices!  Use this to create a map that tells
  // which old cell goes into which slot in the new numbering.
  int *aiCellMapping = new int[iNumCells()];
  multimap<int, int>::iterator iter;
  iC = 0;
  for (iter = mOrderedCells.begin(); iter != mOrderedCells.end(); iter++) {
    vMessage(4, "Old cell: %d  New cell: %d  Sum: %d\n",
	     iter->second, iC, iter->first);
    aiCellMapping[iC] = iter->second;
    iC++;
  }
  assert(iC == iNumCells());

  // Finally, update the cell information.  Fundamentally, re-ordering
  // the cells results in one or more circular permutations of cells
  // (A->B->Q->L->A, etc).  So follow each of these.
  bool *aqUpdated = new bool[iNumCells()];
  for (iC = 0; iC < iNumCells(); iC++) {
    aqUpdated[iC] = false;
  }

  int iNumUpdated = 0;
  int iFirst = 0;
  while (iNumUpdated != iNumCells()) {
    while (aqUpdated[iFirst]) iFirst++;
    assert(iFirst < iNumCells());
    if (aiCellMapping[iFirst] == iFirst) {
      aqUpdated[iFirst] = true;
      iNumUpdated++;
    }
    else {
      assert(pCCell(iFirst)->eType() == Cell::eTriCell);
      TriCell TC;
      TC = *(dynamic_cast<TriCell*>(pCCell(iFirst)));
      int iNext = iFirst;
      while (aiCellMapping[iNext] != iFirst) {
	(*pCCell(iNext)) = (*pCCell(aiCellMapping[iNext]));
	aqUpdated[iNext] = true;
	iNumUpdated++;
	iNext = aiCellMapping[iNext];
      }
      (*pCCell(iNext)) = TC;
      aqUpdated[iNext] = true;
      iNumUpdated++;
    }
  }
  if (!(qValid ())) {
    vMessage(2, "aborted!\n");
    vFoundBug ("Cell re-ordering for 2D meshes");
  }
  vMessage(2, "done.\n");
  delete [] aiCellMapping;
  delete [] aqUpdated;
}

void Mesh2D::vReorderFaces()
  // Reorder faces according to vertex indices, so that faces with the
  // lowest sum of vertex indices are first in order.  The hardest part
  // is actually to get the faces and cells reconnected afterwards.
{
  vMessage(2, "Renumbering faces...");
  assert(qSimplicial());
  multimap<int, int> mOrderedFaces;
  int iF, iNF = iNumFaces();
  for (iF = 0; iF < iNF; iF++) {
    Face *pF = pFFace(iF);

    int iSum = iVertIndex(pF->pVVert(0)) + iVertIndex(pF->pVVert(1));
    mOrderedFaces.insert(pair<int, int>(iSum, iF));
  }
  // Now we have a list of all the faces, sorted in order by smallest
  // average of vertex indices!  Use this to create a map of which -old-
  // face becomes which -new- face.
  int *aiFaceMapping = new int[iNF];
  multimap<int, int>::iterator iter;
  iF = 0;
  for (iter = mOrderedFaces.begin(); iter != mOrderedFaces.end(); iter++) {
    vMessage(4, "Old face: %d  New face: %d  Sum: %d\n",
	     iter->second, iF, iter->first);
    aiFaceMapping[iF] = iter->second;
    iF++;
  }
  assert(iF == iNF);

  // Finally, update the cell information.  Fundamentally, re-ordering
  // the cells results in one or more circular permutations of cells
  // (A->B->Q->L->A, etc).  So follow each of these.
  bool *aqUpdated = new bool[iNF];
  for (iF = 0; iF < iNF; iF++) {
    aqUpdated[iF] = false;
  }

  int iNumUpdated = 0;
  int iFirst = 0;
  while (iNumUpdated != iNF) {
    while (aqUpdated[iFirst]) iFirst++;
    assert(iFirst < iNF);
    if (aiFaceMapping[iFirst] == iFirst) {
      aqUpdated[iFirst] = true;
      iNumUpdated++;
    }
    else {
      assert(pFFace(iFirst)->eType() == Face::eEdgeFace);
      EdgeFace EF;
      EF = *(dynamic_cast<EdgeFace*>(pFFace(iFirst)));
      int iNext = iFirst;
      while (aiFaceMapping[iNext] != iFirst) {
	(*pFFace(iNext)) = (*pFFace(aiFaceMapping[iNext]));
	aqUpdated[iNext] = true;
	iNumUpdated++;
	iNext = aiFaceMapping[iNext];
      }
      (*pFFace(iNext)) = EF;
      aqUpdated[iNext] = true;
      iNumUpdated++;
    }
  }

  if (!(qValid ())) {
    vMessage(2, "aborted!\n");
    vFoundBug ("Face re-ordering for 2D meshes");
  }
  vMessage(2, "done.\n");
  delete [] aqUpdated;
  delete [] aiFaceMapping;
}

//@ Read 2D mesh file 
Mesh2D::Mesh2D (const char strBaseFileName[], const int iQualMeas,
		const Bdry2D* pB2DIn)
  :Mesh(), VAEdgeF(), VAEdgeBF(), VATri(), VAQuad(),
   VAIntEdgeBF(), iOnEdge(0), iOnBdryEdge(0), iOnBdryEdge2(0),
   iInterior(0), iInternalBdry(0), iInternalBdry2(0), B2D()
{
  SUMAA_LOG_EVENT_BEGIN (INIT_TRI);
  SUMAA_LOG_EVENT_BEGIN (READ_MESH_2D);
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(2, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);

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

  vReadFile_Mesh2D (strBaseFileName, iNVerts, iNFaces, iNCells,
		    iNBFaces,  iNIntBFaces, qFaceVert, qCellVert,
		    qFaceCell, qCellFace, qCellRegion, qBFaceFace,
		    qBFaceVert, qBFaceBC, qIntBFaceFace,
		    qIntBFaceVert, qIntBFaceBC, VAVerts,
		    a2iFaceVert, a2iFaceCell, a2iCellVert,
		    a2iCellFace, aiCellRegion, aiBFaceFace,
		    aiBFaceBC, a2iBFaceVert, a2iIntBFaceFace,
		    aiIntBFaceBC, a2iIntBFaceVert);

  if (VAVerts.iSize () == 0)
    vFatalError ("No vertex coordinate data read", "2D 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.",
		 "2D 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.", "2D 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];
	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);
      }
      if (iNCells != 0 && iNCells != iTmpNCells + 1)
	vFatalError ("Number of cells in mesh doesn't match number given.",
		     "2D mesh input");
      if ((iNBFaces != 0 || iNIntBFaces != 0) &&
	  iNBFaces + 2*iNIntBFaces != iTmpNBFaces) // sic
	vFatalError ("Number of bdry faces in mesh doesn't match number given.",
		     "2D mesh input");
      if ((iNVerts != 0 && iNVerts != iTmpNVerts + 1)
	  || iTmpNVerts + 1 != VAVerts.iSize())
	vFatalError ("Number of verts in mesh doesn't match number given.",
		     "2D mesh input");
      iNCells = iTmpNCells + 1;
      iNBFaces = iTmpNBFaces - 2 * iNIntBFaces;
      iNVerts = iTmpNVerts + 1;
    }

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

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

      //@@@ 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.
      delete [] aiCellRegion;
    }

    // Now do assignment for face->cell/vert and cell->face.
    for (int i = 0; i < iNFaces; i++) {
      int iV0 = a2iFaceVert[i][0];
      int iV1 = a2iFaceVert[i][1];

      int iC0 = a2iFaceCell[i][0];
      int iC1 = a2iFaceCell[i][1];

      Face *pF = pFFace (i);
      Cell *pC0, *pC1;
      enum Face::eFaceLoc eFL;

      if (iC0 < 0 || iC1 < 0) {
	eFL = Face::eBdryFace;  // Might be internal bdry, but that's fixed
				// later.
	if (iC0 >= 0) {
	  pC0 = pCCell(iC0);
	  pC0->vAddFace(pF);
	  pC1 = pBFInvalidBFace;
	}
	else {
	  pC0 = pBFInvalidBFace;
	  pC1 = pCCell(iC1);
	  pC1->vAddFace(pF);
	}
      }
      else {
	eFL = Face::eInterior;
	pC0 = pCCell (iC0);
	pC1 = pCCell (iC1);
	pC0->vAddFace(pF);
	pC1->vAddFace(pF);
      }

      // Tell the face everything it needs to know.
      pFFace (i)->vAssign (pC0, pC1, pVVert (iV0), pVVert (iV1));
      if (pC0 != pBFInvalidBFace && pC1 != pBFInvalidBFace) {
	// These verts need to have their hint face set, so that the
	// face can be retrieved by getting the hint when handling
	// BC's.
	pVVert(iV0)->vSetHintFace(pFFace(i));
	pVVert(iV1)->vSetHintFace(pFFace(i));
      }
      pFFace (i)->vSetFaceLoc(eFL);
    } // Done assigning for this face.

    // Set up bdry faces, but not the Bdry2D stuff yet.

    // Now the internal bdry faces, if any.
    if (qIntBFaceFace) {
      for (int iIBF = 0; iIBF < iNIntBFaces; iIBF++) {
	// This array has to have live data, or we wouldn't be here.
	Face *pF0 = pFFace(a2iIntBFaceFace[iIBF][0]);
	Face *pF1 = pFFace(a2iIntBFaceFace[iIBF][1]);
	pF0->vSetFaceLoc(Face::eBdryTwoSide);
	pF1->vSetFaceLoc(Face::eBdryTwoSide);
	
	IntBdryEdge *pIBE = pIBEIntBdryEdge(iIBF);
	pIBE->vAssign(pF0, pF1);
	pF0->vAddCell(pIBE);
	pF1->vAddCell(pIBE);
	
	// Need this both for checking and for setting up patches.
	Vert *pV0 = pIBE->pVVert(0);
	Vert *pV1 = pIBE->pVVert(1);
	
	pV0->vSetType(Vert::eBdryTwoSide);
	pV1->vSetType(Vert::eBdryTwoSide);
	
	// Verts may or may not be given.  If they are given, they must
	// match what the faces tell us.
	if (qIntBFaceVert) {
	  Vert *pV0_Read = pVVert(a2iIntBFaceVert[iIBF][0]);
	  Vert *pV1_Read = pVVert(a2iIntBFaceVert[iIBF][1]);
	  
	  assert(pV0 == pV0_Read);
	  assert(pV1 == pV1_Read);
	}
      } // Done with setting up internal bdry faces w/ face links
    }
    
    // Now turn to regular bdry faces
    for (int iBF = 0; iBF < iNBFaces; iBF++) {
      Face *pF;
      Vert *pV0, *pV1;
      assert(qBFaceFace || qBFaceVert);
      if (qBFaceFace) {
	// This array has to have live data, or we wouldn't be here.
	pF = pFFace(aiBFaceFace[iBF]);
	pV0 = pF->pVVert(0);
	pV1 = pF->pVVert(1);
	
	if (qBFaceVert) {
	  // Do some consistency checking.
	  assert((pV0 == pVVert(a2iBFaceVert[iBF][0]) &&
		  pV1 == pVVert(a2iBFaceVert[iBF][1]))
		 ||
		 (pV0 == pVVert(a2iBFaceVert[iBF][1]) &&
		  pV1 == pVVert(a2iBFaceVert[iBF][0])));
	}
      }
      else {
	assert(qBFaceVert && !qBFaceFace);
	
	pV0 = pVVert(a2iBFaceVert[iBF][0]);
	pV1 = pVVert(a2iBFaceVert[iBF][1]);
	pF = pV0->pFHintFace();
	if (!pF->qHasVert(pV1)) {
	  pF = pV1->pFHintFace();
	  if (!pF->qHasVert(pV0)) {
	    // Now we're in for it: need to walk around the vertex to
	    // find the actually face we're after.
	    Vert *pVTmp = pF->pVVert (0);
	    if (pVTmp == pV1)
	      pVTmp = pF->pVVert (1);
	    Cell *pCTmp = pF->pCCellLeft ();
	    if (pCTmp->eType () == Cell::eBdryEdge)
	      pCTmp = pF->pCCellRight ();
	    while (!pF->qHasVert (pV0)) {
	      assert (pCTmp->eType () == Cell::eTriCell);
	      Vert *pVNew = dynamic_cast<TriCell*>(pCTmp)->pVVertOpposite(pF);
	      pF = dynamic_cast<TriCell*>(pCTmp)->pFFaceOpposite (pVTmp);
	      pVTmp = pVNew;
	      pCTmp = pF->pCCellOpposite (pCTmp);
	    }; // Done finding the right boundary face the hard way
	  }
	}
	assert(pF->qHasVert(pV0));
	assert(pF->qHasVert(pV1));
      }
      pF->vSetFaceLoc(Face::eBdryFace);
      
      BFace *pBF = pBFBFace(iBF);
      pBF->vAssign(pF);
      pF->vAddCell(pBF);
      
      if (pV0->iVertType() == Vert::eBdryTwoSide) {
	pV0->vSetType(Vert::eBdryApex);
      }
      else {
	pV0->vSetType(Vert::eBdry);
      }
      if (pV1->iVertType() == Vert::eBdryTwoSide) {
	pV1->vSetType(Vert::eBdryApex);
      }
      else {
	pV1->vSetType(Vert::eBdry);
      }
    } // Done setting up all BFace stuff

    if (pB2DIn) {
      // Copy boundary data
      B2D = (*pB2DIn);
      vMessage(2, "Number of boundary patches: %d\n", B2D.iNumPatches());

      // Now infer which bdry faces are on which patches.  Start by
      // identifying all bdry verts, and finding which patches, if any,
      // they fall on.  While it's not very efficient to do this by
      // checking all bdry verts against all patches, this approach has
      // the virtue of working, and for meshes that aren't absolutely
      // huge, it should be no problem from the runtime point of view.
      multimap<Vert*, BdryPatch2D*> mmVertPatch;
      for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
	Vert *pV = pVVert(iV);
// 	vMessage(2, "Vert at (%10.6f, %10.6f) lies on patch %10p\n",
// 		 pV->dX(), pV->dY(), pBP2D);
	for (int iP = B2D.iNumPatches() - 1; iP >= 0; iP--) {
	  BdryPatch2D *pBP2D = B2D[iP];
	  if (!pBP2D->qDisjointFrom(pV)) {
	    vMessage(2, "Vert at (%10.6f, %10.6f) lies on patch %10p\n",
		     pV->dX(), pV->dY(), pBP2D);
	    mmVertPatch.insert(pair<Vert*, BdryPatch2D*>(pV, pBP2D));
	  }
	}
      }

      // Definitively identify bdry apexes: any vert on multiple patches
      // is an apex.
      for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
	Vert *pV = pVVert(iV);
	if (mmVertPatch.count(pV) > 1) {
	  pV->vSetType(Vert::eBdryApex);
	}
      }

      // Now for every bdry face, identify a patch common to its verts.
      for (int iBF = iNumBdryFaces() - 1; iBF >= 0; iBF--) {
	BFace *pBF = pBFBFace(iBF);
	Vert *pV0 = pBF->pVVert(0);
	Vert *pV1 = pBF->pVVert(1);
	// How many patches for each one?

	pair<multimap<Vert*, BdryPatch2D*>::iterator,
	  multimap<Vert*, BdryPatch2D*>::iterator> pRange0, pRange1;
	pRange0 = mmVertPatch.equal_range(pV0);
	pRange1 = mmVertPatch.equal_range(pV1);
	multimap<Vert*, BdryPatch2D*>::iterator iter0, iter1;
	iter0 = pRange0.first;
	BdryPatch2D* pBP2D = NULL;
	for (iter0 = pRange0.first; 
	     pBP2D == NULL && iter0 != pRange0.second;
	     iter0++) {
	  for (iter1 = pRange1.first;
	       pBP2D == NULL && iter1 != pRange1.second;
	       iter1++) {
	    if (iter0->second == iter1->second) {
	      pBP2D = iter0->second;
	    }
	  }
	}
	assert(pBP2D != NULL);
	vMessage(2, "BFace %3d lies on patch %10p\n", iBF, pBP2D);
	vMessage(2, "  coords: (%10.6f, %10.6f) to (%10.6f, %10.6f)\n",
		 pV0->dX(), pV0->dY(), pV1->dX(), pV1->dY());
	pBF->vSetPatch(pBP2D);
      }

      // Finally, reset all the periodic data, if any.
      vResetPeriodicData();

      // Check the periodicity
#ifndef NDEBUG
      for (int iBF = 0; iBF < iNumBdryFaces(); iBF++) {
	BFace* pBF = pBFBFace(iBF);
	BdryPatch2D* pBP2D = dynamic_cast<BdryPatch2D*>(pBF->pPatchPointer());
	if (pBP2D->qIsPeriodic()) {
	  Periodic* pP = pBP2D->pPGetPeriodicData();
	  const BdryPatch2D* pBP2DOther = dynamic_cast<const BdryPatch2D*>
	    (pP->pBPOtherPatch(pBP2D));
	  vMessage(2, "BFace %d is on patch %p opposite patch %p\n",
		   iBF, pBP2D, pBP2DOther);
	}
	else {
	  vMessage(2, "BFace %d is on patch %p\n", iBF, pBP2D);
	}
      }
      vMessage(2, "%d bfaces, %d patches, %d periodic patches\n", 
	       iNumBdryFaces(), B2D.iNumPatches(), B2D.iNumPeriodic());
#endif
      qLengthScaleFromCellSizes = false;
    }
    else {
      // Recreate boundary data by inference.

      // 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 iPatch = 0, iBPt = 0;
      for (int ii = 0; ii < iNVerts; aiVertMap[ii++] = -1) {}

      B2D.vSetSize(iNBFaces + iNIntBFaces, iNVerts);

      // Now the internal bdry faces, if any.
      if (qIntBFaceFace) {
	for (int iIBF = 0; iIBF < iNIntBFaces; iIBF++) {
	  // This array has to have live data, or we wouldn't be here.
 	  Face *pF0 = pFFace(a2iIntBFaceFace[iIBF][0]);
 	  Face *pF1 = pFFace(a2iIntBFaceFace[iIBF][1]);

	  IntBdryEdge *pIBE = pIBEIntBdryEdge(iIBF);

	  // Need this both for checking and for setting up patches.
	  Vert *pV0 = pIBE->pVVert(0);
	  Vert *pV1 = pIBE->pVVert(1);

	  // Need to create a patch for this bdry edge
	  int iV0 = iVertIndex(pV0);
	  int iV1 = iVertIndex(pV1);

	  if (aiVertMap[iV0] == -1) {
	    aiVertMap[iV0] = iBPt;
	    B2D.vAddPoint(pVVert(iV0)->adCoords());
	    B2D.vSetBdryPoint(iBPt, false);
	    iBPt++;
	  }

	  if (aiVertMap[iV1] == -1) {
	    aiVertMap[iV1] = iBPt;
	    B2D.vAddPoint(pVVert(iV1)->adCoords());
	    B2D.vSetBdryPoint(iBPt, false);
	    iBPt++;
	  }

	  int iBP0 = aiVertMap[iV0];
	  int iBP1 = aiVertMap[iV1];
	  assert(iBP0 >= 0 && iBP1 >= 0);

	  // Grab both region tags, working from IBF -> face -> cell.
	  Cell *pCL = pF0->pCCellOpposite(pIBE);
	  int iLeftReg = pCL->iRegion();
	  Cell *pCR = pF1->pCCellOpposite(pIBE);
	  int iRightReg = pCR->iRegion();
	  BdryPatch2D *pBP =
	    new BdrySeg(&B2D, iInvalidBC, iInvalidBC, iLeftReg, iRightReg,
			iBP0, iBP1);
	  B2D.vAddPatchToList(pBP);
	  // vSetPatch can handle getting the patch for an internal bdry
	  // fixed up, too.
	  pIBE->vSetPatch(iPatch, &B2D);
	  iPatch++;
	}
      } // Done with setting up internal bdry faces w/ face links

      // Now turn to regular bdry faces
      for (int iBF = 0; iBF < iNBFaces; iBF++) {
	Face *pF;
	Vert *pV0, *pV1;
	assert(qBFaceFace || qBFaceVert);
	if (qBFaceFace) {
	  // This array has to have live data, or we wouldn't be here.
	  pF = pFFace(aiBFaceFace[iBF]);
	  pV0 = pF->pVVert(0);
	  pV1 = pF->pVVert(1);

	}
	else {
	  assert(qBFaceVert && !qBFaceFace);

	  pV0 = pVVert(a2iBFaceVert[iBF][0]);
	  pV1 = pVVert(a2iBFaceVert[iBF][1]);
	  pF = pV0->pFHintFace();
	  if (!pF->qHasVert(pV1)) {
	    pF = pV1->pFHintFace();
	    if (!pF->qHasVert(pV0)) {
	      // Now we're in for it: need to walk around the vertex to
	      // find the actually face we're after.
	      Vert *pVTmp = pF->pVVert (0);
	      if (pVTmp == pV1)
		pVTmp = pF->pVVert (1);
	      Cell *pCTmp = pF->pCCellLeft ();
	      if (pCTmp->eType () == Cell::eBdryEdge)
		pCTmp = pF->pCCellRight ();
	      while (!pF->qHasVert (pV0)) {
		assert (pCTmp->eType () == Cell::eTriCell);
		Vert *pVNew = dynamic_cast<TriCell*>(pCTmp)->pVVertOpposite(pF);
		pF = dynamic_cast<TriCell*>(pCTmp)->pFFaceOpposite (pVTmp);
		pVTmp = pVNew;
		pCTmp = pF->pCCellOpposite (pCTmp);
	      }; // Done finding the right boundary face the hard way
	    }
	  }
	  assert(pF->qHasVert(pV0));
	  assert(pF->qHasVert(pV1));
	}
	pF->vSetFaceLoc(Face::eBdryFace);

	BFace *pBF = pBFBFace(iBF);

	// Need to create a patch for this bdry edge
	int iV0 = iVertIndex(pV0);
	int iV1 = iVertIndex(pV1);

	if (aiVertMap[iV0] == -1) {
	  aiVertMap[iV0] = iBPt;
	  B2D.vAddPoint(pVVert(iV0)->adCoords());
	  B2D.vSetBdryPoint(iBPt, false);
	  iBPt++;
	}

	if (aiVertMap[iV1] == -1) {
	  aiVertMap[iV1] = iBPt;
	  B2D.vAddPoint(pVVert(iV1)->adCoords());
	  B2D.vSetBdryPoint(iBPt, false);
	  iBPt++;
	}

	int iBP0 = aiVertMap[iV0];
	int iBP1 = aiVertMap[iV1];
	assert(iBP0 >= 0 && iBP1 >= 0);

	// Grab the BC from the incoming face data.
	int iF = iFaceIndex(pF);
	int iC, iBC;
	if (a2iFaceCell[iF][0] < 0) {
	  iBC = abs(a2iFaceCell[iF][0]);
	  iC  =     a2iFaceCell[iF][1];
	}
	else {
	  iBC = abs(a2iFaceCell[iF][1]);
	  iC  =     a2iFaceCell[iF][0];
	}
	assert(iBC > 0 && iC >= 0);

	// Find the adjacent region.
	Cell *pC = pCCell(iC);
	int iReg = pC->iRegion();

	BdryPatch2D* pBP;
	if (pF->pCCellLeft() == pC) {
	  // Cell is on the left, bdry on the right.
	  pBP = new BdrySeg(&B2D, iInvalidBC, iBC, iReg, iInvalidRegion,
			    iBP0, iBP1);
	  B2D.vAddPatchToList(pBP);
	  // vSetPatch can handle getting the patch for an internal bdry
	  // fixed up, too.
	  dynamic_cast<BdryEdge*>(pBF)->vSetPatch(iPatch, &B2D);
	  iPatch++;
	}
	else {
	  // Cell is on the right, bdry on the left.
	  pBP = new BdrySeg(&B2D, iBC, iInvalidBC, iInvalidRegion, iReg,
			    iBP0, iBP1);
	  B2D.vAddPatchToList(pBP);
	  // vSetPatch can handle getting the patch for an internal bdry
	  // fixed up, too.
	  dynamic_cast<BdryEdge*>(pBF)->vSetPatch(iPatch, &B2D);
	  iPatch++;
	}
      } // Done setting up all BFace stuff, including building patches.
      delete [] aiVertMap;
      assert (iPatch == iNBFaces + iNIntBFaces);
      qLengthScaleFromCellSizes = true;
    }
    vSetAllHintFaces ();
    if (!(qValid ()))
      vFoundBug ("2D mesh input from user-format face->cell-style file");

    if (!qBFaceFace && !qBFaceVert && !qBFaceBC) {
      vWarning ("Insufficient boundary data read to be useful; ignoring.");
    }
    //@@@ Further consistency checking
    if (qCellFace) {
      vMessage (1, "Checking consistency of cell-face and face-cell data...");
      for (int iC = 0; iC < iNCells; iC++) {
	Cell *pC = pCCell (iC);
	for (int ii = 0; ii < 3; ii++) {
	  int iF = a2iCellFace[iC][ii];
	  Face *pF = pFFace (iF);
	  if (!(pC->qHasFace (pF)))
	    vFatalError ("Cell-face data inconsistent with face-cell data",
			 "2D mesh input");
	}
      }
      vMessage (1, "OK\n");
    }
    if (qCellVert) {
      vMessage (1, "Checking consistency of cell-vert and face-vert data...");
      for (int iC = 0; iC < iNCells; iC++) {
	Cell *pC = pCCell (iC);
	for (int ii = 0; ii < 3; ii++) {
	  int iV = a2iCellVert[iC][ii];
	  Vert *pV = pVVert (iV);
	  if (!(pC->qHasVert (pV)))
	    vFatalError ("Cell-vert data inconsistent with face-vert data",
			 "2D mesh input");
	}
      }
      vMessage (1, "OK\n");
    }

    if (!qValid ())
      vFatalError ("Mesh was invalid.", "2D face-cell connected mesh input");
  } // Done converting face-based data

  else if (qCellVert) {
    if (iNCells == 0)
      vFatalError ("Mesh size unknown.", "2D 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.",
		     "2D 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);
    }

    if (qBFaceVert) {
      // Set up data for all the boundary edges
      for (int iBF = 0; iBF < iNBFaces; iBF++) {
	Vert *apV[2];
	apV[0] = pVVert (a2iBFaceVert[iBF][1]);
	apV[1] = pVVert (a2iBFaceVert[iBF][0]);
	EdgeBFaceCV *pEBFCV = VABdryCV.pTNewEntry ();
	pEBFCV->vSetVerts (2, apV);
	if (qBFaceBC)
	  pEBFCV->vSetBC (aiBFaceBC[iBF]);
      }
    }
    vMessage (1, "Converting triangle list to edge-cell format.\n");
    //@@@@ Convert the triangles to edge-cell data
    vConvertFromCellVert (VACV, VAQuadCV, VABdryCV);
    if (pB2DIn) {
      // Copy boundary data
      B2D = (*pB2DIn);
      vMessage(2, "Number of boundary patches: %d\n", B2D.iNumPatches());

      // Now infer which bdry faces are on which patches.  Start by
      // identifying all bdry verts, and finding which patches, if any,
      // they fall on.  While it's not very efficient to do this by
      // checking all bdry verts against all patches, this approach has
      // the virtue of working, and for meshes that aren't absolutely
      // huge, it should be no problem from the runtime point of view.
      multimap<Vert*, BdryPatch2D*> mmVertPatch;
      for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
	Vert *pV = pVVert(iV);
// 	vMessage(2, "Vert at (%10.6f, %10.6f) lies on patch %10p\n",
// 		 pV->dX(), pV->dY(), pBP2D);
	for (int iP = B2D.iNumPatches() - 1; iP >= 0; iP--) {
	  BdryPatch2D *pBP2D = B2D[iP];
	  if (!pBP2D->qDisjointFrom(pV)) {
	    vMessage(2, "Vert at (%10.6f, %10.6f) lies on patch %10p\n",
		     pV->dX(), pV->dY(), pBP2D);
	    mmVertPatch.insert(pair<Vert*, BdryPatch2D*>(pV, pBP2D));
	  }
	}
      }

      // Definitively identify bdry apexes: any vert on multiple patches
      // is an apex.
      for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
	Vert *pV = pVVert(iV);
	if (mmVertPatch.count(pV) > 1) {
	  pV->vSetType(Vert::eBdryApex);
	}
      }

      // Now for every bdry face, identify a patch common to its verts.
      for (int iBF = iNumBdryFaces() - 1; iBF >= 0; iBF--) {
	BFace *pBF = pBFBFace(iBF);
	Vert *pV0 = pBF->pVVert(0);
	Vert *pV1 = pBF->pVVert(1);
	// How many patches for each one?

	pair<multimap<Vert*, BdryPatch2D*>::iterator,
	  multimap<Vert*, BdryPatch2D*>::iterator> pRange0, pRange1;
	pRange0 = mmVertPatch.equal_range(pV0);
	pRange1 = mmVertPatch.equal_range(pV1);
	multimap<Vert*, BdryPatch2D*>::iterator iter0, iter1;
	iter0 = pRange0.first;
	BdryPatch2D* pBP2D = NULL;
	for (iter0 = pRange0.first; 
	     pBP2D == NULL && iter0 != pRange0.second;
	     iter0++) {
	  for (iter1 = pRange1.first;
	       pBP2D == NULL && iter1 != pRange1.second;
	       iter1++) {
	    if (iter0->second == iter1->second) {
	      pBP2D = iter0->second;
	    }
	  }
	}
	assert(pBP2D != NULL);
	vMessage(2, "BFace %3d lies on patch %10p\n", iBF, pBP2D);
	vMessage(2, "  coords: (%10.6f, %10.6f) to (%10.6f, %10.6f)\n",
		 pV0->dX(), pV0->dY(), pV1->dX(), pV1->dY());
	pBF->vSetPatch(pBP2D);
      }

      // Finally, reset all the periodic data, if any.
      vResetPeriodicData();

      // Check the periodicity
#ifndef NDEBUG
      for (int iBF = 0; iBF < iNumBdryFaces(); iBF++) {
	BFace* pBF = pBFBFace(iBF);
	BdryPatch2D* pBP2D = dynamic_cast<BdryPatch2D*>(pBF->pPatchPointer());
	if (pBP2D->qIsPeriodic()) {
	  Periodic* pP = pBP2D->pPGetPeriodicData();
	  const BdryPatch2D* pBP2DOther = dynamic_cast<const BdryPatch2D*>
	    (pP->pBPOtherPatch(pBP2D));
	  vMessage(2, "BFace %d is on patch %p opposite patch %p\n",
		   iBF, pBP2D, pBP2DOther);
	}
	else {
	  vMessage(2, "BFace %d is on patch %p\n", iBF, pBP2D);
	}
      }
      vMessage(2, "%d bfaces, %d patches, %d periodic patches\n", 
	       iNumBdryFaces(), B2D.iNumPatches(), B2D.iNumPeriodic());
#endif
      qLengthScaleFromCellSizes = false;
    }
    else {
      // Recreate boundary data by inference.

      // 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 iPatch = 0, iBPt = 0;
      for (int ii = 0; ii < iNVerts; aiVertMap[ii++] = -1) {}

      B2D.vSetSize(iNBFaces + iNIntBFaces, iNVerts);

      // Now the internal bdry faces, if any.
      if (qIntBFaceFace) {
	for (int iIBF = 0; iIBF < iNIntBFaces; iIBF++) {
	  // This array has to have live data, or we wouldn't be here.
 	  Face *pF0 = pFFace(a2iIntBFaceFace[iIBF][0]);
 	  Face *pF1 = pFFace(a2iIntBFaceFace[iIBF][1]);

	  IntBdryEdge *pIBE = pIBEIntBdryEdge(iIBF);

	  // Need this both for checking and for setting up patches.
	  Vert *pV0 = pIBE->pVVert(0);
	  Vert *pV1 = pIBE->pVVert(1);

	  // Need to create a patch for this bdry edge
	  int iV0 = iVertIndex(pV0);
	  int iV1 = iVertIndex(pV1);

	  if (aiVertMap[iV0] == -1) {
	    aiVertMap[iV0] = iBPt;
	    B2D.vAddPoint(pVVert(iV0)->adCoords());
	    B2D.vSetBdryPoint(iBPt, false);
	    iBPt++;
	  }

	  if (aiVertMap[iV1] == -1) {
	    aiVertMap[iV1] = iBPt;
	    B2D.vAddPoint(pVVert(iV1)->adCoords());
	    B2D.vSetBdryPoint(iBPt, false);
	    iBPt++;
	  }

	  int iBP0 = aiVertMap[iV0];
	  int iBP1 = aiVertMap[iV1];
	  assert(iBP0 >= 0 && iBP1 >= 0);

	  // Grab both region tags, working from IBF -> face -> cell.
	  Cell *pCL = pF0->pCCellOpposite(pIBE);
	  int iLeftReg = pCL->iRegion();
	  Cell *pCR = pF1->pCCellOpposite(pIBE);
	  int iRightReg = pCR->iRegion();
	  BdryPatch2D *pBP =
	    new BdrySeg(&B2D, iInvalidBC, iInvalidBC, iLeftReg, iRightReg,
			iBP0, iBP1);
	  B2D.vAddPatchToList(pBP);
	  // vSetPatch can handle getting the patch for an internal bdry
	  // fixed up, too.
	  pIBE->vSetPatch(iPatch, &B2D);
	  iPatch++;
	}
      } // Done with setting up internal bdry faces w/ face links

      // Now turn to regular bdry faces
      for (int iBF = 0; iBF < iNBFaces; iBF++) {
	Face *pF;
	Vert *pV0, *pV1;
	assert(qBFaceFace || qBFaceVert);
	if (qBFaceFace) {
	  // This array has to have live data, or we wouldn't be here.
	  pF = pFFace(aiBFaceFace[iBF]);
	  pV0 = pF->pVVert(0);
	  pV1 = pF->pVVert(1);

	}
	else {
	  assert(qBFaceVert && !qBFaceFace);

	  pV0 = pVVert(a2iBFaceVert[iBF][0]);
	  pV1 = pVVert(a2iBFaceVert[iBF][1]);
	  pF = pV0->pFHintFace();
	  if (!pF->qHasVert(pV1)) {
	    pF = pV1->pFHintFace();
	    if (!pF->qHasVert(pV0)) {
	      // Now we're in for it: need to walk around the vertex to
	      // find the actually face we're after.
	      Vert *pVTmp = pF->pVVert (0);
	      if (pVTmp == pV1)
		pVTmp = pF->pVVert (1);
	      Cell *pCTmp = pF->pCCellLeft ();
	      if (pCTmp->eType () == Cell::eBdryEdge)
		pCTmp = pF->pCCellRight ();
	      while (!pF->qHasVert (pV0)) {
		assert (pCTmp->eType () == Cell::eTriCell);
		Vert *pVNew = dynamic_cast<TriCell*>(pCTmp)->pVVertOpposite(pF);
		pF = dynamic_cast<TriCell*>(pCTmp)->pFFaceOpposite (pVTmp);
		pVTmp = pVNew;
		pCTmp = pF->pCCellOpposite (pCTmp);
	      }; // Done finding the right boundary face the hard way
	    }
	  }
	  assert(pF->qHasVert(pV0));
	  assert(pF->qHasVert(pV1));
	}
	pF->vSetFaceLoc(Face::eBdryFace);

	BFace *pBF = pBFBFace(iBF);

	// Need to create a patch for this bdry edge
	int iV0 = iVertIndex(pV0);
	int iV1 = iVertIndex(pV1);

	if (aiVertMap[iV0] == -1) {
	  aiVertMap[iV0] = iBPt;
	  B2D.vAddPoint(pVVert(iV0)->adCoords());
	  B2D.vSetBdryPoint(iBPt, false);
	  iBPt++;
	}

	if (aiVertMap[iV1] == -1) {
	  aiVertMap[iV1] = iBPt;
	  B2D.vAddPoint(pVVert(iV1)->adCoords());
	  B2D.vSetBdryPoint(iBPt, false);
	  iBPt++;
	}

	int iBP0 = aiVertMap[iV0];
	int iBP1 = aiVertMap[iV1];
	assert(iBP0 >= 0 && iBP1 >= 0);

	// Grab the BC from the incoming face data.
	int iF = iFaceIndex(pF);
	int iC, iBC;
	if (a2iFaceCell[iF][0] < 0) {
	  iBC = abs(a2iFaceCell[iF][0]);
	  iC  =     a2iFaceCell[iF][1];
	}
	else {
	  iBC = abs(a2iFaceCell[iF][1]);
	  iC  =     a2iFaceCell[iF][0];
	}
	assert(iBC > 0 && iC >= 0);

	// Find the adjacent region.
	Cell *pC = pCCell(iC);
	int iReg = pC->iRegion();

	BdryPatch2D* pBP;
	if (pF->pCCellLeft() == pC) {
	  // Cell is on the left, bdry on the right.
	  pBP = new BdrySeg(&B2D, iInvalidBC, iBC, iReg, iInvalidRegion,
			    iBP0, iBP1);
	  B2D.vAddPatchToList(pBP);
	  // vSetPatch can handle getting the patch for an internal bdry
	  // fixed up, too.
	  dynamic_cast<BdryEdge*>(pBF)->vSetPatch(iPatch, &B2D);
	  iPatch++;
	}
	else {
	  // Cell is on the right, bdry on the left.
	  pBP = new BdrySeg(&B2D, iBC, iInvalidBC, iInvalidRegion, iReg,
			    iBP0, iBP1);
	  B2D.vAddPatchToList(pBP);
	  // vSetPatch can handle getting the patch for an internal bdry
	  // fixed up, too.
	  dynamic_cast<BdryEdge*>(pBF)->vSetPatch(iPatch, &B2D);
	  iPatch++;
	}
      } // Done setting up all BFace stuff, including building patches.
      delete [] aiVertMap;
      assert (iPatch == iNBFaces + iNIntBFaces);
      qLengthScaleFromCellSizes = true;
    }
    vSetAllHintFaces ();

    if (!(qValid ()))
      vFoundBug ("2D mesh input from user-format FE-style file");
    // For V0.2:
    // Consistency checking is much harder for data that centers
    // around cell-vertex connectivity.  For example, the only
    // meaningful information that can be gleaned from cell-face or
    // face-cell data is that the reciprocity relationships are the
    // same for the data file and the mesh.  This is a lot of work for
    // very little gain, so I'm skipping it for now.
  } // Done translating cell-vert/bface-vert/bface-BC data

  if (qFaceVert)
    delete [] a2iFaceVert;
  if (qCellVert)
    delete [] a2iCellVert;
  if (qFaceCell)
    delete [] a2iFaceCell;
  if (qCellFace)
    delete [] a2iCellFace;
  if (qBFaceFace)
    delete [] aiBFaceFace;
  if (qBFaceVert)
    delete [] a2iBFaceVert;
  if (qBFaceBC)
    delete [] aiBFaceBC;
  vMessage(2, "%d bfaces, %d patches, %d periodic patches\n",
	   iNumBdryFaces(), B2D.iNumPatches(), B2D.iNumPeriodic());
  SUMAA_LOG_EVENT_END (READ_MESH_2D);
  SUMAA_LOG_EVENT_END (INIT_TRI);
}

void Mesh2D::vIdentifyVertexTypesGeometrically() const
  // In addition to some topological deductions, this routine also
  // detects verts at which the boundary has a sharp corner.  This is
  // critical for keeping boundaries from degenerating with repeated
  // coarsening.
  //
  // Note that no effort at all is made to define an apex as a point
  // where two patches meet, because a mesh input from file will have
  // one patch per bdry edge.  This would make every bdry vert a
  // eBdryApex, and that's no good.  However, boundary condition changes
  // are detected, and verts where the BC changes are tagged as
  // eBdryApex.
{
  // Tentatively set the type of every vertex to eInterior.  Then modify
  // that to accomodate boundary verts (both internal and regular).
  for (int iV = 0; iV < iNumVerts(); iV++) {
    Vert *pV = pVVert(iV);
    if (pV->iVertType() != Vert::eBdryApex) 
      pV->vSetType(Vert::eInterior);
  }

  // Mark internal boundaries
  int iBF;
  for (iBF = 0; iBF < iNumBdryFaces(); iBF++) {
    BFace *pBF = pBFBFace(iBF);
    if (pBF->eType() == Cell::eIntBdryEdge) {
      Face *pF = pBF->pFFace();
      if (pF->pVVert(0)->iVertType() != Vert::eBdryApex) 
	pF->pVVert(0)->vSetType(Vert::eBdryTwoSide);
      if (pF->pVVert(1)->iVertType() != Vert::eBdryApex) 
	pF->pVVert(1)->vSetType(Vert::eBdryTwoSide);
    }
  }

  // Now mark regular boundaries.  In the process, tag intersections of
  // regular and internal boundaries as type eBdryApex.  Finally, make a
  // sorted list of bdry normals (sorted by vertex) for both regular and
  // internal boundaries.
  std::multiset<BdryNorm> msBN;
  for (iBF = 0; iBF < iNumBdryFaces(); iBF++) {
    BFace *pBF = pBFBFace(iBF);
    Face *pF = pBF->pFFace();
    Cell *pCL = pF->pCCellLeft();
    Cell *pCR = pF->pCCellRight();

    Vert *pV0 = pF->pVVert(0);
    Vert *pV1 = pF->pVVert(1);

    if (pCL->eType() == Cell::eBdryEdge ||
	pCR->eType() == Cell::eBdryEdge) {
      switch (pV0->iVertType()) {
      case Vert::eBdryApex:
	break;
      case Vert::eBdryTwoSide:
	pV0->vSetType(Vert::eBdryApex);
	break;
      default:
	pV0->vSetType(Vert::eBdry);
	break;
      }

      switch (pV1->iVertType()) {
      case Vert::eBdryApex:
	break;
      case Vert::eBdryTwoSide:
	pV1->vSetType(Vert::eBdryApex);
	break;
      default:
	pV1->vSetType(Vert::eBdry);
	break;
      }
    }

    double adNorm[3];
    pBF->vUnitNormal(adNorm);
    msBN.insert(BdryNorm(pV0, adNorm));
    msBN.insert(BdryNorm(pV1, adNorm));
  }
  assert(int(msBN.size()) == iNumBdryFaces() * 2);

  // At this point, msBN is a sorted set of BdryNorm's, with multiple
  // entries having the same vertex.
  int iRecordsChecked = 0;
  std::multiset<BdryNorm>::iterator itFirst, itLast;
  itFirst = itLast = msBN.begin();
  do {
    Vert *pV = itFirst->pVVert();
    // Should already be tagged as a bdry vert.
    assert(pV->iVertType() != Vert::eInterior);
    do {
      itLast++;
    } while (itLast != msBN.end() && itLast->pVVert() == pV);
    itLast --;
    assert(itLast->pVVert() == pV);

    // If there are three or more distinct BFaces incident on a vertex,
    // it's definitely a eBdryApex, and should already be tagged as
    // such.
    assert(itFirst != itLast);
    int iDist = 1;
    {
      std::multiset<BdryNorm>::iterator itTemp = itFirst;
      itTemp++;
      while (itTemp != itLast) {
	itTemp++; iDist++;
      }
    }

    if (iDist >= 2) {
      assert(pV->iVertType() == Vert::eBdryApex);
      pV->vSetType(Vert::eBdryApex);
    }
    else {
      // There'd better be precisely two bdry edges:
      assert(iDist == 1);

      // Now check how close their normals are to being parallel.
      // Deviations of as much as 20 degrees are permitted.
      double dDot = dDOT2D(itFirst->adNormal(), itLast->adNormal());
      if (dDot < 0.94) {
	// This one is now an apex!
	pV->vSetType(Vert::eBdryApex);
      }
    }
    iRecordsChecked += iDist + 1;
    itFirst = ++itLast;
  } while (iRecordsChecked < 2*iNumBdryFaces());
}


int Mesh2D::iQueueEncroachedBdryEntities(InsertionQueue& IQ,
					 const std::set<BFace*>& spBF) const
{
  // Now queue up all encroached bdry entities
  // Can be used both to add all initially encroached BFaces and
  // subsegments to the insertion queue and to subsequently add
  // newly-created entities that are still encroached.

  // If the list is empty, then we're in the interior and there's
  // nothing to do.
  if (spBF.empty()) return 0;

  // If bdry's can't be changed, there's no point in queueing bdry
  // entities.
  if (!qBdryChangesAllowed()) return 0;

  // First, check all the bdry faces...
  int iAdded = 0;
  vMessage(4, "%d bfaces, %d patches, %d periodic patches\n",
	   iNumBdryFaces(), B2D.iNumPatches(), B2D.iNumPeriodic());
  std::set<BFace*>::iterator iterBF;
  for (iterBF = spBF.begin(); iterBF != spBF.end(); iterBF++) {
    BdryEdgeBase* pBE = (dynamic_cast<BdryEdgeBase*>(*iterBF));
    if (pBE->pVVert(0)->qIsStructured() &&
	pBE->pVVert(1)->qIsStructured()) continue;
    if (pBE->eIsEncroached(eEncroachmentType()) != eClean) {
      if (IQ.qAddEntry(InsertionQueueEntry(pBE))) {
	iAdded ++;
      }
      vMessage(4, "Queueing bdry edge (%10.6f, %10.6f) (%10.6f, %10.6f)\n",
	       pBE->pVVert(0)->dX(), pBE->pVVert(0)->dY(),
	       pBE->pVVert(1)->dX(), pBE->pVVert(1)->dY()); 
      // Now add the periodic face if one exists.
      BdryPatch* pBP = pBE->pPatchPointer();
      if (pBP->qIsPeriodic()) {
	vMessage(4, "Is periodic!\n");
	Periodic *pP = pBP->pPGetPeriodicData();
	Face *pFOpp = const_cast<Face*>(pP->pFFaceOpposite(pBP, pBE->pFFace()));
	if (!pFOpp->qValid()) {
	  // There isn't a periodic face for this one yet; it'll be
	  // added in a subsequent insertion.
	  continue;
	}
	BdryEdgeBase* pBEOpp =
	  dynamic_cast<BdryEdgeBase*>(pFOpp->pCCellLeft());
	if (pBEOpp == NULL) {
	  pBEOpp = dynamic_cast<BdryEdgeBase*>(pFOpp->pCCellRight());
	}
	vMessage(4, "Queueing periodic pair (%10.6f, %10.6f) (%10.6f, %10.6f)\n",
		 pFOpp->pVVert(0)->dX(), pFOpp->pVVert(0)->dY(),
		 pFOpp->pVVert(1)->dX(), pFOpp->pVVert(1)->dY()); 
	assert(pBEOpp != NULL);
	assert(dynamic_cast<const BdryEdge*>(pBEOpp) != NULL);
	if (IQ.qAddEntry(InsertionQueueEntry(pBEOpp))) {
	  iAdded ++;
	}
      } // Done with the periodic case.
    } // Done with queueing an encroached BFace
  } // Done with loop over all candidate BFace's
  vMessage(3, "Added a total of %d encroached bdry edges to the queue.\n",
	   iAdded);
  return iAdded;
}

int Mesh2D::iQueueEncroachedBdryEntities(InsertionQueue& IQ,
					 const std::set<BFace*>& spBF,
					 const double adLoc[],
					 const double dOffset) const
{
  // Now queue up all encroached bdry entities
  // Can be used both to add all initially encroached BFaces and
  // subsegments to the insertion queue and to subsequently add
  // newly-created entities that are still encroached.

  // If the list is empty, then we're in the interior and there's
  // nothing to do.
  if (spBF.empty()) return 0;

  // If bdry's can't be changed, there's no point in queueing bdry
  // entities.
  if (!qBdryChangesAllowed()) return 0;

  // First, check all the bdry faces...
  int iAdded = 0;
  vMessage(4, "%d bfaces, %d patches, %d periodic patches\n",
	   iNumBdryFaces(), B2D.iNumPatches(), B2D.iNumPeriodic());
  std::set<BFace*>::iterator iterBF;
  for (iterBF = spBF.begin(); iterBF != spBF.end(); iterBF++) {
    BdryEdgeBase* pBE = (dynamic_cast<BdryEdgeBase*>(*iterBF));
    if (pBE->eIsPointEncroaching(adLoc, eEncType, true) != eClean) {
      if (IQ.qAddEntry(InsertionQueueEntry(pBE), dOffset)) {
	iAdded ++;
      }
      vMessage(4, "Queueing bdry edge (%10.6f, %10.6f) (%10.6f, %10.6f)\n",
	       pBE->pVVert(0)->dX(), pBE->pVVert(0)->dY(),
	       pBE->pVVert(1)->dX(), pBE->pVVert(1)->dY()); 
      // Now add the periodic face if one exists.
      BdryPatch* pBP = pBE->pPatchPointer();
      if (pBP->qIsPeriodic()) {
	vMessage(4, "Is periodic\n");
	Periodic *pP = pBP->pPGetPeriodicData();
	Face *pFOpp = const_cast<Face*>(pP->pFFaceOpposite(pBP, pBE->pFFace()));
	if (!pFOpp->qValid()) {
	  // There isn't a periodic face for this one yet; it'll be
	  // added in a subsequent insertion.
	  continue;
	}
	BdryEdgeBase* pBEOpp =
	  dynamic_cast<BdryEdgeBase*>(pFOpp->pCCellLeft());
	if (pBEOpp == NULL) {
	  pBEOpp = dynamic_cast<BdryEdgeBase*>(pFOpp->pCCellRight());
	}
	vMessage(4, "Queueing periodic pair (%10.6f, %10.6f) (%10.6f, %10.6f)\n",
		 pFOpp->pVVert(0)->dX(), pFOpp->pVVert(0)->dY(),
		 pFOpp->pVVert(1)->dX(), pFOpp->pVVert(1)->dY()); 
	assert(pBEOpp != NULL);
	assert(dynamic_cast<const BdryEdge*>(pBEOpp) != NULL);
	if (IQ.qAddEntry(InsertionQueueEntry(pBEOpp))) {
	  iAdded ++;
	}
      } // Done with the periodic case.
    }
  }
  vMessage(3, "Added a total of %d encroached bdry edges to the queue.\n",
	   iAdded);
  return iAdded;
}

void Mesh2D::vBuildBdryPatches()
  // Given a Mesh2D (but no knowledge of the underlying geometry), make
  // a new set of patches, one per bdry face.  This is used, in
  // particular, to re-create bdry patch info after mesh coarsening.
{
  int iBF;
  int iNBFaces = iNumTotalBdryFaces();
  int iNVerts = iNumVerts();

  // Record the boundary condition and region info.
  int *aiRightReg = new int[iNBFaces];
  int *aiLeftReg  = new int[iNBFaces];
  int *aiRightBC  = new int[iNBFaces];
  int *aiLeftBC   = new int[iNBFaces];

  for (iBF = 0; iBF < iNBFaces; iBF++) {
    BFace *pBF = pBFBFace(iBF);
    Face *pF = pBF->pFFace(0);
    Cell *pC = pF->pCCellOpposite(pBF);
    assert(pC->eType() == Cell::eTriCell);

    Vert *pV0 = pBF->pVVert(0);
    Vert *pV1 = pBF->pVVert(1);
    Vert *pVOpp = dynamic_cast<TriCell*>(pC)->pVVertOpposite(pF);

    int iOrient = iOrient2D(pV0, pV1, pVOpp);
    bool qIntBFace = (pBF->iNumFaces() == 2);

    if (qIntBFace) {
      // Internal bdry face!

      // FIX ME!  Untested code, because coarsening with internal bdrys
      // hasn't been written yet.

      assert(pBF->eType() == Cell::eIntBdryEdge);
      // Identify the other cell
      Face *pFAlt = pBF->pFFace(1);
      Cell *pCAlt = pFAlt->pCCellOpposite(pBF);

      // Get regions on both sides.
      switch (iOrient) {
      case 1:
	aiLeftReg[iBF]  = pC   ->iRegion();
	aiRightReg[iBF] = pCAlt->iRegion();
	break;
      case -1:
	aiLeftReg[iBF]  = pCAlt->iRegion();
	aiRightReg[iBF] = pC   ->iRegion();
	break;
      default:
	// Unreachable for valid meshes
	assert(0);
	break;
      }
      aiRightBC[iBF] = aiLeftBC[iBF] = iInvalidBC;
    }
    else { // Regular bdry face
      aiRightReg[iBF] = aiLeftReg[iBF] = iInvalidRegion;
      aiRightBC[iBF]  = aiLeftBC[iBF]  = iInvalidBC;

      // Get regions on both sides.
      switch (iOrient) {
      case 1:
	aiLeftReg[iBF]  = pC ->iRegion();
	aiRightBC[iBF]  = pBF->iBdryCond();
	break;
      case -1:
	aiRightReg[iBF] = pC ->iRegion();
	aiLeftBC[iBF]   = pBF->iBdryCond();
	break;
      default:
	// Unreachable for valid meshes
	assert(0);
	break;
      }
    } // Done with real bdry faces
    // Do some post-checking to be sure this matches qualitatively with
    // the data for the patch.
#ifndef NDEBUG
    BdryPatch *pBP = pBF->pPatchPointer();
    assert( ( aiRightReg[iBF] == pBP->iRightRegion() ||
	      aiRightReg[iBF] == pBP->iLeftRegion() )
	    &&
	    ( aiLeftReg[iBF]  == pBP->iRightRegion() ||
	      aiLeftReg[iBF]  == pBP->iLeftRegion() ) );
    assert( ( aiRightBC[iBF] == pBP->iBdryCond() ||
	      aiRightBC[iBF] == iInvalidBC )
	    &&
	    ( aiLeftBC[iBF]  == pBP->iBdryCond() ||
	      aiLeftBC[iBF]  == iInvalidBC ) );
#endif
  } // End of loop over BFaces


  // Throw away all the existing bdry patch info.
  B2D.vSetSize(iNBFaces, iNVerts);
  B2D.vClearPoints();

  // Now begin again, creating a new bdry patch for every BdryEdge.

  // 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) {}

  int iBPt = 0;
  for (iBF = 0; iBF < iNBFaces; iBF++) {
    BFace *pBF = pBFBFace(iBF);
    int iV0 = iVertIndex(pBF->pVVert(0));
    int iV1 = iVertIndex(pBF->pVVert(1));

    if (aiVertMap[iV0] == -1) {
      aiVertMap[iV0] = iBPt;
      B2D.vAddPoint(pVVert(iV0)->adCoords());
      B2D.vSetBdryPoint(iBPt, false);
      iBPt++;
    }

    if (aiVertMap[iV1] == -1) {
      aiVertMap[iV1] = iBPt;
      B2D.vAddPoint(pVVert(iV1)->adCoords());
      B2D.vSetBdryPoint(iBPt, false);
      iBPt++;
    }

    int iBP0 = aiVertMap[iV0];
    int iBP1 = aiVertMap[iV1];
    assert(iBP0 >= 0 && iBP1 >= 0);

    BdryPatch2D *pBP =
      new BdrySeg(&B2D, aiLeftBC[iBF], aiRightBC[iBF], aiLeftReg[iBF],
		  aiRightReg[iBF], iBP0, iBP1);
    B2D.vAddPatchToList(pBP);
    // vSetPatch can handle getting the patch for an internal bdry
    // fixed up, too.
    (dynamic_cast<BdryEdgeBase*>(pBF))->vSetPatch(iBF, &B2D);
  } // Done with boundary cases
  delete [] aiVertMap;

  if (!(qValid ()))
    vFoundBug ("building 2D boundary data from mesh");
}

// Find the cell incident on pV0 that lies astride the direction to pV1.
Cell* Mesh2D::pCBeginPipe(const Vert* const pV0, const Vert* const pV1) const
{
  std::set<Cell*> spCInc;
  std::set<Vert*> spVTmp;
  vNeighborhood(pV0, spCInc, spVTmp);

  // Now check out each cell.
  std::set<Cell*>::iterator iterC;
  for (iterC = spCInc.begin(); iterC != spCInc.end(); iterC++) {
    assert((*iterC)->eType() == Cell::eTriCell);
    TriCell *pTCCand = dynamic_cast<TriCell*>(*iterC);
    assert(pTCCand->qHasVert(pV0));

    Face *pFOpp = pTCCand->pFFaceOpposite(pV0);
    Vert *pVA = pFOpp->pVVert(0);
    Vert *pVB = pFOpp->pVVert(1);

    int iOrientA = iOrient2D(pVA, pV0, pV1);
    int iOrientB = iOrient2D(pVB, pV0, pV1);

    if (iOrientA != iOrientB) {
      // pV1 falls between the lines determined by (pVA, pV0) and (pVB,
      // pV0), or on one of them.  Now, to eliminate a cell that happens
      // to point away from pV1, look at barycentrics.
      double adBary[3];
      pTCCand->vBarycentrics(pV1->adCoords(), adBary);
      int iNegBary = ((adBary[0] < -1.e-10 ? 1 : 0) +
		      (adBary[1] < -1.e-10 ? 1 : 0) +
		      (adBary[2] < -1.e-10 ? 1 : 0));
      int iZeroBary = ((fabs(adBary[0]) < 1.e-10 ? 1 : 0) +
		       (fabs(adBary[1]) < 1.e-10 ? 1 : 0) +
		       (fabs(adBary[2]) < 1.e-10 ? 1 : 0));
      // Combinations of result values:
      // iNeg iZero Case
      //   0    0   pV1 inside pTCCand (Hit!)
      //   0    1   pV1 on pTCCand (Hit!)
      //   0    2   pV1 is a vertex of pTCCand (Hit!)
      //   1    0   pV1 is in the "wedge" of pTCCand but beyond its face
      //              opposite pV0 (Hit! You sunk my battleship!)
      //   1    1   Can't tell yet whether pV1 is on the "proper" side
      //            of pV0 for pTCCand to be the right choice or not. (?)
      //   2    0   pV1 is in the "reverse wedge" of pTCCand (i.e., the
      //            ray from pV0 to pV1 is headed in the opposite
      //            direction to pTCCand) (Miss!)
      if (iNegBary == 0 || (iNegBary == 1 && iZeroBary == 0)) {
	return pTCCand;
      }
      else if (iNegBary == 1 && iZeroBary == 1) {
	// Need to find out whether the edge of pTCCand that is
	// collinear with pV1 is pointing in the same or the opposite
	// direction as pV1.
	assert(iOrientA == 0 || iOrientB == 0);
	Vert *pVTest = (iOrientA == 0) ? pVA : pVB;
	double adVecTest[] = adDIFF2D(pVTest->adCoords(),
				      pV0->adCoords());
	double adVec1[]    = adDIFF2D(pV1->adCoords(),
				      pV0->adCoords());
	double dDot = dDOT2D(adVecTest, adVec1);
	if (dDot > 0)
	  return pTCCand;
	else {
	  // Wrong choice; keep checking
	}
      }
      else {
	// Do nothing; this cell isn't the right one, so keep checking.
      }
    }
  }
  assert(0);
  return pCInvalidCell;
}

void Mesh2D::vWriteBdryPatchFile(const char strBaseFileName[],
				const char strExtraFileSuffix[])
{
  SUMAA_LOG_EVENT_BEGIN(OUTPUT);

  assert(qSimplicial());
  FILE *pFOutFile = NULL;
  char strFileName[1024];
  vPurge();
  sprintf(strFileName, "%s.bpatch%s", strBaseFileName, strExtraFileSuffix);

  /*newfile mesh*/
  vMessage(1,"Opening boundary patch output file %s\n", strFileName);
  if (pFOutFile) fclose(pFOutFile);
  pFOutFile = fopen(strFileName, "w");
  if (NULL == pFOutFile)
    vFatalError("Couldn't open boundary patch output file for writing",
                "2d mesh output");

  int iNumGaussPoints = 18;
  double adS[iNumGaussPoints];
  double adBeginRatio[iNumGaussPoints];
  double adEndRatio[iNumGaussPoints];
  
  double dS = 0.5 * (1 + sqrt(1./3.));

  adS[0] = 0.5; // one gauss point
  adBeginRatio[0] = 0;
  adEndRatio[0] = 1.0;

  adS[1] = 1 - dS; // two gauss points
  adBeginRatio[1] = 0;
  adEndRatio[1] = 0.5;
  adS[2] = dS;
  adBeginRatio[2] = 0.5;
  adEndRatio[2] = 1.0;

  adS[3] = 0.5 - sqrt(0.15); // three gauss points
  adBeginRatio[3] = 0;
  adEndRatio[3] = 5./18.;
  adS[4] = 0.5;
  adBeginRatio[4] = 5./18.;
  adEndRatio[4] = 13./18.;
  adS[5] = 0.5 + sqrt(0.15);
  adBeginRatio[5] = 13./18.;
  adEndRatio[5] = 1.;

  adS[6] = 0.25;  // two times one gauss points
  adBeginRatio[6] = 0;
  adEndRatio[6] = 0.5;
  adS[7] = 0.75;
  adBeginRatio[7] = 0.5;
  adEndRatio[7] = 1.0;

  adS[8] = 0.5 - 0.5*dS; // two times two gauss points
  adBeginRatio[8] = 0;
  adEndRatio[8] = 0.25;
  adS[9] = 0.5*dS;
  adBeginRatio[9] = 0.25;
  adEndRatio[9] = 0.5;
  adS[10] = 1 - 0.5*dS;
  adBeginRatio[10] = 0.5;
  adEndRatio[10] = 0.75;
  adS[11] = 0.5 + 0.5*dS;
  adBeginRatio[11] = 0.75;
  adEndRatio[11] = 1.;

  adS[12] = 0.5*(0.5 - sqrt(0.15));
  adBeginRatio[12] =  0;
  adEndRatio[12] =  5./36.;
  adS[13] =  0.25;
  adBeginRatio[13] =  5./36.;
  adEndRatio[13] =  13./36.;
  adS[14] =  0.5*(0.5 + sqrt(0.15));
  adBeginRatio[14] =  13./36.;
  adEndRatio[14] =  0.5;
  adS[15] = 0.5 + 0.5*(0.5 - sqrt(0.15)); // two times three gauss points
  adBeginRatio[15] = 0.5 + 0;
  adEndRatio[15] = 0.5 + 5./36.;
  adS[16] = 0.5 + 0.25;
  adBeginRatio[16] = 0.5 +5./36.;
  adEndRatio[16] = 0.5 +13./36.;
  adS[17] = 0.5 + 0.5*(0.5 + sqrt(0.15));
  adBeginRatio[17] = 0.5 + 13./36.;
  adEndRatio[17] = 1.;


  for (int iF = 0; iF < iNumFaces(); iF++) {
    Face* pF = pFFace(iF);
    
    if (pF->qIsBdryFace()) {
      Cell *pC = pF->pCCellLeft();
      if ((pC->eType() != Cell::eBdryEdge) && (pC->eType() != Cell::eTriBFace)
	  && (pC->eType() != Cell::eQuadBFace) && (pC->eType() != Cell::eIntBdryEdge)
	  && (pC->eType() != Cell::eIntTriBFace) && (pC->eType() != Cell::eIntQuadBFace)) {
	pC = pF->pCCellRight();
      }

      assert((pC->eType() == Cell::eBdryEdge) || (pC->eType() == Cell::eTriBFace)
	     || (pC->eType() == Cell::eQuadBFace) || (pC->eType() == Cell::eIntBdryEdge)
	     || (pC->eType() == Cell::eIntTriBFace) || (pC->eType() == Cell::eIntQuadBFace));
       

      BdryPatch2D *pBP = dynamic_cast<BdryPatch2D*>((dynamic_cast<BFace*>(pC))->pPatchPointer());
      
      double adGaussLoc[2];
      double adBeginLoc[2];
      double adEndLoc[2];
      double dLength;
      double adNormal[2];
      const double *adVertLocA;
      const double *adVertLocB;
      
      adVertLocA = pC->pVVert(0)->adCoords();
      adVertLocB = pC->pVVert(1)->adCoords();

      fprintf(pFOutFile, "%d", iF);
      
      for (int iGP = 0; iGP < iNumGaussPoints; iGP++) {

 	pBP->vGeodesicDistanceRatioBetweenTwoPoints(adBeginRatio[iGP], adVertLocA, adVertLocB, adBeginLoc);
 	pBP->vGeodesicDistanceRatioBetweenTwoPoints(adEndRatio[iGP], adVertLocA, adVertLocB, adEndLoc);

	dLength = fabs(pBP->dGeodesicDistance(pBP->dParameterAtPoint(adBeginLoc),pBP->dParameterAtPoint(adEndLoc)));

 	pBP->vGeodesicDistanceRatioBetweenTwoPoints(adS[iGP], adVertLocA, adVertLocB, adGaussLoc);

	pBP->vUnitNormal(adGaussLoc, adNormal);
	fprintf(pFOutFile, " %.16g %.16g %.16g %.16g %.16g", adGaussLoc[0], adGaussLoc[1], adNormal[0], adNormal[1], dLength);
      }
      fprintf(pFOutFile, "\n");
    }
  }

  fclose(pFOutFile);
  SUMAA_LOG_EVENT_END(OUTPUT);
}

bool Mesh2D::qCreateQuadFromTris(Cell* const pCA,
				 Cell* const pCB)
{
  assert(pCA->iNumVerts() == 3);
  assert(pCB->iNumVerts() == 3);
  assert(dynamic_cast<const TriCell*>(pCA) != NULL);
  assert(dynamic_cast<const TriCell*>(pCB) != NULL);

  Face *pFCommon = pFCommonFace(pCA, pCB);
  if (!pFCommon->qValid()) return false;

  Cell* pCQuad = pCNewCell(4);

  Face *apFOuter[4];
  int iF, iFOuter = 0;
  for (iF = 0; iF < 3; iF++) {
    Face *pF = pCA->pFFace(iF);
    if (pF != pFCommon) {
      apFOuter[iFOuter++] = pF;
      pF->vReplaceCell(pCA, pCQuad);
    }

    pF = pCB->pFFace(iF);
    if (pF != pFCommon) {
      apFOuter[iFOuter++] = pF;
      pF->vReplaceCell(pCB, pCQuad);
    }
  }
  
  Vert *pVA = pFCommon->pVVert(0);
  Vert *pVB = pFCommon->pVVert(1);

  // Update vertex hint for pVA, if necessary.
  if (pVA->pFHintFace() == pFCommon) {
    bool qGotIt = false;
    iF = 0;
    while (!qGotIt) {
      Face *pF = apFOuter[iF++];
      if (pF->qHasVert(pVA)) {
	qGotIt = true;
	pVA->vSetHintFace(pF);
      }
    }
  }

  // Update vertex hint for pVB, if necessary.
  if (pVB->pFHintFace() == pFCommon) {
    bool qGotIt = false;
    iF = 0;
    while (!qGotIt) {
      Face *pF = apFOuter[iF++];
      if (pF->qHasVert(pVB)) {
	qGotIt = true;
	pVB->vSetHintFace(pF);
      }
    }
  }

  pCA->vMarkDeleted();
  pCB->vMarkDeleted();
  pFCommon->vMarkDeleted();
  pCQuad->vAssign(apFOuter);

  assert(pCQuad->qValid());
  assert(apFOuter[0]->qValid());
  assert(apFOuter[1]->qValid());
  assert(apFOuter[2]->qValid());
  assert(apFOuter[3]->qValid());

  assert(pVA->pFHintFace() != pFCommon && pVA->pFHintFace()->qHasVert(pVA));
  assert(pVB->pFHintFace() != pFCommon && pVB->pFHintFace()->qHasVert(pVB));
  return true;
}


bool Mesh2D::qCreateQuadFromVerts(Vert* const pVA, Vert* const pVB,
				  Vert* const pVC, Vert* const pVD)
{
  // Vertices must be specified in cyclic order.  These checks  verify that.
  assert(iOrient2D(pVA, pVB, pVC) == 1);
  assert(iOrient2D(pVB, pVC, pVD) == 1);
  assert(iOrient2D(pVC, pVD, pVA) == 1);
  assert(iOrient2D(pVD, pVA, pVB) == 1);

  // Step one: ensure that all the perimeter edges exist, or swap for them.
  // Because the edges are known to be non-intersecting (based on the
  // previous assertions), then recovering all four of them should be
  // trivial.  Collect face ID's as you go.
  Face *pFAB, *pFBC, *pFCD, *pFDA;
  bool qRes = qRecoverEdge(pVA, pVB, pFAB) &&
    qRecoverEdge(pVB, pVC, pFBC) &&
    qRecoverEdge(pVC, pVD, pFCD) &&
    qRecoverEdge(pVD, pVA, pFDA);
  if (!(qRes && pFAB->qValid() && pFBC->qValid() &&
	pFCD->qValid() && pFDA->qValid() &&
	!pFAB->qDeleted() && !pFBC->qDeleted() &&
	!pFCD->qDeleted() && !pFDA->qDeleted() && 
	pFAB->qHasVert(pVA) && pFAB->qHasVert(pVB) &&
	pFBC->qHasVert(pVB) && pFBC->qHasVert(pVC) &&
	pFCD->qHasVert(pVC) && pFCD->qHasVert(pVD) &&
	pFDA->qHasVert(pVD) && pFDA->qHasVert(pVA))) {
    // Something bad, and unexpected, just happened!
    return false;
  }  

  // Step two: identify triangles and interior edge for quad creation.
  // Each edge must share a triangle with one of the adjacent edges; we
  // just need to find out which is which.
  Cell* pCQuad;

  Cell *pC1, *pC2;
  pC1 = pCCommonCell(pFAB, pFBC);
  Face *pFCommon = pFInvalidFace;
  if (pC1 == pCInvalidCell) {
    // So pFAB and pFDA have a common tri, as do pFBC and pFCD
    pC1 = pCCommonCell(pFAB, pFDA);
    pC2 = pCCommonCell(pFBC, pFCD);
    if ((pC1 == pCInvalidCell) || (pC2 == pCInvalidCell)) {
      // In this case, there's some vert inside the alleged quad; bail!
      return false;
    }
    pFCommon = pFCommonFace(pC1, pC2);
    assert(pFCommon != pFInvalidFace);
    pCQuad = pCNewCell(4);
    pFAB->vReplaceCell(pC1, pCQuad);
    pFBC->vReplaceCell(pC2, pCQuad);
    pFCD->vReplaceCell(pC2, pCQuad);
    pFDA->vReplaceCell(pC1, pCQuad);
  }
  else {
    // The other tri is shared by pFCD and pFDA
    pC2 = pCCommonCell(pFCD, pFDA);
    if (pC2 == pCInvalidCell) {
      // In this case, there's some vert inside the alleged quad; bail!
      return false;
    }
    pFCommon = pFCommonFace(pC1, pC2);
    assert(pFCommon != pFInvalidFace);
    pCQuad = pCNewCell(4);
    pFAB->vReplaceCell(pC1, pCQuad);
    pFBC->vReplaceCell(pC1, pCQuad);
    pFCD->vReplaceCell(pC2, pCQuad);
    pFDA->vReplaceCell(pC2, pCQuad);
  }


  Vert *pV0 = pFCommon->pVVert(0);
  Vert *pV1 = pFCommon->pVVert(1);

  Face *apFOuter[] = {pFAB, pFBC, pFCD, pFDA};
  // Update vertex hint for pV0, if necessary.
  int iF;
  if (pV0->pFHintFace() == pFCommon) {
    bool qGotIt = false;
    iF = 0;
    while (!qGotIt) {
      Face *pF = apFOuter[iF++];
      if (pF->qHasVert(pV0)) {
	qGotIt = true;
	pV0->vSetHintFace(pF);
      }
    }
  }

  // Update vertex hint for pV1, if necessary.
  if (pV1->pFHintFace() == pFCommon) {
    bool qGotIt = false;
    iF = 0;
    while (!qGotIt) {
      Face *pF = apFOuter[iF++];
      if (pF->qHasVert(pV1)) {
	qGotIt = true;
	pV1->vSetHintFace(pF);
      }
    }
  }

  pC1->vMarkDeleted();
  pC2->vMarkDeleted();
  pFCommon->vMarkDeleted();
  pCQuad->vAssign(apFOuter);

  assert(pCQuad->qValid());
  assert(pFAB->qValid());
  assert(pFBC->qValid());
  assert(pFCD->qValid());
  assert(pFDA->qValid());

  assert(pV0->pFHintFace() != pFCommon && pV0->pFHintFace()->qHasVert(pV0));
  assert(pV1->pFHintFace() != pFCommon && pV1->pFHintFace()->qHasVert(pV1));
  return true;
}

Vert* Mesh2D::createVert(const double dX, const double dY, const double dZ)
{
  assert(iFuzzyComp(dZ, 0) == 0);
  Vert *pV = pVNewVert();
  pV->vSetDefaultFlags();
  double adData[] = {dX, dY};
  pV->vSetCoords(2, adData);
  return pV;
}

Face* Mesh2D::createFace(Vert * const pV0, Vert * const pV1)
{
  assert(pV0->qValid());
  assert(pV1->qValid());
  // Check for prior existence.
  Face *pF = findCommonFace(pV0, pV1);
  if (pF == pFInvalidFace) {
    // Face doesn't exist yet.
    pF = pFNewFace();
    pF->vAssign(pCInvalidCell, pCInvalidCell, pV0, pV1);

#ifndef OMIT_VERTEX_HINTS
    // Make sure verts always have a valid hint face...
    if (!pV0->pFHintFace()->qValid())
      pV0->vSetHintFace(pF);
    if (!pV1->pFHintFace()->qValid())
      pV1->vSetHintFace(pF);
#endif
  }

  assert(pF->iFullCheck());

  return pF;
}

