#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "GR_events.h"
#include "GR_BFace.h"
#include "GR_Cell.h"
#include "GR_Face.h"
#include "GR_Geometry.h"
#include "GR_InsertionQueue.h"
#include "GR_Sort.h"
#include "GR_Vec.h"
#include "GR_Vertex.h"
#include "GR_VolMesh.h"
#include "GR_Watson.h"

#include <algorithm>
#include <map>

using std::multimap;
using std::pair;

extern void vReadUGridBinary(const char * const fileName, VolMesh * const pVM);
extern void vReadVGridBinary(const char * const fileName, VolMesh * const pVM);

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

//@ Construct volume mesh by reading data from a file
VolMesh::VolMesh(const char strBaseFileName[], const int iQualMeas,
		 const double dMaxAngle, Bdry3D * const pB3DIn) :
  Mesh(), qOKToChangeSurfaceEdges(true),
  qDoEdgeSwapping(true), qStrictPatchChecking(false),
  dMaxAngleForSurfSwap(dMaxAngle), VATriF(),
  VAQuadF(), VATriBF(), VAQuadBF(), VATet(), VAPyr(), VAPrism(),
  VAHex(), VAIntTriBF(), VAIntQuadBF(), pB3D(pB3DIn),
  qBdryFromMesh(pB3DIn == NULL)
{
  SUMAA_LOG_EVENT_BEGIN(READ_MESH_3D);
  qLengthScaleFromCellSizes = true;
  qBdryFill = false;
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(3, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);


  bool qIsUGridFile = strstr(strBaseFileName, ".ugrid");
  bool qIsVGridFile = strstr(strBaseFileName, ".cogsg");
//   bool qIsVMeshFile = strstr(strBaseFileName, ".vmesh");
//   bool qIsVTKFile = strstr(strBaseFileName, ".vtk");

  if (qIsUGridFile) {
    vReadUGridBinary(strBaseFileName, this);
    int aiBCList[iNumBdryFaces()];
    for (int i = 0; i < iNumBdryFaces(); i++) {
      aiBCList[i] = iDefaultBC;
    }
    vCreatePatchesFromBdryFaces(aiBCList);
  }
  else if (qIsVGridFile) {
    vReadVGridBinary(strBaseFileName, this);
    int aiBCList[iNumBdryFaces()];
    for (int i = 0; i < iNumBdryFaces(); i++) {
      aiBCList[i] = iDefaultBC;
    }
    vCreatePatchesFromBdryFaces(aiBCList);
  }
  else {
    // Now read all of the mesh data from file
    int iNVerts, iNFaces, iNCells, iNBFaces, iNIntBFaces;
    bool qFaceVert, qCellVert, qFaceCell, qCellFace, qCellRegion;
    bool qBFaceFace, qBFaceVert, qBFaceBC;
    bool qIntBFaceFace, qIntBFaceVert, qIntBFaceBC;
    int (*a2iFaceVert)[3], (*a2iFaceCell)[2];
    int (*a2iCellVert)[4], (*a2iCellFace)[4];
    int *aiCellRegion;
    int *aiBFaceFace, *aiBFaceBC, (*a2iBFaceVert)[3];
    int (*a2iIntBFaceFace)[2], *aiIntBFaceBC, (*a2iIntBFaceVert)[3];
    
    vReadFile_VolMesh(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", "volume 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.",
		  "volume mesh input");
    
    int *aiBCList;
    
    //@@@ Convert face-based data to VolMesh
    // This is the easiest conversion case
    
    if (qFaceVert && qFaceCell) {
      // If the number of cells, faces, and boundary faces have been
      // provided, we can check for mesh validity.
      if (iNFaces == 0)
	vFatalError("Mesh size unknown.", "volume mesh input");
      //@@@@ Verify mesh size
      {
	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];
	  int iV2 = a2iFaceVert[iF][2];
	  iTmpNVerts = max(max(iTmpNVerts, iV0), max(iV1, iV2));
	}
	if (iNCells == 0)
	  iNCells = iTmpNCells+1;
	else if (iNCells != iTmpNCells+1)
	  vFatalError("Number of cells in mesh doesn't match number given.",
		      "volume mesh input");
	
	if (iNBFaces == 0)
	  iNBFaces = iTmpNBFaces;
	else if (iNBFaces != iTmpNBFaces)
	  vFatalError("Number of bdry faces in mesh doesn't match number given.",
		      "volume mesh input");
	
	if (iNVerts == 0)
	  iNVerts = iTmpNVerts+1;
	else if (iNVerts != iTmpNVerts+1)
	  vFatalError("Number of verts in mesh doesn't match number given.",
		      "volume mesh input");
      }
      
      // FIX ME
      // To do Euler's formula, we would need to identify all the edges,
      // which isn't worth doing yet.
      
      vSetupTriFaces(iNFaces);       // Ensure the existence of enough
      vSetupTetCells(iNCells);       // blank faces, bdry faces, and cells.
      vSetupTriBFaces(iNBFaces);
      
      aiBCList = new int[iNBFaces];
      
      
      int iBFace = 0;
      for (int i = 0; i < iNFaces; i ++) {
	// Add entry to cell-face table or bdry face list on the fly
	int iC0 = a2iFaceCell[i][0];
	int iC1 = a2iFaceCell[i][1];
	Cell *pC0, *pC1;
	Face *pF = pFFace(i);
	if (iC0 < 0) {
	  BFace* pBF = pBFBFace(iBFace);
	  pBF->vAssign(pF);
	  pF->vSetFaceLoc(Face::eBdryFace);
	  aiBCList[iBFace] = abs(iC0);
	  pC0 = pBF;
	  iBFace ++;
	}
	else {
	  pC0 = pCCell(iC0);
	  pC0->vAddFace(pF);
	}
	
	if (iC1 < 0) {
	  BFace *pBF = pBFBFace(iBFace);
	  pBF->vAssign(pF);
	  pF->vSetFaceLoc(Face::eBdryFace);
	  aiBCList[iBFace] = abs(iC1);
	  pC1 = pBF;
	  iBFace ++;
	}
	else {
	  pC1 = pCCell(iC1);
	  pC1->vAddFace(pF);
	}
	
	// Tell the face everything it needs to know.
	int iV0 = a2iFaceVert[i][0];
	int iV1 = a2iFaceVert[i][1];
	int iV2 = a2iFaceVert[i][2];
	pFFace(i)->vAssign(pC0, pC1, pVVert(iV0), pVVert(iV1), pVVert(iV2));
      }
      assert(iBFace == iNBFaces);
      
      //@@@ Assign BC's and check boundary consistency
      // If the BC's were read on a boundary-face by boundary-face
      // basis, then these BC's silently override anything read
      // elsewhere. Also, verify that the boundary face and its
      // corresponding face have the same vertices, if data exists to do
      // this check.
      
      //@@@@ The easy case:  real face ID read
      if (qBFaceFace && (qBFaceBC || qBFaceVert)) {
	for (int iBF = 0; iBF < iNBFaces; iBF++) {
	  int iFace = aiBFaceFace[iBF];
	  Face *pF = pFFace(iFace);
	  //@@@@@ Assign BC
	  if (qBFaceBC) {
	    int iBC = aiBFaceBC[iBF];
	    Cell *pC = pF->pCCellLeft();
	    if (pC->eType() == Cell::eTet) {
	      pC = pF->pCCellRight();
	    }
	    assert(pC->eType() == Cell::eTriBFace);
	    aiBCList[iBF] = abs(iBC);
	  } // Done assigning BC
	  //@@@@@ Verify vertex consistency
	  if (qBFaceVert) {
	    int iV0 = iVertIndex(pF->pVVert(0));
	    int iV1 = iVertIndex(pF->pVVert(1));
	    int iV2 = iVertIndex(pF->pVVert(2));
	    if ( (iV0 != a2iBFaceVert[iBF][0] &&
		  iV0 != a2iBFaceVert[iBF][1] &&
		  iV0 != a2iBFaceVert[iBF][2])
		 ||
		 (iV1 != a2iBFaceVert[iBF][0] &&
		  iV1 != a2iBFaceVert[iBF][1] &&
		  iV1 != a2iBFaceVert[iBF][2])
		 ||
		 (iV2 != a2iBFaceVert[iBF][0] &&
		  iV2 != a2iBFaceVert[iBF][1] &&
		  iV2 != a2iBFaceVert[iBF][2]) ) {
	      vFatalError("Verts of bdry face don't match verts of actual face",
			  "volume mesh input");
	    }
	  } // Done consistency checking for vertices on boundary
	} // Done looping over all bdry faces
      } // Done with cases where face ID was read for each bdry face
      
      //@@@@ Assign BC's if bface-vert and bdry cond data have been read
      else if (qBFaceVert && qBFaceBC) {
	for (int iBF = 0; iBF < iNBFaces; iBF++) {
	  int iV0 = a2iBFaceVert[iBF][0];
	  int iV1 = a2iBFaceVert[iBF][1];
	  int iV2 = a2iBFaceVert[iBF][2];
	  Vert *pV0 = pVVert(iV0);
	  Vert *pV1 = pVVert(iV1);
	  Vert *pV2 = pVVert(iV2);
	  // The bface must already exist, but which one is it?
	  // Find all incident faces for pV0.
	  std::set<Cell*> spCIncident;
	  std::set<Vert*> spVNeigh;
	  std::set<BFace*> spBFIncident;
	  BFace *pBF = pBFInvalidBFace; /* Init for picky compilers. */
	  vNeighborhood(pV0, spCIncident, spVNeigh, &spBFIncident);
	  assert(spVNeigh.count(pV1) == 1);
	  assert(spVNeigh.count(pV2) == 1);
	  assert(spBFIncident.size() > 0);
	  for (std::set<BFace*>::iterator iter = spBFIncident.begin();
	       iter != spBFIncident.end(); iter++) {
	    pBF = *iter;
	    assert(pBF->qHasVert(pV0));
	    if (pBF->qHasVert(pV1) && pBF->qHasVert(pV2)) break;
	  }
	  assert(pBF->qHasVert(pV1) && pBF->qHasVert(pV2));
	  aiBCList[iBF] = abs(aiBFaceBC[iBF]);
	} // Done looping over all bdry faces to assign BC's
      } // Done assigning BC's based on bdry verts alone
      else 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 < pC->iNumFaces(); ii++) {
	    int iF = a2iCellFace[iC][ii];
	    Face *pF = pFFace(iF);
	    if (!(pC->qHasFace(pF)))
	      vFatalError("Cell-face data inconsistent with face-cell data",
			  "volume 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 < pC->iNumVerts(); ii++) {
	    int iV = a2iCellVert[iC][ii];
	    Vert *pV = pVVert(iV);
	    if (!(pC->qHasVert(pV)))
	      vFatalError("Cell-vert data inconsistent with face-vert data",
			  "volume mesh input");
	  }
	}
	vMessage(1, "OK\n");
      }
      
    } // Done converting face-based data
    else if (qCellVert) {
      // This is the cell-vert conversion case.
      if (iNCells == 0)
	vFatalError("Mesh size unknown.", "volume 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];
	  int iV3 = a2iCellVert[iC][3];
	  iTmpNVerts = max(iTmpNVerts, max(max(iV0, iV1), max(iV2, iV3)));
	}
	if (iNVerts != 0 && iNVerts != iTmpNVerts+1)
	  vFatalError("Number of verts in mesh doesn't match number given.",
		      "volume mesh input");
      }
      
      assert(iNBFaces >= 0);
      aiBCList = new int[iNBFaces];
      
      // FIX ME 0.2.3
      // All format translation needs to be modularized and pulled out for
      // easy access when converting from and to user meshes for
      // refinement.
      vConvertFromCellVert(iNCells, iNBFaces, a2iCellVert,
			   a2iBFaceVert, aiBFaceBC, aiBCList);
      
      if (!(qValid()))
	vFoundBug("volume mesh input from user-format FE-style file");
      // FIX ME 0.2.3:
      // 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
    
    // FIX ME: Add code to set up bdry data after 3D mesh input conversion.
    
    //@@@ Assign cells to regions
    // The default region, set when the cell was created, is 1.
    if (qCellRegion) {
      for (int iCell = 0; iCell < iNumCells(); iCell++)
	pCCell(iCell)->vSetRegion(aiCellRegion[iCell]);
      // Done with aiCellRegion, so get rid of it.
      free(aiCellRegion);
    }
    
    //@@@ Tag faces and verts
    for (int iVert = 0; iVert < iNumVerts(); iVert++)
      pVVert(iVert)->vSetType(Vert::eInterior);
    for (int iFace = 0; iFace < iNumFaces(); iFace++) {
      Face *pF = pFFace(iFace);
      Cell *pCL = pF->pCCellLeft();
      Cell *pCR = pF->pCCellRight();
      if (pF->qIsBdryFace())
	continue; // Do boundaries separately.
      // Mark internal boundaries
      if (pCL->iRegion() != pCR->iRegion()) {
	pF->vSetFaceLoc(Face::eBdryTwoSide);
	pF->pVVert(0)->vSetType(Vert::eBdryTwoSide);
	pF->pVVert(1)->vSetType(Vert::eBdryTwoSide);
	pF->pVVert(2)->vSetType(Vert::eBdryTwoSide);
      }
      else // Not an internal boundary
	pF->vSetFaceLoc(Face::eInterior);
    }
    
    // Tag ordinary boundary verts, and tag verts that lie at the
    // intersection of internal and regular boundary faces as boundary
    // apexes, so that they can never be smoothed or removed.
    // At this point, geometric tagging of BdryApex and BdryCurve
    // vertices isn't done.
    for (int iBFace = 0; iBFace < iNumBdryFaces(); iBFace++) {
      BFace *pBF = pBFBFace(iBFace);
      for (int iV = 0; iV < pBF->iNumVerts(); iV++) {
	Vert *pV = pBF->pVVert(iV);
	switch (pV->iVertType()) {
	case Vert::eBdryTwoSide:
	  pV->vSetType(Vert::eBdryApex);
	  break;
	case Vert::eInterior:
	  pV->vSetType(Vert::eBdry);
	  break;
	default:
	  // Nothing needs to be done
	  break;
	}
      }
    }

    if (pB3DIn) {
      // 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*, BdryPatch3D*> mmVertPatch;
      for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
	Vert *pV = pVVert(iV);
	for (int iP = pB3D->iNumPatches() - 1; iP >= 0; iP--) {
	  BdryPatch3D *pBP3D = (*pB3D)[iP];
	  if (!pBP3D->qDisjointFrom(pV)) {
	    vMessage(4, "Vert at (%10.6f, %10.6f, %10.6f) lies on patch %10p\n",
		     pV->dX(), pV->dY(), pV->dZ(), pBP3D);
	    mmVertPatch.insert(pair<Vert*, BdryPatch3D*>(pV, pBP3D));
	  }
	}
      }
      
      // 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);
	Vert *pV2 = pBF->pVVert(2);
	pV0->vSetType(Vert::eBdry);
	pV1->vSetType(Vert::eBdry);
	pV2->vSetType(Vert::eBdry);
	
	// How many patches for each one?
	pair<multimap<Vert*, BdryPatch3D*>::iterator,
	  multimap<Vert*, BdryPatch3D*>::iterator> pRange0, pRange1, pRange2;
	pRange0 = mmVertPatch.equal_range(pV0);
	pRange1 = mmVertPatch.equal_range(pV1);
	pRange2 = mmVertPatch.equal_range(pV2);
	multimap<Vert*, BdryPatch3D*>::iterator iter0, iter1, iter2;
	BdryPatch3D* pBP3D = NULL;
	for (iter0 = pRange0.first;
	     pBP3D == NULL && iter0 != pRange0.second;
	     iter0++) {
	  for (iter1 = pRange1.first;
	       pBP3D == NULL && iter1 != pRange1.second;
	       iter1++) {
	    if (iter0->second == iter1->second) {
	      for (iter2 = pRange2.first;
		   pBP3D == NULL && iter2 != pRange2.second;
		   iter2++) {
		if (iter0->second == iter2->second) {
		  pBP3D = iter0->second;
		} // Inner if
	      } // Loop for vert 3
	    } // 1 and 2 have a face in common
	  } // Loop for vert 2
	} // Loop for vert 1
	assert(pBP3D != NULL);
	vMessage(4, "BFace %3d lies on patch %10p\n", iBF, pBP3D);
	vMessage(4, "  coords: (%10.6f, %10.6f, %10.6f),\n",
		 pV0->dX(), pV0->dY(), pV0->dZ());
	vMessage(4, "  coords: (%10.6f, %10.6f, %10.6f),\n",
		 pV1->dX(), pV1->dY(), pV1->dZ());
	vMessage(4, "  coords: (%10.6f, %10.6f, %10.6f),\n",
		 pV2->dX(), pV2->dY(), pV2->dZ());
	pBF->vSetPatch(pBP3D);
      }
      
      // Definitively identify bdry apexes and bdry curve pts: any vert
      // with faceson
      // more than two patches is an apex; those on exactly two are bdry
      // curve pts.  This assumes specifically that there aren't any
      // internal bdrys.
      for (int iV = iNumVerts() - 1; iV >= 0; iV--) {
	Vert *pV = pVVert(iV);
	if (pV->iVertType() == Vert::eBdry) {
	  int iCount = mmVertPatch.count(pV);
	  if (iCount == 2) {
	    // FIX ME:
	    // Could also be a bizarre apex, where the curve separating
	    // two patches has a kink in it.
	    pV->vSetType(Vert::eBdryCurve);
	  }
	  else if (iCount > 2) {
	    // FIX ME:
	    // Could be a curve with one or more internal faces.
	    pV->vSetType(Vert::eBdryApex);
	  }
	  vMessage(4, "Patch count: %d.  Vert type: %d.\n",
		   iCount, pV->iVertType());
	}
      }
    }
    else {
      vMessage(2, "Creating bdry patches directly from mesh data.\n");
      vMessage(2, "You should specify a bdry file if you have one available.\n");
      vCreatePatchesFromBdryFaces(aiBCList);
    }
    delete [] aiBCList;

    if (qFaceVert)
      free(a2iFaceVert);
    if (qCellVert)
      free(a2iCellVert);
    if (qFaceCell)
      free(a2iFaceCell);
    if (qCellFace)
      free(a2iCellFace);
    if (qBFaceFace)
      free(aiBFaceFace);
    if (qBFaceVert)
      free(a2iBFaceVert);
    if (qBFaceBC)
      free(aiBFaceBC);
  }


  // Check to make sure the correct number of bdry faces exist, etc.
  if (qSimplicial()) {
    assert(4*iNumCells() + iNumBdryFaces() == 2*iNumFaces());
    assert(iNumQuadFaces() == 0);
    assert(iNumQuadBdryFaces() == 0);
    assert(iNumPyrCells() == 0);
    assert(iNumPrismCells() == 0);
    assert(iNumHexCells() == 0);
  }
  pQ = new Quality(this, iQualMeas);
  vSetAllHintFaces();
  if (!(qMeshOrientationOK(true)))
    vFatalError("Invalid or improperly-oriented mesh", "3D mesh input");
  SUMAA_LOG_EVENT_END(READ_MESH_3D);
}

VolMesh::VolMesh(const int iQualMeas)
  : Mesh (), qOKToChangeSurfaceEdges(true),
    qDoEdgeSwapping(true), qStrictPatchChecking(true),
    dMaxAngleForSurfSwap(0), VATriF(), VAQuadF(),
    VATriBF(), VAQuadBF(), VATet(), VAPyr(), VAPrism(), VAHex(),
    VAIntTriBF(), VAIntQuadBF(), pB3D(), qBdryFromMesh(true)
{
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(3, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);
  qAllowBdryChanges = true;
  qBdryFill = false;
}

VolMesh::VolMesh(const VolMesh& VM)
  : Mesh (VM),
    qOKToChangeSurfaceEdges(VM.qOKToChangeSurfaceEdges),
    qDoEdgeSwapping(VM.qDoEdgeSwapping),
    qStrictPatchChecking(VM.qStrictPatchChecking),
    dMaxAngleForSurfSwap(VM.dMaxAngleForSurfSwap),
    VATriF(VM.iNumTriFaces()), VAQuadF(VM.iNumQuadFaces()),
    VATriBF(VM.iNumTriBdryFaces()), VAQuadBF(VM.iNumQuadBdryFaces()),
    VATet(VM.iNumTetCells()), VAPyr(VM.iNumPyrCells()),
    VAPrism(VM.iNumPrismCells()), VAHex(VM.iNumHexCells()),
    VAIntTriBF(VM.iNumIntTriBdryFaces()),
    VAIntQuadBF(VM.iNumIntQuadBdryFaces()), pB3D(VM.pB3D),
    qBdryFromMesh(VM.qBdryFromMesh)
{
  qAllowBdryChanges = VM.qAllowBdryChanges;
  qLengthScaleFromCellSizes = VM.qLengthScaleFromCellSizes;
  qBdryFill = VM.qBdryFill;
  iIntSkip = VM.iIntSkip;

  // At this point, don't try copying more than the bare essentials of
  // the smoothing struct.
  {
    SMsmooth_data *pSMSmooth = static_cast<SMsmooth_data*>(VM.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(3, iTech, iFunc, dAcceptValue, &pvSmoothData);
    SMinitGlobalMinValue(pvSmoothData);
  }

  if (qSimplex) {
    assert(VM.iNumQuadFaces() == 0);
    assert(VM.iNumPyrCells() == 0);
    assert(VM.iNumPrismCells() == 0);
    assert(VM.iNumHexCells() == 0);
    assert(VM.iNumQuadBdryFaces() == 0);
    assert(VM.iNumIntQuadBdryFaces() == 0);
  }

  // Clear all face and cell data.  The data not specified is
  int i;

  // invalid by default.
  for (i = iNumFaces() - 1; i >= 0; i--)
    pFFace(i)->vAssign(pCInvalidCell, pCInvalidCell);
  for (i = iNumCells() - 1; i >= 0; i--)
    pCCell(i)->vAssign(pFInvalidFace, pFInvalidFace, pFInvalidFace);

  // For every face, copy the -meaning- of the pointers; the pointers
  // themselves are specifically -not- copied.  In the process, set up
  // the cell-to-face pointers on the fly.
  for (i = iNumFaces() - 1; i >= 0; i--) {
    Face *pFOld = VM.pFFace(i);
    Face *pF = pFFace(i);
    Cell *pCLOld = pFOld->pCCellLeft();
    Cell *pCROld = pFOld->pCCellRight();
    Cell *pCL = pCInvalidCell, *pCR = pCInvalidCell;

    // Translate the cells, taking into account that some "cells" are
    // actually boundaries.
    if (pCLOld->qValid()) {
      switch (pCLOld->eType()) {
      case Cell::eTriBFace:
      case Cell::eQuadBFace:
      case Cell::eIntTriBFace:
      case Cell::eIntQuadBFace:
	{
	  int iBFL = VM.iBFaceIndex(dynamic_cast<BFace*>(pCLOld));
	  BFace *pBFL = pBFBFace(iBFL);
	  pBFL->vSetPatch(dynamic_cast<TriBFaceBase*>(pCLOld)->pPatchPointer() );
	  pCL = pBFL;
	  break;
	}
      case Cell::eTet:
      case Cell::ePyr:
      case Cell::ePrism:
      case Cell::eHex:
	{
	  int iCL = VM.iCellIndex(pCLOld);
	  pCL = pCCell(iCL);
	  break;
	}
      default:
	assert(0);
      }
    }

    if (pCROld->qValid()) {
      switch (pCROld->eType()) {
      case Cell::eTriBFace:
      case Cell::eQuadBFace:
      case Cell::eIntTriBFace:
      case Cell::eIntQuadBFace:
	{
	  int iBFR = VM.iBFaceIndex(dynamic_cast<BFace*>(pCROld));
	  BFace *pBFR = pBFBFace(iBFR);
	  pBFR->vSetPatch(dynamic_cast<TriBFaceBase*>(pCROld)->pPatchPointer() );
	  pCR = pBFR;
	}
      case Cell::eTet:
      case Cell::ePyr:
      case Cell::ePrism:
      case Cell::eHex:
	{
	  int iCR = VM.iCellIndex(pCROld);
	  pCR = pCCell(iCR);
	}
      default:
	assert(0);
      }
    }

    Vert *apV[] = {pVInvalidVert, pVInvalidVert, pVInvalidVert,
		   pVInvalidVert};
    for (int ii = 0; ii < pFOld->iNumVerts(); ii++)
      apV[ii] = pVVert(VM.iVertIndex(pFOld->pVVert(ii)));
    pF->vAssign(pCL, pCR, apV[0], apV[1], apV[2], apV[3]);
    if (pCL->qValid()) pCL->vAddFace(pF);
    if (pCR->qValid()) pCR->vAddFace(pF);
  }

  vSetAllHintFaces();

  assert(qValid());
  if (!(qValid()))
    vFoundBug("VolMesh copying");
}

VolMesh::~VolMesh()
{
  VATriF.~VATriFace();
  VAQuadF.~VAQuadFace();
  VATriBF.~VATriBFace();
  VAQuadBF.~VAQuadBFace();
  VATet.~VATetCell();
  VAPyr.~VAPyrCell();
  VAPrism.~VAPrismCell();
  VAHex.~VAHexCell();
  VAIntTriBF.~VAIntTriBFace();
  VAIntQuadBF.~VAIntQuadBFace();
  if (qBdryFromMesh) {
    delete pB3D;       /// Bdry info
    pB3D = NULL;
  }
}

VolMesh& VolMesh::operator=(const VolMesh& VM)
{
  if (this != &VM) {
    int i;
    pB3D = VM.pB3D;
    qBdryFromMesh = VM.qBdryFromMesh;

    VAVerts = VM.VAVerts;
    pQ = VM.pQ;

    // At this point, don't try copying more than the bare essentials of
    // the smoothing struct.
    {
      SMsmooth_data *pSMSmooth = static_cast<SMsmooth_data*>(VM.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(3, iTech, iFunc, dAcceptValue, &pvSmoothData);
      SMinitGlobalMinValue(pvSmoothData);
    }

    qSimplex = VM.qSimplex;
    qOKToChangeSurfaceEdges = VM.qOKToChangeSurfaceEdges;
    qDoEdgeSwapping = VM.qDoEdgeSwapping;
    qStrictPatchChecking = VM.qStrictPatchChecking;
    dMaxAngleForSurfSwap = VM.dMaxAngleForSurfSwap;
    qAllowBdryChanges = VM.qAllowBdryChanges;
    qLengthScaleFromCellSizes = VM.qLengthScaleFromCellSizes;
    qBdryFill = VM.qBdryFill;
    iIntSkip = VM.iIntSkip;

    // Make sure that the arrays in *this are the right size.
    vSetupTriFaces(VM.iNumTriFaces());
    vSetupTetCells(VM.iNumTetCells());
    vSetupTriBFaces(VM.iNumTriBdryFaces());
    vSetupIntTriBFaces(VM.iNumIntTriBdryFaces());

    if (qSimplex) {
      assert(VM.iNumQuadFaces() == 0);
      assert(VM.iNumPyrCells() == 0);
      assert(VM.iNumPrismCells() == 0);
      assert(VM.iNumHexCells() == 0);
      assert(VM.iNumQuadBdryFaces() == 0);
      assert(VM.iNumIntQuadBdryFaces() == 0);
    }
    vSetupQuadFaces(VM.iNumQuadFaces());
    vSetupPyrCells(VM.iNumPyrCells());
    vSetupPrismCells(VM.iNumPrismCells());
    vSetupHexCells(VM.iNumHexCells());
    vSetupQuadBFaces(VM.iNumQuadBdryFaces());
    vSetupIntQuadBFaces(VM.iNumIntQuadBdryFaces());

    // Clear all face and cell data.  The data not specified is
    // invalid by default.
    for (i = iNumFaces() - 1; i >= 0; i--)
      pFFace(i)->vAssign(pCInvalidCell, pCInvalidCell);
    for (i = iNumCells() - 1; i >= 0; i--)
      pCCell(i)->vAssign(pFInvalidFace, pFInvalidFace);

    // For every face, copy the -meaning- of the pointers; the pointers
    // themselves are specifically -not- copied.  In the process, set up
    // the cell-to-face pointers on the fly.
    for (i = iNumFaces() - 1; i >= 0; i--) {
      Face *pFOld = VM.pFFace(i);
      Face *pF = pFFace(i);
      Cell *pCLOld = pFOld->pCCellLeft();
      Cell *pCROld = pFOld->pCCellRight();
      Cell *pCL = pCInvalidCell, *pCR = pCInvalidCell;

      // Translate the cells, taking into account that some "cells" are
      // actually boundaries.
      switch (pCLOld->eType()) {
	case Cell::eTriBFace:
	case Cell::eQuadBFace:
	case Cell::eIntTriBFace:
	case Cell::eIntQuadBFace:
	{
	  int iBFL = VM.iBFaceIndex(dynamic_cast<BFace*>(pCLOld));
	  BFace *pBFL = pBFBFace(iBFL);
	  pBFL->vSetPatch(dynamic_cast<BFace*>(pCLOld)->pPatchPointer() );
	  pCL = pBFL;
	  break;
	}
	case Cell::eTet:
	case Cell::ePyr:
	case Cell::ePrism:
	case Cell::eHex:
	{
	  int iCL = VM.iCellIndex(pCLOld);
	  pCL = pCCell(iCL);
	  break;
	}
	default:
	assert(0);
      }

      switch (pCROld->eType()) {
	case Cell::eTriBFace:
	case Cell::eQuadBFace:
	case Cell::eIntTriBFace:
	case Cell::eIntQuadBFace:
	{
	  int iBFR = VM.iBFaceIndex(dynamic_cast<BFace*>(pCROld));
	  BFace *pBFR = pBFBFace(iBFR);
	  pBFR->vSetPatch(dynamic_cast<BFace*>(pCROld)->pPatchPointer() );
	  pCR = pBFR;
	}
	case Cell::eTet:
	case Cell::ePyr:
	case Cell::ePrism:
	case Cell::eHex:
	{
	  int iCR = VM.iCellIndex(pCROld);
	  pCR = pCCell(iCR);
	}
	default:
	assert(0);
      }

      Vert *apV[] = {pVInvalidVert, pVInvalidVert, pVInvalidVert,
		    pVInvalidVert};
      for (int ii = 0; ii < pFOld->iNumVerts(); ii++)
	apV[ii] = pVVert(VM.iVertIndex(pFOld->pVVert(ii)));
      pF->vAssign(pCL, pCR, apV[0], apV[1], apV[2], apV[3]);
      pCL->vAddFace(pF);
      pCR->vAddFace(pF);
    }

    vSetAllHintFaces();

    assert(qValid());
    if (!(qValid()))
      vFoundBug("VolMesh copying");
  }
  return (*this);
}

void VolMesh::vConvertFromCellVert(const int iNCells,
				   const int iNBdryFaces,
				   const int a2iCellVert[][4],
				   const int a2iBFaceVert[][3],
				   const int aiBFaceBC[],
				   int*& aiBCList)
{
  qSimplex = true;
  vAllowEdgeSwapping();
  vAllowSurfaceEdgeChanges();

  // Now create list of all faces and sort to pair the two sides of the face
  FaceSort *aFS = new FaceSort[4*iNCells + iNBdryFaces];
  int i, iFace = 0;
  for (i = 0; i < iNBdryFaces; i++) {
    Vert* pV0 = pVVert(a2iBFaceVert[i][0]);
    Vert* pV1 = pVVert(a2iBFaceVert[i][1]);
    Vert* pV2 = pVVert(a2iBFaceVert[i][2]);
    int iBC = aiBFaceBC[i];
    aFS[iFace++] = FaceSort(pV0, pV2, pV1, -abs(iBC));
  }
  for (i = 0; i < iNCells; i++) {
    Vert* pV0 = pVVert(a2iCellVert[i][0]);
    Vert* pV1 = pVVert(a2iCellVert[i][1]);
    Vert* pV2 = pVVert(a2iCellVert[i][2]);
    Vert* pV3 = pVVert(a2iCellVert[i][3]);
    aFS[iFace++] = FaceSort(pV0, pV1, pV2, i);
    aFS[iFace++] = FaceSort(pV0, pV3, pV1, i);
    aFS[iFace++] = FaceSort(pV1, pV3, pV2, i);
    aFS[iFace++] = FaceSort(pV2, pV3, pV0, i);
  }
  int iNFaceSortRecs = iFace;
  // Now sFS contains all the faces in sorted order.
  sort(aFS, aFS + iNFaceSortRecs);

  // Find out how many faces there actually are.
  int iNFaces = 0, iNActualBdryFaces = 0;
  for (i = 0; i < iNFaceSortRecs; i++) {
    iNFaces++;
    if (i == iNFaceSortRecs - 1) {
      iNActualBdryFaces++;
      continue;
    }
    FaceSort FS0 = aFS[i  ];
    FaceSort FS1 = aFS[i+1];

    // Two faces must share the same verts, point in opposite
    // directions, and have different cells on their "left" side
    if ((FS0.pVVert(0) == FS1.pVVert(0) ||
	 FS0.pVVert(0) == FS1.pVVert(1) ||
	 FS0.pVVert(0) == FS1.pVVert(2))
	&&
	(FS0.pVVert(1) == FS1.pVVert(0) ||
	 FS0.pVVert(1) == FS1.pVVert(1) ||
	 FS0.pVVert(1) == FS1.pVVert(2))
	&&
	(FS0.pVVert(2) == FS1.pVVert(0) ||
	 FS0.pVVert(2) == FS1.pVVert(1) ||
	 FS0.pVVert(2) == FS1.pVVert(2))) {
      assert(FS0.iLeftCell() > FS1.iLeftCell());
      if (FS1.iLeftCell() < 0) {
	// A bdry face, with BC data given.
	iNActualBdryFaces++;
      }
      i++; // In any event, need to skip over this pair.
    }
    else {
      iNActualBdryFaces++;
    }
  }

  assert(4*iNCells == 2*iNFaces - iNActualBdryFaces);
  assert(iNBdryFaces == 0 || iNActualBdryFaces == iNBdryFaces);

  if (iNBdryFaces != iNActualBdryFaces) {
    // Need to resize aiBCList
    delete [] aiBCList;
    aiBCList = new int[iNActualBdryFaces];
  }

  // Make certain that the right amount of space is available for
  // connectivity data
  vSetupTriFaces(iNFaces);
  vSetupTriBFaces(iNActualBdryFaces);
  vSetupTetCells(iNCells);

  // Read out the data into the face-vert, face-cell, and cell-face
  // tables
  int iBFace = 0;
  iFace = 0;
  for (i = 0; i < iNFaceSortRecs; i++) {
    FaceSort FS0 = aFS[i  ];
    FaceSort FS1 = aFS[i+1];

    // Two faces must share the same verts, point in opposite
    // directions, and have different cells on their "left" side
    if ((i < iNFaceSortRecs - 1) &&
	((FS0.pVVert(0) == FS1.pVVert(0) ||
	  FS0.pVVert(0) == FS1.pVVert(1) ||
	  FS0.pVVert(0) == FS1.pVVert(2))
	 &&
	 (FS0.pVVert(1) == FS1.pVVert(0) ||
	  FS0.pVVert(1) == FS1.pVVert(1) ||
	  FS0.pVVert(1) == FS1.pVVert(2))
	 &&
	 (FS0.pVVert(2) == FS1.pVVert(0) ||
	  FS0.pVVert(2) == FS1.pVVert(1) ||
	  FS0.pVVert(2) == FS1.pVVert(2)))) {
#ifndef NDEBUG
      {
	double adN0[3], adN1[3];
	vNormal(FS0.pVVert(0)->adCoords(), FS0.pVVert(1)->adCoords(),
		FS0.pVVert(2)->adCoords(), adN0);
	vNormal(FS1.pVVert(0)->adCoords(), FS1.pVVert(1)->adCoords(),
		FS1.pVVert(2)->adCoords(), adN1);
	vNORMALIZE3D(adN0);
	vNORMALIZE3D(adN1);
	assert(0 == iFuzzyComp(-1., dDOT3D(adN0, adN1)));
      }
#endif
      assert(FS0.iLeftCell() > FS1.iLeftCell());
      i++;
      Face* pF = pFFace(iFace++);
      Cell *pC0, *pC1;
      if (FS1.iLeftCell() < 0) { // Boundary face
	BFace* pBF = pBFBFace(iBFace);
	pBF->vAssign(pF);
	aiBCList[iBFace] = abs(FS1.iLeftCell());
	pC0 = pBF;
	iBFace++;
      }
      else {
	pC0 = pCCell(FS1.iLeftCell());
	pC0->vAddFace(pF);
      }
      assert(FS0.iLeftCell() >= 0);
      pC1 = pCCell(FS0.iLeftCell());
      pC1->vAddFace(pF);
      Vert* pVert0 = FS0.pVVert(0);
      Vert* pVert1 = FS0.pVVert(1);
      Vert* pVert2 = FS0.pVVert(2);
      pF->vAssign(pC0, pC1, pVert0, pVert1, pVert2);
    }
    else {
      // Create a BFace out of whole cloth.
      assert(FS0.iLeftCell() >= 0);
      Face *pF = pFFace(iFace++);
      Cell *pC = pCCell(FS0.iLeftCell());
      Vert* pVert0 = FS0.pVVert(0);
      Vert* pVert1 = FS0.pVVert(1);
      Vert* pVert2 = FS0.pVVert(2);
      pC->vAddFace(pF);
      BFace *pBF = pBFBFace(iBFace);
      pBF->vAssign(pF);
      aiBCList[iBFace] = iDefaultBC;
      iBFace++;
      pF->vAssign(pBF, pC, pVert0, pVert1, pVert2);
    }
  }
  assert(iBFace == iNActualBdryFaces);
  assert(iFace == iNFaces);
  vSetAllHintFaces();
  delete [] aFS;
  if (!(qValid()))
    vFoundBug("conversion from cell->vert to face->cell mesh in 3D");
}

// This is the more general version, that will handle all
// GRUMMP-supported 3D cell types, and both tri and quad bdry faces.
void VolMesh::
vConvertFromCellVert(const int nTet, const int nPyr, const int nPrism,
		     const int nHex, const int nBdryTri, const int nBdryQuad,
		     const int tetToVert[][4], const int pyrToVert[][5],
		     const int prismToVert[][6], const int hexToVert[][8],
		     const int bdryTriToVert[][3],
		     const int bdryQuadToVert[][4])
{
  qSimplex = (nPyr + nPrism + nHex + nBdryQuad == 0);
  vAllowEdgeSwapping();
  vAllowSurfaceEdgeChanges();

  assert(nBdryTri + nBdryQuad > 0);

  // Now create list of all faces and sort to pair the two sides of the face
  int iNTriFaces = (4*nTet + 4*nPyr + 2*nPrism + nBdryTri) / 2;
  int iNQuadFaces = (nPyr + 3*nPrism + 6*nHex + nBdryQuad) / 2;
  FaceSort *aFSTri = new FaceSort[iNTriFaces*2];
  FaceSort *aFSQuad = new FaceSort[iNQuadFaces*2];

  int i, iTriFace = 0, iQuadFace = 0;
  // First put the bdry tris in the list
  for (i = 0; i < nBdryTri; i++) {
    Vert* pV0 = pVVert(bdryTriToVert[i][0]);
    Vert* pV1 = pVVert(bdryTriToVert[i][1]);
    Vert* pV2 = pVVert(bdryTriToVert[i][2]);
    aFSTri[iTriFace++] = FaceSort(pV0, pV2, pV1, pVInvalidVert, -iDefaultBC, 0);
  }

  // Now the bdry quads in their list
  for (i = 0; i < nBdryQuad; i++) {
    Vert* pV0 = pVVert(bdryQuadToVert[i][0]);
    Vert* pV1 = pVVert(bdryQuadToVert[i][1]);
    Vert* pV2 = pVVert(bdryQuadToVert[i][2]);
    Vert* pV3 = pVVert(bdryQuadToVert[i][3]);
    aFSQuad[iQuadFace++] = FaceSort(pV0, pV3, pV2, pV1, -iDefaultBC, 0);
  }

  vMessage(2, "Adding faces of tets\n");
  // Now all four faces of the tets.
  for (i = 0; i < nTet; i++) {
    Vert* pV0 = pVVert(tetToVert[i][0]);
    Vert* pV1 = pVVert(tetToVert[i][1]);
    Vert* pV2 = pVVert(tetToVert[i][2]);
    Vert* pV3 = pVVert(tetToVert[i][3]);
    aFSTri[iTriFace++] = FaceSort(pV0, pV1, pV2, pVInvalidVert, i, 0);
    aFSTri[iTriFace++] = FaceSort(pV0, pV3, pV1, pVInvalidVert, i, 1);
    aFSTri[iTriFace++] = FaceSort(pV1, pV3, pV2, pVInvalidVert, i, 2);
    aFSTri[iTriFace++] = FaceSort(pV2, pV3, pV0, pVInvalidVert, i, 3);
  }

  vMessage(2, "Adding faces of pyramids\n");
  // All five faces of the pyramids.
  for (i = 0; i < nPyr; i++) {
    // The following five lines include a translation from the UGrid
    // vertex ordering for a pyramid to the GRUMMP (<- Nastran(?))
    // ordering. 
    Vert* pV0 = pVVert(pyrToVert[i][0]);
    Vert* pV1 = pVVert(pyrToVert[i][1]);
    Vert* pV2 = pVVert(pyrToVert[i][4]);
    Vert* pV3 = pVVert(pyrToVert[i][3]);
    Vert* pV4 = pVVert(pyrToVert[i][2]);
    aFSTri[iTriFace++] = FaceSort(pV0, pV1, pV4, pVInvalidVert, i + nTet, 0);
    aFSTri[iTriFace++] = FaceSort(pV1, pV2, pV4, pVInvalidVert, i + nTet, 1);
    aFSTri[iTriFace++] = FaceSort(pV2, pV3, pV4, pVInvalidVert, i + nTet, 2);
    aFSTri[iTriFace++] = FaceSort(pV3, pV0, pV4, pVInvalidVert, i + nTet, 3);
    aFSQuad[iQuadFace++] = FaceSort(pV0, pV3, pV2, pV1, i + nTet, 4);
  }
  
  vMessage(2, "Adding faces of prisms\n");
  // All five faces of the prisms.
  for (i = 0; i < nPrism; i++) {
    Vert* pV0 = pVVert(prismToVert[i][0]);
    Vert* pV1 = pVVert(prismToVert[i][1]);
    Vert* pV2 = pVVert(prismToVert[i][2]);
    Vert* pV3 = pVVert(prismToVert[i][3]);
    Vert* pV4 = pVVert(prismToVert[i][4]);
    Vert* pV5 = pVVert(prismToVert[i][5]);
    aFSQuad[iQuadFace++] = FaceSort(pV0, pV1, pV4, pV3, i + nTet + nPyr, 0);
    aFSQuad[iQuadFace++] = FaceSort(pV1, pV2, pV5, pV4, i + nTet + nPyr, 1);
    aFSQuad[iQuadFace++] = FaceSort(pV2, pV0, pV3, pV5, i + nTet + nPyr, 2);
    aFSTri[iTriFace++] = FaceSort(pV0, pV2, pV1, pVInvalidVert,
				  i + nTet + nPyr, 3);
    aFSTri[iTriFace++] = FaceSort(pV3, pV4, pV5, pVInvalidVert,
				  i + nTet + nPyr, 4);
  }

  vMessage(2, "Adding faces of hexes\n");
  // All six faces of the hexes.
  for (i = 0; i < nHex; i++) {
    Vert* pV0 = pVVert(hexToVert[i][0]);
    Vert* pV1 = pVVert(hexToVert[i][1]);
    Vert* pV2 = pVVert(hexToVert[i][2]);
    Vert* pV3 = pVVert(hexToVert[i][3]);
    Vert* pV4 = pVVert(hexToVert[i][4]);
    Vert* pV5 = pVVert(hexToVert[i][5]);
    Vert* pV6 = pVVert(hexToVert[i][6]);
    Vert* pV7 = pVVert(hexToVert[i][7]);
    aFSQuad[iQuadFace++] = FaceSort(pV0, pV3, pV2, pV1,
				    i + nTet + nPyr + nPrism, 0);
    aFSQuad[iQuadFace++] = FaceSort(pV0, pV1, pV5, pV4,
				    i + nTet + nPyr + nPrism, 1);
    aFSQuad[iQuadFace++] = FaceSort(pV1, pV2, pV6, pV5,
				    i + nTet + nPyr + nPrism, 2);
    aFSQuad[iQuadFace++] = FaceSort(pV2, pV3, pV7, pV6,
				    i + nTet + nPyr + nPrism, 3);
    aFSQuad[iQuadFace++] = FaceSort(pV3, pV0, pV4, pV7,
				    i + nTet + nPyr + nPrism, 4);
    aFSQuad[iQuadFace++] = FaceSort(pV4, pV5, pV6, pV7,
				    i + nTet + nPyr + nPrism, 5);
  }

  assert(iQuadFace == iNQuadFaces*2);
  assert(iTriFace == iNTriFaces*2);

  // Now sort all this stuff.
  vMessage(2, "Sorting tri faces\n");
  sort(aFSTri, aFSTri + iNTriFaces*2);
  vMessage(2, "Sorting quad faces\n");
  sort(aFSQuad, aFSQuad + iNQuadFaces*2);

  // Make certain that the right amount of space is available for
  // connectivity data
  vSetupTriFaces(iNTriFaces);
  vSetupQuadFaces(iNQuadFaces);

  vSetupTriBFaces(nBdryTri);
  vSetupQuadBFaces(nBdryQuad);

  vSetupTetCells(nTet);
  vSetupPyrCells(nPyr);
  vSetupPrismCells(nPrism);
  vSetupHexCells(nHex);

  // Read out the data into the face-vert, face-cell, and cell-face
  // tables.  This is made easier by the requirement that bdry face data
  // -must- be given, so all faces must match up in pairs.
  int iTriBFace = 0;
  iTriFace = 0;
  vMessage(2, "Setting up tri faces\n");
  for (i = 0; i < iNTriFaces*2; i+=2) {
    if ((i % 500000) == 0 && i > 0) {
      vMessage(2, "Setting tri face %d out of %d\n", i/2, iNTriFaces);
    }
    FaceSort FS0 = aFSTri[i  ];
    FaceSort FS1 = aFSTri[i+1];

    assert(qIdenticalVerts(FS0, FS1));
    assert(FS0.iLeftCell() > FS1.iLeftCell());

    Face* pF = pFFace(iTriFace++);
    Cell *pC0, *pC1;
    if (FS1.iLeftCell() < 0) { // Boundary face
      BFace* pBF = pBFBFace(iTriBFace);
      pBF->vAssign(pF);
      pC0 = pBF;
      iTriBFace++;
    }
    else {
      pC0 = pCCell(FS1.iLeftCell());
      pC0->vAddFace(pF, FS1.iWhichFaceInCell());
    }

    assert(FS0.iLeftCell() >= 0);
    pC1 = pCCell(FS0.iLeftCell());
    pC1->vAddFace(pF, FS0.iWhichFaceInCell());
    
    Vert* pVert0 = FS0.pVVert(0);
    Vert* pVert1 = FS0.pVVert(1);
    Vert* pVert2 = FS0.pVVert(2);
    pF->vAssign(pC0, pC1, pVert0, pVert1, pVert2);
  }
  assert(iTriBFace == nBdryTri);
  assert(iTriFace == iNTriFaces);

  int iQuadBFace = 0;
  iQuadFace = 0;
  vMessage(2, "Setting up quad faces\n");
  for (i = 0; i < iNQuadFaces*2; i+=2) {
    if ((i % 500000) == 0) {
      vMessage(2, "Setting quad face %d out of %d\n", i/2, iNQuadFaces);
    }
    FaceSort FS0 = aFSQuad[i  ];
    FaceSort FS1 = aFSQuad[i+1];

    assert(iQuadFace*2 == i);
    assert(qIdenticalVerts(FS0, FS1));
    assert(FS0.iLeftCell() > FS1.iLeftCell());

    Face* pF = pFFace(iQuadFace + iNTriFaces);
    iQuadFace++;
    Cell *pC0, *pC1;
    if (FS1.iLeftCell() < 0) { // Boundary face
      BFace* pBF = pBFBFace(iQuadBFace + nBdryTri);
      pBF->vAssign(pF);
      pC0 = pBF;
      iQuadBFace++;
    }
    else {
      pC0 = pCCell(FS1.iLeftCell());
      pC0->vAddFace(pF, FS1.iWhichFaceInCell());
    }

    assert(FS0.iLeftCell() >= 0);
    pC1 = pCCell(FS0.iLeftCell());
    pC1->vAddFace(pF, FS0.iWhichFaceInCell());
    
    Vert* pVert0 = FS0.pVVert(0);
    Vert* pVert1 = FS0.pVVert(1);
    Vert* pVert2 = FS0.pVVert(2);
    Vert* pVert3 = FS0.pVVert(3);
    pF->vAssign(pC0, pC1, pVert0, pVert1, pVert2, pVert3);
  }
  assert(iQuadBFace == nBdryQuad);
  assert(iQuadFace == iNQuadFaces);

  vSetAllHintFaces();
  delete [] aFSTri;
  delete [] aFSQuad;
  if (!(qValid()))
    vFoundBug("conversion from cell->vert to face->cell mesh in 3D");
}

//@ Retrieval of entities by index
Face*  VolMesh::pFFace(const int i) const
{
  assert(i >= 0 && i < iNumFaces());
  if (i < iNumTriFaces())
    return VATriF.pTEntry(i);
  else
    return VAQuadF.pTEntry(i - iNumTriFaces());
}

BFace* VolMesh::pBFBFace(const int i) const
{
  assert(i >= 0 && i < iNumTotalBdryFaces());
  int iNTriBF = iNumTriBdryFaces();
  if (i < iNTriBF)
    return VATriBF.pTEntry(i);

  int iNRegBF = iNTriBF + iNumQuadBdryFaces();
  assert(iNRegBF == iNumBdryFaces());
  if (i < iNRegBF)
    return VAQuadBF.pTEntry(i - iNTriBF);

  int iNTQIntTriBF = iNRegBF + iNumIntTriBdryFaces();
  if (i < iNTQIntTriBF)
    return VAIntTriBF.pTEntry(i - iNRegBF);

  else
    return VAIntQuadBF.pTEntry(i - iNTQIntTriBF);
}

BFace* VolMesh::pBFIntBFace(const int i) const
{
  int iNIntTriBF = iNumIntTriBdryFaces();
  if (i < iNIntTriBF)
    return VAIntTriBF.pTEntry(i);

  else
    return VAIntQuadBF.pTEntry(i - iNIntTriBF);
}

Cell*  VolMesh::pCCell(const int i) const
{
  assert(i >= 0 && i < iNumCells());
  int iNTet = iNumTetCells();

  if (i < iNTet)
    return VATet.pTEntry(i);

  int iNTPyr = iNTet + iNumPyrCells();
  if (i < iNTPyr)
    return VAPyr.pTEntry(i - iNTet);

  int iNTPPrism = iNTPyr + iNumPrismCells();
  if (i < iNTPPrism)
    return VAPrism.pTEntry(i - iNTPyr);

  else
    return VAHex.pTEntry(i - iNTPPrism);
}

int VolMesh::iBFaceIndex(const BFace* const pBF) const
{
  switch (pBF->eType()) {
  case Cell::eTriBFace:
    return VATriBF.iIndex(dynamic_cast<const TriBFace*>(pBF));
  case Cell::eQuadBFace:
    assert(!qSimplex);
    return VAQuadBF.iIndex(dynamic_cast<const QuadBFace*>(pBF)) + VATriBF.iSize();
  case Cell::eIntTriBFace:
    return VAIntTriBF.iIndex(dynamic_cast<const IntTriBFace*>(pBF)) + iNumBdryFaces();
  case Cell::eIntQuadBFace:
    assert(!qSimplex);
    return VAIntQuadBF.iIndex(dynamic_cast<const IntQuadBFace*>(pBF)) +
      iNumBdryFaces() + VAIntTriBF.iSize();
  default:
    assert(0);
    return(-1);
  }
}

int VolMesh::iCellIndex(const Cell* const pC) const
{
  assert(pC->eType() == Cell::eTet ||
	 (!qSimplex &&
	  (pC->eType() == Cell::ePyr ||
	   pC->eType() == Cell::ePrism ||
	   pC->eType() == Cell::eHex) ) );
  if (pC->eType() == Cell::eTet)
    return VATet.iIndex(dynamic_cast<const TetCell*>(pC));
  else if (pC->eType() == Cell::ePyr)
    return VAPyr.iIndex(dynamic_cast<const PyrCell*>(pC)) + VATet.iSize();
  else if (pC->eType() == Cell::ePrism)
    return VAPrism.iIndex(dynamic_cast<const PrismCell*>(pC))
	  + VAPyr.iSize() + VATet.iSize();
  else {
    return VAHex.iIndex(dynamic_cast<const HexCell*>(pC)) + VAPrism.iSize()
      + VAPyr.iSize() + VATet.iSize();
  }
}

int VolMesh::iFaceIndex(const Face* const pF) const
{
  assert(pF->eType() == Face::eTriFace ||
	 (pF->eType() == Face::eQuadFace && !qSimplex));
  if (pF->eType() == Face::eTriFace)
    return VATriF.iIndex(dynamic_cast<const TriFace*>(pF));
  else {
    return VAQuadF.iIndex(dynamic_cast<const QuadFace*>(pF)) + VATriF.iSize();
  }
}

bool VolMesh::
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 (VAVerts.qValid( static_cast<Vert*>(pvEnt) )      ||
	  VATriF.qValid( static_cast<TriFace*>(pvEnt) )    ||
	  VATet.qValid( static_cast<TetCell*>(pvEnt) )     ||
	  VAQuadF.qValid( static_cast<QuadFace*>(pvEnt) )  ||
	  VAPyr.qValid( static_cast<PyrCell*>(pvEnt) )     ||
	  VAPrism.qValid( static_cast<PrismCell*>(pvEnt) ) ||
	  VAHex.qValid( static_cast<HexCell*>(pvEnt) )     );
}

//@ Create new entities
Face* VolMesh::pFNewFace(const int iNV)
{
  assert(iNV == 3 || iNV == 4);
  switch (iNV) {
  case 3:
    {
      TriFace *pTF = VATriF.pTNewEntry();
      pTF->vSetDefaultFlags();
      return pTF;
    }
  case 4:
    {
      QuadFace *pQF = VAQuadF.pTNewEntry();
      pQF->vSetDefaultFlags();
      return pQF;
    }
  default:
    assert(0);
    return (pFInvalidFace);
  }
}

BFace* VolMesh::pBFNewBFace(const int iNBV, const bool qIsInternal)
{
  assert(iNBV == 3 || iNBV == 4);
  BFace *pBFRetVal = pBFInvalidBFace;
  switch (iNBV) {
  case 3:
    if (qIsInternal) {
      pBFRetVal = VAIntTriBF.pTNewEntry();
      pBFRetVal->vSetType(Cell::eIntTriBFace);
    }
    else {
      pBFRetVal = VATriBF.pTNewEntry();
      pBFRetVal->vSetType(Cell::eTriBFace);
    }
    break;
  case 4:
    if (qIsInternal) {
      pBFRetVal = VAIntQuadBF.pTNewEntry();
      pBFRetVal->vSetType(Cell::eIntQuadBFace);
    }
    else {
      pBFRetVal = VAQuadBF.pTNewEntry();
      pBFRetVal->vSetType(Cell::eQuadBFace);
    }
    break;
  default:
    assert(0);
    return (pBFInvalidBFace);
  }
  assert(pBFRetVal->qValid());
  pBFRetVal->vClearFaces();
  pBFRetVal->vSetDefaultFlags();
  return pBFRetVal;
}

Cell* VolMesh::pCNewCell(const int iNF, const int iNV)
{
  assert((iNF == 4 && iNV == 4) || // Tetrahedron
	 (iNF == 5 && (iNV == 5 || iNV== 6)) || // Pyramid and prism
	 (iNF == 6 && iNV == 8)); // Hexahedron
  switch (iNV) {
  case 4:
    {
      TetCell *pTet = VATet.pTNewEntry();
      pTet->vSetDefaultFlags();
      pTet->vSetType(Cell::eTet);
      return pTet;
    }
  case 5:
    {
      PyrCell *pPyr = VAPyr.pTNewEntry();
      pPyr->vSetDefaultFlags();
      pPyr->vSetType(Cell::ePyr);
      return pPyr;
    }
  case 6:
    {
      PrismCell *pPrism = VAPrism.pTNewEntry();
      pPrism->vSetDefaultFlags();
      pPrism->vSetType(Cell::ePrism);
      return pPrism;
    }
  case 8:
    {
      HexCell *pHex = VAHex.pTNewEntry();
      pHex->vSetDefaultFlags();
      pHex->vSetType(Cell::eHex);
      return pHex;
    }
  default:
    assert(0);
    return (pCInvalidCell);
  }
}

//@ Weed out empty entries in the connectivity tables
void VolMesh::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
  if (iNumTriFaces() > 0) {
    TriFace *pTFFor  = VATriF.pTFirst();
    TriFace *pTFBack = VATriF.pTLast();
    
    int iNewSize = 0;
    assert(pTFFor != pTFBack);
    
    while (pTFFor != pTFBack) {
      while (pTFFor != pTFBack && !pTFFor->qDeleted()) {
	pTFFor = VATriF.pTNext(pTFFor);
	iNewSize++;
      }
      
      while (pTFBack != pTFFor && pTFBack->qDeleted())
	pTFBack = VATriF.pTPrev(pTFBack);
      
      if (pTFFor != pTFBack) {
	assert(pTFFor->qDeleted());
	iNewSize++;
	(*pTFFor) = (*pTFBack);
	pTFBack->vMarkDeleted();
	pTFFor = VATriF.pTNext(pTFFor);
      }
      else if (!pTFFor->qDeleted()) {
	// This case only arises if there are no deleted faces to be purged.
	iNewSize++;
      }
    } // Loop over all faces
    // This resets the number of entries in use to iNewSize.
    VATriF.vSetup(iNewSize);
  } // Done with tris
  
  if (iNumQuadFaces() > 0) {
    QuadFace *pQFFor  = VAQuadF.pTFirst();
    QuadFace *pQFBack = VAQuadF.pTLast();
    
    int iNewSize = 0;
    assert(pQFFor != pQFBack);
    
    while (pQFFor != pQFBack) {
      while (pQFFor != pQFBack && !pQFFor->qDeleted()) {
	pQFFor = VAQuadF.pTNext(pQFFor);
	iNewSize++;
      }
      
      while (pQFBack != pQFFor && pQFBack->qDeleted())
	pQFBack = VAQuadF.pTPrev(pQFBack);
      
      if (pQFFor != pQFBack) {
	assert(pQFFor->qDeleted());
	iNewSize++;
	(*pQFFor) = (*pQFBack);
	pQFBack->vMarkDeleted();
	pQFFor = VAQuadF.pTNext(pQFFor);
      }
      else if (!pQFFor->qDeleted()) {
	// This case only arises if there are no deleted faces to be purged.
	iNewSize++;
      }
    } // Loop over all faces
    // This resets the number of entries in use to iNewSize.
    VAQuadF.vSetup(iNewSize);
  } // Done with quads
    
  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

void VolMesh::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
  // Purge tetrahedra
  {
    int iFor = 0, iBack = iNumTetCells() - 1;

    int iNewSize = 0;
    // The case of exactly one valid cell
    if ((iFor == iBack) &&
	pCCell(iFor)->qValid() &&
	!pCCell(iFor)->qDeleted()) {
      iNewSize = 1;
      iFor = 1; // Necessary because of an assertion later.
    }

    while (iFor < iBack) {
      while (!pCCell(iFor)->qDeleted() && iFor < iBack){
	iFor++;
	iNewSize++;
      }

      while (pCCell(iBack)->qDeleted() && iFor < iBack)
	iBack--;

      if (iFor != iBack) {
	Cell *pCFor = pCCell(iFor);
	Cell *pCBack = pCCell(iBack);

	assert(pCFor->qDeleted());
	assert(!pCBack->qDeleted());
	(*pCFor) = (*pCBack);
	pCBack->vMarkDeleted();
// #ifndef NDEBUG
//	for (int i = 0; i < pCFor->iNumFaces(); i++)
//	  assert(pCFor->pFFace(i)->iFullCheck());
// #endif
	iNewSize++;
	iFor++;
      }
      else if (!pCCell(iFor)->qDeleted()) {
	iNewSize++;
	iFor++;
      }
    }
    assert(iFor == iNewSize);
    // Reset the number of tet cells in use
    VATet.vSetup(iNewSize);
  } // Tetrahedron iteration

  // Purge pyramids
  {
    int iFor = iNumTetCells(), iBack = iFor + iNumPyrCells() - 1;

    int iNewSize = 0;
    // The case of exactly one valid cell
    if ((iFor == iBack) &&
	pCCell(iFor)->qValid() &&
	!pCCell(iFor)->qDeleted()) {
      iNewSize = 1;
      iFor = 1; // Necessary because of an assertion later.
    }

    while (iFor < iBack) {
      while (!pCCell(iFor)->qDeleted() && iFor < iBack){
	iFor++;
	iNewSize++;
      }

      while (pCCell(iBack)->qDeleted() && iFor < iBack)
	iBack--;

      if (iFor != iBack) {
	Cell *pCFor = pCCell(iFor);
	Cell *pCBack = pCCell(iBack);

	assert(pCFor->qDeleted());
	assert(!pCBack->qDeleted());
	(*pCFor) = (*pCBack);
	pCBack->vMarkDeleted();
// #ifndef NDEBUG
//	for (int i = 0; i < pCFor->iNumFaces(); i++)
//	  assert(pCFor->pFFace(i)->iFullCheck());
// #endif
	iNewSize++;
	iFor++;
      }
      else if (!pCCell(iFor)->qDeleted()) {
	iNewSize++;
	iFor++;
      }
    }
    assert(iFor == iNewSize + iNumTetCells());
    // Reset the number of tet cells in use
    VAPyr.vSetup(iNewSize);
  } // Pyramid iteration

  // Purge prisms
  {
    int iFor = iNumTetCells() + iNumPyrCells(),
      iBack = iFor + iNumPrismCells() - 1;

    int iNewSize = 0;
    // The case of exactly one valid cell
    if ((iFor == iBack) &&
	pCCell(iFor)->qValid() &&
	!pCCell(iFor)->qDeleted()) {
      iNewSize = 1;
      iFor = 1; // Necessary because of an assertion later.
    }

    while (iFor < iBack) {
      while (!pCCell(iFor)->qDeleted() && iFor < iBack){
	iFor++;
	iNewSize++;
      }

      while (pCCell(iBack)->qDeleted() && iFor < iBack)
	iBack--;

      if (iFor != iBack) {
	Cell *pCFor = pCCell(iFor);
	Cell *pCBack = pCCell(iBack);

	assert(pCFor->qDeleted());
	assert(!pCBack->qDeleted());
	(*pCFor) = (*pCBack);
	pCBack->vMarkDeleted();
// #ifndef NDEBUG
//	for (int i = 0; i < pCFor->iNumFaces(); i++)
//	  assert(pCFor->pFFace(i)->iFullCheck());
// #endif
	iNewSize++;
	iFor++;
      }
      else if (!pCCell(iFor)->qDeleted()) {
	iNewSize++;
	iFor++;
      }
    }
    assert(iFor == iNewSize + iNumTetCells() + iNumPyrCells());
    // Reset the number of prism cells in use
    VAPrism.vSetup(iNewSize);
  } // Prism iteration

  // Purge hexahedra
  {
    int iFor = iNumTetCells() + iNumPyrCells() + iNumPrismCells(),
      iBack = iFor + iNumHexCells() - 1;

    int iNewSize = 0;
    // The case of exactly one valid cell
    if ((iFor == iBack) &&
	pCCell(iFor)->qValid() &&
	!pCCell(iFor)->qDeleted()) {
      iNewSize = 1;
      iFor = 1; // Necessary because of an assertion later.
    }

    while (iFor < iBack) {
      while (!pCCell(iFor)->qDeleted() && iFor < iBack){
	iFor++;
	iNewSize++;
      }

      while (pCCell(iBack)->qDeleted() && iFor < iBack)
	iBack--;

      if (iFor != iBack) {
	Cell *pCFor = pCCell(iFor);
	Cell *pCBack = pCCell(iBack);

	assert(pCFor->qDeleted());
	assert(!pCBack->qDeleted());
	(*pCFor) = (*pCBack);
	pCBack->vMarkDeleted();
// #ifndef NDEBUG
//	for (int i = 0; i < pCFor->iNumFaces(); i++)
//	  assert(pCFor->pFFace(i)->iFullCheck());
// #endif
	iNewSize++;
	iFor++;
      }
      else if (!pCCell(iFor)->qDeleted()) {
	iNewSize++;
	iFor++;
      }
    }
    assert(iFor ==
	   iNewSize + iNumTetCells() + iNumPyrCells() + iNumPrismCells());
    // Reset the number of hex cells in use
    VAHex.vSetup(iNewSize);
  } // Hexahedron iteration

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

void VolMesh::vPurgeBdryFaces()
  // 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
  {
    TriBFace *pTBFFor  = VATriBF.pTFirst();
    TriBFace *pTBFBack = VATriBF.pTLast ();

    int iNewSize = 0;

    while (pTBFFor != pTBFBack) {
      while (!pTBFFor->qDeleted() && pTBFFor != pTBFBack){
	pTBFFor = VATriBF.pTNext(pTBFFor);
	iNewSize++;
      }

      while (pTBFBack->qDeleted() && pTBFFor != pTBFBack)
	pTBFBack = VATriBF.pTPrev(pTBFBack);

      if (pTBFFor != pTBFBack) {
	assert(pTBFFor->qDeleted());
	iNewSize++;
	(*pTBFFor) = (*pTBFBack);
	pTBFBack->vMarkDeleted();
#ifndef NDEBUG
	for (int i = 0; i < pTBFFor->iNumFaces(); i++)
	  assert(pTBFFor->pFFace(i)->iFullCheck());
#endif
	pTBFFor = VATriBF.pTNext(pTBFFor);
      }
      else if (!pTBFFor->qDeleted())
	iNewSize++;
    }
    // Reset the number of tri bdry faces in use
    VATriBF.vSetup(iNewSize);
  } // Triangular boundary face iteration
  {
    IntTriBFace *pITBFFor  = VAIntTriBF.pTFirst();
    IntTriBFace *pITBFBack = VAIntTriBF.pTLast ();

    int iNewSize = 0;

    while (pITBFFor != pITBFBack) {
      while (!pITBFFor->qDeleted() && pITBFFor != pITBFBack){
	pITBFFor = VAIntTriBF.pTNext(pITBFFor);
	iNewSize++;
      }

      while (pITBFBack->qDeleted() && pITBFFor != pITBFBack)
	pITBFBack = VAIntTriBF.pTPrev(pITBFBack);

      if (pITBFFor != pITBFBack) {
	assert(pITBFFor->qDeleted());
	iNewSize++;
	(*pITBFFor) = (*pITBFBack);
	pITBFBack->vMarkDeleted();
#ifndef NDEBUG
	for (int i = 0; i < pITBFFor->iNumFaces(); i++)
	  assert(pITBFFor->pFFace(i)->iFullCheck());
#endif
	pITBFFor = VAIntTriBF.pTNext(pITBFFor);
      }
      else if (!pITBFFor->qDeleted())
	iNewSize++;
    }
    // Reset the number of internal tri bdry faces in use
    VAIntTriBF.vSetup(iNewSize);
  } // Internal triangular boundary face iteration
  {
    QuadBFace *pQBFFor  = VAQuadBF.pTFirst();
    QuadBFace *pQBFBack = VAQuadBF.pTLast ();

    int iNewSize = 0;

    while (pQBFFor != pQBFBack) {
      while (!pQBFFor->qDeleted() && pQBFFor != pQBFBack){
	pQBFFor = VAQuadBF.pTNext(pQBFFor);
	iNewSize++;
      }

      while (pQBFBack->qDeleted() && pQBFFor != pQBFBack)
	pQBFBack = VAQuadBF.pTPrev(pQBFBack);

      if (pQBFFor != pQBFBack) {
	assert(pQBFFor->qDeleted());
	iNewSize++;
	(*pQBFFor) = (*pQBFBack);
	pQBFBack->vMarkDeleted();
#ifndef NDEBUG
	for (int i = 0; i < pQBFFor->iNumFaces(); i++)
	  assert(pQBFFor->pFFace(i)->iFullCheck());
#endif
	pQBFFor = VAQuadBF.pTNext(pQBFFor);
      }
      else if (!pQBFFor->qDeleted())
	iNewSize++;
    }
    // Reset the number of quad bdry faces in use
    VAQuadBF.vSetup(iNewSize);
  } // Quadangular boundary face iteration
  {
    IntQuadBFace *pITBFFor  = VAIntQuadBF.pTFirst();
    IntQuadBFace *pITBFBack = VAIntQuadBF.pTLast ();

    int iNewSize = 0;

    while (pITBFFor != pITBFBack) {
      while (!pITBFFor->qDeleted() && pITBFFor != pITBFBack){
	pITBFFor = VAIntQuadBF.pTNext(pITBFFor);
	iNewSize++;
      }

      while (pITBFBack->qDeleted() && pITBFFor != pITBFBack)
	pITBFBack = VAIntQuadBF.pTPrev(pITBFBack);

      if (pITBFFor != pITBFBack) {
	assert(pITBFFor->qDeleted());
	iNewSize++;
	(*pITBFFor) = (*pITBFBack);
	pITBFBack->vMarkDeleted();
#ifndef NDEBUG
	for (int i = 0; i < pITBFFor->iNumFaces(); i++)
	  assert(pITBFFor->pFFace(i)->iFullCheck());
#endif
	pITBFFor = VAIntQuadBF.pTNext(pITBFFor);
      }
      else if (!pITBFFor->qDeleted())
	iNewSize++;
    }
    // Reset the number of internal quad bdry faces in use
    VAIntQuadBF.vSetup(iNewSize);
  } // Internal quadrilateral boundary face iteration
  if (qLocalLog) {
    SUMAA_LOG_EVENT_END(PURGE);
  }
}

void VolMesh::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;
  bool qAnyChanges = false;

  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->pFHintFace()->qValid()) {
      pV->vMarkDeleted();
    }
    if (pV->qDeleted()) {// If vert is deleted, copy over it.
      aiVertInd[iV] = -1;
      qAnyChanges = true;
    }
    else {
      aiVertInd[iV] = iVActive;
      Vert *pVA = pVVert(iVActive);
      // Copy the data from its old to its new location
      (*pVA) = (*pV);
      iVActive++;
    }
  }

  // If there are no verts deleted, then skip all the hard work.
  if (!qAnyChanges) {
    delete [] aiVertInd;
    if (qLocalLog) {
      SUMAA_LOG_EVENT_END(PURGE);
    }
    return;
  }

  // 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->qDeleted()) 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]);
  }
  delete [] aiVertInd;

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

void VolMesh::vCheckForITAPSBdryErrors()
{
  // This routine is a hack to make sure that all faces on the bdry have
  // a BFace attached to them.  Any faces that have a missing cell are
  // automatically assigned a BFace with a default BC instead.  If
  // called at the wrong time, this routine could easily break a mesh in
  // a way that ITAPS can't fix, so it should only be called from vPurge,
  // which in turn is called from the ITAPS implementation only when
  // allowing entity handles to change.
  for (int iF = 0; iF < iNumFaces(); iF++) {
    Face *pF = pFFace(iF);
    if (pF->qDeleted()) continue;
    Cell *pCL = pF->pCCellLeft();
    Cell *pCR = pF->pCCellRight();
    if (!pCL->qValid()) {
      // Add a BFace on this side
      assert(pCR->qValid());
      BFace *pBF = pBFNewBFace(pF->iNumVerts());
      pBF->vAddFace(pF);
      pF->vAddCell(pBF);
      assert(pBF == pF->pCCellLeft());
    }
    else if (!pCR->qValid()) {
      BFace *pBF = pBFNewBFace(pF->iNumVerts());
      pBF->vAddFace(pF);
      pF->vAddCell(pBF);
      assert(pBF == pF->pCCellRight());
    }
  }
}

void VolMesh::vPurge()
{
  SUMAA_LOG_EVENT_BEGIN(PURGE);
#ifdef ITAPS
  vCheckForITAPSBdryErrors();
#endif
  vPurgeCells();
  vPurgeFaces();
  vPurgeBdryFaces();
  vPurgeVerts();
  SUMAA_LOG_EVENT_END(PURGE);
}

//@ Validity checking of entire mesh
bool VolMesh::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()) {
      Face *pF = pV->pFHintFace();
      if (pF->qValid())
	qRetVal = qRetVal && pF->qHasVert(pV);
    }
  }

  for (i = iNumFaces() - 1; i >= 0 && qRetVal; i--) {
    Face *pF = pFFace(i);
    if (pF->qDeleted()) 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);
    if (!pC->qValid() || pC->qDeleted()) continue;
    int iFull = pC->iFullCheck();
    bool qRegOK = (pC->iRegion() != iInvalidRegion);
    qRetVal = qRetVal && (iFull && qRegOK);
    // Check that the cell is actually closed
    qRetVal = qRetVal && pC->qIsClosed();
    // Provided that the cell actually exists, be sure that it has a
    // positive volume.
    // Test temporarily disabled to confirm topo validity without geom check.
//     if (iFull == 1 && pC->eType() == CellSkel::eTet)
//       qRetVal = qRetVal && (pC->dSize() > 0);
  }

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

double VolMesh::dBoundarySize() const
{
  double dRetVal = 0;
  for (int i = iNumBdryFaces() - 1; i >= 0; i--)
    dRetVal += pBFBFace(i)->dSize();
  return dRetVal;
}

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

bool VolMesh::qWatertight() const
{
  double adSum[] = {0., 0., 0.};
  double adTemp[3];
  for (int i = iNumBdryFaces() - 1; i >= 0; i--) {
    BFace *pBF = pBFBFace(i);
    if (pBF->iFullCheck()) {
      pBF->vVecSize(adTemp);
      adSum[0] += adTemp[0];
      adSum[1] += adTemp[1];
      adSum[2] += adTemp[2];
    }
  }
  double dCompSize = dBoundarySize();
  double dInvCompSize = 0.001/dCompSize;
  vSCALE3D(adSum, dInvCompSize);

  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:\n");
    vMessage(2, "  x: %12.5g\n  y: %12.5g\n  z: %12.5g\n", adSum[0],
	     adSum[1], adSum[2]);
    return false;
  }
}

bool VolMesh::qMeshOrientationOK(const bool qFix)
  // Check and correct some bad orientation situations.  Specifically,
  // if a mesh is specified with all tets left-handed, this will be
  // detected and corrected.  More general cases, including isolated
  // inverted tetrahedra in the interior of the mesh, are not yet
  // handled.  Also, meshes with degenerate tetrahedra are not fixed.
{
  // The status of the two tets on either side of a face is summarized
  // by the formula 5*(left-status) + right-status, where *-status is -1
  // for reverse-oriented tets, 0 for flat tets, 1 for
  // correctly-oriented tets, 2 for boundaries, and 3 for non-tets.  The
  // following enum gives symbolic values to the results.
  enum {eInvInv   = -6, eInvFlat  = -5, eInvGood  = -4, eInvBdry  = -3,
	eInvNontet = -2,
	eFlatInv  = -1, eFlatFlat =  0, eFlatGood =  1, eFlatBdry =  2,
	eFlatNontet = 3,
	eGoodInv  =  4, eGoodFlat =  5, eGoodGood =  6, eGoodBdry =  7,
	eGoodNontet = 8,
	eBdryInv  =  9, eBdryFlat =  10, eBdryGood =  11, eBdryBdry = 12,
	eBdryNontet = 13,
	eNontetInv  =  14, eNontetFlat =  15, eNontetGood =  16,
	eNontetBdry = 17, eNontetNontet = 18};
  int i, iGood, iBad, iFixed, iUnfixed;
  // Check is performed on a face-by-face basis.  Performing this check
  // on a cell-by-cell basis for tets is not necessarily thorough
  // enough, as some faces will never be used to check an orientation.
  // A more general solution to this problem may have to be found for
  // non-tetrahedral cells.
  iGood = iBad = iFixed = iUnfixed = 0;
  for (i = 0; i < iNumFaces(); i++) {
    Face *pF = pFFace(i);
    Cell *pCL = pF->pCCellLeft();
    Cell *pCR = pF->pCCellRight();
    Vert *pV0 = pF->pVVert(0);
    Vert *pV1 = pF->pVVert(1);
    Vert *pV2 = pF->pVVert(2);
    Vert *pVL = pVInvalidVert;
    Vert *pVR = pVInvalidVert;
    int iStatL, iStatR, iStatus;

    // Check orientation on the left.
    switch (pCL->eType()) {
    case Cell::eTet:
      pVL = dynamic_cast<TetCell*>(pCL)->pVVertOpposite(pF);
      iStatL = iOrient3D(pV0, pV2, pV1, pVL);
      break;
    case Cell::eTriBFace:
    case Cell::eQuadBFace:
    case Cell::eIntTriBFace:
    case Cell::eIntQuadBFace:
      iStatL = 2;
      break;
    default:
      iStatL = 3;
    }

    // Check orientation on the right.
    switch (pCR->eType()) {
    case Cell::eTet:
      pVR = dynamic_cast<TetCell*>(pCR)->pVVertOpposite(pF);
      iStatR = iOrient3D(pV0, pV1, pV2, pVR);
      break;
    case Cell::eTriBFace:
    case Cell::eQuadBFace:
    case Cell::eIntTriBFace:
    case Cell::eIntQuadBFace:
      iStatR = 2;
      break;
    default:
      iStatR = 3;
    }

    assert(!pF->qHasVert(pVR));
    assert(!pF->qHasVert(pVL));

    iStatus = iStatL*5 + iStatR;

    switch (iStatus) {
    case eGoodGood:
    case eBdryGood:
    case eGoodBdry:
    case eGoodNontet:
    case eNontetGood:
    case eNontetBdry:
    case eBdryNontet:
    case eNontetNontet:
      iGood++;
      break;
    case eInvInv:
    case eInvBdry:
    case eBdryInv:
    case eInvNontet:
    case eNontetInv:
      iBad++;
      if (qFix) {
	pF->vInterchangeCells();
	assert(!pF->pCCellRight()->qHasVert(pVR));
	assert(!pF->pCCellLeft()->qHasVert(pVL));
	iFixed++;
      }
      else iUnfixed++;
      break;
    case eFlatInv:
      if (qFix) {
	pF->vInterchangeCells();
	Vert *pVTemp = pVL;
	pVL = pVR;
	pVR = pVTemp;
	assert(pF->pCCellRight()->qHasVert(pVR));
	assert(pF->pCCellLeft()->qHasVert(pVL));
      }
      // This case has now been converted, so fall through.
    case eFlatGood:
      iUnfixed++;
      iBad++;
      break;
    case eInvFlat:
      if (qFix) {
	pF->vInterchangeCells();
	Vert *pVTemp = pVL;
	pVL = pVR;
	pVR = pVTemp;
	assert(pF->pCCellRight()->qHasVert(pVR));
	assert(pF->pCCellLeft()->qHasVert(pVL));
      }
      // This case has now been converted, so fall through.
    case eGoodFlat:
      iUnfixed++;
      iBad++;
      break;
    case eFlatFlat:
    case eFlatBdry:
    case eBdryFlat:
    case eNontetFlat:
    case eFlatNontet:
      iUnfixed++;
      iBad++;
      break;
    case eInvGood:
      if (qFix) {
	pF->vInterchangeCells();
	Vert *pVTemp = pVL;
	pVL = pVR;
	pVR = pVTemp;
	assert(pF->pCCellRight()->qHasVert(pVR));
	assert(pF->pCCellLeft()->qHasVert(pVL));
      }
      // This case has now been converted, so fall through.
    case eGoodInv:
      iUnfixed++;
      iBad++;
      break;
    case eBdryBdry:
      // This is a really weird case that should never come up.
      // The only fix is to remove the face and the corresponding bdry
      // faces.
      iUnfixed++;
      iBad++;
      break;
    default:
      // Should never be able to get here at all.
      assert(0);
      break;
    }
  }
  assert(iGood + iBad == iNumFaces());
  assert(iBad == iFixed + iUnfixed);
  vMessage(2, "Checked %d faces in mesh; %s\n",
	   iNumFaces(), "bounding tets classified as follows:\n");
  vMessage(2, "  Correct for a right-handed mesh:   %9d\n", iGood);
  vMessage(2, "  Incorrect for a right-handed mesh: %9d\n", iBad);
  if (iBad != 0) {
    if (qFix) {
      vMessage(2, "    Of these, %5d were corrected.\n", iFixed);
      vMessage(2, "              %5d could not be corrected.\n", iUnfixed);
    }
    else {
      vMessage(2, "No attempt was made to correct these errors.\n");
    }
  }
  if (iUnfixed != 0)
    vFatalError("Mesh can not be easily made into a right-handed mesh",
		"processing volume mesh data file");
  else {
    // Make sure nothing got screwed up.
    vMessage(1, "Volume mesh is properly oriented.\n");
  }
  return (qValid());
}

void VolMesh::vClassifyTets() const
     // Classify tets according to the scheme of Bern, Chew, Eppstein,
     // and Ruppert.
{
  static const double dLilDihed = 20., dBigDihed = 150.;
  static const double dLilSolid =  3., dBigSolid = 240.;
  typedef enum {eNonTet = 1, eRound = 2, eNeedle = 4, eWedge = 8,
		eSpindle = 16, eSliver = 32, eCap = 64} eTetType;

  SUMAA_LOG_EVENT_BEGIN(CLASSIFY_TETS);
  int iRound = 0, iNeedle = 0, iWedge = 0, iSpindle = 0, iSliver = 0,
    iCap = 0;

  int i, *aiFaceType = new int[iNumFaces()];
  for (i = 0; i < iNumFaces(); i++) aiFaceType[i] = 0;

  vMessage(2, "\nClassifying tets according to which pathologies they have...\n");
  for (int iC = iNumCells() - 1; iC >= 0; iC--) {
    Cell *pC = pCCell(iC);
    // Don't do anything for cells which have been deleted or are
    // otherwise bogus.
    if (pC->iFullCheck() != 1) continue;
    if (pC->eType() == Cell::eTet) {
      double adDihed[6];
      int iNDihed;
      dynamic_cast<TetCell*>(pC)->vAllDihed(adDihed, &iNDihed);

      int iSmall = 0, iLarge = 0;
      for (i = 0; i < iNDihed; i++)
	if      (adDihed[i] < dLilDihed) iSmall++;
	else if (adDihed[i] > dBigDihed) iLarge++;

      double dS1 = adDihed[0] + adDihed[1] + adDihed[2] - 180;
      double dS2 = adDihed[0] + adDihed[3] + adDihed[4] - 180;
      double dS3 = adDihed[1] + adDihed[3] + adDihed[5] - 180;
      double dS4 = adDihed[2] + adDihed[4] + adDihed[5] - 180;
      int iSmallSolid = ( (dS1 < dLilSolid ? 1 : 0) +
			  (dS2 < dLilSolid ? 1 : 0) +
			  (dS3 < dLilSolid ? 1 : 0) +
			  (dS4 < dLilSolid ? 1 : 0) );
      int iLargeSolid = ( (dS1 > dBigSolid ? 1 : 0) +
			  (dS2 > dBigSolid ? 1 : 0) +
			  (dS3 > dBigSolid ? 1 : 0) +
			  (dS4 > dBigSolid ? 1 : 0) );
      assert(iLarge <= 3);
      assert(iSmall <= 4);
      assert(iLargeSolid <= 1);

      int iThisType;
      if      (iLargeSolid == 1)         {iThisType = eCap; iCap ++;}
      else if (iLarge > 0 && iSmall > 0) {iThisType = eSliver; iSliver ++;}
      else if (iLarge > 0)               {iThisType = eSpindle; iSpindle ++;}
      else if (iSmall > 0)               {iThisType = eWedge; iWedge ++;}
      else if (iSmallSolid > 0)          {iThisType = eNeedle; iNeedle ++;}
      else                               {iThisType = eRound; iRound ++;}

      double adRes[1];
      int iNumRes;
      vShortestToRadius(pC, &iNumRes, adRes);
      assert(iNumRes == 1);
      vMessage(3, "[DETAIL] Tet: %7d  Type: %1d  B: %.5f\n",
	       iC, iThisType, adRes[0]);

      for (int ii = 0; ii < 4; ii++) {
	int iF = iFaceIndex(pC->pFFace(ii));
	aiFaceType[iF] += iThisType;
      }
    } // if (tet)
    else {
      for (int ii = 0; ii < pC->iNumFaces(); ii++) {
	int iF = iFaceIndex(pC->pFFace(ii));
	aiFaceType[iF] += eNonTet;
      }
    }
  } // End of loop over cells

  // Loop over the boundary faces and make due note of them
  for (i = 0; i < iNumBdryFaces(); i++) {
    BFace *pBF = pBFBFace(i);
    if (!pBF->qDeleted()) {
      int iF = iFaceIndex(pBF->pFFace());
      aiFaceType[iF] += eNonTet;
    }
  }

  int iTotal = iRound + iNeedle + iWedge + iSpindle + iSliver + iCap;
  vMessage(2, "Analyzed a total of %d tetrahedra.\n", iTotal);
  vMessage(2, "Number of well-shaped tets:%8d (%5.2f%%)\n", iRound,
	   100*double(iRound)/double(iTotal));
  vMessage(2, "                   needles:%8d (%5.2f%%)\n", iNeedle,
	   100*double(iNeedle)/double(iTotal));
  vMessage(2, "                    wedges:%8d (%5.2f%%)\n", iWedge,
	   100*double(iWedge)/double(iTotal));
  vMessage(2, "                  spindles:%8d (%5.2f%%)\n", iSpindle,
	   100*double(iSpindle)/double(iTotal));
  vMessage(2, "                   slivers:%8d (%5.2f%%)\n", iSliver,
	   100*double(iSliver)/double(iTotal));
  vMessage(2, "                      caps:%8d (%5.2f%%)\n", iCap,
	   100*double(iCap)/double(iTotal));

  // Now tabulate, for each face, the type of its bounding cells
  int a2iFaces[7][7];
  for (i = 0; i < 7; i++)
    for (int ii = 0; ii < 7; ii++)
      a2iFaces[i][ii] = 0;

  for (i = 0; i < iNumFaces(); i++) {
    if (pFFace(i)->qDeleted()) continue;
    switch (aiFaceType[i]) {
    case eNonTet + eNonTet: // Neither cell is a tet of any kind (bdry,
			    // pyr, prism, hex)
      a2iFaces[0][0] += 1;
      break;
    case eNonTet + eRound: // One round, one not a tet
      a2iFaces[0][1] += 1;
      break;
    case eNonTet + eNeedle: // One needle, one not a tet
      a2iFaces[0][2] += 1;
      break;
    case eNonTet + eWedge: // One wedge, one not a tet
      a2iFaces[0][3] += 1;
      break;
    case eNonTet + eSpindle: // One spindle, one not a tet
      a2iFaces[0][4] += 1;
      break;
    case eNonTet + eSliver: // One sliver, one not a tet
      a2iFaces[0][5] += 1;
      break;
    case eNonTet + eCap: // One cap, one not a tet
      a2iFaces[0][6] += 1;
      break;
    case eRound + eRound: // Both round
      a2iFaces[1][1] += 1;
      break;
    case eRound + eNeedle: // One needle, one round
      a2iFaces[1][2] += 1;
      break;
    case eRound + eWedge: // One wedge, one round
      a2iFaces[1][3] += 1;
      break;
    case eRound + eSpindle: // One spindle, one round
      a2iFaces[1][4] += 1;
      break;
    case eRound + eSliver: // One sliver, one round
      a2iFaces[1][5] += 1;
      break;
    case eRound + eCap: // One cap, one round
      a2iFaces[1][6] += 1;
      break;
    case eNeedle + eNeedle: // Both needles
      a2iFaces[2][2] += 1;
      break;
    case eNeedle + eWedge: // One wedge, one needle
      a2iFaces[2][3] += 1;
      break;
    case eNeedle + eSpindle: // One spindle, one needle
      a2iFaces[2][4] += 1;
      break;
    case eNeedle + eSliver: // One sliver, one needle
      a2iFaces[2][5] += 1;
      break;
    case eNeedle + eCap: // One cap, one needle
      a2iFaces[2][6] += 1;
      break;
    case eWedge + eWedge: // Both wedges
      a2iFaces[3][3] += 1;
      break;
    case eWedge + eSpindle: // One spindle, one wedge
      a2iFaces[3][4] += 1;
      break;
    case eWedge + eSliver: // One sliver, one wedge
      a2iFaces[3][5] += 1;
      break;
    case eWedge + eCap: // One cap, one wedge
      a2iFaces[3][6] += 1;
      break;
    case eSpindle + eSpindle: // Both spindles
      a2iFaces[4][4] += 1;
      break;
    case eSpindle + eSliver: // One sliver, one spindle
      a2iFaces[4][5] += 1;
      break;
    case eSpindle + eCap: // One cap, one spindle
      a2iFaces[4][6] += 1;
      break;
    case eSliver + eSliver: // Both slivers
      a2iFaces[5][5] += 1;
      break;
    case eSliver + eCap: // One cap, one sliver
      a2iFaces[5][6] += 1;
      break;
    case eCap + eCap: // Both caps
      a2iFaces[6][6] += 1;
      break;
    default: // Should never get here
      assert(0);
      break;
    }
  }

  vMessage(2, "\nBreakdown of faces by the types of cells on either side\n\n");
  vMessage(2, "         Non-tet   Round  Needle   Wedge Spindle  Sliver     Cap\n");
  vMessage(2, "Non-tet %8d %7d %7d %7d %7d %7d %7d\n",
	   a2iFaces[0][0], a2iFaces[0][1], a2iFaces[0][2],
	   a2iFaces[0][3], a2iFaces[0][4], a2iFaces[0][5],
	   a2iFaces[0][6]);
  vMessage(2, "Round   %16d %7d %7d %7d %7d %7d\n",
	   a2iFaces[1][1], a2iFaces[1][2], a2iFaces[1][3],
	   a2iFaces[1][4], a2iFaces[1][5], a2iFaces[1][6]);
  vMessage(2, "Needle  %24d %7d %7d %7d %7d\n",
	   a2iFaces[2][2], a2iFaces[2][3], a2iFaces[2][4],
	   a2iFaces[2][5], a2iFaces[2][6]);
  vMessage(2, "Wedge   %32d %7d %7d %7d\n",
	   a2iFaces[3][3], a2iFaces[3][4], a2iFaces[3][5],
	   a2iFaces[3][6]);
  vMessage(2, "Spindle %40d %7d %7d\n",
	   a2iFaces[4][4], a2iFaces[4][5], a2iFaces[4][6]);
  vMessage(2, "Sliver  %48d %7d\n",
	   a2iFaces[5][5], a2iFaces[5][6]);
  vMessage(2, "Cap     %56d\n\n", a2iFaces[6][6]);

  delete [] aiFaceType;
  SUMAA_LOG_EVENT_END(CLASSIFY_TETS);
}

void
VolMesh::vMeshInBrick(const double dXMin, const double dXMax,
		      const double dYMin, const double dYMax,
		      const double dZMin, const double dZMax)
{
  // x runs left to right
  // y runs bottom to top
  // z runs back to front
  Vert *pVLoLeftBack = pVNewVert();
  Vert *pVLoLeftFront = pVNewVert();
  Vert *pVHiLeftFront = pVNewVert();
  Vert *pVHiLeftBack = pVNewVert();
  Vert *pVLoRightBack = pVNewVert();
  Vert *pVLoRightFront = pVNewVert();
  Vert *pVHiRightFront = pVNewVert();
  Vert *pVHiRightBack = pVNewVert();

  double adLoc[3];
  adLoc[0] = dXMin;
  adLoc[1] = dYMin;
  adLoc[2] = dZMin;
  pVLoLeftBack->vSetCoords(3, adLoc);
  adLoc[2] = dZMax;
  pVLoLeftFront->vSetCoords(3, adLoc);
  adLoc[1] = dYMax;
  pVHiLeftFront->vSetCoords(3, adLoc);
  adLoc[2] = dZMin;
  pVHiLeftBack->vSetCoords(3, adLoc);

  adLoc[0] = dXMax;
  adLoc[1] = dYMin;
  adLoc[2] = dZMin;
  pVLoRightBack->vSetCoords(3, adLoc);
  adLoc[2] = dZMax;
  pVLoRightFront->vSetCoords(3, adLoc);
  adLoc[1] = dYMax;
  pVHiRightFront->vSetCoords(3, adLoc);
  adLoc[2] = dZMin;
  pVHiRightBack->vSetCoords(3, adLoc);

  // Set up initial mesh, with five tets filling the big hex.

  vSetupTriFaces(16);
  vSetupTriBFaces(12);
  vSetupTetCells(5);

  // Begin with the faces on the outside of the hex
  // Those that meet at the LoLeftFront corner
  pFFace(0)->vAssign(pCCell(0), pBFBFace(0), pVLoLeftFront,
		     pVLoRightFront, pVHiLeftFront);
  pFFace(1)->vAssign(pCCell(0), pBFBFace(1), pVLoLeftFront,
		     pVLoLeftBack, pVLoRightFront);
  pFFace(2)->vAssign(pCCell(0), pBFBFace(2), pVLoLeftFront,
		     pVHiLeftFront, pVLoLeftBack);

  // Those that meet at the HiRightFront corner
  pFFace(3)->vAssign(pCCell(1), pBFBFace(3), pVHiRightFront,
		     pVHiLeftFront, pVLoRightFront);
  pFFace(4)->vAssign(pCCell(1), pBFBFace(4), pVHiRightFront,
		     pVLoRightFront, pVHiRightBack);
  pFFace(5)->vAssign(pCCell(1), pBFBFace(5), pVHiRightFront,
		     pVHiRightBack, pVHiLeftFront);

  // Those that meet at the HiLeftBack corner
  pFFace(6)->vAssign(pCCell(2), pBFBFace(6), pVHiLeftBack,
		     pVLoLeftBack, pVHiLeftFront);
  pFFace(7)->vAssign(pCCell(2), pBFBFace(7), pVHiLeftBack,
		     pVHiLeftFront, pVHiRightBack);
  pFFace(8)->vAssign(pCCell(2), pBFBFace(8), pVHiLeftBack,
		     pVHiRightBack, pVLoLeftBack);

  // Those that meet at the LoRightBack corner
  pFFace(9)->vAssign(pCCell(3), pBFBFace(9), pVLoRightBack,
		     pVHiRightBack, pVLoRightFront);
  pFFace(10)->vAssign(pCCell(3), pBFBFace(10), pVLoRightBack,
		      pVLoLeftBack, pVHiRightBack);
  pFFace(11)->vAssign(pCCell(3), pBFBFace(11), pVLoRightBack,
		      pVLoRightFront, pVLoLeftBack);

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

  // Now the four internal faces.
  pFFace(12)->vAssign(pCCell(0), pCCell(4), pVLoLeftBack,
		      pVHiLeftFront, pVLoRightFront);
  pFFace(13)->vAssign(pCCell(1), pCCell(4), pVLoRightFront,
		      pVHiLeftFront, pVHiRightBack);
  pFFace(14)->vAssign(pCCell(2), pCCell(4), pVHiLeftFront,
		      pVLoLeftBack, pVHiRightBack);
  pFFace(15)->vAssign(pCCell(3), pCCell(4), pVLoRightFront,
		      pVHiRightBack, pVLoLeftBack);

  // Set up the bdry faces
  for (int iBF = 0; iBF < 12; iBF++)
    pBFBFace(iBF)->vAssign(pFFace(iBF));

  // Set up the cells
  pCCell(0)->vAssign(pFFace( 0), pFFace( 1), pFFace( 2), pFFace(12));
  pCCell(1)->vAssign(pFFace( 3), pFFace( 4), pFFace( 5), pFFace(13));
  pCCell(2)->vAssign(pFFace( 6), pFFace( 7), pFFace( 8), pFFace(14));
  pCCell(3)->vAssign(pFFace( 9), pFFace(10), pFFace(11), pFFace(15));
  pCCell(4)->vAssign(pFFace(12), pFFace(13), pFFace(14), pFFace(15));

  vSetAllHintFaces();
  assert(qValid());
}


VolMesh::VolMesh(const VAVertex& VA, const int iQualMeas,
		 const double dMaxAngle)
  : Mesh(VA), qOKToChangeSurfaceEdges(true),
    qDoEdgeSwapping(true), qStrictPatchChecking(true),
    dMaxAngleForSurfSwap(dMaxAngle), VATriF(),
    VAQuadF(), VATriBF(), VAQuadBF(), VATet(), VAPyr(), VAPrism(),
    VAHex(), VAIntTriBF(), VAIntQuadBF(), pB3D(), qBdryFromMesh(true)
  // Create a volume mesh of a point cloud, placed inside a big box.
{
  SUMAA_LOG_EVENT_BEGIN(INIT_TETRA);
  vMessage(0, "Generating volume mesh from point cloud.\n");
  qBdryFill = false;
  pQ = new Quality(this, iQualMeas);
  SMinitSmoothing(3, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);

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

  int iVert, iNVerts = VA.iSize();
  for (iVert = 0; iVert < iNVerts; iVert++) {
    double dX = VA[iVert].dX();
    double dY = VA[iVert].dY();
    double dZ = VA[iVert].dZ();
    dXMin = min(dXMin, dX);
    dXMax = max(dXMax, dX);
    dYMin = min(dYMin, dY);
    dYMax = max(dYMax, dY);
    dZMin = min(dZMin, dZ);
    dZMax = max(dZMax, dZ);
  }

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

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

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

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

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

  vMeshInBrick(dXLo, dXHi, dYLo, dYHi, dZLo, dZHi);

  //@@ Insert all points from the point set
  // Currently these are inserted lexicographically and swapping is
  // performed after each insertion.  Don't insert the verts forming the
  // big box.
  vMessage(1, "Adding remaining vertices to the mesh...\n");
  int iSwaps = 0;
  int iTenth = 1 + iNumVerts() / 10;
  if (iTenth < 25) iTenth = 25;
  SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
  int iV;
  for (iV = 0; iV < iNumVerts() - 8; iV++) {
    // Don't insert vertices that are already in the mesh
    Vert *pVInsert = pVVert(iV);

    int iC = iNumCells() - 1;
    Cell *pCGuess;
    while ((pCGuess = pCCell(iC))->iFullCheck() != 1) iC --;
    int iSwapsThisTime = 0;

    // Allow swapping after insertion; force insertion to take place
    // without snap-to-edge/face.  NOTE:  This vertex already exists in
    // the point set.
    qInsertPoint(pVInsert->adCoords(), pCGuess, &iSwapsThisTime,
		 1, 1, pVInsert);
    iSwaps += iSwapsThisTime;
    if ((iV+1) % iTenth == 0)
      vMessage(1, "%5d %s  %6d swaps performed.\n",
	       iV, "points inserted from the point cloud.", iSwaps);
  }

  vMessage(1, "%5d %s  %6d swaps performed.\n",
	   iV, "points inserted from the point cloud.", iSwaps);
  SUMAA_LOG_EVENT_END(INSERT_VERTEX);

  vMessage(1, "All verts added; improving connectivity.\n");

  // That's it.  The entire mesh is relevant and no exterior
  // cleanup is required.

  vAllowSwapRecursion();
  vAllowSurfaceEdgeChanges();
  vAllowEdgeSwapping();
  vSetSwapType(eDelaunay);
  iSwap();

  if (!qValid())
    vFoundBug("creation of constrained 3D mesh from point cloud");
  vMessage(1, "Done creating constrained volume mesh.\n");
  SUMAA_LOG_EVENT_END(INIT_TETRA);
}

int VolMesh::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;
  std::set<BFace*>::iterator iterBF;
  for (iterBF = spBF.begin(); iterBF != spBF.end(); iterBF++) {
    TriBFaceBase* pTBF = (dynamic_cast<TriBFaceBase*>(*iterBF));
    if (pTBF->eIsEncroached(eEncType) != eClean) {
      if (IQ.qAddEntry(InsertionQueueEntry(pTBF))) {
	iAdded ++;
      }
    }
  }

  // Now check all the bdry subsegs bounded by the BFaces given.  The
  // assumption is that bdry subsegs that are encroached will be
  // encroached by points that they share a tet with them.  This
  // assumption is certainly valid for CDT's, but may not be generally
  // valid.
  std::set<FaceSort> sFSEncSubseg;
  vCheckBdryFragmentForSubsegEncroachment(spBF, sFSEncSubseg);
  std::set<FaceSort>::iterator iter;
  for (iter = sFSEncSubseg.begin(); iter != sFSEncSubseg.end();
       iter++) {
    FaceSort FS = *iter;
    if (IQ.qAddEntry(InsertionQueueEntry(FS.pVVert(0), FS.pVVert(1)))) {
      iAdded ++;
    }
  }
  vMessage(4, "Added a total of %d encroached bdry entities to the queue.\n",
	   iAdded);
  return iAdded;
}

int VolMesh::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;
  std::set<BFace*>::iterator iterBF;
  for (iterBF = spBF.begin(); iterBF != spBF.end(); iterBF++) {
    TriBFaceBase* pTBF = (dynamic_cast<TriBFaceBase*>(*iterBF));
    if (pTBF->eIsPointEncroaching(adLoc) != eClean) {
      if (IQ.qAddEntry(InsertionQueueEntry(pTBF), dOffset)) {
	iAdded ++;
      }
    }
  }

  // Now check all the bdry subsegs bounded by the BFaces given.  The
  // assumption is that bdry subsegs that are encroached will be
  // encroached by points that they share a tet with them.  This
  // assumption is certainly valid for CDT's, but may not be generally
  // valid.
  std::set<FaceSort> sFSEncSubseg;
  vCheckBdryFragmentForSubsegEncroachment(spBF, sFSEncSubseg, adLoc);
  std::set<FaceSort>::iterator iter;
  for (iter = sFSEncSubseg.begin(); iter != sFSEncSubseg.end(); iter++) {
    FaceSort FS = *iter;
    if (IQ.qAddEntry(InsertionQueueEntry(FS.pVVert(0), FS.pVVert(1)), dOffset))
      {
	iAdded ++;
      }
  }
  vMessage(3, "Added a total of %d encroached bdry entities to the queue.\n",
	   iAdded);
  return iAdded;
}

void VolMesh::vProtectSmallAngles()
{
  // For every vertex tagged as a small angle vertex (== has incident
  // segments with a small angle between them), cut all incident edges
  // at a distance of about the length scale at the vertex.
  for (int iV = iNumVerts() - 1; iV >= 0; iV --) {
    Vert *pV = pVVert(iV);
    if (pV->qIsSmallAngleCorner() && pV->iVertType() == Vert::eBdryApex) {
      // Now we find the neighbors that lie on bdry segments and split
      // edges connecting pV to those neighbors.

      // The usual call to vNeighborhood.
      std::set<Vert*> spVNeigh;
      std::set<BFace*> spBFInc;
      {
	bool qBdryVert;
	std::set<Cell*> spCTmp;
	vNeighborhood(pV, spCTmp, spVNeigh, &spBFInc, &qBdryVert);
      }

      // Given the bdry faces incident on pV, find all boundary segments
      // incident on it.

      vMessage(4, "  Finding all bdry subsegs incident on vertex %d...\n",
	       iV);
      std::set<Vert*>::iterator iterV;
      for (iterV = spVNeigh.begin(); iterV != spVNeigh.end(); iterV++) {
	bool qIsBdrySeg = false;
	Vert *pVNeigh = *iterV;
	// Now check whether these two are connected by a bdry subseg.
	if (pVNeigh->iVertType() == Vert::eBdryApex ||
	    pVNeigh->iVertType() == Vert::eBdryCurve) {
	  vMessage(4, "  Checking for bdry vertex %d\n",
		   iVertIndex(pVNeigh));
	  // pVNeigh is definitely on -some- bdry subseg.  Now check
	  // whether the two share (at least) two bdry faces.
	  std::set<BFace*> spBFCommon;
	  std::set<BFace*>::iterator iterBF;
	  for (iterBF = spBFInc.begin(); iterBF != spBFInc.end(); iterBF++) {
	    BFace *pBF = *iterBF;
	    assert(pBF->qHasVert(pV));
	    if (pBF->qHasVert(pVNeigh)) spBFCommon.insert(pBF);
	  }
	  switch (spBFCommon.size()) {
	  case 0:
	    // No common BFaces; connected across the interior of the mesh.
	    vMessage(4, "  Connecting edge not even on the bdry.\n");
	    break;
	  case 1:
	    // Should never happen.
	    vFatalError("Exactly one common bdry face between two points?!?",
			"bdry subseg list during small angle protection");
	    break;
	  case 2:
	    // Are the two BFace's on the same or different BdryPatch's?
	    {
	      iterBF = spBFCommon.begin();
	      BFace *pBF0 = *iterBF;
	      iterBF++;
	      BFace *pBF1 = *iterBF;
	      assert(pBF0->eType() == Cell::eTriBFace);
	      assert(pBF1->eType() == Cell::eTriBFace);
	      BdryPatch* pPatch0 = pBF0->pPatchPointer();
	      BdryPatch* pPatch1 = pBF1->pPatchPointer();
	      if (pPatch0 == pPatch1) {
		vMessage(4, "  Connected by edge inside patch.\n");
		qIsBdrySeg = false;
	      }
	      else {
		qIsBdrySeg = true;
	      }
	    }
	    break;
	  default:
	    // For more than two shared BFace's, this edge must be the
	    // meeting of more than two facets.
	    qIsBdrySeg = true;
	    break;
	  } // Done with switch over # of common BFace's
	} // Done checking whether it's really a bdry subseg
	if (!qIsBdrySeg) {
	  std::set<Vert*>::iterator iterVTmp = iterV--;
	  spVNeigh.erase(iterVTmp);
	}
      } // Loop over all neighboring verts

      // Calculate the distance from pV to the new ring of verts that
      // will be inserted to protect the corner.  To prevent clashes
      // with other nearby corners, this can never be more than half the
      // distance to the nearest corner.  Using half the length scale
      // also handles cases of nearby non-incident segments or facets.
      double dDistance = pV->dLS() / 3.;

      // For each incident segment, insert the point.
      for (iterV = spVNeigh.begin(); iterV != spVNeigh.end(); iterV++) {
	// Split this edge.
	Vert *pVNeigh = *iterV;
	double adVec[] = adDIFF3D(pVNeigh->adCoords(),
				  pV->adCoords());
	vNORMALIZE3D(adVec);
	adVec[XDIR] = pV->dX() + adVec[XDIR] * dDistance;
	adVec[YDIR] = pV->dY() + adVec[YDIR] * dDistance;
	adVec[ZDIR] = pV->dZ() + adVec[ZDIR] * dDistance;

	Vert *pVNew = pVNewVert();
	pVNew->vSetCoords(3, adVec);

	// Find a BFace that contains both pV and pVNeigh.  Everything
	// in LpBFInc contains pV, so that's already taken care of.
	{
	  std::set<Vert*> spVTmp;
	  bool qBdryVert;
	  std::set<Cell*> spCTmp;
	  vNeighborhood(pV, spCTmp, spVTmp, &spBFInc, &qBdryVert);
	}
	BFace *pBF = pBFInvalidBFace;
	std::set<BFace*>::iterator iterBF;
	for (iterBF = spBFInc.begin(); iterBF != spBFInc.end(); iterBF++) {
	  pBF = *iterBF;
	  if (pBF->qHasVert(pVNeigh)) break;
	}
	assert(pBF->qHasVert(pV));
	assert(pBF->qHasVert(pVNeigh));
	SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
	iInsertOnBdryEdge(pVNew, pV, pVNeigh, pBF->pFFace(), true);
	pVNew->vSetType(Vert::eBdryCurve);
	vUpdateLengthScale(pVNew);
	SUMAA_LOG_EVENT_END(INSERT_VERTEX);
	vMessage(4, "Inserted between verts at (%12g, %12g, %12g)",
		 pV->dX(), pV->dY(), pV->dZ());
	vMessage(4, "and (%12g, %12g, %12g)\n",
		 pVNeigh->dX(), pVNeigh->dY(), pVNeigh->dZ());
	vMessage(4, "  New point located at (%12g, %12g, %12g)\n",
		 pVNew->dX(), pVNew->dY(), pVNew->dZ());
      } // Done splitting this subsegment
    } // Done splitting all subsegs incident on this small angle apex
  } // Done with all small angle apexes
}

void VolMesh::vCreatePatchesFromBdryFaces(const int aiBCList[])
  // The array aiBCList tells what the bdry conditions are for each
  // BFace.  This data makes it trivial to assemble the Bdry3D object.
{
  if (pB3D) delete pB3D;
  qBdryFromMesh = true;
  pB3D = new Bdry3D(VAVerts);

  for (int iBF = 0; iBF < iNumBdryFaces(); iBF++) {
    Vert *apV[3];

    BFace *pBF = pBFBFace(iBF);
    Face *pF = pBF->pFFace();
    Cell *pC0 = pF->pCCellLeft();
    Cell *pC1 = pF->pCCellRight();

    apV[0] = pF->pVVert(0);
    apV[1] = pF->pVVert(1);
    apV[2] = pF->pVVert(2);

    int iBCL = iInvalidBC, iBCR = iInvalidBC;
    int iRegR = iInvalidRegion, iRegL = iInvalidRegion;

    if (pC0 == pBF) {
      // Boundary on the left!
      iBCL = aiBCList[iBF];
      iRegR = pC1->iRegion();
    }
    else {
      assert(pC1 == pBF);
      // Boundary on the right!
      iBCR = aiBCList[iBF];
      iRegL = pC0->iRegion();
    }
    BdryPolygon *pBPoly = new BdryPolygon(iBCL, iBCR, iRegL, iRegR, 3,
					  apV);
    pB3D->vAddPatch(pBPoly);
    pBF->vSetPatch(pBPoly);
  } // Done with loop over all BFaces to create Bdry3D object.
  assert(pB3D->iNumBdryPatches() == this->iNumBdryFaces());
}

void VolMesh::vIdentifyVertexTypesGeometrically() const
{
  int iV;
  for (iV = 0; iV < iNumVerts(); iV++) {
    pVVert(iV)->vSetType(Vert::eInterior);
  }

  // Identify and count all boundary vertices and make a list of all the
  // normals for boundary faces incident on a given vertex.  Once this
  // list is sorted, it will be used to divide boundary verts into
  // eBdryApex, eBdryCurve, and eBdry verts.

  int iBF;
  assert(qSimplicial() && iNumIntBdryFaces() == 0);
  BdryNorm *aBN = new BdryNorm[iNumBdryFaces() * 3];
  for (iBF = 0; iBF < iNumBdryFaces(); iBF++) {
    BFace *pBF = pBFBFace(iBF);
    double adNorm[3];
    pBF->vUnitNormal(adNorm);
    for (iV = 0; iV < pBF->iNumVerts(); iV++) {
      aBN[3*iBF + iV] = BdryNorm(pBF->pVVert(iV), adNorm);
    }
  }

  // For each run of identical vertices, identify the number of
  // distinct normals and tag the vertex accordingly:
  //  1 -> eBdry
  //  2 -> eBdryCurve
  //  3 or more -> eBdryApex
  // Mark boundary apexes to NOT be deleted; that is, clear
  // DELETE_REQUEST for these vertices.

  sort(aBN, aBN + 3*iNumBdryFaces());
  int iFirst, iLast;
  int iBdry = 0, iBdryCurve = 0, iBdryApex = 0, iTotal = 0;
  for (iFirst = 0; iFirst < 3*iNumBdryFaces(); iFirst = iLast+1) {
    Vert *pV = aBN[iFirst].pVVert();
    iTotal++;

    // Make sure this vert hasn't been touched before.
    assert(pV->iVertType() == Vert::eInterior);

    // Identify the first and last record with this vertex in it.
    iLast = iFirst;
    while (iLast <= 3*iNumBdryFaces() && aBN[iLast].pVVert() == pV) {
      iLast++;
    };
    iLast --;
    assert(aBN[iLast].pVVert() == pV);

    // For each boundary vertex, find the number of distinct normals
    // for bdry faces incident on it.  If there are more than two,
    // this is a eBdryApex.  If there are zero, this is just an
    // ordinary eBdry vert.  Otherwise, there should be exactly two.
    //   In this case, find the edges that divide these normals into
    // two groups.  If the edges are roughly collinear, then the
    // vertex is on a eBdryCurve.  If the edges aren't collinear, then
    // the vertex is a eBdryApex.  This is a knife case, which may be
    // uncommon in practice but should still be handled correctly.
    int iNumNorms = 0;
    int aiNormTag[1000]; // Should be far fewer bdry faces incident on
    // a single bdry vert.
    int iNorm, iTag = 0;

    assert(iFirst != iLast);
    int iDist = iLast - iFirst;

    vMessage(4, "Vertex %d has %d incident bdry faces.\n",
	     iVertIndex(pV), iDist+1);

    for (iNorm = 0; iNorm <= iDist; iNorm++)
      aiNormTag[iNorm] = -1;

    for (iNorm = 0; iNorm <= iDist; iNorm++) {
      if (aiNormTag[iNorm] == -1) {
	aiNormTag[iNorm] = iTag;

	const double *adNormN = aBN[iFirst + iNorm].adNormal();
	for (int iiF = iNorm+1; iiF <= iDist; iiF++) {
	  if (aiNormTag[iiF] == -1) {
	    const double *adNormF = aBN[iFirst + iiF].adNormal();
	    double dDot = dDOT3D(adNormF, adNormN);
	    // All normals are directed into the domain, so there's no
	    // need to take the absolute value of dDot.  And taking the
	    // absolute value can cause faces that have normals nearly
	    // opposite each other to test wrong here.
	    if (dDot > cos(dMaxAngleForSurfSwap * M_PI/180)) {
	      aiNormTag[iiF] = iTag;
	    }
	  }
	} // Loop over rest of bdry faces for this vert
	iTag++;
      } // If this bdry faces is untagged, process it
    } // Check all bdry faces incident on the vertex.
    iNumNorms = iTag;
    assert(iNumNorms >= 1);

    // Tag the vertex appropriately.
    // Fix for 0.2.1:  There is a knife case in which there are only
    // two distinct normals at a eBdryApex.
    switch (iNumNorms) {
    case 1:
      pV->vSetType(Vert::eBdry);
      iBdry++;
      break;
    case 2:
      pV->vSetType(Vert::eBdryCurve);
      iBdryCurve++;
      break;
    default:
      pV->vSetType(Vert::eBdryApex);
      iBdryApex++;
      pV->vMarkToKeep();
      break;
    } // end switch
  }
  vMessage(2, "Out of %d verts on the bdry, %d were apexes, %d on curves.\n",
	   iTotal, iBdryApex, iBdryCurve);
  assert(iTotal == iBdry + iBdryApex + iBdryCurve);
  delete [] aBN;
}

bool VolMesh::qBdryFacesOKToSwap(const TriBFaceBase * const pBF0,
				 const TriBFaceBase * const pBF1) const
{
  if (qStrictPatchChecking) {
    return (pBF0->pPatchPointer() == pBF1->pPatchPointer());
  }
  else {
    // Check that the BC's are the same and that the normals are more or
    // less aligned (within the tolerances of dMaxAngleForSurfSwap).
    if (pBF0->iBdryCond() == pBF1->iBdryCond()) {
      double adNorm0[3], adNorm1[3];
      pBF0->vUnitNormal(adNorm0);
      pBF1->vUnitNormal(adNorm1);
      double dDot = dDOT3D(adNorm0, adNorm1);
      return (dDot >= cos(dMaxAngleForSurfSwap * M_PI/180));
    }
    else {
      return false;
    }
  }
}

typedef struct _vert_vert_ {
  int iV0, iV1;
  _vert_vert_ (const int i0 = -1, const int i1 = -1) :
    iV0(i0), iV1(i1) {}
  bool operator<(const struct _vert_vert_ v2) const {
    return (iV0 < v2.iV0 ||
	    (iV0 == v2.iV0 && iV1 < v2.iV1));
  }
}
VertVert;

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

extern "C"
{
  static int qRenumComp3D (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 VolMesh::
vReorderVerts_RCM (void)
{
  vMessage(2, "Renumbering vertices...");
  assert(qSimplicial());
  int iVert, iFace, iNVerts = iNumVerts (), iNFaces = iNumFaces ();
  //@@ Construct the vert-vert connectivity
  // Make a list containing each edge (once in each direction) and sort
  int iNEdges;
  VertVert *aVV;
  {
    std::set<VertVert> sVV;
    for (iFace = iNFaces - 1; iFace >= 0; iFace--) {
      Face *pF = pFFace (iFace);
      int iV0 = iVertIndex (pF->pVVert (0));
      int iV1 = iVertIndex (pF->pVVert (1));
      int iV2 = iVertIndex (pF->pVVert (2));
      sVV.insert(VertVert (iV1, iV0));
      sVV.insert(VertVert (iV0, iV1));
      sVV.insert(VertVert (iV1, iV2));
      sVV.insert(VertVert (iV2, iV1));
      sVV.insert(VertVert (iV2, iV0));
      sVV.insert(VertVert (iV0, iV2));
    }
    // This is actually twice the number of unique edges, because both
    // directions appear here.  iNEdges shouldn't be roughly the same as
    // iNFaces.
    iNEdges = sVV.size();
    aVV = new VertVert[iNEdges];
    VertVert *pVV = copy(sVV.begin(), sVV.end(), aVV);
    assert(pVV - aVV == iNEdges);
  }

  // Tabulate starting locations and degree for each vertex.
  std::vector<int> viStart(iNVerts + 1), viDegree(iNVerts);
  std::vector<bool> vqUsed(iNVerts);
  int iStart = 0, iNAttachedVerts = 0;
  for (iVert = 0; iVert < iNVerts; iVert++) {
    viStart[iVert] = iStart;
    for (; aVV[iStart].iV0 == iVert && iStart < iNEdges; iStart++) {}
    if (iStart != viStart[iVert])
      iNAttachedVerts++;
    viDegree[iVert] = iStart - viStart[iVert];
    vqUsed[iVert] = false;
  }
  viStart[iNVerts] = iNEdges;

  // 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 = viStart[iVert + 1] - viStart[iVert];
      if (iDeg < iMinDeg && iDeg > 0 && !vqUsed[iVert]) {
	iMinDeg = iDeg;
	iMinVert = iVert;
	if (iMinDeg == 2)
	  break;	// Can never be lower than this, so why go on?
      }
    }
    assert (iMinVert != -1);

    // Initialize
    vqUsed[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 = viStart[iV]; iNeigh < viStart[iV + 1]; iNeigh++) {
	  // Grab a candidate
	  assert (aVV[iNeigh].iV0 == iV);
	  int iCand = aVV[iNeigh].iV1;
	  if (!vqUsed[iCand]) {
	    // Add it if it isn't already in the list
	    aRReorder[iNumRenumbered].iV = iCand;
	    aRReorder[iNumRenumbered].iDeg = viDegree[iCand];
	    vqUsed[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), qRenumComp3D);
    } while (iEnd != iBegin);
  }

  for (iVert = 0; iNumRenumbered != iNVerts; iVert++) {
    if (!vqUsed[iVert]) {
      aRReorder[iNumRenumbered].iV = iVert;
      vqUsed[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[3];		// Only two verts to an edge

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

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

void VolMesh::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());
  std::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.
  std::vector<int> aiCellMapping(iNumCells());
  std::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.
  std::vector<bool> aqUpdated(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::eTet);
      TetCell TC;
      TC = *(dynamic_cast<TetCell*>(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 3D meshes");
  }
  vMessage(2, "done.\n");
}

void VolMesh::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());
  std::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.
  std::vector<int> aiFaceMapping(iNF);
  std::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.
  std::vector<bool> aqUpdated(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::eTriFace);
      TriFace TF;
      TF = *(dynamic_cast<TriFace*>(pFFace(iFirst)));
      int iNext = iFirst;
      while (aiFaceMapping[iNext] != iFirst) {
	(*pFFace(iNext)) = (*pFFace(aiFaceMapping[iNext]));
	aqUpdated[iNext] = true;
	iNumUpdated++;
	iNext = aiFaceMapping[iNext];
      }
      (*pFFace(iNext)) = TF;
      aqUpdated[iNext] = true;
      iNumUpdated++;
    }
  }

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

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

