#ifndef GR_Mesh
#define GR_Mesh 1

#include <vector>
#include <set>

#include "GR_config.h"
#include "GR_Face.h"
#include "GR_BFace.h"
#include "GR_Cell.h"
#include "GR_Quality.h"
#include "GR_VarArrayBase.h"
#include "GR_Vertex.h"
#include "GR_VertTree.h"
#include "OptMS.h"
#include "SMsmooth.h"

enum eSwap {eNone, eDelaunay, eMaxDihed, eMinSine, eFromQual};
typedef double (*GR_LSfunc) (const double *adLoc);  // VC7

/**\brief This is the base class for all mesh data.
 *
 * Vertex data is the primary thing stored internally; connectivity data is
 * declared in the derived classes.
 */
class Mesh {
private:
  /// Operator= disallowed.
  Mesh& operator=(const Mesh&) {assert(0); return(*this);}
protected:
  std::set<Vert*> spVUpdateLS;
  ///
  void *pvSmoothData;
  /// Variable-sized array for storing vertex data.
  VAVertex VAVerts;
  ///
  double dLengthFactor;
  GR_LSfunc dfLengthScale;    // VC7
  /// Grading constant
  double dLipschitzAlpha;
  /// Quality measure.
  mutable Quality *pQ;
public:
  void vSetAlpha(const double dA) {dLipschitzAlpha = dA;}
  double dAlpha() const {return dLipschitzAlpha;}
protected:
  /// Resolution constant
  double dMeshResolution;
protected:
  int iIntSkip;
  bool qBdryFill:1;
  ///
  bool qSwapRecur:1;
  ///
  bool qAllowBdryChanges:1;
  ///
  bool qLengthScaleInitialized:1;
  bool qLengthScaleFromCellSizes:1;
  ///
  bool qSimplex:1;
  ///
  eSwap eSwapMeasure:3;
  ///
  eEncroachType eEncType:3;
/// Constructor and destructor
public:
  enum eMeshType {eMesh2D, eSurfMesh, eVolMesh};
  enum eInsertType {eBadCell, eBdryFaceEncroach, eBdrySegEncroach};
  ///
  Mesh() :
    spVUpdateLS(), pvSmoothData(NULL), VAVerts(),
    dLengthFactor(1), dfLengthScale(NULL), dLipschitzAlpha(1), pQ(NULL),
    dMeshResolution(1), iIntSkip(3), qBdryFill(true),  qSwapRecur(true),
    qAllowBdryChanges(true), qLengthScaleInitialized(false),
    qLengthScaleFromCellSizes(false),
    qSimplex(true), eSwapMeasure(eDelaunay), eEncType(eLens)
    {
      // The value of iIntSkip is the correct default for all meshes.
      // The value of qBdryFill is the proper default for 2D and surface;
      // the 3D constructors override this.
      // Can't initialize Quality at the Mesh level.
      SMinitSmoothing(2, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);
    }
  ///
  Mesh(const VAVertex& VAPts) :
    spVUpdateLS(), pvSmoothData(NULL), VAVerts(VAPts),
    dLengthFactor(1), dfLengthScale(NULL), dLipschitzAlpha(1), pQ(NULL),
    dMeshResolution(1), iIntSkip(3), qBdryFill(true),  qSwapRecur(true),
    qAllowBdryChanges(true), qLengthScaleInitialized(false),
    qLengthScaleFromCellSizes(false),
    qSimplex(true), eSwapMeasure(eDelaunay), eEncType(eLens)
    {
      // The value of iIntSkip is the correct default for all meshes.
      // The value of qBdryFill is the proper default for 2D and surface;
      // the 3D constructors override this.
      SMinitSmoothing(2, 'f', OPTMS_DEFAULT, OPTMS_DEFAULT, &pvSmoothData);
    }
  ///
  Mesh(const Mesh& M) :
    spVUpdateLS(M.spVUpdateLS), pvSmoothData(NULL), VAVerts(M.VAVerts),
    dLengthFactor(M.dLengthFactor), dfLengthScale(M.dfLengthScale),
    dLipschitzAlpha(M.dLipschitzAlpha), pQ(NULL),
    dMeshResolution(M.dMeshResolution), iIntSkip(M.iIntSkip),
    qBdryFill(M.qBdryFill), qSwapRecur(M.qSwapRecur),
    qAllowBdryChanges(M.qAllowBdryChanges),
    qLengthScaleInitialized(M.qLengthScaleInitialized),
    qLengthScaleFromCellSizes(M.qLengthScaleFromCellSizes),
    qSimplex(M.qSimplex),
    eSwapMeasure(M.eSwapMeasure), eEncType(M.eEncType)

    {
      SMsmooth_data *pSMSmooth = static_cast<SMsmooth_data*>(M.pvSmoothData);
      int iFunc = pSMSmooth->smooth_param->function_id;
      int iTech = pSMSmooth->smooth_param->smooth_technique;
      double dAcceptValue = pSMSmooth->smooth_param->lap_accept_value;
      SMconvertToDegrees(static_cast<SMsmooth_data *>(pSMSmooth)->smooth_param
			 ->function_id, &dAcceptValue);
      int iDim = pSMSmooth->dimension;
      SMinitSmoothing(iDim, iTech, iFunc, dAcceptValue, &pvSmoothData);
      SMinitGlobalMinValue(pvSmoothData);

      pQ = new Quality(*M.pQ);
      pQ->vSetMesh(this);
    }
  ///
  virtual ~Mesh()
    {
      if (pvSmoothData) SMfinalizeSmoothing(pvSmoothData);
      if (pQ) delete pQ;
    }
/// Inquiry functions
public:
  void vSetVerts(VAVertex* const pVA) {VAVerts = *pVA;}
  ///
  const VAVertex* pVAGetVerts() const {return (&VAVerts);}
  ///
  int iVertIndex(const Vert* const pV) const { return (VAVerts.iIndex(pV)); }
  ///
  virtual int iFaceIndex(const Face* const pF) const = 0;
  ///
  virtual int iBFaceIndex(const BFace* const pBF) const = 0;
  ///
  virtual int iCellIndex(const Cell* const pC) const = 0;
  ///
  virtual bool qIsValidEnt(void* pvEnt) const = 0;
    ///
  eSwap eSwapType() const {return eSwapMeasure;}
  ///
  void vSetSwapType(const eSwap eSwapIn) {eSwapMeasure = eSwapIn;}
  ///
  bool qBdryChangesAllowed() const {return qAllowBdryChanges;}
  ///
  void vAllowBdryChanges() {qAllowBdryChanges = true;}
  ///
  void vDisallowBdryChanges() {qAllowBdryChanges = false;}
  ///
  bool qSwapRecursionAllowed() const {return qSwapRecur;}
  ///
  void vAllowSwapRecursion() {qSwapRecur = true;}
  ///
  void vDisallowSwapRecursion() {qSwapRecur = false;}
  ///
//void vSetLengthFunc(double (*dfLengthFuncIn)(const double *adLoc))
  void vSetLengthFunc(GR_LSfunc dfLengthFuncIn)   // VC7
    { dfLengthScale = dfLengthFuncIn; dLengthFactor = 1.; }
  ///
  double dLengthScale(const double *adLoc) const
    { return dLengthFactor*dfLengthScale(adLoc); }
  ///
  void vInitLengthScale(const double dRes = 1, const double dGrade = 1);
  void vRecomputeLengthScale();
  ///
  void vUpdateLengthScale(Vert * const pV = pVInvalidVert);
  /// FIX ME:  The following interface is a kludge so that 2D meshing
  //  can easily use the new length scale stuff.  This should go away
  //  once 2D is also using the new insertion code.
  void vUpdateLengthScale(std::set<Vert*>& spV);
  ///
  void vSetLengthScale(Vert * const pV);
  ///
  void vSetLengthScale(const double dLen);
  ///
  void vMarkForLengthScaleCheck(Vert * const pV) {
    spVUpdateLS.insert(pV);
  }
  ///
  void vUpdateInfluence();
  void vCheckEdgeLengths() const;
  virtual eMeshType eType() const = 0;
  /// Return the number of vertices in the database; some may be deleted.
  int iNumVerts() const {return VAVerts.iSize();}

  /// Return the number of faces in the database; some may be deleted.
  virtual int iNumFaces() const = 0;
  ///
  virtual int iNumBdryFaces() const = 0;
  ///
  virtual int iNumIntBdryFaces() const = 0;
  ///
  int iNumTotalBdryFaces() const {return iNumBdryFaces() + iNumIntBdryFaces();}

  /// Return the number of cells in the database; some may be deleted.
  virtual int iNumCells() const = 0;

  /// Return a cell by index.
  virtual Cell* pCCell(const int i) const = 0;
  ///
  const VAVertex* pVAVertexData() const {return &VAVerts;}

  /// Return a vertex by index.
  Vert* pVVert(const int i) const {return VAVerts.pTEntry(i);}

  /// Return a face by index.
  virtual Face* pFFace(const int i) const = 0;

  /// Return the handle of an unused vertex in the database.
  Vert* pVNewVert()
    {
      Vert *pV = VAVerts.pTNewEntry();
      pV->vSetDefaultFlags();
      return (pV);
    }
  /// Create a vertex at the given location.
  virtual Vert* createVert(const double dX, const double dY,
			   const double dZ) = 0;

  ///
  virtual BFace* pBFBFace(const int i) const = 0;
  /// perimeter length or  area:
  virtual double dBoundarySize() const {assert(0); return 0;}
  /// interior   area  or volume:
  virtual double dInteriorSize() const  {assert(0); return 0;}
  ///
  virtual bool qWatertight() const  {assert(0); return 0;}
  /// Conversion to/from simplicial form
  bool qSimplicial() const {return qSimplex;}
  ///
  virtual void vMakeSimplicial() {assert(0);}
  ///
  void vAllowNonSimplicial() {qSimplex = false;}
  ///
  void vSetQualityMeasure(const int iQual) {
    delete pQ;
    pQ = new Quality(this, iQual);
  }
  void vEvaluateQuality() const {pQ->vEvaluate();}
  double dEvaluateQuality(const CellSkel* const pCS) const
    {return pQ->dEvaluate(pCS);}
  void vImportQualityData(const double adData[]) const
    {pQ->vImportQualityData(adData);}
  void vWriteQualityToFile(char strQualFileName[]) const
    {pQ->vWriteToFile(strQualFileName);}
public:
  /// Validity checking
  virtual bool qValid() const = 0;
  /// Mesh optimization routines
  virtual bool qDoSwap(const Vert* const pVVertA,
		       const Vert* const pVVertB,
		       const Vert* const pVVertC,
		       const Vert* const pVVertD,
		       const Vert* const pVVertE = pVInvalidVert) const = 0;
  ///
  virtual int iSwap(const int iMaxPasses = 2,
		    const bool qAlreadyMarked = false) = 0;
  ///
  virtual int iFaceSwap(Face* const pF) = 0;
  ///
  virtual void vPurgeBdryFaces() = 0;
  ///
  virtual void vPurgeFaces() = 0;
  ///
  virtual void vPurgeCells() = 0;
  ///
  virtual void vPurgeVerts() = 0;
  ///
  virtual void vPurge() = 0;
  ///
public:
  ///
  void vReorder(void) {
    vReorderVerts_RCM();
    vReorderCells();
    vReorderFaces();
    for (int iF = 0; iF < iNumFaces(); iF++) {
      Face *pF = pFFace(iF);
      if (pF->qIsBdryFace() && pF->pCCellLeft()->qIsBdryCell()) {
	pF->vInterchangeCellsAndVerts();
      }
    }
  }
protected:
  virtual void vReorderVerts_RCM(void) = 0;
  virtual void vReorderCells(void) = 0;
  virtual void vReorderFaces(void) = 0;
public:
  inline void vSetAllHintFaces() const;
  ///
  inline void vClearAllHintFaces() const;
  // Generic insertion
  virtual bool qInsertPoint(const double adPoint[], Cell* const pCGuess,
			    int * const piSwaps, const bool qSwap = true,
			    const bool qForce = false,
			    Vert* pVNew = pVInvalidVert) = 0;
  virtual int iInsertPointWatson(WatsonInfo&) = 0;
  // Generic Watson hull cleanup, probably only needs to be implemented in 3D.
  virtual void vCleanupWatson(WatsonInfo&) {}
  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 = 0;
  // Use this next function to clear out all points lying within a given
  // ball when doing insertion on / near boundaries.  Returns true iff
  // the ball has no interior points left.
  virtual bool qDeleteInteriorPointsInBall(const Cell *pCSeed,
					   const double adLoc[],
					   const double dRadius,
					   int& iNDeleted) = 0;
  int iQueueEncroachedBdryEntities(InsertionQueue&) const;
  virtual int iQueueEncroachedBdryEntities(InsertionQueue&,
					   const std::set<BFace*>&) const = 0;
  virtual int iQueueEncroachedBdryEntities(InsertionQueue& IQ,
					   const std::set<BFace*>& spBF,
					   const double adLoc[],
					   const double dOffset = 0) const = 0;
  ///
  void vSetEncroachmentType(const eEncroachType eET) { eEncType = eET; }
  eEncroachType eEncroachmentType() const {return eEncType;}

  // The worst shape quality measure (shortest edge to circumradius) allowed.
  virtual double dWorstAllowedCellShape() const = 0;
  // Mesh coarsening routines
private:
  int iMarkChain(Vert * const pVSurf, Vert * const pVFirst,
		 const VertConnect aVC[],
		 const double dMaxDist, const int iSkipNum);
  void vDeactivateChain(Vert * const pVSurf, Vert * const pVFirst,
			const VertConnect aVC[], const double dMaxDist) const;
  void vAddNormalVerts(const VertConnect aVC[],
		       const Vert::VertType VTVertTypeToStartFrom);
  void vAddToCoarseVerts(const VertConnect aVC[]) const;
  void vSelectCoarseVerts(VertConnect aVC[]);
  void vDeactivateAll() const;
  int iActivate(const Vert::VertType VT) const;
  int iInteriorSkip() const {return iIntSkip;}
  bool qBoundaryFill() const {return qBdryFill;}
  virtual Cell* pCBeginPipe(const Vert* const pVLast,
			    const Vert* const pVNew) const = 0;
public:
  void vSetInteriorSkip(const int iNewSkip)
    {
      assert(iNewSkip > 0);
      iIntSkip = iNewSkip;
    }
  void vSetBoundaryFill(const bool qNewFill) {qBdryFill = qNewFill;}
  virtual double dCoarsenConstant() const = 0;
  void vCoarsen(const double dLengthRatio = 2);
  void vCoarsenToLengthScale(const bool qHeal = false);
  virtual void vIdentifyVertexTypesGeometrically() const = 0;
  void vInitConflictGraph(VertConnect aVC[]) const;
  void vUpdateConflictGraph(VertConnect aVC[]) const;
  void vRemoveTaggedVerts();
  virtual int iNumBdryPatches() = 0;
  virtual void vBuildBdryPatches() {}
  virtual bool qRemoveVertByContraction(Vert * const pV, int& iNewSwaps,
					Vert * apVPreferred[] = NULL,
					const int iNPref = 0) = 0;
  virtual bool qRemoveVert(Vert * const pV, int& iNewSwaps) = 0;
  bool qIsLocallyDelaunay() const;
// Mesh smoothing routines
  virtual void vSetSmoothingGoal(const int iGoal) = 0;
  void vSetSmoothingTechnique(const int iTech)
    {
      SMsetSmoothTechnique(pvSmoothData, iTech);
    }
  void vSetSmoothingThreshold(const double dThresh)
    {
      SMsetSmoothThreshold(pvSmoothData, dThresh);
    }
  void vResetThreshold()
    {
      // For floating threshold after at least one smoothing pass, this
      // will reset the threshold using the global minimum value from
      // the last pass.  In other cases, the threshold value remains
      // unchanged.
      double dThresh =
	(static_cast<SMsmooth_data*>(pvSmoothData))->
	smooth_param->lap_accept_value;
      SMconvertToDegrees((static_cast<SMsmooth_data *>(pvSmoothData))
			 ->smooth_param->function_id, &dThresh);
      SMsetSmoothThreshold(pvSmoothData, dThresh);
    }
  virtual int iSmooth(const int iMaxPasses = 10) = 0;
// Mesh adaptation and related routines
  virtual bool qDelaunayize() = 0;
  virtual void vAdaptToLengthScale(const eEncroachType eET);
};

/// Return info about entities in the one-cell ball around the vertex.
void vNeighborhood(const Vert* const pVert,
		   std::set<Cell*>& spCInc,
		   std::set<Vert*>& spVNeigh,
		   std::set<BFace*>* pspBFInc = NULL,
		   bool *pqBdryVert = NULL,
		   std::set<Face*>* pspFNearby = NULL);

#ifndef BROKEN_INLINE

//@@ Compute hint faces for all vertices
inline void Mesh::vClearAllHintFaces() const
{
  int i;
  for (i = 0; i < iNumVerts(); i++)
    pVVert(i)->vClearHintFace();
}

//@@ Compute hint faces for all vertices
inline void Mesh::vSetAllHintFaces() const
{
  // Have to clear vertex hints -before- checking validity, else bad
  // hints will cause an irrelevant assertion.
#ifndef NDEBUG
  vClearAllHintFaces();
  assert(qValid());
#endif

  int i;
  for (i = 0; i < iNumFaces(); i++) {
    Face* pF = pFFace(i);
    if (pF->qDeleted() || !(pF->iFullCheck())) continue;
    for (int ii = pF->iNumVerts() - 1; ii >= 0; ii --) {
      Vert* pV = pF->pVVert(ii);
      if (pV->qValid()) {
	  pV->vSetHintFace(pF);
      }
    }
  }

  // Bdry faces are prefered as vertex hints, so overwrite as needed.
  for (i = 0; i < iNumBdryFaces(); i++) {
    BFace* pBF = pBFBFace(i);
    if (pBF->qDeleted()) continue;
    Face* pF = pBF->pFFace(0);
    if (pF->qDeleted() || !(pF->iFullCheck())) continue;
    for (int ii = pF->iNumVerts() - 1; ii >= 0; ii --) {
      Vert* pV = pF->pVVert(ii);
      if (pV->qValid()) {
	pV->vSetHintFace(pF);
      }
    }
  }
}
#endif // Done with inlines
#endif // Done with class Mesh
