#ifndef GR_VolMesh
#define GR_VolMesh 1

#include "GR_config.h"
#include "GR_List.h"
#include "GR_Mesh.h"
#include "GR_VarArray3D.h"

///
class VolMesh : public Mesh {
private:
  ///
  bool qOKToChangeSurfaceEdges:1, qDoEdgeSwapping:1,
    qStrictPatchChecking:1; 
  ///
  // This angle is used to determine when swapping of surface edges is
  // allowed, and also to decide when a surface edge is actually on a
  // bdry curve or not (in the absence of other info and of surface
  // patch recovery code)
  double dMaxAngleForSurfSwap; // In degrees
  /// Three-dimensional topology
  VATriFace       VATriF;      /// Faces  
  VAQuadFace      VAQuadF;     /// Faces
  VATriBFace      VATriBF;     /// BFaces
  VAQuadBFace     VAQuadBF;    /// BFaces 
  VATetCell       VATet;       /// Cells  
  VAPyrCell       VAPyr;       /// Cells
  VAPrismCell     VAPrism;     /// Cells
  VAHexCell       VAHex;       /// Cells
  VAIntTriBFace   VAIntTriBF;  /// Internal bdry faces
  VAIntQuadBFace  VAIntQuadBF; /// Internal bdry faces
  Bdry3D          *pB3D;       /// Bdry info
  bool            qBdryFromMesh; 
public:
  /// This should be used with care
  VolMesh& operator=(const VolMesh& VM);
  /// 
  VolMesh(const VolMesh&);
  ///
  VolMesh(const char strBaseFileName[], const int iQualMeas = 2,
	  const double dMaxAngle = 0, Bdry3D * const pB3DIn = NULL);
  ///
  VolMesh(const VAVertex& VA, const int iQualMeas = 2,
	  const double dMaxAngle = 0);
  ///
  VolMesh(SurfMesh &SM, const int iQualMeas = 2,
	  const double dRes = 1, const double dGrade = 1,
	  const double dMaxAngle = 0);
  ///
  VolMesh(const int iQualMeas = 2);
  ///
  virtual ~VolMesh();
  ///
private:
  void vCreatePatchesFromBdryFaces(const int aiBCList[]);
public:
  void vMeshInBrick(const double dXMin, const double dXMax,
		    const double dYMin, const double dYMax,
		    const double dZMin, const double dZMax);
  ///
  void vConvertFromCellVert(const int iNCells, const int iNBdryFaces,
			    const int a2iCellVert[][4],
			    const int a2iBFaceVert[][3],
			    const int aiBFaceBC[], int*& aiBCList);
  ///
  void vConvertFromCellVert(const int nTet, const int nPyr, const int nPrism,
			    const int nHex, const int nBdryTri, const int nBdryQuad,
			    const int tetToVert[][4], const int pyrToVert[][5],
			    const int prismToVert[][6], const int hexToVert[][8],
			    const int bdryTriToVert[][3],
			    const int bdryQuadToVert[][4]);
  ///
  eMeshType eType() const {return eVolMesh;}
  ///
  int iNumFaces() const {return VATriF.iSize() + VAQuadF.iSize();}
  ///
  int iNumTriFaces() const {return VATriF.iSize();}
  ///
  int iNumQuadFaces() const {return VAQuadF.iSize();}
  ///
  int iNumBdryFaces() const {return VATriBF.iSize() + VAQuadBF.iSize();}
  ///
  int iNumTriBdryFaces() const {return VATriBF.iSize();}
  ///
  int iNumQuadBdryFaces() const {return VAQuadBF.iSize();}
  ///
  int iNumIntBdryFaces() const
    {return VAIntTriBF.iSize() + VAIntQuadBF.iSize();}
  ///
  int iNumIntTriBdryFaces() const {return VAIntTriBF.iSize();}
  ///
  int iNumIntQuadBdryFaces() const {return VAIntQuadBF.iSize();}
  ///
  int iNumCells() const {return VATet.iSize() + VAPyr.iSize() +
			   VAPrism.iSize() + VAHex.iSize();}
  ///
  int iNumTetCells() const {return VATet.iSize();}
  ///
  int iNumPyrCells() const {return VAPyr.iSize();}
  ///
  int iNumPrismCells() const {return VAPrism.iSize();}
  ///
  int iNumHexCells() const {return VAHex.iSize();}
  ///
  Face*    pFFace(const int i) const;
  ///
  BFace* pBFBFace(const int i) const;
  ///
  BFace* pBFIntBFace(const int i) const;
  ///
  Cell*    pCCell(const int i) const;
  /// Index retrieval
  virtual int iFaceIndex(const Face* const pF) const;
  ///
  virtual int iBFaceIndex(const BFace* const pBF) const;
  ///
  virtual int iCellIndex(const Cell* const pC) const;
  ///
  virtual bool qIsValidEnt(void* pvEnt) const;
  /// Create a flock of null entries
  void vSetupTriFaces(const int i)  {VATriF  .vSetup(i);}
  ///
  void vSetupQuadFaces(const int i) {VAQuadF .vSetup(i);}
  ///
  void vSetupTriBFaces(const int i) {VATriBF .vSetup(i);}
  ///
  void vSetupQuadBFaces(const int i){VAQuadBF.vSetup(i);}
  ///
  void vSetupIntTriBFaces(const int i) {VAIntTriBF .vSetup(i);}
  ///
  void vSetupIntQuadBFaces(const int i){VAIntQuadBF.vSetup(i);}
  ///
  void vSetupTetCells(const int i)  {VATet   .vSetup(i);}
  ///
  void vSetupPyrCells(const int i)  {VAPyr   .vSetup(i);}
  ///
  void vSetupPrismCells(const int i){VAPrism .vSetup(i);}
  ///
  void vSetupHexCells(const int i)  {VAHex   .vSetup(i);}
  /// Add entries to connectivity tables
  Face*    pFNewFace(const int iNV = 3);
  ///
  BFace* pBFNewBFace(const int iNBV = 3, const bool qIsInternal = false);
  ///
  Cell*    pCNewCell(const int iNF = 4, const int iNV = 4);
  ///
  Cell*    pCNewCell(const Cell* const pC); // Returned cell has
						  // same type and
						  // data as old one. 
  /// Create a vertex at the given location.
  virtual Vert* createVert(const double dX, const double dY,
			   const double dZ);
  /// Clean up connectivity tables
  void vPurgeCells();
  ///
  void vPurgeBdryFaces();
  ///
  void vPurgeFaces();
  ///
  void vPurgeVerts();
  ///
  void vPurge();
  /// The following routine adds bdry faces where needed after ITAPS
  /// modification.  Doesn't put BFaces on the right patches, though...
  void vCheckForITAPSBdryErrors();
  /// Ensure that the correct number of verts exist
  void vSetupVerts(const int i) {VAVerts.vSetup(i);}
  ///
  void vAssignInitialLengthScales(double dRes, double dAlpha);
protected:
  virtual void vReorderVerts_RCM(void);
  virtual void vReorderCells(void);
  virtual void vReorderFaces(void);
public:
  /// perimeter area:
  double dBoundarySize() const;
  /// interior volume:
  double dInteriorSize() const;
  ///
  bool qWatertight() const;
  ///
  bool qMeshOrientationOK(const bool qFix = true);
/// Conversion to/from simplicial form
  ///
  void vMakeSimplicial() {
    assert(iNumCells() == iNumTetCells());
    qSimplex = true;
  }

/// Validity checking
public:
  bool qValid() const;
/// Insert points into the mesh
  int iInsertOnEdge(Vert* const pVNew, Vert* const pV0,
			    Vert* const pV1, Face* const pF,
			    const int qSwap);
  ///
  int iInsertOnBdryEdge(Vert* const pVNew, Vert* const pV0,
			Vert* const pV1, Face* const pF,
			const int qSwap);
  ///
  int iInsertOnFace(Vert* const pVNew, Face* const pF,
		    Cell* const pC, const int qSwap);
  ///
  int iInsertInInterior(Vert* const pVNew, Cell* const pC,
			const int qSwap);
  void vInsertWatsonInterior(Vert* const pVNew,
			     const std::set<Cell*>& spCToRemove,
			     const std::set<Face*>& spFToRemove,
			     const std::set<Face*>& spFHull);
  void vInsertWatsonBoundary(Vert* const pVNew,
			     std::set<Cell*>& spCToRemove,
			     std::set<Face*>& spFToRemove,
			     std::set<Face*>& spFHull,
			     const std::set<BFace*>& spBFEncroached);
  // Bdry encroachment check.
  virtual int iQueueEncroachedBdryEntities(InsertionQueue& IQ,
					   const std::set<BFace*>& spBF) const;
  virtual int iQueueEncroachedBdryEntities(InsertionQueue& IQ,
					   const std::set<BFace*>& spBF,
					   const double adLoc[],
					   const double dOffset = 0) const;
public:
  ///
  bool qInsertPoint(const double adPoint[], Cell* const pCGuess,
		    int *const piSwaps, const bool qSwap = true,
		    const bool qForce = false,
		    Vert* pVNew = pVInvalidVert);
  int iInsertPointWatson(WatsonInfo& WI);
  void vCleanupWatson(WatsonInfo& WI);
  // The worst shape quality measure allowed.
  virtual double dWorstAllowedCellShape() const
    {
      switch (eEncType) {
      case eBall:
	return sqrt(0.09375); // sqrt(3) sqrt(2) / 8
      case eLens:
	return sqrt(0.1875);  // sqrt(3) / 4
      case eNewLens:
	return sqrt(0.375);   // sqrt(3) sqrt(2) / 4
	// The mythical best-case scenario would be to use sqrt(3) / 2,
	// because then slivers would be impossible.
      default:
	assert(0);
	return 1;
      }
    }
  bool qInsertPointWatson(const double adPoint[], Cell* const pCSeed);
  bool qInsertPointWatson(Vert* const pVNew, Cell* const pCSeed);
  virtual bool qDeleteInteriorPointsInBall(const Cell *pCSeed,
					   const double adLoc[],
					   const double dRadius,
					   int& iNDeleted);
  virtual void vGetWatsonData(const double adPoint[3],
			      Cell* const pCSeed,
			      std::set<Cell*>& spCSearch,
			      std::set<Face*>& spFHull,
			      std::set<Face*>& spFToRemove,
			      std::set<BFace*>& spBFEncroached) const;
// Mesh optimization routines
/// Sliver fix
public:
  int iRepairBadCells(bool qAllowInsertion = false);
  ///
  int iBreakBadBdryCells(bool qAllowInsertion = false);
  ///
  int iSplitLargeAngle(Cell *pC,
                       const double dSplitAngle,
		       const double * const adDihed = 0,
		       const bool qOnlyInterior = false);
// // Mesh coarsening routines
/// // Mesh smoothing routines  
public:
  void vProtectSmallAngles();
  ///
  virtual bool qDoSwap(const Vert* const pVVertA,
		       const Vert* const pVVertB,
		       const Vert* const pVVertC,
		       const Vert* const pVVertD,
		       const Vert* const pVVertE) const;
  ///
  virtual int iSwap(const int iMaxPasses = 2,
		    const bool qAlreadyMarked = false);
  ///
  int iFaceSwap(Face* const pF);
  ///
  bool qSurfaceEdgeChangesAllowed() const {return qOKToChangeSurfaceEdges;}
  ///
  void vAllowSurfaceEdgeChanges(const double dMaxAngle = 0)
    {
      assert(dMaxAngle >= 0 && dMaxAngle < 180.);
      qOKToChangeSurfaceEdges = true;
      dMaxAngleForSurfSwap = dMaxAngle;
    }
  ///
  double dMaxAngleForSurfaceEdgeChanges() const {return dMaxAngleForSurfSwap;}
  ///
  void vDisallowSurfaceEdgeChanges() {qOKToChangeSurfaceEdges = false;}
  ///
  bool qStrictPatchCheckingUsed() const {return qStrictPatchChecking;}
  ///
  void vSetStrictPatchChecking() {qStrictPatchChecking = true;}
  ///
  void vSetLenientPatchChecking() {qStrictPatchChecking = false;}
  virtual int iNumBdryPatches() {return pB3D->iNumPatches();}
  ///
  bool qBdryFacesOKToSwap(const TriBFaceBase * const pBF0,
			  const TriBFaceBase * const pBF1) const;
  ///
  bool qEdgeSwappingAllowed() const {return qDoEdgeSwapping;}
  ///
  void vAllowEdgeSwapping() {qDoEdgeSwapping = true;}
  ///
  void vDisallowEdgeSwapping() {qDoEdgeSwapping = false;}

  ///
  void vRefineToLengthScale(const bool qSameSize,
			    const int iMaxPasses = 40);
  ///
  virtual bool qRemoveVertByContraction(Vert * const pV, int& iNewSwaps,
			   Vert * apVPreferred[] = NULL,
			   const int iNPref = 0);
  ///
  virtual bool qRemoveVert(Vert * const pV, int& iNewSwaps);
  ///
  virtual void vIdentifyVertexTypesGeometrically() const;
  virtual double dCoarsenConstant() const {return 0.36;}
  ///
  int iSmooth(const int iMaxPasses = 10);
  void vSetSmoothingGoal(const int iGoal)
    {
      SMsetSmoothFunction(pvSmoothData, iGoal+20);
    }
  ///
  int iSmoothVertex(Vert * const pV);
  ///
  bool qDelaunayize();
  ///
  int iReconfigure(TetCell* apTCTets[], TriFace* const pFFace,
		   Vert* const pVVertA, Vert* const pVVertB,
		   Vert* const pVVertC, Vert* const pVVertD,
		   Vert* const pVVertE);
  ///
  int iReconfTet23(TetCell* apTCTets[], TriFace* const pFFace,
		   Vert* const pVVertA, Vert* const pVVertB,
		   Vert* const pVVertC, Vert* const pVVertD,
		   Vert* const pVVertE);
  ///
  int iReconfTet32(TetCell* apTCTets[], TriFace* const pFFace,
		   Vert* const pVVertA, Vert* const pVVertB,
		   Vert* const pVVertC, Vert* const pVVertD,
		   Vert* const pVVertE);
  ///
  int iReconfTet22(TetCell* apTCTets[], TriFace* const pFFace,
		   Vert* const pVVertA, Vert* const pVVertB,
		   Vert* const pVVertC, Vert* const pVVertD,
		   Vert* const pVVertE);
  ///
  int iReconfTet44(TetCell* apTCTets[], TriFace* const pFFace,
		   Vert* const pVVertA, Vert* const pVVertB,
		   Vert* const pVVertC, Vert* const pVVertD,
		   Vert* const pVVertE);
  ///
  int iEdgeSwap3D(Face* const pF, Vert* pVNorth,
		  Vert* pVSouth, Vert* const pVOther,
		  const bool qAllowAnyValid = false,
		  const Vert* const pVWant0 = pVInvalidVert,
		  const Vert* const pVWant1 = pVInvalidVert);
  ///
  int iBdryEdgeSwap3D(Face* const pF, Vert* pVNorth,
                      Vert* pVSouth, Vert* const pVOther);
private:
  double dQualFunc(Vert* const pV0, Vert* const pV1, 
		   Vert* const pV2, Vert* const pV3) const;
 public:
  // Boundary recovery routines
  void vMarkCleanNeighbors(Cell* pC);
  int iRecoverFace(Mesh2D* pM2D, const int iVerts[3],
		   int * const piPointsInserted = NULL,
		   const bool qAllowEdgeSwap = true);
  int iRecoverEdge(const Vert* const pV0,
		   const Vert* const pV1,
		   const bool qInsertionOK = false);
  int iRecoverEdge(const Vert* const pV0,
		   const Vert* const pV1,
		   SurfMesh& SM,
		   const bool qInsertionOK = false);
  void vSplitEdge(const Vert* const pV0,
		  const Vert* const pV1,
		  SurfMesh& SM);
  // Find the cell incident on pV0 that covers the ray towards pV1 (may be
  // a multiple choice question if pV0 and pV1 are adjacent, or if there
  // is a face containing one but not the other and coplanar with both. 
  Cell* pCBeginPipe(const Vert* const pV0, const Vert* const pV1) const;
  int iReamPipe(const Vert* const pVInit,
		const Vert* const pVEnd,
		bool& qHitEdge, Face* apFSkipped[], int& iNSkip);
  void vDeleteExternal(const SurfMesh* pSM);
  void vFindFaceContaining(const Vert * const pV0,
			   const Vert * const pV1,
			   const Vert * const pV2,
			   Face *&pF, bool &qSuccess);
  void vIdentifyBadEdgesAndFaces(const Vert* const pV0,
				 const Vert* const pV1,
				 const Vert* const pV2,
				 std::set<BadEdge>& sBEBadEdges,
				 std::set<Face*>& spFBadFaces);
public:
  /// Mesh quality assessment
  void vClassifyTets() const;
  ///
  friend class VolQual;
};

// Check to see whether a VolMesh correctly conforms to a surface mesh.
bool qIsBdryConstrained(VolMesh* pVM, SurfMesh* pSM, 
			std::set<ConstrainEdge>& sCESubSegs,
			std::set<ConstrainEdge>& sCEInFacets,
			std::set<ConstrainFace>& sCFSubFacets);


/** Unlikely ever to be called by user code. */
void vWriteFile_VolMesh(VolMesh & Mesh,
		        const char strBaseFileName[],
		        const char strExtraFileSuffix[] = "");
///
void vReadFile_VolMesh(const char * const strBaseFileName,
		       int& iNumVerts,
		       int& iNumFaces,
		       int& iNumCells,
		       int& iNumBdryFaces,
                       int& iNumIntBdryFaces,
		       bool& qFaceVert,
		       bool& qCellVert,
		       bool& qFaceCell,
		       bool& qCellFace,
		       bool& qCellRegion,
		       bool& qBFaceFace,
		       bool& qBFaceVert,
		       bool& qBFaceBC,
		       bool& qIntBFaceFace,
		       bool& qIntBFaceVert,
                       bool& qIntBFaceBC,
		       VAVertex& VA,
		       int (*&a2iFaceVert)[3],
		       int (*&a2iFaceCell)[2],
		       int (*&a2iCellVert)[4],
		       int (*&a2iCellFace)[4],
		       int *&aiCellRegion,
		       int *&aiBFaceFace,
		       int *&aiBFaceBC,
		       int (*&a2iBFaceVert)[3],
                       int (*&a2iIntBFaceFace)[2],
                       int *&aiIntBFaceBC,
                       int (*&a2iIntBFaceVert)[3]);
#endif
