#include <queue>

#include "GR_config.h"
#include "GR_events.h"
#include "GR_misc.h"
#include "GR_Constrain.h"
#include "GR_Geometry.h"
#include "GR_SurfMesh.h"
#include "GR_VolMesh.h"

bool qIsBdryConstrained(VolMesh* pVM, SurfMesh* pSM,
			std::set<ConstrainEdge>& sCESubSegs,
			std::set<ConstrainEdge>& sCEInFacets,
			std::set<ConstrainFace>& sCFSubFacets)
{
  assert(pVM->qSimplicial());
  assert(pSM->qSimplicial());

  //@@ Make sure that are the vertex hints are valid in both meshes
  pVM->vSetAllHintFaces();
  pSM->vSetAllHintFaces();

  //@@ Clear the lists of missing edges and faces
  sCESubSegs.clear();
  sCEInFacets.clear();
  sCFSubFacets.clear();

  // The old code for this went through a vertex at a time, if you can
  // imagine.  Looking at the problem now, the obvious choice seems to
  // be to make sets of subsements, in-facet edges, and subfacets from
  // the surface mesh.  Then make sets of edges and faces in the volume
  // mesh.  Finally, subtract these sets to get the real list of
  // unconstrained edges and faces.
  std::set<ConstrainEdge> sCESurfSegs, sCESurfInFacets, sCEVolEdges;
  std::set<ConstrainFace> sCFSurfTris, sCFVolTris;

  for (int iEdge = pSM->iNumFaces() - 1; iEdge >= 0; iEdge--) {
    Face *pF = pSM->pFFace(iEdge);
    if (pF->qDeleted()) continue;
    Cell *pCR = pF->pCCellRight();
    Cell *pCL = pF->pCCellLeft();

    BdryPatch3D *pBPLeft = pSM->pBdryPatch(pSM->iCellIndex(pCL));
    BdryPatch3D *pBPRight = pSM->pBdryPatch(pSM->iCellIndex(pCR));

    int iV0 = pSM->iVertIndex(pF->pVVert(0));
    int iV1 = pSM->iVertIndex(pF->pVVert(1));
    
    if (pBPLeft == pBPRight) {
      // Interior to a facet
      sCESurfInFacets.insert(ConstrainEdge(iV0, iV1));
    }
    else {
      // On a segment
      sCESurfSegs.insert(ConstrainEdge(iV0, iV1));
    }
  }

  // Now for the triangles in the surface list.
  for (int iCell = pSM->iNumCells() - 1; iCell >= 0; iCell--) {
    Cell *pC = pSM->pCCell(iCell);
    if (pC->qDeleted()) continue;
    
    int iV0 = pSM->iVertIndex(pC->pVVert(0));
    int iV1 = pSM->iVertIndex(pC->pVVert(1));
    int iV2 = pSM->iVertIndex(pC->pVVert(2));

    sCFSurfTris.insert(ConstrainFace(iV0, iV1, iV2));
  }

  // Build comparable sets from the volume mesh.
  for (int iFace = pVM->iNumFaces() - 1; iFace >= 0; iFace--) {
    Face *pF = pVM->pFFace(iFace);
    if (pF->qDeleted()) continue;
    
    int iV0 = pVM->iVertIndex(pF->pVVert(0));
    int iV1 = pVM->iVertIndex(pF->pVVert(1));
    int iV2 = pVM->iVertIndex(pF->pVVert(2));

    sCFVolTris.insert(ConstrainFace(iV0, iV1, iV2));
    sCEVolEdges.insert(ConstrainEdge(iV0, iV1));
    sCEVolEdges.insert(ConstrainEdge(iV1, iV2));
    sCEVolEdges.insert(ConstrainEdge(iV2, iV0));
  }
  
  // Now use set_difference to find out what's missing.  Yes, there will
  // be lots of entries from the volume mesh that don't match anything,
  // but this construction is still only order of the size of the volume
  // mesh....
  std::insert_iterator<std::set<ConstrainEdge> >
    iiSubSegs(sCESubSegs, sCESubSegs.begin()),
    iiInFacets(sCEInFacets, sCEInFacets.begin());
  std::set_difference(sCESurfSegs.begin(), sCESurfSegs.end(),
		      sCEVolEdges.begin(), sCEVolEdges.end(),
		      iiSubSegs);
  std::set_difference(sCESurfInFacets.begin(), sCESurfInFacets.end(),
		      sCEVolEdges.begin(), sCEVolEdges.end(),
		      iiInFacets);

  std::insert_iterator<std::set<ConstrainFace> > iiFaces(sCFSubFacets,
							 sCFSubFacets.begin());
  std::set_difference(sCFSurfTris.begin(), sCFSurfTris.end(),
		      sCFVolTris.begin(), sCFVolTris.end(),
		      iiFaces);

  //@@ Issue a report on the results
  if (sCESubSegs.empty() && sCEInFacets.empty() && sCFSubFacets.empty()) {
    vMessage(1, "All edges and faces are properly constrained.\n");
    return true;
  }
  else {
    vMessage(1, "Missing %d input edges, %d/%d total edges, %d/%d faces.\n",
	     static_cast<int>(sCESubSegs.size()),
	     static_cast<int>(sCESubSegs.size()+sCEInFacets.size()),
	     pSM->iNumFaces(), static_cast<int>(sCFSubFacets.size()),
	     pSM->iNumCells());
    return false;
  }
}

int VolMesh::iRecoverEdge(const Vert* const pV0,
			  const Vert* const pV1,
			  const bool qInsertionOK)
{
  assert(qSimplicial());

  if (qInsertionOK) {
    // FIX ME: Shouldn't bother with this unless the segment is
    // encroached.  Why?  Because if no segments are encroached, then
    // we can always recover the bdry segs.  (If segments are present
    // but encroached, that can always be fixed later.)

    // Curved bdry: replace with code that splits the underlying patch
    // boundary properly, instead of simply taking the midpoint of a
    // segmetn.
    double adNewCoords[] = {0.5 * (pV0->dX() + pV1->dX()),
			    0.5 * (pV0->dY() + pV1->dY()),
			    0.5 * (pV0->dZ() + pV1->dZ())};
    Vert *pVBdry0 = pB3D->pVVert(iVertIndex(pV0));
    Vert *pVBdry1 = pB3D->pVVert(iVertIndex(pV1));
    pB3D->pVAddVertex(pVBdry0, pVBdry1, adNewCoords);
    {
      int iNSwaps = 0, iCell = 0;
      Cell *pC = pCCell(iCell);
      while (!pC->iFullCheck()) {
	pC = pCCell(++iCell);
      }
      // Any live tet will do, since the domain is concave.
      SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
#ifndef NDEBUG
      bool qOK =
#endif
	qInsertPoint(adNewCoords, pC, &iNSwaps, true, true);
      SUMAA_LOG_EVENT_END(INSERT_VERTEX);
      // Mark this vert correctly
      pVVert(iNumVerts()-1)->vSetType(Vert::eBdryCurve);
      assert(qOK);
    }
    // That particular edge certainly doesn't need to be recovered; a
    // return value of 1 makes sure we never try to again.
    return 1;
  }

  int iChange, iNSkip, iNPasses = 0, iRetVal = 0;
  Face* apFSkipped[20];
  do {
    iChange = iNSkip = 0;
    bool qHitEdge = false;
    iNPasses ++;
    vMessage(4, "Pass %d to remove this edge\n", iNPasses);

    // Walk from pV0 towards pV1 until either (a) you hit an edge
    // or (b) you get to pV1.  Keep track of all faces you had to skip
    // over.
    int iRes = iReamPipe(pV0, pV1, qHitEdge, apFSkipped, iNSkip);
    if (iRes == 1) // Pipe is unclogged
      return 1;
    if (iRes == 0) { // Changes made
      iChange = 1;
      iRetVal = -1;
    }
    if (qHitEdge) { // Try from the other end
      qHitEdge = false;
      iNSkip = 0;
      iRes = iReamPipe(pV1, pV0, qHitEdge, apFSkipped, iNSkip);
      if (iRes == 1) // Pipe is unclogged; not bloody likely
	return 1;
      if (iRes == 0) { // Something has changed
	iChange = 1;
	iRetVal = -1;
      }
    }

    // If an edge was hit, try to remove it
    if (qHitEdge) {
      vMessage(3, "Unresolved edge hit!  Edge not recovered.\n");
    }
  } while (iChange && iNPasses < 5);
  return iRetVal;
}

void VolMesh::vIdentifyBadEdgesAndFaces(const Vert* const pV0,
					const Vert* const pV1,
					const Vert* const pV2,
					std::set<BadEdge>& sBEBadEdges,
					std::set<Face*>& spFBadFaces)
{
  // Now proceed to recover the face.  First, find a cell containing an
  // edge of the required face.
  Cell *pC = pCBeginPipe(pV0, pV1);
  assert(pC->qHasVert(pV0));
  assert(pC->qHasVert(pV1));

  sBEBadEdges.clear();
  spFBadFaces.clear();

  // If the face already exists within this cell, great.
  if (pC->qHasVert(pV2)) return;

  // Walk around edge to find cell that crosses the surface face.
  int qDone;
  do {
    qDone = 1;
    for (int iF = 0; iF < 4; iF++) {
      Face *pF = pC->pFFace(iF);
      if ( ! (pF->qHasVert(pV0) && pF->qHasVert(pV1)) ) continue;
      assert(pF->eType() == Face::eTriFace);
      // Walk through this face if the search vertex is on the opposite
      // side of the plane of the face from the other vertex in pC.
      int iOrient0 = iOrient3D(pF->pVVert(0), pF->pVVert(1),
			       pF->pVVert(2), pV2);
      int iOrient1 = iOrient3D(pF->pVVert(0), pF->pVVert(1),
			       pF->pVVert(2),
			       dynamic_cast<TetCell*>(pC)->pVVertOpposite(pF));
      if (iOrient0 * iOrient1 == -1) { // Don't walk in case of ties.
	qDone = 0;
	pC = pF->pCCellOpposite(pC);
	assert(pC->eType() == Cell::eTet);
	break;
      }
    } // end loop over faces of a tet
  } while (!qDone);
  // Now we have a cell containing both pV0 and pV1 and also spanning
  // the direction from pV0 towards pV2.  Again, we may already be
  // done.
  if (pC->qHasVert(pV2)) return;

  //@@@ Identify the bad edge and faces
  Vert *pVNorth = pVInvalidVert, *pVSouth = pVInvalidVert,
    *pVOther = pVInvalidVert;
  Face *pF;
  for (int iF = 0; iF < 4; iF++) {
    // Get a face which has only one of pV0 and pV1
    pF = pC->pFFace(iF);
    if (pF->qHasVert(pV0) && pF->qHasVert(pV1)) continue;
    // Identify the verts at either end of the bad edge
    if (pF->pVVert(0) == pV0 || pF->pVVert(0) == pV1) {
      pVOther = pF->pVVert(0);
      if (iOrient3D(pV0, pV1, pV2, pF->pVVert(1)) == 1) {
	pVNorth = pF->pVVert(1);
	pVSouth = pF->pVVert(2);
      }
      else {
	pVNorth = pF->pVVert(2);
	pVSouth = pF->pVVert(1);
      }
    }
    else if (pF->pVVert(1) == pV0 || pF->pVVert(1) == pV1) {
      pVOther = pF->pVVert(1);
      if (iOrient3D(pV0, pV1, pV2, pF->pVVert(0)) == 1) {
	pVNorth = pF->pVVert(0);
	pVSouth = pF->pVVert(2);
      }
      else {
	pVNorth = pF->pVVert(2);
	pVSouth = pF->pVVert(0);
      }
    }
    else {
      pVOther = pF->pVVert(2);
      assert(pF->pVVert(2) == pV0 || pF->pVVert(2) == pV1);
      if (iOrient3D(pV0, pV1, pV2, pF->pVVert(0)) == 1) {
	pVNorth = pF->pVVert(0);
	pVSouth = pF->pVVert(1);
      }
      else {
	pVNorth = pF->pVVert(1);
	pVSouth = pF->pVVert(0);
      }
    }
    break;
  }

  sBEBadEdges.insert(BadEdge(pVNorth, pVSouth, pVOther, pF));
  spFBadFaces.insert(pF);

  std::queue<BadEdge> qBE;
  qBE.push(BadEdge(pVNorth, pVSouth, pVOther, pF));
  //@@@ Walk around the edge, identifying all bad faces and edges
  while (!qBE.empty()) {
    // Set up starting face, cell, and other vert of starting face
    BadEdge BE = qBE.front();
    qBE.pop();
    Face *pFStart;
    pFStart = pF = BE.pF;
    pVNorth = BE.pVN;
    pVSouth = BE.pVS;
    Vert *pV = BE.pVO;
    pC = pF->pCCellLeft();
    do {
      // Grab new face incident on edge and add to bad face list
      assert(pC->eType() == Cell::eTet);
      TetCell* pTC = dynamic_cast<TetCell*>(pC);
      Face *pFPrev = pF;
      pF = pTC->pFFaceOpposite(pV);
      spFBadFaces.insert(pF);
      // Grab other vertex of that face
      pV = pTC->pVVertOpposite(pFPrev);
      pC = pF->pCCellOpposite(pC);
      // If vertex is not one of pV0, pV1, pV2, there's another bad edge
      // here. Identify it and store it.
      if (pV != pV0 && pV != pV1 && pV != pV2) {
	int iOrient = iOrient3D(pV0, pV1, pV2, pV);
	assert(iOrient != 0);
	if (iOrient == 1) {
	  // pV is on the "north" side of the surface face
	  BadEdge BETmp(pV, pVSouth, pVNorth, pF);
	  if (sBEBadEdges.count(BETmp) == 0) {
	    sBEBadEdges.insert(BETmp);
	    qBE.push(BETmp);
	  }
	}
	else {
	  // pV is on the "south" side of the surface face
	  BadEdge BETmp(pVNorth, pV, pVSouth, pF);
	  if (sBEBadEdges.count(BETmp) == 0) {
	    sBEBadEdges.insert(BETmp);
	    qBE.push(BETmp);
	  }
	}
      } // Done adding new bad edge
    } while(pF != pFStart);
  } // Done adding all bad edges
  vMessage(4, "Found %d edges intersecting desired face.\n",
	   static_cast<int>(sBEBadEdges.size()));
  vMessage(4, "Found %d faces intersecting desired face.\n",
	   static_cast<int>(spFBadFaces.size()));
}

//@ Recover a surface face that's missing from the volume mesh
// This is done by swapping if possible and by point insertion if necessary
int VolMesh::iRecoverFace(Mesh2D* pM2D, const int aiVerts[3],
			  int * const piPointsInserted,
			  const bool qAllowEdgeSwap)
{
  //   The following assertion chokes when called from
  //   VolMesh::VolMesh(Bdry3D&)
  //   assert(qValid());
  assert(pM2D->qValid());
  vMessage(3, "%s %d, %d, and %d.\n",
	   "Trying to recover surface face composed of vertices ",
	   aiVerts[0], aiVerts[1], aiVerts[2]);
  //@@ Find all faces and edges violating the given surface face.
  Vert *pV0 = pVVert(aiVerts[0]);
  Vert *pV1 = pVVert(aiVerts[1]);
  Vert *pV2 = pVVert(aiVerts[2]);
  // First, make sure that all the edges of the surface face are
  // present.  If they aren't, this means that they are interior edges
  // to a surface facet (and so weren't recovered earlier when doing
  // bdry subsegments); and that there is a vertex encroaching on them
  // (which makes it so that the Delaunay tetrahedralization doesn't
  // match the constrained DT).  If the edge is missing, recover it
  // without point insertion (this should always be possible, but note
  // that there's still a check to be sure that it works).
  int iEdgeAOK = iRecoverEdge(pV0, pV1);
  int iEdgeBOK = iRecoverEdge(pV1, pV2);
  int iEdgeCOK = iRecoverEdge(pV2, pV0);
  if ((iEdgeAOK != 1) || (iEdgeBOK != 1) || (iEdgeCOK != 1))
    vFatalError("Failed to recover a non-subsegment bdry edge.",
		"3D constrained Delaunay tetrahedralization.");

  std::set<BadEdge> sBEBadEdges;
  std::set<Face*> spFBadFaces;
  vIdentifyBadEdgesAndFaces(pV0, pV1, pV2, sBEBadEdges, spFBadFaces);

  //@@ For each face in the list, swap it if possible.
  // Swap only if the face is of type T32 and all of its equatorial
  // points are on the same side of or in the plane of the surface face.
  bool qChange, qGlobalChange = false;
  do {
    qChange = false;
    std::queue<Face*> qpFBad;
    for (std::set<Face*>::iterator iterF = spFBadFaces.begin();
	 iterF != spFBadFaces.end(); iterF++) {
      qpFBad.push(*iterF);
    }

    while (!qpFBad.empty()) {
      Face* pF = qpFBad.front();
      qpFBad.pop();
      assert(pF->eType() == Face::eTriFace);
      TriFace *pTF = dynamic_cast<TriFace*>(pF);
      if (pTF->qDeleted()) continue;
      Vert *pVVertA, *pVVertB, *pVVertC, *pVVertD, *pVVertE;
      TetCell* apTCTets[4];
      int iNTets;
      Vert *pVJunk0, *pVJunk1, *pVJunk2;
      eFaceCat eFC = pTF->eCategorizeFace(pVVertA, pVVertB, pVVertC,
					  pVVertD, pVVertE, apTCTets, iNTets,
					  pVJunk0, pVJunk1, pVJunk2);

      if (eFC == eT32) {
	assert(iNTets == 3);
	Face *apFDel[3];
	apFDel[0] = pFCommonFace(apTCTets[0], apTCTets[1]);
	apFDel[1] = pFCommonFace(apTCTets[1], apTCTets[2]);
	apFDel[2] = pFCommonFace(apTCTets[2], apTCTets[0]);
	assert(apFDel[0]->qValid());
	assert(apFDel[1]->qValid());
	assert(apFDel[2]->qValid());
	int iOrientD = iOrient3D(pV0, pV1, pV2, pVVertD);
	int iOrientE = iOrient3D(pV0, pV1, pV2, pVVertE);
	if (iOrientD * iOrientE != -1) continue;
	assert(pTF == apFDel[0] || pTF == apFDel[1] || pTF == apFDel[2]);

	qChange = true;
	qGlobalChange = true;
#ifndef NDEBUG
	int iSwaps =
#endif
	  iReconfTet32(apTCTets, pTF, pVVertA, pVVertB, pVVertC,
		       pVVertD, pVVertE);
	assert(iSwaps == 1);
	// Check each of the three faces to be sure that either
	// they've been deleted or all of their verts lie on the
	// same side of the target face (including ties).
	for (int ii = 0; ii < 3; ii++) {
	  Face *pFTest = apFDel[ii];
	  if (pFTest->qDeleted()) {
	    spFBadFaces.erase(pFTest);
	  }
	  else {
	    Vert *pVTest0 = pFTest->pVVert(0);
	    Vert *pVTest1 = pFTest->pVVert(1);
	    Vert *pVTest2 = pFTest->pVVert(2);
	    int iOrient0 = iOrient3D(pV0, pV1, pV2, pVTest0);
	    int iOrient1 = iOrient3D(pV0, pV1, pV2, pVTest1);
	    int iOrient2 = iOrient3D(pV0, pV1, pV2, pVTest2);
	    if ((iOrient0*iOrient1 != -1) &&
		(iOrient1*iOrient2 != -1) &&
		(iOrient2*iOrient0 != -1)) {
	      spFBadFaces.erase(pFTest);
	    }
	  }
	}
      } // T32 face
    } // End loop over all bad faces
    vIdentifyBadEdgesAndFaces(pV0, pV1, pV2, sBEBadEdges, spFBadFaces);
    // Do this until no more progress can be made by simple swapping or
    // until all offending edges are removed
  } while (qChange && !spFBadFaces.empty());

  // If no swapping has been done at all, then we can do edge swapping
  // based on the original bad edge data.  Also, don't allow too much
  // explosion of bad edges before we give up entirely.

  // I don't think the following implementation has any hope of working.
  if (!spFBadFaces.empty() && !qGlobalChange &&
      sBEBadEdges.size() < 40 && qAllowEdgeSwap) {
     vMessage(4, "Trying edge swapping...\n");
     bool qSuccess = false;
     for (std::set<BadEdge>::iterator iter = sBEBadEdges.begin();
	  iter != sBEBadEdges.end() && !qSuccess; iter++) {
       BadEdge BE = *iter;
       qSuccess = (0 != iEdgeSwap3D(BE.pF, BE.pVN, BE.pVS, BE.pVO,
 				   true));
     }
     // If an edge has been removed, then recurse.
     if (qSuccess)
       iRecoverFace(pM2D, aiVerts, piPointsInserted, qAllowEdgeSwap);
  }

  if (spFBadFaces.empty()) {
    vMessage(4, "All bad faces removed by swapping\n");
    return 1;
  }
  else
    vMessage(4, "After swapping, %d bad faces remain\n",
	     static_cast<int>(spFBadFaces.size()));
  if (piPointsInserted == NULL)
    return 0;

  //@@ If swapping fails, insert a point to split an offending face.
  // The bad face should be split at its circumcenter, so long as this
  // doesn't violate encroachment rules.  Not yet implemented,
  // obviously.
  vMessage(2, "Would need to try face splitting to recover faces; not implemented.\n");
  return 0;
}

void VolMesh::vMarkCleanNeighbors(Cell* pC)
{
  int iReg = pC->iRegion();
  for (int iF = 0; iF < pC->iNumFaces(); iF++) {
    Face *pF = pC->pFFace(iF);
    // If the face has been deleted, then we need not do anything.
    if (pF->qValid() && !pF->qDeleted()) {
      Cell *pCCand = pF->pCCellOpposite(pC);
      assert(pCCand->qValid());

      // Be sure never to deal with this face again
      if (pCCand->iRegion() == iReg) continue;
      else if (pCCand->iRegion() == iInvalidRegion) {
	// If the opposite cell hasn't been marked, mark it the same as
	// this one.
	pCCand->vSetRegion(iReg);
	vMarkCleanNeighbors(pCCand);
      } // Candidate cell needs to be tagged.
    } // Face is OK
  } // Loop over faces
}

// Find the cell incident on pV0 that lies astride the direction to pV1.
Cell* VolMesh::pCBeginPipe(const Vert* const pV0, const Vert* const pV1) const
{
  Face *pF = pV0->pFHintFace();
  Cell *pC = pF->pCCellLeft();
  if (pC->eType() == Cell::eTriBFace) {
    pC = pF->pCCellRight();
    assert(pC->eType() == Cell::eTet);
  }
  bool qFoundIt = false;
#ifndef NDEBUG
  int iTrials = 0;
#endif
  while (!qFoundIt) {
    qFoundIt = true;
#ifndef NDEBUG
    iTrials ++;
    assert(iTrials < 1000);
#endif
    Face *pFOpp = dynamic_cast<TetCell*>(pC)->pFFaceOpposite(pV0);
    for (int iF = 0; iF < 4; iF++) {
      Face *pFTrial = pC->pFFace(iF);
      if (pFTrial == pFOpp) continue;
      int iSign = pFTrial->pCCellLeft() == pC ? -1 : 1;
      if (iOrient3D(pFTrial->pVVert(0), pFTrial->pVVert(1),
		    pFTrial->pVVert(2), pV1) == - iSign) {
	// Incorrect orientation; walk through this face and try again
	pC = pFTrial->pCCellOpposite(pC);
	qFoundIt = false;
	break;
      }
    }
  }
  return pC;
}

// Find the face (other than pFLast) through which the line from pV0 to
// pV1 crosses the surface of Cell pC.
static void vWhichFace(const Cell* const pC, const Face* const pFLast,
		       const Vert* const pV0, const Vert* const pV1,
		       const Face*& pF, bool& qHitEdge)
{
  qHitEdge = false;
  for (int iF = 0; iF < 4; iF++) {
    pF = pC->pFFace(iF);
    // Don't give back the same face that was previously found
    if (pF == pFLast) continue;
    // Don't give back a face containing one of the endpoints
    if (pF->qHasVert(pV0) || pF->qHasVert(pV1)) continue;
    const Vert *pVA = pF->pVVert(0);
    const Vert *pVB = pF->pVVert(1);
    const Vert *pVC = pF->pVVert(2);
    // If pV0 and pV1 aren't on opposite sides of the face, then you'll
    // have to try another.
    if (iOrient3D(pVA, pVB, pVC, pV0) != -iOrient3D(pVA, pVB, pVC, pV1))
      continue;
    int iOrientAB = iOrient3D(pVA, pVB, pV0, pV1);
    if (iOrientAB == 0) {qHitEdge = true; return;}  // Hit an edge
    int iOrientBC = iOrient3D(pVB, pVC, pV0, pV1);
    if (iOrientBC == 0) {qHitEdge = true; return;}  // Hit an edge
    int iOrientCA = iOrient3D(pVC, pVA, pV0, pV1);
    if (iOrientCA == 0) {qHitEdge = true; return;}  // Hit an edge
    if (iOrientAB == iOrientBC && iOrientAB == iOrientCA) return;
  }
}

//@ Work your way through a pipe trying to remove faces.
int VolMesh::iReamPipe(const Vert* const pVInit,
		       const Vert* const pVEnd,
		       bool& qHitEdge, Face* apFSkipped[],
		       int& iNSkip)
     // Return 1 if the pipe is cleared, 0 if progress is made but the
     // pipe still has internal faces, and -1 if no changes were made.
{
  // A new version of this routine, including edge swapping internally.

  // Find the first cell at the pVInit end of the pipe from pVInit to pVEnd
  Cell* pCStart0 = pCBeginPipe(pVInit, pVEnd);

  vMessage(4, "Pipe begins: %5d  Pipe ends: %5d\n",
	   iVertIndex(pVInit), iVertIndex(pVEnd));

  // Walk from pVInit towards pVEnd until either (a) you hit an edge or (b)
  // you get to pVEnd.  Keep track of how many faces you had to skip
  // over. Any time you swap a face, you must begin the ream again,
  // because the identities of cells and faces change with swapping, and
  // you could end up holding an invalid cell or face pointer, and
  // cratering the code.

  // Find which face the edge (pVInit,pVEnd) passes through
  Face *pFCand = pFInvalidFace, *pFLast = pFInvalidFace;
  Cell *pC = pCStart0;
  iNSkip = 0;
  TetCell* apTCTets[4];
  int iNTets;
  bool qChange = false, qReset = false, qEndOfPipe = false;
  int iCount = 0;
  eFaceCat ePrevFC = eOther;
  Vert *pVPrevD = pVInvalidVert, *pVPrevE = pVInvalidVert;

  do {
    if (pC->qHasVert(pVEnd)) {
      vMessage(4, "  Pipe already reamed.\n");
      return 1;
    }
    {
      const Face* pFTmp = pFInvalidFace;
      vWhichFace(pC, pFLast, pVInit, pVEnd, pFTmp, qHitEdge);
      pFCand = const_cast<Face*>(pFTmp);
    }

    if (qHitEdge) {
      // Try to swap away this edge in the volume mesh.  First
      // identify the edge.
      Vert *pVA = pFCand->pVVert(0);
      Vert *pVB = pFCand->pVVert(1);
      Vert *pVC = pFCand->pVVert(2);

      Vert *pVNorth, *pVSouth, *pVOther;

      if (iOrient3D(pVA, pVB, pVInit, pVEnd) == 0) {
	pVNorth = pVA;
	pVSouth = pVB;
	pVOther = pVC;
      }
      else if (iOrient3D(pVB, pVC, pVInit, pVEnd) == 0) {
	pVNorth = pVB;
	pVSouth = pVC;
	pVOther = pVA;
      }
      else {
	assert(iOrient3D(pVC, pVA, pVInit, pVEnd) == 0);
	pVNorth = pVC;
	pVSouth = pVA;
	pVOther = pVB;
      }
      int iSwaps = iEdgeSwap3D(pFCand, pVNorth, pVSouth,
			       pVOther, true, pVInit, pVEnd);
      if (iSwaps > 0) {
	// Success!  Well, at least we got rid of the offending edge.
	// So restart through the pipe.
	vMessage(4, "  Successfully flipped an edge.\n");
	qReset = qChange = true;
      }
      else {
	vMessage(4, "  Failed to flip an edge.\n");
      }
      qHitEdge = false;
    }
    else { // The normal case, where you haven't hit an edge
      assert(pFCand->eType() == Face::eTriFace);
      TriFace *pTF = dynamic_cast<TriFace*>(pFCand);
      Vert *pVJunk0, *pVJunk1, *pVJunk2;
      Vert *pVVertA, *pVVertB, *pVVertC, *pVVertD, *pVVertE;
      eFaceCat eFC = pTF->eCategorizeFace(pVVertA, pVVertB, pVVertC,
					  pVVertD, pVVertE,
					  apTCTets, iNTets, pVJunk0, pVJunk1,
					  pVJunk2);
      vMessage(4, " FaceCat = %d %8d %8d %8d %8d %8d\n", eFC,
	       iVertIndex(pVVertA), iVertIndex(pVVertB),
	       iVertIndex(pVVertC), iVertIndex(pVVertD),
	       iVertIndex(pVVertE));
      if ((eFC == eT44)
	  && (pVInit == pVVertD || pVEnd == pVVertE ||
	      pVInit == pVVertE || pVEnd == pVVertD)) {
	// Do T44 swaps whenever the T44 face is first or last in the
	// pipe.
	iReconfTet44(apTCTets, pTF, pVVertA, pVVertB, pVVertC,
		     pVVertD, pVVertE);
	qChange = qReset = true;
      }
      else if (eFC == eT23) {
	// Do T23 swaps at the beginning or end of the pipe
	// OR
	// If the previous and next faces in the pipe share an edge.
	bool qSwap = (pVInit == pVVertD || pVInit == pVVertE ||
		      pVEnd  == pVVertD || pVEnd  == pVVertE);
	if (!qSwap) {
	  // For interior faces of the pipe, check to see whether the
	  // preceding and following pipe faces share an edge.
	  assert(pFLast->qValid());
	  Cell *pCTemp = pFCand->pCCellOpposite(pC);
	  assert(!pCTemp->qHasVert(pVEnd));
	  const Face *pFNext;
	  bool qHitTemp = false;
	  vWhichFace(pCTemp, pFCand, pVInit, pVEnd, pFNext, qHitTemp);
	  assert(pFNext->qValid());
	  if (!qHitTemp) {
	    int iV, iMatch = 0;
	    for (iV = 0; iV < 3; iV++)
	      if (pFNext->qHasVert(pFLast->pVVert(iV))) iMatch++;
	    assert (iMatch == 1 || iMatch == 2);
	    if (iMatch == 2)
	      qSwap = true;
	  }
	}
	if (qSwap) {
	  iReconfTet23(apTCTets, pTF, pVVertA, pVVertB, pVVertC,
		       pVVertD, pVVertE);
	  qChange = qReset = true;
	}
      } // Done with T23 case
      else if (eFC == eT32) {
	// Swap if this is the first face at one end or the other.
	// OR
	// If the previous face was also T32 about the same edge.
	// OR
	// If the previous and subsequent faces share an edge.
	bool qSwap =
	  (pVInit == pVVertA || pVInit == pVVertB ||
	   pVInit == pVVertC || pVEnd  == pVVertA ||
	   pVEnd  == pVVertB || pVEnd  == pVVertC)
	  ||
	  (ePrevFC == eT32 &&
	   ((pVPrevD == pVVertD && pVPrevE == pVVertE) ||
	    (pVPrevD == pVVertE && pVPrevE == pVVertD))) ;
	if (!qSwap) {
	  // For interior faces of the pipe, check to see whether the
	  // preceding and following pipe faces share an edge.
	  assert(pFLast->qValid());
	  Cell *pCTemp = pFCand->pCCellOpposite(pC);
	  assert(!pCTemp->qHasVert(pVEnd));
	  const Face *pFNext;
	  bool qHitTemp = false;
	  vWhichFace(pCTemp, pFCand, pVInit, pVEnd, pFNext, qHitTemp);
	  assert(pFNext->qValid());
	  if (!qHitTemp) {
	    int iV, iMatch = 0;
	    for (iV = 0; iV < 3; iV++)
	      if (pFNext->qHasVert(pFLast->pVVert(iV))) iMatch++;
	    assert (iMatch == 1 || iMatch == 2);
	    if (iMatch == 2)
	      qSwap = true;
	  }
	} // Done checking for previous/next faces sharing an edge
	if (qSwap) {
	  iReconfTet32(apTCTets, pTF, pVVertA, pVVertB, pVVertC,
		       pVVertD, pVVertE);
	  qChange = qReset = true;
	}
      } // Done with T32 case
      else {
	qReset = false;
      }

      ePrevFC = eFC;
      pVPrevD = pVVertD;
      pVPrevE = pVVertE;
    }
    if (qReset) {
      pC = pCBeginPipe(pVInit, pVEnd);
      pFLast = pFInvalidFace;
      qReset = false;
      iNSkip = 0;
      iCount ++;
    }
    else {
      // No swapping, just go on
      apFSkipped[iNSkip] = pFCand;
      iNSkip = iNSkip + 1;
      // Grab the next cell in the pipe
      pFLast = pFCand;
      pC = pFCand->pCCellOpposite(pC);
    }
    // If you made it through and no faces were skipped, you're done

    qEndOfPipe = pC->qHasVert(pVEnd);
    if (qEndOfPipe && iNSkip == 0) {
      vMessage(4, "   Successfully reamed pipe; %d swaps.\n", iCount);
      return 1;
    }
  } while ((!qEndOfPipe) && (iCount < 10) && (iNSkip < 10));
  if (qChange)
    return 0;
  else {
    vMessage(4, "   Failed to ream pipe; %d swaps.\n", iCount);
    return -1;
  }
}

int VolMesh::iRecoverEdge(const Vert* const pV0In,
			  const Vert* const pV1In,
			  SurfMesh& SM,
			  const bool qInsertionOK)
{
  assert(qSimplicial());

  int iNSkip = 0, iRetVal = 0;
  Face* apFSkipped[20];
  bool qHitEdge = false;

  // A modified approach to edge recovery:  volume mesh edges
  // Walk from pV0 towards pV1 until either (a) you hit an edge
  // or (b) you get to pV1.  Keep track of all faces you had to skip
  // over.
  int iRes, iLaps = 0;
  const Vert *pV0 = pV0In;
  const Vert *pV1 = pV1In;
  do {
    iRes = iReamPipe(pV0, pV1, qHitEdge, apFSkipped, iNSkip);
    if (iRes == 1) // Pipe is unclogged
      return 1;
    const Vert *pVTmp = pV0;
    pV0 = pV1;
    pV1 = pVTmp;
    iLaps++;
  } while (iRes == 0 && iLaps <= 6); // Changes made; reverse and try again

  // If an edge was hit, this (almost certainly) means that an edge in
  // the surface mesh and an edge in the volume mesh both lie in the
  // same coplanar facet and that the four points in question are
  // co-circular.  In this situation, the easiest fix is to swap away
  // the offending edge in the -surface- mesh, leaving the volume mesh
  // edge alone.  (While in sufficiently complex degenerate cases, a
  // single swap isn't enough to fix this, multiple swaps (perhaps for
  // different recover-edge targets) will in fact suffice.)

  // After making this forced swap, return with value -1, to indicate
  // that something useful was done, even though an edge is not
  // certain to have been recovered.  Actually, you can do even
  // better: when the edge -has- been recovered, you know that, so you
  // can return 1.

  if (qHitEdge) {
    vMessage(4, "  Hit an edge! Trying to remove it...\n");
    int iV0 = iVertIndex(pV0);
    int iV1 = iVertIndex(pV1);
    Vert *pVS0 = SM.pVVert(iV0);
    Vert *pVS1 = SM.pVVert(iV1);
    assert(dDIST3D(pVS0->adCoords(), pV0->adCoords()) < 1.e-12);
    assert(dDIST3D(pVS1->adCoords(), pV1->adCoords()) < 1.e-12);
    Face *pFS = pFInvalidFace;
    {
      std::set<Vert*> spVNear;
      std::set<Face*> spFNear;
      std::set<Cell*> spCTmp;
      std::set<BFace*> spBFTmp;
      bool qTmp;
      vNeighborhood(pVS0, spCTmp, spVNear, &spBFTmp, &qTmp, &spFNear);
      assert(spVNear.count(pVS1) == 1);
      std::set<Face*>::iterator iterF;
      for (iterF = spFNear.begin(); iterF != spFNear.end(); iterF++) {
	pFS = *iterF;
	if (pFS->qHasVert(pVS0) && pFS->qHasVert(pVS1)) break;
      }
      assert(pFS->qHasVert(pVS0));
      assert(pFS->qHasVert(pVS1));
    }

    // Identify the cells and patches on either side (patches should
    // be the same).
    Cell *pCLeft = pFS->pCCellLeft();
    Cell *pCRight = pFS->pCCellRight();
    BdryPatch3D *pBPLeft = SM.pBdryPatch(SM.iCellIndex(pCLeft));
    BdryPatch3D *pBPRight = SM.pBdryPatch(SM.iCellIndex(pCRight));
    if (pBPLeft != pBPRight) {
      // A surface edge not internal to a facet.

      // Try to swap away this edge in the volume mesh.  First
      // identify the edge.
      Vert *pVA = apFSkipped[0]->pVVert(0);
      Vert *pVB = apFSkipped[0]->pVVert(1);
      Vert *pVC = apFSkipped[0]->pVVert(2);

      Vert *pVNorth, *pVSouth, *pVOther;
      Vert *pVV0 = pVVert(iV0);
      Vert *pVV1 = pVVert(iV1);

      if (iOrient3D(pVA, pVB, pVS0, pVS1) == 0) {
	pVNorth = pVA;
	pVSouth = pVB;
	pVOther = pVC;
      }
      else if (iOrient3D(pVB, pVC, pVS0, pVS1) == 0) {
	pVNorth = pVB;
	pVSouth = pVC;
	pVOther = pVA;
      }
      else {
	assert(iOrient3D(pVC, pVA, pVS0, pVS1) == 0);
	pVNorth = pVC;
	pVSouth = pVA;
	pVOther = pVB;
      }
      int iSwaps = iEdgeSwap3D(apFSkipped[0], pVNorth, pVSouth,
			       pVOther, true, pVV0, pVV1);
      if (iSwaps > 0) {
	// Success!  Well, at least we got rid of the offending edge.
	// So return -1; if the correct edge showed up in the result,
	// we'll figure it out.
	vMessage(4, "   Swapped away the edge in the volume mesh.\n");
	return -1;
      }
      else if (qInsertionOK) {
	// Edge swap failed; have to solve this by splitting.
	double adNewLoc[] = {0.5*(pVS0->dX() + pVS1->dX()),
			     0.5*(pVS0->dY() + pVS1->dY()),
			     0.5*(pVS0->dZ() + pVS1->dZ())};

	Vert *pVNewSurf = SM.pVNewVert();
	pVNewSurf->vSetCoords(3, adNewLoc);
	SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
	iSwaps =  SM.iInsertOnFace(pVNewSurf, pFS, pCLeft, true);
	SUMAA_LOG_EVENT_END(INSERT_VERTEX);

	SUMAA_LOG_EVENT_BEGIN(INSERT_VERTEX);
	qInsertPoint(adNewLoc, apFSkipped[0]->pCCellLeft(),
		     &iSwaps, 1, 1);
	Vert *pVNew = pVVert(iNumVerts() - 1);
	pVNew->vSetType(Vert::eBdryCurve);
	// We may eventually work at trying to remove verts that were
	// added to constrain the boundary.  If/when that is done, we'll
	// be glad that this vert is marked for deletion.
	pVNew->vMarkToDelete();
	SUMAA_LOG_EVENT_END(INSERT_VERTEX);

	vMessage(4, "   Split the edge in both surface and volume meshes.\n");
	// This edge definitely won't need to be recovered now...
	return 1;
      }
      else {
	vMessage(4, "   Not splitting; did nothing.\n");
      }
    } // Different patches on the two sides of the surface face.
    else {
      // Internal facet edge (most common case; four or more
      // cocircular points in a facet boundary)

      // Are the four verts forming those two cells (a) coplanar and (b)
      // cocircular?
      assert(pCLeft->eType() == Cell::eTriCell);
      assert(pCRight->eType() == Cell::eTriCell);
      Vert *pVLeft  = dynamic_cast<TriCell*>(pCLeft )->pVVertOpposite(pFS);
      Vert *pVRight = dynamic_cast<TriCell*>(pCRight)->pVVertOpposite(pFS);

      if (iOrient3D(pV0, pV1, pVLeft, pVRight) != 0) {
	vMessage(4, "    Peculiar.  I was expecting coplanar triangles.\n");
	vMessage(4, "    Nothing I can do here.\n");
	return 0;
      }

      Vert VDummy;
      double adNorm[3];
      vNormal3D(pV0->adCoords(), pV1->adCoords(), pVLeft->adCoords(),
		adNorm);
      vNORMALIZE3D(adNorm);
      double adDummyLoc[] = {pV0->dX() + adNorm[0],
			     pV0->dY() + adNorm[1],
			     pV0->dZ() + adNorm[2]};
      VDummy.vSetCoords(3, adDummyLoc);

      int iInside = iInsphere(pV0, pV1, pVLeft, &VDummy, pVRight);
      if (iInside != 0) {
	vMessage(4, "    Peculiar.  I was expecting cocircular verts.\n");
	vMessage(4, "    Nothing I can do here.\n");
	return 0;
      }

      // Find the offending edge in the VolMesh.  (apFSkipped[0] is a
      // Face that contains that edge, but which edge is it?)
      Vert *pVA = apFSkipped[0]->pVVert(0);
      Vert *pVB = apFSkipped[0]->pVVert(1);
      Vert *pVC = apFSkipped[0]->pVVert(2);
      int iVA = iVertIndex(pVA);
      int iVB = iVertIndex(pVB);
      int iVC = iVertIndex(pVC);

      int iVLeft = SM.iVertIndex(pVLeft);
      int iVRight = SM.iVertIndex(pVRight);

      // First, the easy way: is one of the edges of ABC the same as
      // edge LR in the surface mesh?

      bool qChangeSurf = false, qExactRecovery;
      if ( ((iVLeft == iVA && iVRight == iVB) ||
	    (iVLeft == iVB && iVRight == iVA))
	   ||
	   ((iVLeft == iVB && iVRight == iVC) ||
	    (iVLeft == iVC && iVRight == iVB))
	   ||
	   ((iVLeft == iVC && iVRight == iVA) ||
	    (iVLeft == iVA && iVRight == iVC)) ) {
	// This swap will produce in the surface mesh the edge that
	// exists in the volume mesh.
	qExactRecovery = true;
	qChangeSurf = true;
      }
      else {
	// This swap will eliminate an edge of the surface mesh that
	// intersects the desired edge of the volume mesh, but won't
	// recover the volume edge.
	qExactRecovery = false;

	// Check each possible edge in turn
	bool qABGuilty = ((iOrient3D(pV0, pV1, pVA, pVB) == 0) &&
			  (iInsphere(pV0, pV1, pVLeft, &VDummy, pVA) ==
			   0) &&
			  (iInsphere(pV0, pV1, pVLeft, &VDummy, pVB) ==
			   0));
	bool qBCGuilty = ((iOrient3D(pV0, pV1, pVB, pVC) == 0) &&
			  (iInsphere(pV0, pV1, pVLeft, &VDummy, pVB) ==
			   0) &&
			  (iInsphere(pV0, pV1, pVLeft, &VDummy, pVC) ==
			   0));
	bool qCAGuilty = ((iOrient3D(pV0, pV1, pVC, pVA) == 0) &&
			  (iInsphere(pV0, pV1, pVLeft, &VDummy, pVC) ==
			   0) &&
			  (iInsphere(pV0, pV1, pVLeft, &VDummy, pVA) ==
			   0));
	if (qABGuilty || qBCGuilty || qCAGuilty) {
	  qChangeSurf = true;
	}
      }
      if (qChangeSurf) {
	// Reconfigure the SurfMesh to remove the edge and return.
	SM.vDisallowSwapRecursion();
#ifndef NDEBUG
	int iSwaps =
#endif
	  SM.iReconfigure(dynamic_cast<EdgeFace*>(pFS));
	assert(iSwaps == 1);
	if (qExactRecovery) {
	  vMessage(4, "   Made a co-circular surface swap; should have gotten the edge I wanted.\n");
	  return  1;
	}
	else {
	  vMessage(4, "   Made a co-circular surface swap; didn't get the edge I wanted.\n");
	  return -1;
	}
      }
      else {
	vMessage(4, "   Failed to do anything to the edge.\n");
      }
    } // Edge internal to facet
  } // Hit an edge square on

  return iRetVal;
}

//@ Delete all entities lieing outside the surface mesh
void VolMesh::vDeleteExternal(const SurfMesh* pSM)
{
  assert(qSimplicial());
  assert(pSM->qSimplicial());
  int iCell;
  // Clear the domain tag for each cell in the volume mesh.  At this
  // point, cells are being tagged only as internal or external;
  // multiple subdomains will come later.
  for (iCell = 0; iCell < iNumCells(); iCell++)
    pCCell(iCell)->vSetRegion(iInvalidRegion);

  // Likewise for BdryFaces, so that they will be cleaned up properly.
  // Also, mark the faces attached to these as Face::eBdryFace, so that
  // convex hull faces that are part of the mesh are treated properly.
  for (int iBFace = 0; iBFace < iNumBdryFaces(); iBFace++) {
    BFace *pBF = pBFBFace(iBFace);
    pBF->vSetRegion(iInvalidRegion);
    pBF->vMarkDeleted();
    pBF->pFFace()->vSetFaceLoc(Face::eBdryFace);
  }

  Cell **apCDetach = new Cell*[pSM->iNumCells()];
  Face **apFSurface = new Face*[pSM->iNumCells()];

  //@@ For each face in the surface mesh, find the matching volume face
  for (iCell = 0; iCell < pSM->iNumCells(); iCell++) {
    Cell *pCS = pSM->pCCell(iCell);
    // Note: These are triangular surface mesh cells
    // Find the identities of the vertices in the volume mesh
    Vert* pVV0 = pVVert(pSM->iVertIndex(pCS->pVVert(0)));
    Vert* pVV1 = pVVert(pSM->iVertIndex(pCS->pVVert(1)));
    Vert* pVV2 = pVVert(pSM->iVertIndex(pCS->pVVert(2)));

    //@@@ Find the face that contains all of these verts

    // First, find an edge containing two of them and the info needed to
    // walk around that edge
    Cell *pC = pCBeginPipe(pVV0, pVV1);
    assert(pC->qHasVert(pVV0));
    assert(pC->qHasVert(pVV1));
    Face *pF = pFInvalidFace;
    Vert *pV = pVInvalidVert;
    for (int iF = 0; iF < 4; iF++) {
      pF = pC->pFFace(iF);
      if (pF->qHasVert(pVV0) && pF->qHasVert(pVV1)) break;
    }
    for (int iV = 0; iV < 3; iV++) {
      pV = pF->pVVert(iV);
      if (pV != pVV0 && pV != pVV1) break;
    }

    // Walk around the edge until you hit the correct face
    Face *pFStart = pFInvalidFace;
#ifndef NDEBUG
    bool qTurned = false;
#endif
    while (!pF->qHasVert(pVV2) && pF != pFStart) {
      if (pFStart == pFInvalidFace) pFStart = pF;
      assert(pC->eType() == Cell::eTet);
      TetCell* pTC = dynamic_cast<TetCell*>(pC);
      Face *pFPrev = pF;
      pF = pTC->pFFaceOpposite(pV);
      pV = pTC->pVVertOpposite(pFPrev);
      pC = pF->pCCellOpposite(pC);
      if (pC->eType() == Cell::eTriBFace && !pF->qHasVert(pVV2)) {
	// If you've hit the boundary without finding the vertex you
	// were after, reverse direction and go back the other way.
#ifndef NDEBUG
	assert(!qTurned);
	qTurned = true;
#endif
	pFPrev = pFStart = pF;
	pF = pTC->pFFaceOpposite(pV);
	pV = pTC->pVVertOpposite(pFPrev);
	pC = pF->pCCellOpposite(pTC);
      }
    }

    assert(pF->qHasVert(pVV0));
    assert(pF->qHasVert(pVV1));
    assert(pF->qHasVert(pVV2));

    // Is the cyclic ordering for the volume face and the surface cell
    // the same, or opposite?
    double adNormVol[3];
    pF->vNormal(adNormVol);
    vNORMALIZE3D(adNormVol);
    assert(pCS->eType() == Cell::eTriCell);
    double adVec01[] = adDIFF3D(pCS->pVVert(1)->adCoords(),
				pCS->pVVert(0)->adCoords());
    double adVec02[] = adDIFF3D(pCS->pVVert(2)->adCoords(),
				pCS->pVVert(0)->adCoords());
    double adNormSurf[3];
    vCROSS3D(adVec01, adVec02, adNormSurf);
    vNORMALIZE3D(adNormSurf);

    double adNormPatch[3];
    BdryPatch3D *pBP3D = pSM->pBdryPatch(iCell);
    pBP3D->vUnitNormal(pCS->pVVert(0)->adCoords(), adNormPatch);

    double dPatchDot = -dDOT3D(adNormVol, adNormPatch);
    double dSurfDot  = -dDOT3D(adNormVol, adNormSurf);

    // Mark the interior and exterior cells adjacent to this face
    int iRightHandPatch = iFuzzyComp(dPatchDot,0.);
    int iRightHandSurf  = iFuzzyComp(dSurfDot,0.);
    assert(abs(iRightHandPatch) == 1);
    assert(abs(iRightHandSurf ) == 1);

    // These next few lines are relevant eventually for getting
    // multi-region support up and going again.
    int iRightRegion = pBP3D->iRightRegion();
    int iLeftRegion  = pBP3D->iLeftRegion ();

    bool qIsInteriorBoundary = false;
    if (iRightRegion != iOutsideRegion &&
	iLeftRegion != iOutsideRegion)
      qIsInteriorBoundary = true;

    apFSurface[iCell] = pF;

    // Topology-based tagging isn't working as well as I'd like, so I'm
    // going to try geometry-based tagging.

    double adFaceCent[3], adCellCent[3];
    pF->vCentroid(adFaceCent);
    pF->pCCellLeft()->vCentroid(adCellCent);
    double adLeftVec[] = adDIFF3D(adCellCent, adFaceCent);
    vNORMALIZE3D(adLeftVec);
    double dDot = dDOT3D(adLeftVec, adNormSurf);
    if (qIsInteriorBoundary) {
      // For now, use the old topology-based scheme, which -seemed- to
      // work okay.
      assert(iRightRegion != iOutsideRegion);
      assert(iLeftRegion  != iOutsideRegion);
      apCDetach[iCell] = pCInvalidCell;
      if (iRightHandPatch == 1) {
	pF->pCCellRight()->vSetRegion(iRightRegion);
	pF->pCCellLeft() ->vSetRegion(iLeftRegion);
      }
      else {
	pF->pCCellRight()->vSetRegion(iLeftRegion);
	pF->pCCellLeft() ->vSetRegion(iRightRegion);
      }
    }
    else {
      // If dDot < 0, then for sure this cell lies inside, because the
      // surface mesh is oriented this way on purpose.
      assert(iFuzzyComp(dDot, 0) != 0);
      int iIntRegion = ((iRightRegion == iOutsideRegion) ?
			iLeftRegion : iRightRegion);
      vMessage(4, "Face: %3d.  Verts %3d, %3d, %3d.\n",
	       iFaceIndex(pF), iVertIndex(pF->pVVert(0)),
	       iVertIndex(pF->pVVert(1)), iVertIndex(pF->pVVert(2)));
      vMessage(4, "  Vert %3d @ (%10.6G, %10.6G, %10.6G)\n",
	       iVertIndex(pF->pVVert(0)), pF->pVVert(0)->dX(),
	       pF->pVVert(0)->dY(), pF->pVVert(0)->dZ());
      vMessage(4, "  Vert %3d @ (%10.6G, %10.6G, %10.6G)\n",
	       iVertIndex(pF->pVVert(1)), pF->pVVert(1)->dX(),
	       pF->pVVert(1)->dY(), pF->pVVert(1)->dZ());
      vMessage(4, "  Vert %3d @ (%10.6G, %10.6G, %10.6G)\n",
	       iVertIndex(pF->pVVert(2)), pF->pVVert(2)->dX(),
	       pF->pVVert(2)->dY(), pF->pVVert(2)->dZ());
      vMessage(4, "Surface mesh normal: (%10.6G, %10.6G, %10.6G)\n",
	       adNormSurf[0], adNormSurf[1], adNormSurf[2]);
      vMessage(4, "Volume mesh normal:  (%10.6G, %10.6G, %10.6G) Dot: %8.3G\n",
	       adNormVol[0], adNormVol[1], adNormVol[2],
	       (dDOT3D(adNormSurf, adNormVol)));
      vMessage(4, "Patch normal:        (%10.6G, %10.6G, %10.6G) Dot: %8.3G\n",
	       adNormPatch[0], adNormPatch[1], adNormPatch[2],
	       (dDOT3D(adNormSurf, adNormPatch)));
      if (dDot < 0) {
	pC = pF->pCCellRight();
	vMessage(4, "Detaching cell %3d. Verts %3d, %3d, %3d, %3d.\n",
		 iCellIndex(pC), iVertIndex(pC->pVVert(0)),
		 iVertIndex(pC->pVVert(1)), iVertIndex(pC->pVVert(2)),
		 iVertIndex(pC->pVVert(3)));
	apCDetach[iCell] = pF->pCCellRight();
	pF->pCCellRight()->vSetRegion(iOutsideRegion);
	pF->pCCellLeft() ->vSetRegion(iIntRegion);
      }
      else {
	pC = pF->pCCellLeft();
	vMessage(4, "Detaching cell %3d. Verts %3d, %3d, %3d, %3d.\n",
		 iCellIndex(pC), iVertIndex(pC->pVVert(0)),
		 iVertIndex(pC->pVVert(1)), iVertIndex(pC->pVVert(2)),
		 iVertIndex(pC->pVVert(3)));
	apCDetach[iCell] = pF->pCCellLeft();
	pF->pCCellLeft() ->vSetRegion(iOutsideRegion);
	pF->pCCellRight()->vSetRegion(iIntRegion);
      }
    }

       // Old topology-based tagging; fails mysteriously for at least
       // the hex-multi case; more if you try to fix it. ;-)
//     // Tag cells as inside/outside based on the patch orientation
//     if (iRightHandPatch == 1) {
//       pF->pCCellRight()->vSetRegion(iRightRegion);
//       pF->pCCellLeft() ->vSetRegion(iLeftRegion);
//     }
//     else {
//       pF->pCCellRight()->vSetRegion(iLeftRegion);
//       pF->pCCellLeft() ->vSetRegion(iRightRegion);
//     }

//     // Tag cell to be detached from the mesh based on surface mesh
//     // orientation
//     if (qIsInteriorBoundary) {
//       assert(iRightRegion != iOutsideRegion);
//       assert(iLeftRegion  != iOutsideRegion);
//       apCDetach[iCell] = pCInvalidCell;
//     }
//     else {
//       // There's something wrong with the following condition.  It works
//       // nearly all the time, but for the hex-multi case, it causes
//       // problems.
//       if (iFuzzyComp(dPatchDot, dSurfDot) == 0) {
//	// When the surface mesh and the patch are oriented the same
//	// way, there are no problems.
//	if (iRightHandPatch == 1) {
//	  apCDetach[iCell] = pF->pCCellLeft();
//	}
//	else {
//	  apCDetach[iCell] = pF->pCCellRight();
//	}
//       }
//       else {
//	// In this case, do a geometric check.
//	double adFaceCent[3], adCellCent[3];
//	pF->vCentroid(adFaceCent);
//	pF->pCCellLeft()->vCentroid(adCellCent);
//	double adLeftVec[] = adDIFF3D(adCellCent, adFaceCent);
//	vNORMALIZE3D(adLeftVec);
//	double dDot = dDOT3D(adLeftVec, adNormSurf);
//	// If dDot > 0, then for sure this cell lies inside, because the
//	// surface mesh is oriented this way on purpose.
//	assert(iFuzzyComp(dDot, 0) != 0);
//	if (dDot > 0) {
//	  apCDetach[iCell] = pF->pCCellRight();
//	}
//	else {
//	  apCDetach[iCell] = pF->pCCellLeft();
//	}
//       }
//     }
  } // End of loop associating surface faces with volume faces.

#ifndef NDEBUG
  vMessage(4, "Initial markings...\n");
  double dSum = 0;
  double adSum[iOutsideRegion+1];
  for (int iiR = 0; iiR <= iOutsideRegion; adSum[iiR++] = 0) {}
  for (iCell = 0; iCell < iNumCells(); iCell++) {
    Cell *pC = pCCell(iCell);
    if (pC->qDeleted()) continue;
    double dSize = pC->dSize();
    int iReg = pC->iRegion();
    if (iReg != iOutsideRegion) dSum += dSize;
    adSum[iReg] += dSize;
    vMessage(4, "Cell %3d. Verts %3d, %3d, %3d, %3d.  Size: %10.4g (%10.4g) Region %d (%10.4g).\n",
	     iCell, iVertIndex(pC->pVVert(0)),
	     iVertIndex(pC->pVVert(1)), iVertIndex(pC->pVVert(2)),
	     iVertIndex(pC->pVVert(3)), dSize, dSum, iReg, adSum[iReg]);
  }
#endif

  // For each marked cell in the volume mesh, recursively mark unmarked
  // neighbors
  for (iCell = 0; iCell < iNumCells(); iCell++) {
    Cell *pC = pCCell(iCell);
    if (pC->iRegion() != iInvalidRegion) vMarkCleanNeighbors(pC);
  }

#ifndef NDEBUG
  vMessage(4, "Intermediate markings...\n");
  for (int iiR = 0; iiR <= iOutsideRegion; adSum[iiR++] = 0) {}
  for (iCell = 0; iCell < iNumCells(); iCell++) {
    Cell *pC = pCCell(iCell);
    if (pC->qDeleted()) continue;
    double dSize = pC->dSize();
    int iReg = pC->iRegion();
    if (iReg != iOutsideRegion) dSum += dSize;
    adSum[iReg] += dSize;
    vMessage(4, "Cell %3d. Verts %3d, %3d, %3d, %3d.  Size: %10.4g (%10.4g) Region %d (%10.4g).\n",
	     iCell, iVertIndex(pC->pVVert(0)),
	     iVertIndex(pC->pVVert(1)), iVertIndex(pC->pVVert(2)),
	     iVertIndex(pC->pVVert(3)), dSize, dSum, iReg, adSum[iReg]);
  }
#endif

  // Create the actual boundary that you'll need.
  for (iCell = 0; iCell < pSM->iNumCells(); iCell++) {
    Face *pF = apFSurface[iCell];
    Cell *pC = apCDetach [iCell];

    if (pC->qValid()) {
      // This is the regular boundary case.
      if (pC->eType() != Cell::eTriBFace &&
	  pC->eType() != Cell::eQuadBFace) {
	// Detach the outside cell from the face on the boundary
	pC->vRemoveFace(pF);
	pC->vMarkDeleted();
	// Create a new boundary face and attach it to the face
	BFace *pBF = pBFNewBFace();
	pBF->vAssign(pF);
	assert(pBF->eType() == Cell::eTriBFace);
	pBF->vSetPatch(pSM->pBdryPatch(iCell));
	pF->vReplaceCell(pC, pBF);
	pF->vSetFaceLoc(Face::eBdryFace);
      }
      else {
	// Already have a BFace as the cell.
	assert(pF->iFaceLoc() == Face::eBdryFace);
	pC->vMarkNotDeleted();
	pC->vSetRegion(iInsideRegion);
      }
    }
    else {
      // Internal boundary.  Need to duplicate the face and add an
      // internal bdry face between the two copies.
      assert(pF->eType() == Face::eTriFace);
      Cell *pCL = pF->pCCellLeft();
      Cell *pCR = pF->pCCellRight();
      assert(pCL->iFullCheck() == 1);
      assert(pCR->iFullCheck() == 1);

      Face *pFNew = pFNewFace();
      pFNew->vAssign(pCL, pCR, pF->pVVert(0), pF->pVVert(1),
		     pF->pVVert(2));
      pFNew->vCopyAllFlags(*pF);
      pFNew->vSetFaceLoc(Face::eBdryTwoSide);
      pF->vSetFaceLoc(Face::eBdryTwoSide);

      // Create a new internal tri boundary face
      BFace *pBFInt = pBFNewBFace(3, true);
      pBFInt->vAddFace(pF);
      pBFInt->vAddFace(pFNew);
      pBFInt->vSetPatch(pSM->pBdryPatch(iCell));

      // Update cell connectivity
      pCR->vReplaceFace(pF, pFNew);
      pF->vReplaceCell(pCR, pBFInt);
      pFNew->vReplaceCell(pCL, pBFInt);

      assert(pCL->iFullCheck() == 1);
      assert(pCR->iFullCheck() == 1);
      assert(pBFInt->iFullCheck() == 1);
      assert(pF->iFullCheck());
      assert(pFNew->iFullCheck());
    }
  }
  delete [] apFSurface;

  delete [] apCDetach;

#ifndef NDEBUG
  vMessage(4, "Final markings...\n");
  dSum = 0;
  for (int iiR = 0; iiR <= iOutsideRegion; adSum[iiR++] = 0) {}
  for (iCell = 0; iCell < iNumCells(); iCell++) {
    Cell *pC = pCCell(iCell);
    if (pC->qDeleted()) continue;
    double dSize = pC->dSize();
    int iReg = pC->iRegion();
    if (iReg != iOutsideRegion) dSum += dSize;
    adSum[iReg] += dSize;
    vMessage(4, "Cell %3d. Verts %3d, %3d, %3d, %3d.  Size: %10.4g (%10.4g) Region %d (%10.4g).\n",
	     iCell, iVertIndex(pC->pVVert(0)),
	     iVertIndex(pC->pVVert(1)), iVertIndex(pC->pVVert(2)),
	     iVertIndex(pC->pVVert(3)), dSize, dSum, iReg, adSum[iReg]);
  }
#endif

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

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

  int iFace;
  for (iFace = 0; iFace < iNumFaces(); iFace++) {
    Face *pF = pFFace(iFace);
    if (pF->pCCellLeft()->qDeleted() &&
	pF->pCCellRight()->qDeleted())
      pF->vMarkDeleted();
  }

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

  vPurge();
  assert(qValid());

#ifndef NDEBUG
  for (iCell = 0; iCell < iNumCells(); iCell++)
    assert(pCCell(iCell)->iRegion() != 0);
#endif
}

void VolMesh::vFindFaceContaining(const Vert * const pV0,
				  const Vert * const pV1,
				  const Vert * const pV2,
				  Face *&pF, bool &qSuccess)
{
  qSuccess = false;

  //@ First, find an edge connecting two of the verts.
  Cell *pC = pCBeginPipe(pV0, pV1);

  // A failure exit, because pV0 and pV1 do not share an edge in the
  // volume mesh.
  if (! (pC->qHasVert(pV0) && pC->qHasVert(pV1)) ) return;
  assert(pC->qHasVert(pV0));
  assert(pC->qHasVert(pV1));

  if (pC->qHasVert(pV2)) {
    qSuccess = true;
    for (int iF = 0; iF < pC->iNumFaces(); iF++) {
      pF = pC->pFFace(iF);
      if (pF->qHasVert(pV0) && pF->qHasVert(pV1) && pF->qHasVert(pV2))
    return;
  }
    assert(0);
  }

  // Two of the three verts are in pC.  Walk around this edge to find a
  // cell that contains pV2.
  pF = pFInvalidFace;
  Vert *pV = pVInvalidVert;
  for (int iF = 0; iF < 4; iF++) {
    pF = pC->pFFace(iF);
    if (pF->qHasVert(pV0) && pF->qHasVert(pV1)) break;
  }
  for (int iV = 0; iV < 3; iV++) {
    pV = pF->pVVert(iV);
    if (pV != pV0 && pV != pV1) break;
  }

  // Walk around the edge until you hit the correct face
  Face *pFStart = pFInvalidFace;
#ifndef NDEBUG
  bool qTurned = false;
#endif
  while (!pF->qHasVert(pV2) && pF != pFStart) {
    if (pFStart == pFInvalidFace) pFStart = pF;
    assert(pC->eType() == Cell::eTet);
    TetCell* pTC = dynamic_cast<TetCell*>(pC);
    Face *pFPrev = pF;
    pF = pTC->pFFaceOpposite(pV);
    pV = pTC->pVVertOpposite(pFPrev);
    pC = pF->pCCellOpposite(pC);
    if (pC->eType() == Cell::eTriBFace && !pF->qHasVert(pV2)) {
      // If you've hit the boundary without finding the vertex you
      // were after, reverse direction and go back the other way.
#ifndef NDEBUG
      assert(!qTurned);
      qTurned = true;
#endif
      pFPrev = pFStart = pF;
      pF = pTC->pFFaceOpposite(pV);
      pV = pTC->pVVertOpposite(pFPrev);
      pC = pF->pCCellOpposite(pTC);
    }
  }

  assert(pF->qHasVert(pV0));
  assert(pF->qHasVert(pV1));
  assert(pF->qHasVert(pV2));
  assert(pC->qHasFace(pF));
  qSuccess = true;
}
