#include "iMesh.h"
#include "iMesh_GRUMMP.hh"
using namespace ITAPS_GRUMMP;

#include <iostream>
#include <vector>
#include "GR_Entity.h"
#include "GR_Geometry.h"
#include "GR_VolMesh.h"
#include "iMesh_GRUMMP_misc.hh"
#include "iMesh_GRUMMP_EntIter.hh"
#include "iMesh_GRUMMP_Workset.hh"

// Unless otherwise indicated, these funcs should be 3D-specific, in
// iMesh_Volume

void iMesh_Volume::vCreateMesh()
{
  VolMesh *pVM = new VolMesh();
  pM = pVM;
}

void iMesh_Volume::load(const char *name, const char *, int* err)
{
  *err = iBase_SUCCESS;
  VolMesh *pVM = new VolMesh(name);
  pM = pVM;
  // Set up the initial (GRUMMP) face hash table.
  for (int iF = pM->iNumFaces() - 1; iF >= 0; iF--) {
    vAddFaceToHash(pM->pFFace(iF));
  }
  iNVert = pM->iNumVerts();
  iNTri = pVM->iNumTriFaces();
  iNQuad = pVM->iNumQuadFaces();
  iNTet = pVM->iNumTetCells();
  iNPyr = pVM->iNumPyrCells();
  iNPrism = pVM->iNumPrismCells();
  iNHex = pVM->iNumHexCells();
  iNLineSeg = iNPolygon = iNSept = iNPolyhedron = 0;
  assert(iNTri + iNQuad == (mmFaceHash.size()));
}

void iMesh_Volume::save (const char *name, const char *, int* err)
{
  *err = iBase_SUCCESS;
  VolMesh *pVM = dynamic_cast<VolMesh*>(pM);
  vWriteFile_VolMesh(*pVM, name);
}

void iMesh_Volume::getAdjTable (int* adjacency_table)
{
  // Array is known to be okay before this call.

  // The contents of this adjacency table should be:
  //    V E F R
  //  V I U L L
  //  E U U U U
  //  F I U I I
  //  R I U I I
  // where I = immediate
  //       L = local traversal
  //       U = unavailable

  // The adjacency info available in GRUMMP is always the same.

  adjacency_table[0] = iBase_ALL_ORDER_1;
  adjacency_table[1] = iBase_UNAVAILABLE;
  adjacency_table[2] = iBase_ALL_ORDER_1;
  adjacency_table[3] = iBase_ALL_ORDER_1;

  adjacency_table[4] = iBase_UNAVAILABLE;
  adjacency_table[5] = iBase_UNAVAILABLE;
  adjacency_table[6] = iBase_UNAVAILABLE;
  adjacency_table[7] = iBase_UNAVAILABLE;

  adjacency_table[8] = iBase_ALL_ORDER_1;
  adjacency_table[9] = iBase_UNAVAILABLE;
  adjacency_table[10]= iBase_ALL_ORDER_1;
  adjacency_table[11]= iBase_ALL_ORDER_1;

  adjacency_table[12]= iBase_ALL_ORDER_1;
  adjacency_table[13]= iBase_UNAVAILABLE;
  adjacency_table[14]= iBase_ALL_ORDER_1;
  adjacency_table[15]= iBase_ALL_ORDER_1;
}

void iMesh_Volume::getEntities(const int entity_type,
				      const int entity_topology,
				      iBase_EntityHandle** entity_handles,
				      int* entity_handles_allocated,
				      int* entity_handles_size,
				      int *err)

{
  *err = iBase_SUCCESS;

  VolMesh *pVM = dynamic_cast<VolMesh*>(pM);

  if (entity_topology == iMesh_POINT ||
      (entity_topology == iMesh_ALL_TOPOLOGIES &&
       entity_type     == iBase_VERTEX) ) {
    ////////////////////////////////////////////////
    // Returning all vertices.
    ////////////////////////////////////////////////
    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNVert,
		   iBase_EntityHandle);

    unsigned int iNV = pM->iNumVerts(); // This is >= iNVert
    unsigned int iNOrigV = iNV;
    unsigned int iV, iEH;
    for (iV = 0, iEH = 0;
	 iV < iNOrigV; iV++) {
      Vert *pV = pM->pVVert(iV);
      if (pV->qDeleted()) {
	iNV--;
      }
      else {
	(*entity_handles)[iEH] = pV;
	iEH++;
      }
    }
    assert(iEH == iNVert);
    assert(iNV == iNVert);
    assert(iNV <= iNOrigV);
    *entity_handles_size = iEH;
  }
  else if (entity_topology == iMesh_POLYGON) {
    ////////////////////////////////////////////////
    // Returning all polygonal faces, of which there are zero.
    ////////////////////////////////////////////////
    int iNPoly = 0;
    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNPoly,
		   iBase_EntityHandle);
    *entity_handles_size = 0;
  }
  else if (entity_topology == iMesh_TRIANGLE) {
    ////////////////////////////////////////////////
    // Returning all triangular faces.
    ////////////////////////////////////////////////
    int iNTriTemp = 0;
    iNTriTemp = pVM->iNumTriFaces();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNTriTemp,
		   iBase_EntityHandle);

    int iTri, iEH;
    for (iTri = 0, iEH = 0;
	 iTri < iNTriTemp; iTri++) {
      Face *pF = pM->pFFace(iTri);
      if (pF->qDeleted()) {
	iNTriTemp--;
      }
      else {
	(*entity_handles)[iEH] = pF;
	iEH++;
      }
    }
    *entity_handles_size = iNTriTemp;
    assert(iEH  == *entity_handles_size);
  } // Done with triangular faces.
  else if (entity_topology == iMesh_QUADRILATERAL) {
    ////////////////////////////////////////////////
    // Returning all quadrilateral faces.
    ////////////////////////////////////////////////
    int iNQuadTemp = 0, iNTriTemp = 0;
    iNTriTemp = pVM->iNumTriFaces();
    iNQuadTemp = pVM->iNumQuadFaces();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNQuadTemp,
		   iBase_EntityHandle);

    int iQuad, iEH;
    for (iQuad = iNTriTemp, iEH = 0;
	 iQuad < iNQuadTemp + iNTriTemp; iQuad++) {
      Face *pF = pVM->pFFace(iQuad);
      if (pF->qDeleted()) {
	iNQuadTemp--;
      }
      else {
	(*entity_handles)[iEH] = pF;
	iEH++;
      }
    }
    *entity_handles_size = iNQuadTemp;
    assert(iEH == *entity_handles_size);
  } // Done with quadrilateral cells.
  else if (entity_topology == iMesh_ALL_TOPOLOGIES &&
	   entity_type == iBase_FACE) {
    ////////////////////////////////////////////////
    // Returning all faces.
    ////////////////////////////////////////////////
    int iNFace = 0;
    iNFace = pVM->iNumFaces();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNFace, iBase_EntityHandle);

    int iFace, iEH;
    for (iFace = 0, iEH = 0;
	 iFace < iNFace; iFace++) {
      Face *pF = pVM->pFFace(iFace);
      if (pF->qDeleted()) {
	iNFace--;
      }
      else {
	(*entity_handles)[iEH] = pF;
	iEH++;
      }
    }
    *entity_handles_size = iNFace;
    assert(iEH == *entity_handles_size);
  } // Done with all faces.
  else if (entity_topology == iMesh_TETRAHEDRON) {
    ////////////////////////////////////////////////
    // Returning all tets.
    ////////////////////////////////////////////////
    int iNTetTemp = pVM->iNumTetCells();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNTetTemp,
		   iBase_EntityHandle);

    int iTet, iEH;
    for (iTet = 0, iEH = 0;
	 iTet < iNTetTemp; iTet++) {
      Cell *pC = pVM->pCCell(iTet);
      if (pC->qDeleted()) {
	iNTetTemp--;
      }
      else {
	(*entity_handles)[iEH] = pC;
	iEH++;
      }
    }
    *entity_handles_size = iNTetTemp;
    assert(iEH == *entity_handles_size);
  }
  else if (entity_topology == iMesh_PYRAMID) {
    ////////////////////////////////////////////////
    // Returning all pyramids.
    ////////////////////////////////////////////////
    int iNPyrTemp = pVM->iNumPyrCells();
    int iBegin = pVM->iNumTetCells();
    int iEnd = iBegin + iNPyrTemp;

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNPyrTemp, iBase_EntityHandle);

    int iPyr, iEH;
    for (iPyr = iBegin, iEH = 0;
	 iPyr < iEnd; iPyr++) {
      Cell *pC = pVM->pCCell(iPyr);
      if (pC->qDeleted()) {
	iNPyrTemp--;
      }
      else {
	(*entity_handles)[iEH] = pC;
	iEH++;
      }
    }
    *entity_handles_size = iNPyrTemp;
    assert(iEH == *entity_handles_size);
  }
  else if (entity_topology == iMesh_PRISM) {
    ////////////////////////////////////////////////
    // Returning all prisms.
    ////////////////////////////////////////////////
    int iNPrismTemp = pVM->iNumPrismCells();
    int iBegin = pVM->iNumTetCells() + pVM->iNumPyrCells();
    int iEnd = iBegin + iNPrismTemp;

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNPrismTemp, iBase_EntityHandle);

    int iPrism, iEH;
    for (iPrism = iBegin, iEH = 0;
	 iPrism < iEnd; iPrism++) {
      Cell *pC = pVM->pCCell(iPrism);
      if (pC->qDeleted()) {
	iNPrismTemp--;
      }
      else {
	(*entity_handles)[iEH] = pC;
	iEH++;
      }
    }
    *entity_handles_size = iNPrismTemp;
    assert(iEH == *entity_handles_size);
  }
  else if (entity_topology == iMesh_HEXAHEDRON) {
    ////////////////////////////////////////////////
    // Returning all hexs.
    ////////////////////////////////////////////////
    int iNHexTemp = pVM->iNumHexCells();
    int iBegin = pVM->iNumTetCells() + pVM->iNumPyrCells()
      + pVM->iNumPrismCells();
    int iEnd = iBegin + iNHexTemp;

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNHexTemp, iBase_EntityHandle);

    int iHex, iEH;
    for (iHex = iBegin, iEH = 0;
	 iHex < iEnd; iHex++) {
      Cell *pC = pVM->pCCell(iHex);
      if (pC->qDeleted()) {
	iNHexTemp--;
      }
      else {
	(*entity_handles)[iEH] = pC;
	iEH++;
      }
    }
    *entity_handles_size = iNHexTemp;
    assert(iEH == *entity_handles_size);
  }
  else if (entity_topology == iMesh_ALL_TOPOLOGIES &&
	   entity_type == iBase_REGION) {
    ////////////////////////////////////////////////
    // Returning all cells!
    ////////////////////////////////////////////////
    int iNCell = pVM->iNumCells();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNCell, iBase_EntityHandle);

    int iCell, iEH;
    for (iCell = 0, iEH = 0;
	 iCell < iNCell; iCell++) {
      Cell *pC = pVM->pCCell(iCell);
      if (pC->qDeleted()) {
	iNCell--;
      }
      else {
	(*entity_handles)[iEH] = pC;
	iEH++;
      }
    }
    *entity_handles_size = iNCell;
    assert(iEH == *entity_handles_size);
  } // Done with all cells.
  else if (entity_topology == iMesh_ALL_TOPOLOGIES &&
	   entity_type == iBase_ALL_TYPES) {
    ////////////////////////////////////////////////
    // Returning all entities!
    ////////////////////////////////////////////////

    vCopyRootSet( entity_handles, entity_handles_allocated,
		       entity_handles_size);
    if (*err != iBase_SUCCESS) return;
  }
  else if (entity_topology == iMesh_POLYHEDRON ||
	   entity_topology == iMesh_SEPTAHEDRON ||
	   entity_topology == iMesh_LINE_SEGMENT ||
	   (entity_topology == iMesh_ALL_TOPOLOGIES &&
	    entity_type     == iBase_EDGE) ) {
    ////////////////////////////////////////////////
    // There aren't any of these.
    ////////////////////////////////////////////////
    int iN = 0;
    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iN, iBase_EntityHandle);
    // Nothing else to do, except return iN when I can.
    *entity_handles_size = 0;
  }
}

void iMesh_Volume::createEntArr
(/*in*/ const int new_entity_topology,
 /*in*/ const iBase_EntityHandle* lower_order_entity_handles,
 /*in*/ const int lower_order_entity_handles_size,
 /*out*/ iBase_EntityHandle** new_entity_handles,
 /*out*/ int* new_entity_handles_allocated,
 /*out*/ int* new_entity_handles_size,
 /*inout*/ int** status,
 /*inout*/ int* status_allocated,
 /*out*/ int* status_size, int *err)
{
  *err = iBase_SUCCESS;

  *new_entity_handles_size = 0;
  *status_size = 0;
  int vertsPerEnt;
  switch (new_entity_topology) {
  case iMesh_POINT:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    return;
    break;
  case iMesh_TRIANGLE:
    vertsPerEnt = 3;
    break;
  case iMesh_QUADRILATERAL:
    vertsPerEnt = 4;
    break;
  case iMesh_TETRAHEDRON:
    vertsPerEnt = 4;
    break;
  case iMesh_PRISM:
    vertsPerEnt = 6;
    break;
  case iMesh_PYRAMID:
    vertsPerEnt = 5;
    break;
  case iMesh_HEXAHEDRON:
    vertsPerEnt = 8;
    break;
  case iMesh_SEPTAHEDRON:
  case iMesh_POLYHEDRON:
  case iMesh_LINE_SEGMENT:
  case iMesh_POLYGON:
    *err = iBase_NOT_SUPPORTED;
    break;
  case iMesh_ALL_TOPOLOGIES:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    break;
  } // End switch over topologies

  // Now set up arrays and call createEnt to do the actual work
  if (lower_order_entity_handles_size % vertsPerEnt != 0) {
    *err = iBase_INVALID_ENTITY_COUNT;
    return;
  }

  (*status_size) = (*new_entity_handles_size) =
    lower_order_entity_handles_size/vertsPerEnt;
  TRY_ARRAY_SIZE(new_entity_handles, new_entity_handles_allocated,
		 *new_entity_handles_size, iBase_EntityHandle);
  TRY_ARRAY_SIZE(status, status_allocated, *status_size, int);

  if (*err != iBase_SUCCESS) return;

  for (int iLO = 0, iEH = 0, iStat = 0;
       iLO <= lower_order_entity_handles_size - 1;
       iEH++, iStat++, iLO += vertsPerEnt) {
    createEnt(new_entity_topology,
	      lower_order_entity_handles + iLO, vertsPerEnt,
	      (*new_entity_handles)+iEH, (*status) + iStat, err);
  }
}

void iMesh_Volume::createEnt
  (/*in*/ const int new_entity_topology,
   /*in*/ const iBase_EntityHandle* lower_order_entity_handles,
   /*in*/ const int lower_order_entity_handles_size,
   /*out*/ iBase_EntityHandle* new_entity_handle,
   /*out*/ int* status, int *err)
{
  VolMesh *pVM = dynamic_cast<VolMesh*>(pM);

  switch (new_entity_topology) {
  case iMesh_POINT:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    *new_entity_handle = NULL;
    *status = iBase_CREATION_FAILED;
    break;
  case iMesh_TRIANGLE:
    {
      // We're going to create a new TriFace!  That means three
      // vertices.
      if (lower_order_entity_handles_size != 3) {
	*err = iBase_INVALID_ENTITY_COUNT;
	return;
      }

      Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[0]));
      Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[1]));
      Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[2]));

      if (!pV0->qValid() || !pV1->qValid() || !pV2->qValid()) {
	*err = iBase_INVALID_ENTITY_HANDLE;
	*new_entity_handle = NULL;
	*status = iBase_CREATION_FAILED;
	break;
      }

      Vert *apV[] = {pV0, pV1, pV2};
      Face *pF = NULL;
      if (qDoesFaceExist(3, apV, pF)) {
	// Return this face as a handle, with status ALREADY_EXISTED
	assert(pF->qValid());
	*new_entity_handle = pF;
	*status = iBase_ALREADY_EXISTED;
      }
      else {
	// Create it.
	TriFace* pTF = dynamic_cast<TriFace*>(pVM->pFNewFace());
	if (!pTF->qValid()) {
	  *err = iBase_ENTITY_CREATION_ERROR;
	  *new_entity_handle = NULL;
	  *status = iBase_CREATION_FAILED;
	  break;
	}

	pTF->vAssign(pCInvalidCell, pCInvalidCell, pV0, pV1, pV2);
	pV0->vSetHintFace(pTF);
	pV1->vSetHintFace(pTF);
	pV2->vSetHintFace(pTF);
	*new_entity_handle = pTF;
	*status =  iBase_NEW;

	vAddFaceToHash(pTF);
	iNTri++;
	assert(iNTri <= static_cast<unsigned int>(pVM->iNumTriFaces()));
	assert(iNTri + iNQuad == (mmFaceHash.size()));
      }
    }
    break;
  case iMesh_QUADRILATERAL:
    {
      // We're going to create a new QuadFace!  That means four
      // vertices
      if (lower_order_entity_handles_size != 4) {
	*err = iBase_INVALID_ENTITY_COUNT;
	break;
      }

      Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[0]));
      Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[1]));
      Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[2]));
      Vert *pV3 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[3]));

      if (!pV0->qValid() || !pV1->qValid() || !pV2->qValid()
	  || !pV3->qValid()) {
	*err = iBase_INVALID_ENTITY_HANDLE;
	*new_entity_handle = NULL;
	*status = iBase_CREATION_FAILED;
	break;
      }

      Vert *apV[] = {pV0, pV1, pV2, pV3};
      Face *pF = NULL;
      if (qDoesFaceExist(4, apV, pF)) {
	// Return this face as a handle, with status ALREADY_EXISTED
	assert(pF->qValid());
	*new_entity_handle = pF;
	*status = iBase_ALREADY_EXISTED;
      }
      else {
	// Create it.
	QuadFace* pQF = dynamic_cast<QuadFace*>(pVM->pFNewFace(4));
	if (!pQF->qValid()) {
	  *err = iBase_ENTITY_CREATION_ERROR;
	  *new_entity_handle = NULL;
	  *status = iBase_CREATION_FAILED;
	  break;
	}

	pQF->vAssign(pCInvalidCell, pCInvalidCell, pV0, pV1, pV2, pV3);
	pV0->vSetHintFace(pQF);
	pV1->vSetHintFace(pQF);
	pV2->vSetHintFace(pQF);
	pV3->vSetHintFace(pQF);
	*new_entity_handle = pQF;
	*status = iBase_NEW;

	vAddFaceToHash(pQF);
	iNQuad++;
	assert(iNQuad <= static_cast<unsigned int>(pVM->iNumQuadFaces()));
	assert(iNTri + iNQuad == (mmFaceHash.size()));
      }
    }
    break;
  case iMesh_TETRAHEDRON:
    {
      // New tets to be built from faces, never from vertices or edges.
      if (lower_order_entity_handles_size != 4) {
	*err = iBase_INVALID_ENTITY_COUNT;
	break;
      }

      Face *apF[4];
      apF[0] = dynamic_cast<Face*>(static_cast<Entity*>
				   (lower_order_entity_handles[0]));
      apF[1] = dynamic_cast<Face*>(static_cast<Entity*>
				   (lower_order_entity_handles[1]));
      apF[2] = dynamic_cast<Face*>(static_cast<Entity*>
				   (lower_order_entity_handles[2]));
      apF[3] = dynamic_cast<Face*>(static_cast<Entity*>
				   (lower_order_entity_handles[3]));

      if (!apF[0]->qValid() || !apF[1]->qValid() || !apF[2]->qValid() ||
	  !apF[3]->qValid()) {
	Vert *apV[4];
	apV[0] = dynamic_cast<Vert*>(static_cast<Entity*>
				     (lower_order_entity_handles[0]));
	apV[1] = dynamic_cast<Vert*>(static_cast<Entity*>
				     (lower_order_entity_handles[1]));
	apV[2] = dynamic_cast<Vert*>(static_cast<Entity*>
				     (lower_order_entity_handles[2]));
	apV[3] = dynamic_cast<Vert*>(static_cast<Entity*>
				     (lower_order_entity_handles[3]));

	if (!apV[0]->qValid() || !apV[1]->qValid() || !apV[2]->qValid() ||
	    !apV[3]->qValid()) {
	  *err = iBase_INVALID_ENTITY_TYPE;
	  *new_entity_handle = NULL;
	  *status = iBase_CREATION_FAILED;
	}
	// Now set up an array to hold the verts for all four faces.
	iBase_EntityHandle aehVerts[] =
	  {apV[0], apV[1], apV[2],
	   apV[0], apV[3], apV[1],
	   apV[1], apV[3], apV[2],
	   apV[2], apV[3], apV[0]};
	int faces_allocated = 4, faces_size = 4;
	int tmp_status[4], tmp_status_allocated = 4, tmp_status_size;
	// These are required for the silly thing to work
	iBase_EntityHandle *pEH = reinterpret_cast<iBase_EntityHandle*>(apF);
	int *ppStat = tmp_status;
	createEntArr(iMesh_TRIANGLE, aehVerts, 12, 
		     &pEH, &faces_allocated, &faces_size,
		     &ppStat, &tmp_status_allocated, &tmp_status_size,
		     err);
      }

      // Check whether the cell already exists.
      Cell *pCCom01 = pCCommonCell(apF[0], apF[1]);
      Cell *pCCom12 = pCCommonCell(apF[1], apF[2]);
      Cell *pCCom23 = pCCommonCell(apF[2], apF[3]);
      Cell *pCCom30 = pCCommonCell(apF[3], apF[0]);

      if (pCCom01->qValid() && pCCom12->qValid() && pCCom23->qValid() &&
	  pCCom30->qValid() && pCCom01 == pCCom12 && pCCom12 == pCCom23 &&
	  pCCom23 == pCCom30) {
	*new_entity_handle = pCCom01;
	*status = iBase_ALREADY_EXISTED;
      }
      else {
	TetCell *pTCNew = dynamic_cast<TetCell*>(pVM->pCNewCell(4, 4));
	if (!pTCNew->qValid()) {
	  *err = iBase_ENTITY_CREATION_ERROR;
	  *new_entity_handle = NULL;
	  *status = iBase_CREATION_FAILED;
	  break;
	}
	pTCNew->vAssign(apF);
	for (int i = 0; i < 4; i++) {
	  Vert *pVA = apF[i]->pVVert(0);
	  Vert *pVB = apF[i]->pVVert(1);
	  Vert *pVC = apF[i]->pVVert(2);
	  Vert *pVD;
	  {
	    int ii = (i+1)%4;
	    pVD = apF[ii]->pVVert(0);
	    if ((pVD == pVA) || (pVD == pVB) || (pVD == pVC)) {
	      pVD = apF[ii]->pVVert(1);
	      if ((pVD == pVA) || (pVD == pVB) || (pVD == pVC)) {
		pVD = apF[ii]->pVVert(2);
		assert((pVD != pVA) && (pVD != pVB) && (pVD != pVC));
	      }
	    }
	  }

	  if (iOrient3D(pVA, pVB, pVC, pVD) == 1) {
	    assert(!apF[i]->pCCellRight()->qValid());
	    apF[i]->vSetRightCell(pTCNew);
	  }
	  else {
	    assert(!apF[i]->pCCellLeft()->qValid());
	    apF[i]->vSetLeftCell(pTCNew);
	  }
	} // Done connecting faces to cells
	assert(pTCNew->iFullCheck());
	*new_entity_handle = pTCNew;
	*status = iBase_NEW;
	iNTet++;
      } // Done actually creating the new cell.
    } // Done creating tets
    break;
  case iMesh_PRISM:
    // This stuff will work once the GRUMMP mesh mod overhaul is merged
    // into the trunk.
//     {
//       if (lower_order_entity_handles_size != 6) {
// 	*err = iBase_INVALID_ENTITY_COUNT;
// 	break;
//       }
//       Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[0]));
//       Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[1]));
//       Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[2]));
//       Vert *pV3 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[3]));
//       Vert *pV4 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[4]));
//       Vert *pV5 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[5]));
//       assert(pV0 && pV1 && pV2 && pV3 && pV4 && pV5);
//       *new_entity_handle = pVM->createPrismCell(pV0, pV1, pV2, pV3, pV4, pV5);
//       *status = iBase_NEW;
//       iNPrism++;
//     }
//     break;    
  case iMesh_PYRAMID:
    // This stuff will work once the GRUMMP mesh mod overhaul is merged
    // into the trunk.
//     {
//       if (lower_order_entity_handles_size != 5) {
// 	*err = iBase_INVALID_ENTITY_COUNT;
// 	break;
//       }
//       Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[0]));
//       Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[1]));
//       Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[2]));
//       Vert *pV3 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[3]));
//       Vert *pV4 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[4]));
//       assert(pV0 && pV1 && pV2 && pV3 && pV4);
//       *new_entity_handle = pVM->createPyrCell(pV0, pV1, pV2, pV3, pV4);
//       *status = iBase_NEW;
//       iNPyr++;
//     }
//     break;    
  case iMesh_HEXAHEDRON:
    // This stuff will work once the GRUMMP mesh mod overhaul is merged
    // into the trunk.
//     {
//       if (lower_order_entity_handles_size != 8) {
// 	*err = iBase_INVALID_ENTITY_COUNT;
// 	break;
//       }
//       Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[0]));
//       Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[1]));
//       Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[2]));
//       Vert *pV3 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[3]));
//       Vert *pV4 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[4]));
//       Vert *pV5 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[5]));
//       Vert *pV6 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[6]));
//       Vert *pV7 = dynamic_cast<Vert*>(static_cast<Entity*>
// 				      (lower_order_entity_handles[7]));
//       assert(pV0 && pV1 && pV2 && pV3 && pV4 && pV5 && pV6 && pV7);
//       *new_entity_handle = pVM->createHexCell(pV0, pV1, pV2, pV3, pV4,
// 					      pV5, pV6, pV7);
//       *status = iBase_NEW;
//       iNHex++;
//     }
//     break;    
  case iMesh_SEPTAHEDRON:
  case iMesh_POLYHEDRON:
  case iMesh_LINE_SEGMENT:
  case iMesh_POLYGON:
    *err = iBase_NOT_SUPPORTED;
    break;
  case iMesh_ALL_TOPOLOGIES:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    break;
  }
}

void iMesh_Volume::deleteEnt(iBase_EntityHandle entity_handle,
				    int *err)
{
  *err = iBase_SUCCESS;

  Entity *pE = static_cast<Entity*>(entity_handle);
  switch (pE->eEntTopology()) {
  case iMesh_TETRAHEDRON:
  case iMesh_PYRAMID:
  case iMesh_PRISM:
  case iMesh_HEXAHEDRON:
    {
      // In this case, there are no upward adjacencies, so it's always
      // safe to get rid of these.
      vRemoveAllReferencesTo(pE);
      // Remove connections from its edges.
      Cell *pC = dynamic_cast<Cell*>(pE);
      assert(pC);
      for (int ii = pC->iNumFaces() - 1; ii >= 0; ii--) {
	pC->pFFace(ii)->vRemoveCell(pC);
      }
      pC->vMarkDeleted();
      switch (pC->iNumVerts()) {
      case 4:
	iNTet--;
	break;
      case 5:
	iNPyr--;
	break;
      case 6:
	iNPrism--;
	break;
      case 8:
	iNHex--;
	break;
      default:
	// How did someone create an entity like this?
	assert(0);
	break;
      }
      break;
    }
  case iMesh_QUADRILATERAL:
  case iMesh_TRIANGLE:
    // These could have cells that depend on them.
    {
      Face *pF = dynamic_cast<Face*>(pE);
      assert(pF);
      bool qLeftOkay = !pF->pCCellLeft()->qValid() ||
	pF->pCCellLeft()->qDeleted();
      bool qRightOkay = !pF->pCCellRight()->qValid() ||
	pF->pCCellRight()->qDeleted();
      if (pF->qIsBdryFace()) {
	BFace *pBF = pBFInvalidBFace;
	// Identify the bdry face, if it's still attached.
	if (qLeftOkay) {
	  // Only chance is on the right; if it's a BFace, the dynamic
	  // cast will put it in pBF, otherwise NULL.
	  pBF = dynamic_cast<BFace*>(pF->pCCellRight());
	}
	else if (qRightOkay) {
	  // Same thing, other side.
	  pBF = dynamic_cast<BFace*>(pF->pCCellLeft());
	}
	if (pBF) {
	  pF->vRemoveCell(pBF);
	  pBF->vMarkDeleted();
	}
	qLeftOkay = qRightOkay = true;
      }

      if (qLeftOkay && qRightOkay) {
	vRemoveAllReferencesTo(pE);
	pF->vMarkDeleted();
	vRemoveFaceFromHash(pF);
	if (pF->iNumVerts() == 3) iNTri--;
	else                      iNQuad--;
      }
      break;
    }
  case iMesh_POINT:
    // Every vertex that's connected to anything must have a valid
    // hint face.
    {
      Vert *pV = dynamic_cast<Vert*>(pE);
      assert(pV);
      if (!pV->qDeleted()) {
	Face *pF = pV->pFHintFace();
	if (!pF->qValid() || !pF->qHasVert(pV) || pF->qDeleted()) {
	  vRemoveAllReferencesTo(pE);
	  pV->vMarkDeleted();
	  iNVert--;
	}
      }
      break;
    }
  case iMesh_ALL_TOPOLOGIES:
  case iMesh_SEPTAHEDRON:
  case iMesh_POLYHEDRON:
  case iMesh_POLYGON:
  case iMesh_LINE_SEGMENT:
    // These are either illegal or never expected to be supported by GRUMMP.
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    break;
  }
}

struct qNotHasVert : public std::binary_function<Face*, Vert*, bool> {
  bool operator()(const Face* pF, const Vert* pV) const
  { return !pF->qHasVert(pV); }
};

struct qHasVert : public std::binary_function<Face*, Vert*, bool> {
  bool operator()(const Face* pF, const Vert* pV) const
  { return pF->qHasVert(pV); }
};

void iMesh_Volume::getEntAdj(const iBase_EntityHandle entity_handle,
				    const int entity_type_requested,
				    iBase_EntityHandle** adj_entity_handles,
				    int* adj_entity_handles_allocated,
				    int* adj_entity_handles_size, int *err)
{
  Entity *pE = static_cast<Entity*>(entity_handle);
  int eET = pE->eEntType();
  if (eET == entity_type_requested ||
      iBase_EDGE == entity_type_requested) {
    // Not supposed to return any data in this case.
    *adj_entity_handles_size = 0;
    TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated, 0, iBase_EntityHandle);
  }
  else if (eET == iBase_VERTEX &&
	   entity_type_requested != iBase_VERTEX) {
    // Now get fancy with a call that computes the neighborhood.  This
    // is the only expensive adjacency request, requiring a local mesh
    // traversal.
    Vert *pV = dynamic_cast<Vert*>(pE);
    assert(pV->qValid());
    std::set<Cell*> spCInc;
    std::set<Vert*> spVNeigh;
    std::set<BFace*> spBFInc;
    bool qBdryVert;
    std::set<Face*> spFNearby;
    vNeighborhood(pV, spCInc, spVNeigh, &spBFInc, &qBdryVert, &spFNearby);
    if (entity_type_requested == iBase_FACE) {
      // Add faces incident on pV.
      int iCount = std::count_if(spFNearby.begin(), spFNearby.end(),
				 bind2nd(qHasVert(), pV));
      *adj_entity_handles_size = iCount;
      TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		     *adj_entity_handles_size, iBase_EntityHandle);
      std::remove_copy_if(spFNearby.begin(), spFNearby.end(),
			  *adj_entity_handles, bind2nd(qNotHasVert(), pV));
    }
    else if (entity_type_requested == iBase_REGION) {
      *adj_entity_handles_size = spCInc.size();
      TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		     *adj_entity_handles_size, iBase_EntityHandle);
      std::copy(spCInc.begin(), spCInc.end(), *adj_entity_handles);
    }
    else {
      assert(entity_type_requested == iBase_ALL_TYPES);
      // Add faces incident on pV.
      int iCount = std::count_if(spFNearby.begin(), spFNearby.end(),
				 bind2nd(qHasVert(), pV));
      *adj_entity_handles_size = iCount + spCInc.size();
      TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		     *adj_entity_handles_size, iBase_EntityHandle);
      iBase_EntityHandle* tmp =
	std::remove_copy_if(spFNearby.begin(), spFNearby.end(),
			  *adj_entity_handles, bind2nd(qNotHasVert(), pV));
      iBase_EntityHandle* tmp2 =
	std::copy(spCInc.begin(), spCInc.end(), tmp);
      assert(tmp2 - *adj_entity_handles == *adj_entity_handles_size);
    }
  }
  else {
    switch (entity_type_requested) {
    case iBase_VERTEX:
      {
	*adj_entity_handles_size = pE->iNumVerts();
	TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		       *adj_entity_handles_size, iBase_EntityHandle);
	pE->vAllVertHandles(*adj_entity_handles);
	break;
      }
    case iBase_FACE:
      {
	*adj_entity_handles_size = pE->iNumFaces();
	TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		       *adj_entity_handles_size, iBase_EntityHandle);
	pE->vAllFaceHandles(*adj_entity_handles);
	break;
      }
    case iBase_REGION:
      {
	*adj_entity_handles_size = pE->iNumCells();
	TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		       *adj_entity_handles_size, iBase_EntityHandle);
	pE->vAllCellHandles(*adj_entity_handles);
	if (*adj_entity_handles[0] == NULL) {
	  *adj_entity_handles[0] = *adj_entity_handles[1];
	  *adj_entity_handles_size = *adj_entity_handles_size - 1;
	}
	break;
      }
    case iBase_ALL_TYPES:
      {
	*adj_entity_handles_size = pE->iNumCells() + pE->iNumVerts() +
	  pE->iNumFaces();
	TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		       *adj_entity_handles_size, iBase_EntityHandle);
	iBase_EntityHandle *tmp = *adj_entity_handles;
	pE->vAllVertHandles(tmp);
	tmp += pE->iNumVerts();
	pE->vAllFaceHandles(tmp);
	tmp += pE->iNumFaces();
	pE->vAllCellHandles(tmp);
	tmp += pE->iNumCells();
	assert(tmp - *adj_entity_handles == *adj_entity_handles_size);
      }
    case iBase_EDGE:
      {
	// Do nothing, because there are no edges in the 3D mesh.
	// In fact, applications with any sense have multiple ways to
	// decide not to get here, but don't throw an error anyway,
	// because the spec implies we shouldn't.
	*adj_entity_handles_size = 0;
	TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		       *adj_entity_handles_size, iBase_EntityHandle);
	break;
      }
    default:
      // Can't get here, because you'll have already thrown an exception.
      assert(0);
    }
  } // The general case (not vertex nor requesting type of current entity
}

void iMesh_Volume::getEntArrAdj(const iBase_EntityHandle* entity_handles,
				       const int entity_handles_size,
				       const int entity_type_requested,
				       iBase_EntityHandle** adj_entity_handles,
				       int* adj_entity_handles_allocated,
				       int* adj_entity_handles_size,
				       int** offset,
				       int* offset_allocated,
				       int* offset_size, int *err)
{
  // Even though it may be a bit awkward, put the switch on requested type
  // on the -outside-, with separate loops inside each case statement.  This
  // will make the code longer, but presumably faster...
  *err = iBase_SUCCESS;

  *offset_size = entity_handles_size+1;
  TRY_ARRAY_SIZE(offset, offset_allocated, *offset_size, int);

  switch (entity_type_requested) {
  case iBase_VERTEX: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE, iEOut = 0;
    for (iE = 0; iE < entity_handles_size; iE++) {
      (*offset)[iE] = iEOut;
      Entity *pE = static_cast<Entity*>(entity_handles[iE]);
      int ET = pE->eEntType();
      if (ET != entity_type_requested) {
	iEOut += pE->iNumVerts();
      }
    }

    (*offset)[iE] = iEOut;

    // Now set up the output array and transcribe the data.
    TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated, iEOut, iBase_EntityHandle);
    *adj_entity_handles_size = iEOut;

    for (iE = 0; iE < entity_handles_size; iE++) {
      Entity *pE = static_cast<Entity*>(entity_handles[iE]);
      int ET = pE->eEntType();
      if (ET != entity_type_requested) {
	pE->vAllVertHandles(*adj_entity_handles + (*offset)[iE]);
      }
    } // The general case (not vertex nor requesting type of current entity
    break;
  }
  case iBase_EDGE: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE, iEOut = 0;
    for (iE = 0; iE < entity_handles_size; iE++) {
      (*offset)[iE] = iEOut;
      // Since GRUMMP currently stores no edges, there's nothing to do here.
    }
    (*offset)[iE] = iEOut;

    // Now set up the output array and transcribe the data.
    TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated, iEOut, iBase_EntityHandle);
    *adj_entity_handles_size = iEOut;

    for (iE = 0; iE < entity_handles_size; iE++) {
      // Do nothing, because there are no edges represented in the 3D mesh.
      // In fact, applications with any sense have multiple ways to
      // decide not to get here, but don't throw an error anyway,
      // because the spec implies we shouldn't.
    }
    break;
  }
  case iBase_FACE: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE, iEOut = 0;
    for (iE = 0; iE < entity_handles_size; iE++) {
      (*offset)[iE] = iEOut;
      Entity *pE = static_cast<Entity*>(entity_handles[iE]);
      int ET = pE->eEntType();
      if (ET == entity_type_requested) {
	// Not supposed to return any data in this case.
      }
      else if (ET == iBase_VERTEX &&
	       entity_type_requested != iBase_VERTEX) {
	// Now get fancy with a call that computes the neighborhood.  This
	// is the only expensive adjacency request, requiring a local mesh
	// traversal.
	Vert *pV = dynamic_cast<Vert*>(pE);
	assert(pV->qValid());
	std::set<Cell*> spCInc;
	std::set<Vert*> spVNeigh;
	std::set<BFace*> spBFInc;
	bool qBdryVert;
	std::set<Face*> spFNearby;
	vNeighborhood(pV, spCInc, spVNeigh, &spBFInc, &qBdryVert, &spFNearby);
	iEOut += std::count_if(spFNearby.begin(), spFNearby.end(),
			       bind2nd(qHasVert(), pV));
      }
      else {
	iEOut += pE->iNumFaces();
      }
    }

    (*offset)[iE] = iEOut;

    // Now set up the output array and transcribe the data.
    TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated, iEOut, iBase_EntityHandle);
    *adj_entity_handles_size = iEOut;

    iBase_EntityHandle* pEH = *adj_entity_handles;

    for (iE = 0; iE < entity_handles_size; iE++) {
      Entity *pE = static_cast<Entity*>(entity_handles[iE]);
      int ET = pE->eEntType();
      if (ET == entity_type_requested) {
	// Not supposed to return any data in this case.
      }
      else if (ET == iBase_VERTEX &&
	       entity_type_requested != iBase_VERTEX) {
	// Now get fancy with a call that computes the neighborhood.  This
	// is the only expensive adjacency request, requiring a local mesh
	// traversal.
	Vert *pV = dynamic_cast<Vert*>(pE);
	assert(pV->qValid());
	std::set<Cell*> spCInc;
	std::set<Vert*> spVNeigh;
	std::set<BFace*> spBFInc;
	bool qBdryVert;
	std::set<Face*> spFNearby;
	vNeighborhood(pV, spCInc, spVNeigh, &spBFInc, &qBdryVert, &spFNearby);

	// Add faces incident on pV.
	pEH = std::remove_copy_if(spFNearby.begin(), spFNearby.end(),
				  (*adj_entity_handles) + (*offset)[iE],
				  bind2nd(qNotHasVert(), pV));
	assert(pEH == (*adj_entity_handles) + (*offset)[iE+1]);
      }
      else {
	pE->vAllFaceHandles(*adj_entity_handles + (*offset)[iE]);
      }
    }
    break;
  }

  case iBase_REGION: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE, iEOut = 0;
    for (iE = 0; iE < entity_handles_size; iE++) {
      (*offset)[iE] = iEOut;
      Entity *pE = static_cast<Entity*>(entity_handles[iE]);
      int ET = pE->eEntType();
      if (ET == entity_type_requested) {
	// Not supposed to return any data in this case.
      }
      else if (ET == iBase_VERTEX) {
	// Now get fancy with a call that computes the neighborhood.  This
	// is the only expensive adjacency request, requiring a local mesh
	// traversal.
	Vert *pV = dynamic_cast<Vert*>(pE);
	assert(pV->qValid());
	std::set<Cell*> spCInc;
	std::set<Vert*> spVNeigh;
	std::set<BFace*> spBFInc;
	bool qBdryVert;
	std::set<Face*> spFNearby;
	vNeighborhood(pV, spCInc, spVNeigh, &spBFInc, &qBdryVert, &spFNearby);
	iEOut += spCInc.size();
      }
      else {
	iEOut += pE->iNumCells();
      }
    }
    (*offset)[iE] = iEOut;

    // Now set up the output array and transcribe the data.
    TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated, iEOut, iBase_EntityHandle);
    *adj_entity_handles_size = iEOut;

    iBase_EntityHandle* pEH = *adj_entity_handles;

    for (iE = 0; iE < entity_handles_size; iE++) {
      Entity *pE = static_cast<Entity*>(entity_handles[iE]);
      int ET = pE->eEntType();
      if (ET == entity_type_requested) {
	// Not supposed to return any data in this case.
      }
      else if (ET == iBase_VERTEX) {
	// Now get fancy with a call that computes the neighborhood.  This
	// is the only expensive adjacency request, requiring a local mesh
	// traversal.
	Vert *pV = dynamic_cast<Vert*>(pE);
	assert(pV->qValid());
	std::set<Cell*> spCInc;
	std::set<Vert*> spVNeigh;
	std::set<BFace*> spBFInc;
	bool qBdryVert;
	std::set<Face*> spFNearby;
	vNeighborhood(pV, spCInc, spVNeigh, &spBFInc, &qBdryVert, &spFNearby);
	pEH = std::copy(spCInc.begin(), spCInc.end(), pEH);
	assert(pEH == (*adj_entity_handles) + (*offset)[iE+1]);
      }
      else {
	pE->vAllCellHandles(*adj_entity_handles + (*offset)[iE]);
      }
    }
    break;
  }
  case iBase_ALL_TYPES:
    {
      // Pre-compute the output size, and put the correct values in the
      // offset array.
      int iE, iEOut = 0;
      for (iE = 0; iE < entity_handles_size; iE++) {
	(*offset)[iE] = iEOut;
	Entity *pE = static_cast<Entity*>(entity_handles[iE]);
	int ET = pE->eEntType();
	if (ET == iBase_VERTEX) {
	  // Now get fancy with a call that computes the neighborhood.  This
	  // is the only expensive adjacency request, requiring a local mesh
	  // traversal.
	  Vert *pV = dynamic_cast<Vert*>(pE);
	  assert(pV->qValid());
	  std::set<Cell*> spCInc;
	  std::set<Vert*> spVNeigh;
	  std::set<Face*> spFNearby;
	  vNeighborhood(pV, spCInc, spVNeigh, NULL, NULL, &spFNearby);
	  iEOut += spCInc.size();
	  iEOut += std::count_if(spFNearby.begin(), spFNearby.end(),
				 bind2nd(qHasVert(), pV));
	}
	else {
	  iEOut += pE->iNumCells() + pE->iNumFaces() + pE->iNumVerts();
	}
      }
      (*offset)[iE] = iEOut;
      
      // Now set up the output array and transcribe the data.
      TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated, iEOut, iBase_EntityHandle);
      *adj_entity_handles_size = iEOut;
      
      iBase_EntityHandle* pEH = *adj_entity_handles;
      
      for (iE = 0; iE < entity_handles_size; iE++) {
	Entity *pE = static_cast<Entity*>(entity_handles[iE]);
	int ET = pE->eEntType();
	if (ET == entity_type_requested) {
	  // Not supposed to return any data in this case.
	}
	else if (ET == iBase_VERTEX) {
	  // Now get fancy with a call that computes the neighborhood.  This
	  // is the only expensive adjacency request, requiring a local mesh
	  // traversal.
	  Vert *pV = dynamic_cast<Vert*>(pE);
	  assert(pV->qValid());
	  std::set<Cell*> spCInc;
	  std::set<Vert*> spVNeigh;
	  std::set<Face*> spFNearby;
	  vNeighborhood(pV, spCInc, spVNeigh, NULL, NULL, &spFNearby);
	  pEH = std::copy(spCInc.begin(), spCInc.end(), pEH);
	  pEH = std::remove_copy_if(spFNearby.begin(), spFNearby.end(),
				    pEH, bind2nd(qNotHasVert(), pV));
	  
	  assert(pEH == (*adj_entity_handles) + (*offset)[iE+1]);
	}
	else {
	  pE->vAllCellHandles(pEH); pEH += pE->iNumCells();
	  pE->vAllFaceHandles(pEH); pEH += pE->iNumFaces();
	  pE->vAllVertHandles(pEH); pEH += pE->iNumVerts();
	}
      }
      assert(pEH == *adj_entity_handles + *adj_entity_handles_size);
      break;
    }
  } // End of the switch
} // End of iMesh_getEntArrAdj

void 
iMesh_Volume::getEnt2ndAdj(const iBase_EntityHandle entity_handle,
				  const int bridge_entity_type,
				  const int requested_entity_type,
				  iBase_EntityHandle** adjacent_entities,
				  int* adjacent_entities_allocated,
				  int* adjacent_entities_size,
				  int* err)
{
  *err = iBase_SUCCESS;
  Entity *pE = static_cast<Entity*>(entity_handle);
  int input_type = pE->eEntType();

  // If the requested type is EDGE, there's no data to be returned.
  if (requested_entity_type == iBase_EDGE) {
    *adjacent_entities_size = 0;
    return;
  }

  // Look, the key here is pretty simple:  identify short-circuit cases
  // that can be handled by a single ITAPS adjacency request; this
  // includes things starting from a vertex, for which GRUMMP has a
  // built-in call that makes a second adjacency call always
  // unnecessary.  The short-circuit cases (in addition to vertex) are
  // those in which the bridge type matches one of the other two and
  // isn't iBase_ALL_TYPES; these are handled by a first-adjacency call
  // to requested_type.  Otherwise, call getEntAdj, then getEntArrAdj
  // and call it a day.

  if (input_type == iBase_VERTEX) {
    // Huge possibilities for optimization here compared with two
    // adjacency calls, because the second one is never needed.  The
    // vertex adjacency call would require a call to vNeighborhood
    // anyway, and that contains all the info we need.

    // Possible cases:
    // FACE and ALL_TYPES as bridge type turn out to be identical.

    //  ->all->all   = everything from neighborhood calc
    //  ->all->regi  = 1st adj regions
    //  ->all->face  = all faces from neighborhood calc
    //  ->all->edge  = nothing
    //  ->all->vert  = all verts from neighborhood calc

    //  ->regi->all  = everything from neighborhood calc except regions
    //  ->regi->regi = nothing
    //  ->regi->face = all faces from neighborhood calc
    //  ->regi->edge = nothing
    //  ->regi->vert = all verts from neighborhood calc

    //  ->face->all  = all of neighborhood calc except opposite faces 
    //  ->face->regi = 1st adj regions
    //  ->face->face = nothing
    //  ->face->edge = nothing
    //  ->face->vert = most verts from neighborhood calc

    //  ->edge->all  = all of neighborhood calc except opposite faces
    //                    and verts not connected by an edge
    //  ->edge->regi = 1st adj regions
    //  ->edge->face = 1st adj faces
    //  ->edge->edge = nothing
    //  ->edge->vert = most verts from neighborhood calc

    //  ->vert->all  = nothing
    //  ->vert->regi = nothing
    //  ->vert->face = nothing
    //  ->vert->edge = nothing
    //  ->vert->vert = nothing

    if (bridge_entity_type == iBase_VERTEX ||
	(bridge_entity_type == requested_entity_type &&
	 bridge_entity_type != iBase_ALL_TYPES) ||
	requested_entity_type == iBase_EDGE) {
      // This handles eleven of twenty-five cases by returning nothing!
      *adjacent_entities_size = 0;
      TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		     *adjacent_entities_size, iBase_EntityHandle);
    }
    else if ((requested_entity_type == iBase_REGION) ||
	(requested_entity_type == iBase_FACE &&
	 bridge_entity_type == iBase_EDGE)) {
      // This handles another four cases (->e->r, ->f->r, ->a->r, ->e->f).
      getEntAdj(entity_handle, requested_entity_type, adjacent_entities,
		adjacent_entities_allocated, adjacent_entities_size, err);
    }
    else {
      // In all of these cases, we need to compute the neighborhood of
      // the vertex.
      Vert *pV = dynamic_cast<Vert*>(entity_handle);
      std::set<Cell*> spCInc;
      std::set<Vert*> spVNeigh;
      std::set<Face*> spFNearby;
      vNeighborhood(pV, spCInc, spVNeigh, NULL, NULL, &spFNearby);
      int iNCells = spCInc.size();
      int iNFaces = spFNearby.size();
      int iNVerts = spVNeigh.size();

      if (bridge_entity_type >= iBase_REGION) {
	// A total of six cases: bridge = region or all, with requested
	// = vertex, face, or all
	if (requested_entity_type == iBase_ALL_TYPES) {
	  // Two more cases (regi->all and ->all->all).  These are easy:
	  // they take everything

	  if (bridge_entity_type == iBase_REGION) {
	    *adjacent_entities_size = iNFaces + iNVerts;
	  }
	  else {
	    *adjacent_entities_size = iNCells + iNFaces + iNVerts;
	  }
	  
	  TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
			 *adjacent_entities_size, iBase_EntityHandle);
	  int ii = 0;
	  if (bridge_entity_type != iBase_REGION) {
	    std::copy(spCInc.begin(), spCInc.end(), (*adjacent_entities)+ii);
	    ii += iNCells;
	  }
	  
	  std::copy(spVNeigh.begin(), spVNeigh.end(), (*adjacent_entities)+ii);
	  ii += iNVerts;
	  
	  std::copy(spFNearby.begin(), spFNearby.end(),
		    (*adjacent_entities)+ii);
	  ii += iNFaces;
	  
	  assert(ii == *adjacent_entities_size);
	} // Done with requested_type = ALL_TYPES for region and all
	  // bridge
	else if (requested_entity_type == iBase_FACE) {
	  // Two cases (->all->face, ->regi->face); either way, we need
	  // all the faces.
	  *adjacent_entities_size = iNFaces;
	  TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
			 *adjacent_entities_size, iBase_EntityHandle);
	  iBase_EntityHandle* ehResult = std::copy(spFNearby.begin(),
						   spFNearby.end(),
						   *adjacent_entities);
	  assert(ehResult - *adjacent_entities == *adjacent_entities_size);
	}
	else {
	  assert(requested_entity_type == iBase_VERTEX);
	  // Two cases (->all->vert, ->face->vert); either way, we need
	  // all the verts.
	  *adjacent_entities_size = iNVerts;
	  TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
			 *adjacent_entities_size, iBase_EntityHandle);
	  iBase_EntityHandle* ehResult = std::copy(spVNeigh.begin(),
						   spVNeigh.end(),
						   *adjacent_entities);
	  assert(ehResult - *adjacent_entities == *adjacent_entities_size);
	}
      } // Done with bridging via region or all
      else {
	// The other combinations (->face->all, ->face->vert, edge->all,
	// edge->vert) are tricky, or at least tedious.  For now (read:
	// until it turns out to be an issue), these will get done as
	// two single adjacency calls.  For edges, this implies an empty
	// result, but that's fine.
	if (bridge_entity_type == iBase_EDGE) {
	  *adjacent_entities_size = 0;
	}
	else {
	  assert(bridge_entity_type == iBase_FACE &&
		 (requested_entity_type == iBase_ALL_TYPES ||
		  requested_entity_type == iBase_VERTEX));
	  iBase_EntityHandle *intermediate_ents = NULL;
	  int intermediate_ents_allocated = 0, intermediate_ents_size = 0;
	  getEntAdj(entity_handle, bridge_entity_type, &intermediate_ents,
		    &intermediate_ents_allocated, &intermediate_ents_size, err);
	  
	  iBase_EntityHandle *temp_ents = NULL;
	  int temp_ents_allocated = 0, temp_ents_size = 0;
	  int *offset, offset_allocated = 0, offset_size;
	  getEntArrAdj(intermediate_ents, intermediate_ents_size,
		       requested_entity_type, &temp_ents,
		       &temp_ents_allocated, &temp_ents_size,
		       &offset, &offset_allocated, &offset_size, err);
	  std::set<iBase_EntityHandle> setOutput(temp_ents,
						 temp_ents + temp_ents_size);
	  setOutput.erase(entity_handle);
	  *adjacent_entities_size = setOutput.size();
	  TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
			 *adjacent_entities_size, iBase_EntityHandle);
	  std::copy(setOutput.begin(), setOutput.end(),*adjacent_entities);
	  if (intermediate_ents) free(intermediate_ents);
	  if (temp_ents) free(temp_ents);
	  if (offset) free(offset);
	}
      } // Done with the face/edge->all/vert cases
    } // Done with vert->A->B cases that aren't just first adjacencies
  } // Done with vert->A->B
  else if ((bridge_entity_type == input_type ||
	    bridge_entity_type == requested_entity_type) &&
	   bridge_entity_type != iBase_ALL_TYPES) {
    *adjacent_entities_size = 0;
    TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		   *adjacent_entities_size, iBase_EntityHandle);
  }
  else {
    // No better alternative algorithmically than to make two
    // consecutive adjacency calls (one entity, one array).
    iBase_EntityHandle *intermediate_ents = NULL;
    int intermediate_ents_allocated = 0, intermediate_ents_size;
    getEntAdj(entity_handle, bridge_entity_type, &intermediate_ents,
	      &intermediate_ents_allocated, &intermediate_ents_size, err);

    iBase_EntityHandle *temp_ents = NULL;
    int temp_ents_allocated = 0, temp_ents_size = 0;
    int *offset = NULL, offset_allocated = 0, offset_size;
    getEntArrAdj(intermediate_ents, intermediate_ents_size,
		 requested_entity_type, &temp_ents,
		 &temp_ents_allocated, &temp_ents_size,
		 &offset, &offset_allocated, &offset_size, err);

    std::set<iBase_EntityHandle> setOutput(temp_ents,
					   temp_ents + temp_ents_size);
    setOutput.erase(entity_handle);
    *adjacent_entities_size = setOutput.size();
    TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		   *adjacent_entities_size, iBase_EntityHandle);
    std::copy(setOutput.begin(), setOutput.end(),*adjacent_entities);

    if (intermediate_ents) free(intermediate_ents);
    if (temp_ents) free(temp_ents);
    if (offset) free(offset);
  }  
}

