#ifndef _TRIANGLES_H_
#define _TRIANGLES_H_

#include "./include.h"
#include "./BamgOpts.h"
#include "./BamgMesh.h"
#include "./BamgGeom.h"
#include "./Triangle.h"
#include "./VertexOnGeom.h"
#include "./VertexOnVertex.h"
#include "./VertexOnEdge.h"
#include "./ListofIntersectionTriangles.h"


namespace bamg {

	//classes
	class Geometry;
	class CrackedEdge;
	class QuadTree;
	class SubDomain;

	class Mesh {

		public:

			Geometry                    & Gh;                    // Geometry
			Mesh                        & BTh;                   // Background Mesh Bth== *this =>no background
			long                          NbRef;                 // counter of ref on the this class if 0 we can delete
			long                          maxnbv,nbtx;           // nombre max de sommets , de triangles
			long                          nt,nbv,nbt,nbe;        // nb of legal triangles, nb of vertex, of triangles, of edges with reference,
			long                          NbOfQuad;              // nb of quadrangle
			long                          NbSubDomains;
			long                          NbOutT;                // Nb of oudeside triangle
			long                          NbOfTriangleSearchFind;
			long                          NbOfSwapTriangle;
			BamgVertex                   *vertices;
			long                          NbVerticesOnGeomVertex;
			VertexOnGeom                 *VerticesOnGeomVertex;
			long                          NbVerticesOnGeomEdge;
			VertexOnGeom                 *VerticesOnGeomEdge;
			long                          NbVertexOnBThVertex;
			VertexOnVertex               *VertexOnBThVertex;
			long                          NbVertexOnBThEdge;
			VertexOnEdge                 *VertexOnBThEdge;
			long                          NbCrackedVertices;
			long                         *CrackedVertices;
			long                          NbCrackedEdges;
			CrackedEdge                  *CrackedEdges;
			R2                            pmin,pmax;             // extrema
			double                        coefIcoor;             // coef to integer Icoor1;
			Triangle                     *triangles;
			Edge                         *edges;
			QuadTree                     *quadtree;
			BamgVertex                  **ordre;
			SubDomain                    *subdomains;
			ListofIntersectionTriangles   lIntTria;

			//Constructors/Destructors
			Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts);
			Mesh(double* index,double* x,double* y,int nods,int nels);
			Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator
			Mesh(const Mesh &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
			Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
				try {TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);}
				catch(...) { this->~Mesh(); throw; }
			}
			Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
				try { TriangulateFromGeom0(maxnbv,bamgopts);}
				catch(...) { this->~Mesh(); throw; }
			}
			~Mesh(); 

			//Operators
			const BamgVertex & operator[]  (long i) const { return vertices[i];};
			BamgVertex & operator[](long i) { return vertices[i];};
			const Triangle & operator()  (long i) const { return triangles[i];};
			Triangle & operator()(long i) { return triangles[i];};

			//Methods
			void SetIntCoor(const char * from =0);
			double MinimalHmin() {return 2.0/coefIcoor;}
			double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
			I2 toI2(const R2 & P) const {
				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
			R2 toR2(const I2 & P) const {
				return  R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
			}
			void AddVertex(BamgVertex & s,Triangle * t,Icoor2 *  =0) ;
			void Insert();
			void ForceBoundary();
			void FindSubDomain(int OutSide=0);
			long TriangleReferenceList(long*) const;
			void TriangleIntNumbering(long* renumbering);
			void ShowHistogram() const;
			void CrackMesh(BamgOpts* bamgopts);
			void ShowRegulaty() const;
			void SmoothMetric(double raisonmax) ;
			void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ;
			void MaxSubDivision(double maxsubdiv);
			Edge** MakeGeometricalEdgeToEdge();
			long SplitInternalEdgeWithBorderVertices();
			void MakeQuadrangles(double costheta);
			int  SplitElement(int choice);
			void MakeQuadTree();
			void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1);
			long InsertNewPoints(long nbvold,long & NbTSwap) ; 
			void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
			void ReNumberingVertex(long * renu);
			void SmoothingVertex(int =3,double=0.3);
			Metric MetricAt (const R2 &) const;
			GeometricalEdge* ProjectOnCurve( Edge & AB, BamgVertex &  A, BamgVertex & B,double theta, BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
			long Number(const Triangle & t) const  { return &t - triangles;}
			long Number(const Triangle * t) const  { return t - triangles;}
			long Number(const BamgVertex & t) const  { return &t - vertices;}
			long Number(const BamgVertex * t) const  { return t - vertices;}
			long Number(const Edge & t) const  { return &t - edges;}
			long Number(const Edge * t) const  { return t - edges;}
			long Number2(const Triangle * t) const  { return t - triangles; }
			BamgVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
			Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
			void ReadMesh(double* index,double* x,double* y,int nods,int nels);
			void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
			void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
			void ReadMetric(const BamgOpts* bamgopts);
			void WriteMetric(BamgOpts* bamgopts);
			void AddMetric(BamgOpts* bamgopts);
			void BuildMetric0(BamgOpts* bamgopts);
			void BuildMetric1(BamgOpts* bamgopts);
			void AddGeometryMetric(BamgOpts* bamgopts);
			int  isCracked() const {return NbCrackedVertices != 0;}
			void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
			void ReconstructExistingMesh();

			//Inline methods
			inline  void ReMakeTriangleContainingTheVertex(){
				for (int i=0;i<nbv;i++) vertices[i].vint=0, vertices[i].t=NULL;
				for (int i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
			}
			inline  void  UnMarkUnSwapTriangle(){
				for (int i=0;i<nbt;i++)
				 for(int j=0;j<3;j++)
				  triangles[i].SetUnMarkUnSwap(j);
			  }
			inline  void  SetVertexFieldOn(){
				for (int i=0;i<nbv;i++)                    vertices[i].onGeometry=NULL;
				for (int j=0;j<NbVerticesOnGeomVertex;j++) VerticesOnGeomVertex[j].SetOn();
				for (int k=0;k<NbVerticesOnGeomEdge;k++ )  VerticesOnGeomEdge[k].SetOn();
			}	       
			inline  void   SetVertexFieldOnBTh(){
				for (int i=0;i<nbv;i++)                 vertices[i].onGeometry=NULL;
				for (int j=0;j<NbVertexOnBThVertex;j++) VertexOnBThVertex[j].SetOnBTh();
				for (int k=0;k<NbVertexOnBThEdge;k++ )  VertexOnBThEdge[k].SetOnBTh();
			}

		private:
			void TriangulateFromGeom1(long maxnbv,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
			void TriangulateFromGeom0(long maxnbv,BamgOpts* bamgopts);// the real constructor mesh generator
			void Init(long);

	};

	/*Intermediary*/
	TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ;
	TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b);
	void  swap(Triangle *t1,short a1,
				Triangle *t2,short a2,
				BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2);
	int SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,
				TriangleAdjacent & tt1,Icoor2 & dets1,
				Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
	int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret) ;
	inline TriangleAdjacent Previous(const TriangleAdjacent & ta){
		return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);
	}
	inline TriangleAdjacent Next(const TriangleAdjacent & ta){
		return TriangleAdjacent(ta.t,NextEdge[ta.a]);
	}
	inline  TriangleAdjacent Adj(const TriangleAdjacent & a){
		return  a.Adj();
	}
	inline void Adj(GeometricalEdge * & on,int &i){
		int j=i;i=on->DirAdj[i];on=on->Adj[j];
	}
	inline double qualite(const BamgVertex &va,const BamgVertex &vb,const BamgVertex &vc){
		double ret; 
		I2 ia=va,ib=vb,ic=vc;
		I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;
		Icoor2 deta=Det(ab,ac);
		if (deta <=0) ret = -1;
		else {
			double a = sqrt((double) (ac,ac)),
					 b = sqrt((double) (bc,bc)),
					 c = sqrt((double) (ab,ab)),
					 p = a+b+c;
			double h= Max(Max(a,b),c),ro=deta/p;
			ret = ro/h;
		}
		return ret;
	}

}
#endif