#include "iMesh.h"
#include "iMesh_GRUMMP.hh"

#include <vector>
#include "GR_Entity.h"
#include "GR_Geometry.h"
#include "GR_Mesh2D.h"
#include "iMesh_GRUMMP.hh"
#include "iMesh_GRUMMP_misc.hh"
#include "iMesh_GRUMMP_EntIter.hh"
#include "iMesh_GRUMMP_Workset.hh"

#undef SLOW_BUT_SURE

using namespace ITAPS_GRUMMP;

void iMesh_Plane::vCreateMesh()
{
  Mesh2D *pM2D = new Mesh2D();
  pM = pM2D;
}

void iMesh_Plane::load(const char *name, const char* /*options*/,
		       int* err)
{
  *err = iBase_SUCCESS;
  Mesh2D *pM2D = new Mesh2D(name);
  pM = pM2D;
  // Set up the initial (GRUMMP) face hash table.
  for (int iF = pM->iNumFaces() - 1; iF >= 0; iF--) {
    vAddFaceToHash(pM->pFFace(iF));
  }
  iNVert = pM->iNumVerts();
  iNLineSeg = pM2D->iNumFaces();
  iNTri = pM2D->iNumTriCells();
  iNQuad = pM2D->iNumQuadCells();
  iNTet = iNPyr = iNPrism = iNHex = iNPolygon = iNSept = iNPolyhedron = 0;
  assert((iNLineSeg) == mmFaceHash.size());
}

void iMesh_Plane::save (const char *name, const char* /*options*/,
			int *err)
{
  *err = iBase_SUCCESS;
  Mesh2D *pM2D = dynamic_cast<Mesh2D*>(pM);
  vWriteFile_Mesh2D(*pM2D, name);
}

void iMesh_Plane::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 L L U
  //  E I I I U
  //  F I I I U
  //  R U U U U
  // 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_ALL_ORDER_1;
  adjacency_table[2] = iBase_ALL_ORDER_1;
  adjacency_table[3] = iBase_UNAVAILABLE;
  
  adjacency_table[4] = iBase_ALL_ORDER_1;
  adjacency_table[5] = iBase_ALL_ORDER_1;
  adjacency_table[6] = iBase_ALL_ORDER_1;
  adjacency_table[7] = iBase_UNAVAILABLE;
  
  adjacency_table[8] = iBase_ALL_ORDER_1;
  adjacency_table[9] = iBase_ALL_ORDER_1;
  adjacency_table[10]= iBase_ALL_ORDER_1;
  adjacency_table[11]= iBase_UNAVAILABLE;
  
  adjacency_table[12]= iBase_UNAVAILABLE;
  adjacency_table[13]= iBase_UNAVAILABLE;
  adjacency_table[14]= iBase_UNAVAILABLE;
  adjacency_table[15]= iBase_UNAVAILABLE;
}

void iMesh_Plane::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;
  Mesh2D *pM2D = dynamic_cast<Mesh2D*>(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);
    *entity_handles_size = iNVert;
  }
  else if (entity_topology == iMesh_LINE_SEGMENT ||
	   (entity_topology == iMesh_ALL_TOPOLOGIES &&
	    entity_type     == iBase_EDGE) ) {
    ////////////////////////////////////////////////
    // Returning all edges
    ////////////////////////////////////////////////
    int iNFace_db = pM2D->iNumFaces();
    
    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNFace_db,
		   iBase_EntityHandle);
    
    int iF = 0, iEH = 0;
    // Skip to the first non-deleted face
    while (iF < iNFace_db && pM2D->pFFace(iF)->qDeleted()) iF++;
    
    while (iF < iNFace_db) {
      (*entity_handles)[iEH] = pM2D->pFFace(iF);
      iEH++;
      iF++;
      while (iF < iNFace_db && pM2D->pFFace(iF)->qDeleted()) iF++;
    }
    *entity_handles_size = iEH;
    assert(iEH  <= (*entity_handles_allocated));
  } // Done with edges.
  else if (entity_topology == iMesh_TRIANGLE) {
    ////////////////////////////////////////////////
    // Returning all triangular cells.
    ////////////////////////////////////////////////
    int iNTri_db = pM2D->iNumTriCells();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNTri_db,
		   iBase_EntityHandle);

    int iTri = 0, iEH = 0;
    // Skip to the first non-deleted tri cell
    while (iTri < iNTri_db && pM2D->pCCell(iTri)->qDeleted()) iTri++;

    while (iTri < iNTri_db) {
      (*entity_handles)[iEH] = pM2D->pCCell(iTri);
      iEH++;
      iTri++;
      while (iTri < iNTri_db && pM2D->pCCell(iTri)->qDeleted()) iTri++;
    }
    *entity_handles_size = iEH;
    assert(iEH  <= (*entity_handles_allocated));
  } // Done with triangular cells.
  else if (entity_topology == iMesh_QUADRILATERAL) {
    ////////////////////////////////////////////////
    // Returning all quadrilateral cells.
    ////////////////////////////////////////////////
    int iNTri_db = pM2D->iNumTriCells();
    int iNQuad_db = pM2D->iNumQuadCells();
    int iNCell_db = iNTri_db + iNQuad_db;

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNQuad_db,
		   iBase_EntityHandle);

    int iQuad = iNTri_db, iEH = 0;
    // Skip to the first non-deleted quad cell
    while (iQuad < iNCell_db && pM2D->pCCell(iQuad)->qDeleted()) iQuad++;

    while (iQuad < iNCell_db) {
      (*entity_handles)[iEH] = pM2D->pCCell(iQuad);
      iEH++;
      iQuad++;
      while (iQuad < iNCell_db && pM2D->pCCell(iQuad)->qDeleted()) iQuad++;
    }
    *entity_handles_size = iEH;
    assert(iEH  <= (*entity_handles_allocated));
  } // Done with quadrilateral cells.
  else if (entity_topology == iMesh_ALL_TOPOLOGIES &&
	   entity_type == iBase_FACE) {
    ////////////////////////////////////////////////
    // Returning all cells.
    ////////////////////////////////////////////////
    int iNCell_db = pM2D->iNumCells();

    TRY_ARRAY_SIZE(entity_handles, entity_handles_allocated, iNCell_db,
		   iBase_EntityHandle);

    int iCell = 0, iEH = 0;
    // Skip to the first non-deleted cell
    while (iCell < iNCell_db && pM2D->pCCell(iCell)->qDeleted()) iCell++;

    while (iCell < iNCell_db) {
      (*entity_handles)[iEH] = pM2D->pCCell(iCell);
      iEH++;
      iCell++;
      while (iCell < iNCell_db && pM2D->pCCell(iCell)->qDeleted()) iCell++;
    }
    *entity_handles_size = iEH;
    assert(iEH  <= (*entity_handles_allocated));
  } // 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 );
  }
  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 {
    ////////////////////////////////////////////////
    // 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_Plane::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)
{
  Mesh2D *pM2D = dynamic_cast<Mesh2D*>(pM);

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

#ifdef NDEBUG
      Vert *pV0 = static_cast<Vert*>(lower_order_entity_handles[0]);
      Vert *pV1 = static_cast<Vert*>(lower_order_entity_handles[1]);
#else
      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]));

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

      Face *pF = pFInvalidFace;

      Vert *apV[] = {pV0, pV1};
      if (qDoesFaceExist(2, apV, pF)) {
	// Return this face as a handle, with status ALREADY_EXISTED
	assert(pF->qValid());
	*new_entity_handle = pF;
	*status = iBase_ALREADY_EXISTED;
      }
      else {
	pF = pM2D->createFace(pV0, pV1);
 	if (!pF->qValid()) {
 	  *err = iBase_ENTITY_CREATION_ERROR;
 	  *status = iBase_CREATION_FAILED;
 	  *new_entity_handle = NULL;
 	}

	else {
	  *new_entity_handle = pF;
	  *status = iBase_NEW;
	  
	  vAddFaceToHash(pF);
	  iNLineSeg++;
	  assert((iNLineSeg) == mmFaceHash.size());
	}
      }
      break;
    }
  case iMesh_TRIANGLE:
    {
      // New triangles to be built from faces or vertices.
      if (lower_order_entity_handles_size != 3) {
	*err = iBase_INVALID_ENTITY_COUNT;
	*status = iBase_CREATION_FAILED;
	*new_entity_handle = NULL;
	break;
      }
      
      Entity *pE = static_cast<Entity*>(static_cast<Entity*>(lower_order_entity_handles[0]));
      int ET = pE->eEntType();
      Face *apF[3];
      
      if (ET == iBase_VERTEX) {
	// Build edges first, then use the same code as the edge case
	// from there.
#ifdef NDEBUG
	Vert *pV0 = static_cast<Vert*>(lower_order_entity_handles[0]);
	Vert *pV1 = static_cast<Vert*>(lower_order_entity_handles[1]);
	Vert *pV2 = static_cast<Vert*>(lower_order_entity_handles[2]);
#else
	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;
	  *status = iBase_CREATION_FAILED;
	  *new_entity_handle = NULL;
	  break;
	}
#endif
	
	// This is a sneaky cheat on identifying verts for edges that
	// works because the verts are cyclically arranged around the
	// triangle.
	Vert *apV[] = {pV0, pV1, pV2, pV0};
	for (int iF = 0; iF < 3; iF++) {
	  if (!qDoesFaceExist(2, apV + iF, apF[iF])) {
	    apF[iF] = pM2D->pFNewFace();
	    if (!apF[iF]->qValid()) {
	      *err = iBase_ENTITY_CREATION_ERROR;
	      *status = iBase_CREATION_FAILED;
	      *new_entity_handle = NULL;
	      break;
	    }
	    
	    apF[iF]->vAssign(pCInvalidCell, pCInvalidCell, apV[iF], apV[iF+1]);
	    apV[iF]->vSetHintFace(apF[iF]);
	    apV[iF+1]->vSetHintFace(apF[iF]);
	    
	    vAddFaceToHash(apF[iF]);
	    iNLineSeg++;
	    assert((iNLineSeg) == mmFaceHash.size());
	  } // Done creating this new edge
	} 
      } // Done creating/identifying edges from verts.
      else {
#ifdef NDEBUG
	apF[0] = static_cast<Face*>(lower_order_entity_handles[0]);
	apF[1] = static_cast<Face*>(lower_order_entity_handles[1]);
	apF[2] = static_cast<Face*>(lower_order_entity_handles[2]);
#else
	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]));
#endif
      } // Done identifying edges.
      
      // The following also handles the case of getting handed weird
      // stuff instead of legitimate edge faces.
      if (!apF[0]->qValid() || !apF[1]->qValid() || !apF[2]->qValid()) {
	*err = iBase_INVALID_ENTITY_HANDLE;
	*status = iBase_CREATION_FAILED;
	*new_entity_handle = NULL;
	break;
      }
      
      // Verify that the cell doesn't already exist.
      Cell *pCCom01 = pCCommonCell(apF[0], apF[1]);
      Cell *pCCom12 = pCCommonCell(apF[1], apF[2]);
      Cell *pCCom20 = pCCommonCell(apF[2], apF[0]);
      
      if (pCCom01->qValid() && pCCom12->qValid() && pCCom20->qValid() &&
	  pCCom01 == pCCom12 && pCCom12 == pCCom20) {
	*new_entity_handle = pCCom01;
	*status = iBase_ALREADY_EXISTED;
      }
      else {
	// pCNewCell(3) will -always- return a TriCell* cast to Cell*,
	// or die trying.
	TriCell *pTCNew = static_cast<TriCell*>(pM2D->pCNewCell(3));
	if (!pTCNew->qValid()) {
	  *err = iBase_ENTITY_CREATION_ERROR;
	  *status = iBase_CREATION_FAILED;
	  *new_entity_handle = NULL;
	  break;
	}
	
	pTCNew->vAssign(apF);
	Vert *apV[] = {apF[0]->pVVert(0), apF[0]->pVVert(1),
		       apF[1]->pVVert(0)};
	if (apV[0] == apV[2] || apV[1] == apV[2])
	  apV[2] = apF[1]->pVVert(1);
	
	for (int i = 0; i < 3; i++) {
	  Vert *pVA = apF[i]->pVVert(0);
	  Vert *pVB = apF[i]->pVVert(1);
	  Vert *pVC;
	  if (apV[0] == pVA || apV[0] == pVB) {
	    if (apV[1] == pVA || apV[1] == pVB) {
	      pVC = apV[2];
	    }
	    else {
	      pVC = apV[1];
	    }
	  }
	  else {
	    pVC = apV[0];
	  }
	  assert(pVA != pVC && pVB != pVC);
	  
	  if (iOrient2D(pVA, pVB, pVC) == 1) {
	    assert(!apF[i]->pCCellLeft()->qValid());
	    apF[i]->vSetLeftCell(pTCNew);
	  }
	  else {
	    assert(!apF[i]->pCCellRight()->qValid());
	    apF[i]->vSetRightCell(pTCNew);
	  }
	} // Done connecting faces to cells
	pTCNew->vCanonicalizeFaceOrder();
	assert(pTCNew->iFullCheck());
	*new_entity_handle = pTCNew;
	*status = iBase_NEW;
	iNTri++;
      } // Done actually creating the new cell.
      break;
    } // Done creating a triangle
  case iMesh_QUADRILATERAL:
    {
      pM2D->vAllowNonSimplicial();
      // New quadrilaterals to be built from faces or vertices.
      if (lower_order_entity_handles_size != 4) {
	*err = iBase_INVALID_ENTITY_COUNT;
	*status = iBase_CREATION_FAILED;
	*new_entity_handle = NULL;
	break;
      }

      Entity *pE = static_cast<Entity*>(static_cast<Entity*>(lower_order_entity_handles[0]));
      int ET = pE->eEntType();
      Face *apF[4];

      if (ET == iBase_VERTEX) {
	// Build edges first, then use the same code as the edge case
	// from here.
#ifdef NDEBUG
	Vert *pV0 = static_cast<Vert*>(lower_order_entity_handles[0]);
	Vert *pV1 = static_cast<Vert*>(lower_order_entity_handles[1]);
	Vert *pV2 = static_cast<Vert*>(lower_order_entity_handles[2]);
	Vert *pV3 = static_cast<Vert*>(lower_order_entity_handles[3]);
#else
	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;
	  *status = iBase_CREATION_FAILED;
	  *new_entity_handle = NULL;
	  break;
	}
#endif

	// This is a sneaky cheat on identifying verts for edges that
	// works because the verts are cyclically arranged around the
	// quadrilateral.
	Vert *apV[] = {pV0, pV1, pV2, pV3, pV0};
	for (int iF = 0; iF < 4; iF++) {
	  if (!qDoesFaceExist(2, apV + iF, apF[iF])) {
	    apF[iF] = pM2D->pFNewFace();
	  if (!apF[iF]->qValid()) {
	    *err = iBase_ENTITY_CREATION_ERROR;
	    *status = iBase_CREATION_FAILED;
	    *new_entity_handle = NULL;
	    break;
	  }

	  apF[iF]->vAssign(pCInvalidCell, pCInvalidCell, apV[iF], apV[iF+1]);
	  apV[iF]->vSetHintFace(apF[iF]);
	  apV[iF+1]->vSetHintFace(apF[iF]);

	  vAddFaceToHash(apF[iF]);
	  iNLineSeg++;
	  assert((iNLineSeg) == mmFaceHash.size());
	  }
	}
      } // Done creating edges from verts.
      else {
#ifdef NDEBUG
	apF[0] = static_cast<Face*>(lower_order_entity_handles[0]);
	apF[1] = static_cast<Face*>(lower_order_entity_handles[1]);
	apF[2] = static_cast<Face*>(lower_order_entity_handles[2]);
	apF[3] = static_cast<Face*>(lower_order_entity_handles[3]);
#else
	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]));
#endif
      } // Done identifying edges.

      // The following also handles the case of getting handed weird
      // stuff instead of legitimate edge faces.
      if (!apF[0]->qValid() || !apF[1]->qValid() ||
	  !apF[2]->qValid() || !apF[3]->qValid()) {
	*err = iBase_INVALID_ENTITY_HANDLE;
	*status = iBase_CREATION_FAILED;
	*new_entity_handle = NULL;
	break;
      }

      // Verify that the cell doesn't already exist.
      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 {
	// pCNewCell(4) will -always- return a QuadCell* cast to Cell*,
	// or die trying.
	QuadCell *pQCNew = static_cast<QuadCell*>(pM2D->pCNewCell(4));
	if (!pQCNew->qValid()) {
	  *err = iBase_ENTITY_CREATION_ERROR;
	  *status = iBase_CREATION_FAILED;
	  *new_entity_handle = NULL;
	  break;
	}

	pQCNew->vAssign(apF);
	double adApproxCent[] = {0,0};
	// This loop will double count everything.
	for (int i = 0; i < 4; i++) {
	  adApproxCent[0] += apF[i]->pVVert(0)->dX() + apF[i]->pVVert(1)->dX();
	  adApproxCent[1] += apF[i]->pVVert(0)->dY() + apF[i]->pVVert(1)->dY();
	}
	adApproxCent[0] /= 8;
	adApproxCent[1] /= 8;

	Vert VTmp;
	VTmp.vSetCoords(2, adApproxCent);

	for (int i = 0; i < 4; i++) {
	  Vert *pVA = apF[i]->pVVert(0);
	  Vert *pVB = apF[i]->pVVert(1);

	  if (iOrient2D(pVA, pVB, &VTmp) == 1) {
	    assert(!apF[i]->pCCellLeft()->qValid());
	    apF[i]->vSetLeftCell(pQCNew);
	  }
	  else {
	    assert(!apF[i]->pCCellRight()->qValid());
	    apF[i]->vSetRightCell(pQCNew);
	  }
	} // Done connecting faces to cells
	pQCNew->vCanonicalizeFaceOrder();
	assert(pQCNew->iFullCheck());
	*new_entity_handle = pQCNew;
	*status = iBase_NEW;
	iNQuad++;
      } // Done actually creating the new cell.
      break;
    } // Done creating quadrilaterals
    
  case iMesh_TETRAHEDRON:
  case iMesh_PYRAMID:
  case iMesh_PRISM:
  case iMesh_HEXAHEDRON:
  case iMesh_SEPTAHEDRON:
  case iMesh_POLYHEDRON:
  case iMesh_POLYGON:
    *err = iBase_NOT_SUPPORTED;
    break;
  case iMesh_ALL_TOPOLOGIES:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    break;
  }
}

void iMesh_Plane::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;

  Mesh2D *pM2D = dynamic_cast<Mesh2D*>(pM);

  *new_entity_handles_size = 0;
  *status_size = 0;
  switch (new_entity_topology) {
  case iMesh_POINT:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    return;
    break;
  case iMesh_LINE_SEGMENT:
    {
      // We're going to create a new EdgeFace!  That means two vertices.
      if (lower_order_entity_handles_size % 2 != 0) {
	*err = iBase_INVALID_ENTITY_COUNT;
	break;
      }
      
      (*status_size) = (*new_entity_handles_size) =
	lower_order_entity_handles_size/2;
      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);

      for (int iLO = 0, iEH = 0; iLO <= lower_order_entity_handles_size - 1; iEH++)
	{
#ifdef NDEBUG
	  Vert *pV0 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
	  Vert *pV1 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
#else
	  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]));
      
	  if (!pV0->qValid() || !pV1->qValid()) {
	    *err = iBase_INVALID_ENTITY_HANDLE;
	    (*status)[iEH] = iBase_CREATION_FAILED;
	    (*new_entity_handles)[iEH] = NULL;
	    continue;
	  }
#endif
	  Face *pF = pFInvalidFace;
	  Vert *apV[] = {pV0, pV1};
	  if (qDoesFaceExist(2, apV, pF)) {
	    // Return this face as a handle, with status ALREADY_EXISTED
	    assert(pF->qValid());
	    (*new_entity_handles)[iEH] = pF;
	    (*status)[iEH] = iBase_ALREADY_EXISTED;
	  }
	  else {
	    // Create it.
	    pF = pM2D->pFNewFace();
	    if (!pF->qValid()) {
	      *err = iBase_ENTITY_CREATION_ERROR;
	      (*status)[iEH] = iBase_CREATION_FAILED;
	      (*new_entity_handles)[iEH] = NULL;
	      continue;
	    }
	    
	    pF->vAssign(pCInvalidCell, pCInvalidCell, pV0, pV1);
	    pV0->vSetHintFace(pF);
	    pV1->vSetHintFace(pF);
	    (*new_entity_handles)[iEH] = pF;
	    (*status)[iEH] = iBase_NEW;
	    
	    vAddFaceToHash(pF);
	    iNLineSeg++;
	    assert((iNLineSeg) == mmFaceHash.size());
	  }
	}
      break;
    }
  case iMesh_TRIANGLE:
    {
      // New triangles to be built from faces or vertices.
      if (lower_order_entity_handles_size % 3 != 0) {
	*err = iBase_INVALID_ENTITY_COUNT;
	break;
      }
      
      (*status_size) = (*new_entity_handles_size) =
	lower_order_entity_handles_size/3;
      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);

      Entity *pE = static_cast<Entity*>(static_cast<Entity*>(lower_order_entity_handles[0]));
      int ET = pE->eEntType();
      for (int iLO = 0, iEH = 0; iLO <= lower_order_entity_handles_size - 1; iEH++) {
	Face *apF[3];
	
	if (ET == iBase_VERTEX) {
	  // Build edges first, then use the same code as the edge case
	  // from there.
#ifdef NDEBUG
	  Vert *pV0 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
	  Vert *pV1 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
	  Vert *pV2 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
#else
	  Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  
	  if (!pV0->qValid() || !pV1->qValid() || !pV2->qValid()) {
	    *err = iBase_INVALID_ENTITY_HANDLE;
	    (*status)[iEH] = iBase_CREATION_FAILED;
	    (*new_entity_handles)[iEH] = NULL;
	    continue;
	  }
#endif
	  
	  // This is a sneaky cheat on identifying verts for edges that
	  // works because the verts are cyclically arranged around the
	  // triangle.
	  Vert *apV[] = {pV0, pV1, pV2, pV0};
	  for (int iF = 0; iF < 3; iF++) {
	    if (!qDoesFaceExist(2, apV + iF, apF[iF])) {
	      apF[iF] = pM2D->pFNewFace();
	      if (!apF[iF]->qValid()) {
		*err = iBase_ENTITY_CREATION_ERROR;
		(*status)[iEH] = iBase_CREATION_FAILED;
		(*new_entity_handles)[iEH] = NULL;
		continue;
	      }
	      
	      apF[iF]->vAssign(pCInvalidCell, pCInvalidCell, apV[iF], apV[iF+1]);
	      apV[iF]->vSetHintFace(apF[iF]);
	      apV[iF+1]->vSetHintFace(apF[iF]);
	      
	      vAddFaceToHash(apF[iF]);
	      iNLineSeg++;
	      assert((iNLineSeg) == mmFaceHash.size());
	    } // Done creating this new edge
	  }
	} // Done creating/identifying edges from verts.
	else {
#ifdef NDEBUG
	  apF[0] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
	  apF[1] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
	  apF[2] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
#else
	  apF[0] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  apF[1] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  apF[2] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
#endif
	} // Done identifying edges.
	
	// The following also handles the case of getting handed weird
	// stuff instead of legitimate edge faces.
	if (!apF[0]->qValid() || !apF[1]->qValid() || !apF[2]->qValid()) {
	  *err = iBase_INVALID_ENTITY_HANDLE;
	  (*status)[iEH] = iBase_CREATION_FAILED;
	  (*new_entity_handles)[iEH] = NULL;
	  continue;
	}
	
	// Verify that the cell doesn't already exist.
	Cell *pCCom01 = pCCommonCell(apF[0], apF[1]);
	Cell *pCCom12 = pCCommonCell(apF[1], apF[2]);
	Cell *pCCom20 = pCCommonCell(apF[2], apF[0]);
	
	if (pCCom01->qValid() && pCCom12->qValid() && pCCom20->qValid() &&
	    pCCom01 == pCCom12 && pCCom12 == pCCom20) {
	  (*new_entity_handles)[iEH] = pCCom01;
	  (*status)[iEH] = iBase_ALREADY_EXISTED;
	}
	else {
	  // pCNewCell(3) will -always- return a TriCell* cast to Cell*,
	  // or die trying.
	  TriCell *pTCNew = static_cast<TriCell*>(pM2D->pCNewCell(3));
	  if (!pTCNew->qValid()) {
	    *err = iBase_ENTITY_CREATION_ERROR;
	    (*status)[iEH] = iBase_CREATION_FAILED;
	    (*new_entity_handles)[iEH] = NULL;
	    continue;
	  }
	  
	  pTCNew->vAssign(apF);
	  Vert *apV[] = {apF[0]->pVVert(0), apF[0]->pVVert(1),
			 apF[1]->pVVert(0)};
	  if (apV[0] == apV[2] || apV[1] == apV[2])
	    apV[2] = apF[1]->pVVert(1);
	  
	  for (int i = 0; i < 3; i++) {
	    Vert *pVA = apF[i]->pVVert(0);
	    Vert *pVB = apF[i]->pVVert(1);
	    Vert *pVC;
	    if (apV[0] == pVA || apV[0] == pVB) {
	      if (apV[1] == pVA || apV[1] == pVB) {
		pVC = apV[2];
	      }
	      else {
		pVC = apV[1];
	      }
	    }
	    else {
	      pVC = apV[0];
	    }
	    assert(pVA != pVC && pVB != pVC);
	    
	    if (iOrient2D(pVA, pVB, pVC) == 1) {
	      assert(!apF[i]->pCCellLeft()->qValid());
	      apF[i]->vSetLeftCell(pTCNew);
	    }
	    else {
	      assert(!apF[i]->pCCellRight()->qValid());
	      apF[i]->vSetRightCell(pTCNew);
	    }
	  } // Done connecting faces to cells
	  assert(pTCNew->iFullCheck());
	  (*new_entity_handles)[iEH] = pTCNew;
	  (*status)[iEH] = iBase_NEW;
	  iNTri++;
	} // Done actually creating the new cell.
	break;
      } // Done creating a triangle
    }
  case iMesh_QUADRILATERAL:
    {
      pM2D->vAllowNonSimplicial();
      // New quadrilaterals to be built from faces or vertices.
      if (lower_order_entity_handles_size % 4 != 0) {
	*err = iBase_INVALID_ENTITY_COUNT;
	break;
      }
      
      (*status_size) = (*new_entity_handles_size) =
	lower_order_entity_handles_size/4;
      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);
      
      Entity *pE = static_cast<Entity*>(static_cast<Entity*>(lower_order_entity_handles[0]));
      int ET = pE->eEntType();
      
      for (int iLO = 0, iEH = 0; iLO <= lower_order_entity_handles_size - 1; iEH++) {
	Face *apF[4];
	
	if (ET == iBase_VERTEX) {
	  // Build edges first, then use the same code as the edge case
	  // from here.
#ifdef NDEBUG
	  Vert *pV0 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
	  Vert *pV1 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
	  Vert *pV2 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
	  Vert *pV3 = static_cast<Vert*>(lower_order_entity_handles[iLO++]);
#else
	  Vert *pV0 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  Vert *pV1 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  Vert *pV2 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  Vert *pV3 = dynamic_cast<Vert*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  
	  if (!pV0->qValid() || !pV1->qValid() ||
	      !pV2->qValid() || !pV3->qValid()) {
	    *err = iBase_INVALID_ENTITY_HANDLE;
	    (*status)[iEH] = iBase_CREATION_FAILED;
	    (*new_entity_handles)[iEH] = NULL;
	    continue;
	  }
#endif
	  
	  // This is a sneaky cheat on identifying verts for edges that
	  // works because the verts are cyclically arranged around the
	  // quadrilateral.
	  Vert *apV[] = {pV0, pV1, pV2, pV3, pV0};
	  for (int iF = 0; iF < 4; iF++) {
	    if (!qDoesFaceExist(2, apV + iF, apF[iF])) {
	      apF[iF] = pM2D->pFNewFace();
	      if (!apF[iF]->qValid()) {
		*err = iBase_ENTITY_CREATION_ERROR;
		(*status)[iEH] = iBase_CREATION_FAILED;
		(*new_entity_handles)[iEH] = NULL;
		continue;
	      }
	      
	      apF[iF]->vAssign(pCInvalidCell, pCInvalidCell, apV[iF], apV[iF+1]);
	      apV[iF]->vSetHintFace(apF[iF]);
	      apV[iF+1]->vSetHintFace(apF[iF]);
	      
	      vAddFaceToHash(apF[iF]);
	      iNLineSeg++;
	      assert(iNLineSeg == (mmFaceHash.size()));
	    }
	  }
	} // Done creating edges from verts.
	else {
#ifdef NDEBUG
	  apF[0] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
	  apF[1] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
	  apF[2] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
	  apF[3] = static_cast<Face*>(lower_order_entity_handles[iLO++]);
#else
	  apF[0] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  apF[1] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  apF[2] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
	  apF[3] = dynamic_cast<Face*>(static_cast<Entity*>(lower_order_entity_handles[iLO++]));
#endif
	} // Done identifying edges.
	
	  // The following also handles the case of getting handed weird
	  // stuff instead of legitimate edge faces.
	if (!apF[0]->qValid() || !apF[1]->qValid() ||
	    !apF[2]->qValid() || !apF[3]->qValid()) {
	  *err = iBase_INVALID_ENTITY_HANDLE;
	  (*status)[iEH] = iBase_CREATION_FAILED;
	  (*new_entity_handles)[iEH] = NULL;
	  continue;
	}
	
	// Verify that the cell doesn't already exist.
	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_handles)[iEH] = pCCom01;
	  (*status)[iEH] = iBase_ALREADY_EXISTED;
	}
	else {
	  // pCNewCell(4) will -always- return a QuadCell* cast to Cell*,
	  // or die trying.
	  QuadCell *pQCNew = static_cast<QuadCell*>(pM2D->pCNewCell(4));
	  if (!pQCNew->qValid()) {
	    *err = iBase_ENTITY_CREATION_ERROR;
	    (*status)[iEH] = iBase_CREATION_FAILED;
	    (*new_entity_handles)[iEH] = NULL;
	    continue;
	  }
	  
	  pQCNew->vAssign(apF);
	  double adApproxCent[] = {0,0};
	  // This loop will double count everything.
	  for (int i = 0; i < 4; i++) {
	    adApproxCent[0] += apF[i]->pVVert(0)->dX() + apF[i]->pVVert(1)->dX();
	    adApproxCent[1] += apF[i]->pVVert(0)->dY() + apF[i]->pVVert(1)->dY();
	  }
	  adApproxCent[0] /= 8;
	  adApproxCent[1] /= 8;
	  
	  Vert VTmp;
	  VTmp.vSetCoords(2, adApproxCent);
	  
	  for (int i = 0; i < 4; i++) {
	    Vert *pVA = apF[i]->pVVert(0);
	    Vert *pVB = apF[i]->pVVert(1);
	    
	    if (iOrient2D(pVA, pVB, &VTmp) == 1) {
	      assert(!apF[i]->pCCellLeft()->qValid());
	      apF[i]->vSetLeftCell(pQCNew);
	    }
	    else {
	      assert(!apF[i]->pCCellRight()->qValid());
	      apF[i]->vSetRightCell(pQCNew);
	    }
	  } // Done connecting faces to cells
	  assert(pQCNew->iFullCheck());
	  (*new_entity_handles)[iEH] = pQCNew;
	  (*status)[iEH] = iBase_NEW;
	  iNQuad++;
	} // Done actually creating the new cell.
	break;
      } // End of loop
    } // Done creating quadrilaterals
  case iMesh_TETRAHEDRON:
  case iMesh_PYRAMID:
  case iMesh_PRISM:
  case iMesh_HEXAHEDRON:
  case iMesh_SEPTAHEDRON:
  case iMesh_POLYHEDRON:
  case iMesh_POLYGON:
    *err = iBase_NOT_SUPPORTED;
    break;
  case iMesh_ALL_TOPOLOGIES:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
    break;
  }
}

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

  Entity *pE = static_cast<Entity*>(entity_handle);
  switch (pE->eEntTopology()) {
  case iMesh_QUADRILATERAL:
    {
      // 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.
#ifdef NDEBUG
      Cell *pC = static_cast<Cell*>(pE);
#else
      Cell *pC = dynamic_cast<Cell*>(pE);
#endif
      for (int ii = pC->iNumFaces() - 1; ii >= 0; ii--) {
	pC->pFFace(ii)->vRemoveCell(pC);
      }
      pC->vMarkDeleted();
      iNQuad --;
    }
    break;
  case iMesh_TRIANGLE:
    {
      // 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.
#ifdef NDEBUG
      Cell *pC = static_cast<Cell*>(pE);
#else
      Cell *pC = dynamic_cast<Cell*>(pE);
#endif
      for (int ii = pC->iNumFaces() - 1; ii >= 0; ii--) {
	pC->pFFace(ii)->vRemoveCell(pC);
      }
      pC->vMarkDeleted();
      iNTri --;
    }
    break;
  case iMesh_LINE_SEGMENT:
    // These could have cells that depend on them.
    {
#ifdef NDEBUG
      EdgeFace *pEF = static_cast<EdgeFace*>(pE);
#else
      EdgeFace *pEF = dynamic_cast<EdgeFace*>(pE);
      assert(pEF);
#endif

      bool qLeftOkay = !pEF->pCCellLeft()->qValid() ||
	pEF->pCCellLeft()->qDeleted();
      bool qRightOkay = !pEF->pCCellRight()->qValid() ||
	pEF->pCCellRight()->qDeleted();
      if (pEF->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.
#ifdef NDEBUG
	  pBF = static_cast<BFace*>(pEF->pCCellRight());
#else
	  pBF = dynamic_cast<BFace*>(pEF->pCCellRight());
#endif
	}
	else if (qRightOkay) {
	  // Same thing, other side.
#ifdef NDEBUG
	  pBF = static_cast<BFace*>(pEF->pCCellLeft());
#else
	  pBF = dynamic_cast<BFace*>(pEF->pCCellLeft());
#endif
	}
	pEF->vRemoveCell(pBF);
	pBF->vMarkDeleted();
	qLeftOkay = qRightOkay = true;
      }

      if (qLeftOkay && qRightOkay) {
	vRemoveAllReferencesTo(pE);
	pEF->vMarkDeleted();
	vRemoveFaceFromHash(pEF);
	iNLineSeg--;
	assert(iNLineSeg == (mmFaceHash.size()));
      }
      else {
	*err = iBase_INVALID_ENTITY_HANDLE;
      }
      break;
    }
  case iMesh_POINT:
    // Every vertex that's connected to anything must have a valid
    // hint face.
    {
#ifdef NDEBUG
      Vert *pV = static_cast<Vert*>(pE);
#else
      Vert *pV = dynamic_cast<Vert*>(pE);
      assert(pV);
#endif
      if (!pV->qDeleted()) {
	Face *pF = pV->pFHintFace();
	if (!pF->qValid() || !pF->qHasVert(pV) || pF->qDeleted()) {
	  vRemoveAllReferencesTo(pE);
	  pV->vMarkDeleted();
	  iNVert--;
	}
	else {
	  *err = iBase_INVALID_ENTITY_HANDLE;
	}
      }
      break;
    }
  default:
    *err = iBase_INVALID_ENTITY_TOPOLOGY;
  }
}

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_Plane::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) {
    // 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);
    switch (entity_type_requested) {
    case iBase_EDGE:
      // Add faces incident on pV.
      *adj_entity_handles_size = std::count_if(spFNearby.begin(),
					       spFNearby.end(),
					       bind2nd(qHasVert(), pV));
      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));
      break;
    case iBase_FACE:
      *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);
      break;
    case iBase_REGION:
      *adj_entity_handles_size = 0;
      TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		     *adj_entity_handles_size, iBase_EntityHandle);
      break;
    case iBase_ALL_TYPES:
      *adj_entity_handles_size = spCInc.size() +
	std::count_if(spFNearby.begin(), spFNearby.end(),
		      bind2nd(qHasVert(), pV));
      TRY_ARRAY_SIZE(adj_entity_handles, adj_entity_handles_allocated,
		     *adj_entity_handles_size, iBase_EntityHandle);
      {
	iBase_EntityHandle* tmpPtr =
	  std::remove_copy_if(spFNearby.begin(), spFNearby.end(),
			      *adj_entity_handles, bind2nd(qNotHasVert(), pV));
	std::copy(spCInc.begin(), spCInc.end(), tmpPtr);
      }
      break;
    } // End switch
  }
  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_EDGE:
      {
	*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_FACE:
      {
	*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_REGION:
      {
	// 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;
      }
    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);
      }
      break;
    } // End of switch over requested type
  } // The general case (not vertex nor requesting type of current entity
}

void iMesh_Plane::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...

  *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 = 0, 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_REGION: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE = 0, iEOut = 0;
    for (iE = 0; iE < entity_handles_size; iE++) {
      (*offset)[iE] = iEOut;
      // Since there are no regions in 2D, do nothing 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 regions in the 2D 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_EDGE: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE = 0, 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);
	int iCount = std::count_if(spFNearby.begin(), spFNearby.end(),
				   bind2nd(qHasVert(), pV));
	iEOut += iCount;
      }
      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(),
				  pEH, bind2nd(qNotHasVert(), pV));
	assert(pEH <= (*adj_entity_handles + iEOut));
      }
      else {
	pE->vAllFaceHandles(*adj_entity_handles + (*offset)[iE]);
      }
    }
    break;
  }

  case iBase_FACE: {
    // Pre-compute the output size, and put the correct values in the
    // offset array.
    int iE = 0, 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 = 0, 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);
	int iCount = std::count_if(spFNearby.begin(), spFNearby.end(),
				   bind2nd(qHasVert(), pV));
	iEOut += iCount + spCInc.size();
      }
      else {
	iEOut += pE->iNumFaces() + pE->iNumCells() + 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 == 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);

	// Add faces incident on pV.
	pEH = std::remove_copy_if(spFNearby.begin(), spFNearby.end(),
				  pEH, bind2nd(qNotHasVert(), pV));
	pEH = std::copy(spCInc.begin(), spCInc.end(), pEH);
	assert(pEH <= (*adj_entity_handles + iEOut));
      }
      else {
	pE->vAllFaceHandles(pEH);
	pEH += pE->iNumFaces();

	pE->vAllCellHandles(pEH);
	pEH += pE->iNumCells();

	pE->vAllVertHandles(pEH);
	pEH += pE->iNumVerts();
      }
    }
    break;
  } // End of the switch
} // End of iMesh_getEntArrAdj

void 
iMesh_Plane::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 regions are involved at all, return no data.
  if (bridge_entity_type == iBase_REGION ||
      requested_entity_type == iBase_REGION) {
    *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.
    
    // Possible cases:
    // FACE and ALL_TYPES as bridge type turn out to be identical.
    
    //  ->all->all   = everything from neighborhood calc
    //  ->all->face  = 1st adj faces
    //  ->all->edge  = all edges from neighborhood calc
    //  ->all->vert  = all verts from neighborhood calc
    
    //  ->face->all  = all of neighborhood calc except faces
    //  ->face->face = nothing
    //  ->face->edge = all edges from neighborhood calc
    //  ->face->vert = all verts from neighborhood calc
    
    //  ->edge->all  = all of neighborhood calc except edges
    //  ->edge->face = 1st adj faces
    //  ->edge->edge = nothing
    //  ->edge->vert = most verts from neighborhood calc
    
    //  ->vert->all  = 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)) {
      // These six cases do nothing at all.
      *adjacent_entities_size = 0;
    }
    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 iNEdges = spFNearby.size();
      int iNVerts = spVNeigh.size();
      // The following is true both for interior and bdry verts.
      assert(iNCells + iNVerts == iNEdges);
      
      if (requested_entity_type == iBase_FACE) {
	// Two cases (->e->f, ->a->f): return all neighborhood cells
	*adjacent_entities_size = iNCells;
	TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		       *adjacent_entities_size, iBase_EntityHandle);
	std::copy(spCInc.begin(), spCInc.end(), (*adjacent_entities));
      }
      else if (requested_entity_type == iBase_EDGE) {
	// Two cases (->f->e, ->a->e): return all neighborhood edges
	*adjacent_entities_size = iNEdges;
	TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		       *adjacent_entities_size, iBase_EntityHandle);
	std::copy(spFNearby.begin(), spFNearby.end(), (*adjacent_entities));
      }
      else if (requested_entity_type == iBase_VERTEX) {
	if (bridge_entity_type == iBase_EDGE) {
	  // One case (->e->v): return some of the neighborhood verts
	  *adjacent_entities_size = std::count_if(spFNearby.begin(),
						  spFNearby.end(),
						  bind2nd(qHasVert(), pV));
	  TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
			 *adjacent_entities_size, iBase_EntityHandle);
	  std::set<Face*>::iterator iter = spFNearby.begin(),
	    iterEnd = spFNearby.end();
	  int entIndex = 0;
	  for ( ; iter != iterEnd; iter++) {
	    Face *pF = *iter;
	    if (pF->qHasVert(pV)) {
	      if (pF->pVVert(0) == pV) {
		(*adjacent_entities)[entIndex++] = pF->pVVert(1);
	      }
	      else {
		assert(pF->pVVert(1) == pV);
		(*adjacent_entities)[entIndex++] = pF->pVVert(0);
	      }
	    }
	  }
	  assert(entIndex == *adjacent_entities_size);
	}
	else {
	  // Two cases (->f->v, ->a->v): return all neighborhood verts
	  *adjacent_entities_size = iNVerts;
	  TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
			 *adjacent_entities_size, iBase_EntityHandle);
	  std::copy(spVNeigh.begin(), spVNeigh.end(), (*adjacent_entities));
	}
      } // Done with VERTEX output
      else if (requested_entity_type == iBase_ALL_TYPES) {
	// Three more cases (->edge->all, ->face->all, and ->all->all)
	bool qReturnEdges = (bridge_entity_type != iBase_EDGE);
	bool qReturnCells = (bridge_entity_type != iBase_FACE);
	*adjacent_entities_size =
	  iNVerts + (qReturnEdges ? iNEdges : 0)
	  + (qReturnCells ? iNCells : 0);
	
	TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		       *adjacent_entities_size, iBase_EntityHandle);
	int ii = 0;
	if (qReturnCells) {
	  std::copy(spCInc.begin(), spCInc.end(), (*adjacent_entities)+ii);
	  ii += iNCells;
	}
	
	if (qReturnEdges) {
	  std::copy(spFNearby.begin(), spFNearby.end(),
		    (*adjacent_entities)+ii);
	  ii += iNEdges;
	}
	
	// Include the verts regardless
	std::copy(spVNeigh.begin(), spVNeigh.end(), (*adjacent_entities)+ii);
	ii += iNVerts;
	
	assert(ii == *adjacent_entities_size);
      } // Done with requested_type = ALL_TYPES
    } // Done with all cases that require neighborhood calculation
  } // Done with all cases that have input_type == VERTEX
  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 *tmp_ents = NULL;
    int tmp_ents_allocated = 0, tmp_ents_size = 0;
    int *offset = NULL, offset_allocated = 0, offset_size;
    getEntArrAdj(intermediate_ents, intermediate_ents_size,
		 requested_entity_type, &tmp_ents,
		 &tmp_ents_allocated, &tmp_ents_size,
		 &offset, &offset_allocated, &offset_size, err);
    // Need to uniquify this data.
    std::set<iBase_EntityHandle> tmpSet;
    tmpSet.insert(tmp_ents, tmp_ents + tmp_ents_size);
    tmpSet.erase(entity_handle);
    *adjacent_entities_size = tmpSet.size();
    TRY_ARRAY_SIZE(adjacent_entities, adjacent_entities_allocated,
		   *adjacent_entities_size, iBase_EntityHandle);
    std::copy(tmpSet.begin(), tmpSet.end(), (*adjacent_entities));

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

