#ifndef _MESH2_H_
#define _MESH2_H_

#include "../objects/objects.h"
#include "../shared/shared.h"
#include "../include/macros.h"
#include "../toolkits/toolkits.h"

#include "meshtype.h"

#include "objects/Metric.h"
#include "objects/DoubleAndInt4.h"
#include "objects/Direction.h"
#include "objects/Vertex.h"
#include "objects/TriangleAdjacent.h"
#include "objects/Edge.h"
#include "objects/GeometricalVertex.h"
#include "objects/GeometricalEdge.h"
#include "objects/Curve.h"
#include "objects/Triangle.h"
#include "objects/ListofIntersectionTriangles.h"
#include "objects/GeometricalSubDomain.h"
#include "objects/SubDomain.h"
#include "objects/VertexOnGeom.h"
#include "objects/VertexOnVertex.h"
#include "objects/VertexOnEdge.h"
#include "objects/CrackedEdge.h"
#include "objects/Triangles.h"
#include "objects/Geometry.h"

namespace bamg {
	
	/*INLINE functions{{{1*/
	// to sort in descending order
	template<class T>  inline void  HeapSort(T *c,long n){
		/*Intermediary*/
		int l,j,r,i;
		T   crit;
		c--;                    //the array must starts at 1 and not 0 
		if(n<=1) return;        //return if size <=1
		l=n/2+1;                //initialize l and r
		r=n;
		for(;;){
			if(l<=1){
				crit  =c[r];
				c[r--]=c[1];
				if (r==1){c[1]=crit; return;}
			}
			else  crit = c[--l]; 
			j=l;
			for(;;){
				i=j;
				j=2*j;
				if  (j>r) {c[i]=crit;break;}
				if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
				if (crit < c[j]) c[i]=c[j];       //c[j]  > crit -> stock this large value in i(<j)
				else{c[i]=crit;break;}            //c[j]  < crit -> stock crit in i (<j)
			}
		}
	}
	// to sort in descending order and return ordering
	template<class T>  inline void  HeapSort(int** porder,T* c,int n){
		/*Intermediary*/
		int  l,j,r,i;
		T    crit;
		int  pos;
		int* order=NULL;
		order=(int*)xmalloc(n*sizeof(int));
		for(i=0;i<n;i++) order[i]=i+1;
		c--;                    //the array must starts at 1 and not 0 
		order--;
		if(n<=1) return;        //return if size <=1
		l=n/2+1;                //initialize l and r
		r=n;
		for(;;){
			if(l<=1){
				crit  =c[r]; pos=order[r];
				c[r--]=c[1]; order[r+1]=order[1];
				if (r==1){
					c[1]=crit; order[1]=pos;
					order++;
					*porder=order;
					return;
				}
			}
			else  {crit=c[--l]; pos=order[l];}
			j=l;
			for(;;){
				i=j;
				j=2*j;
				if  (j>r) {c[i]=crit;order[i]=pos;break;}
				if ((j<r) && (c[j] < c[j+1]))j++;
				if (crit < c[j]) {c[i]=c[j];order[i]=order[j];} 
				else{c[i]=crit;order[i]=pos;break;}
			}
		}
	}
	inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){
		return A[0]*( B[1]*C[2]-B[2]*C[1])
		  - A[1]*( B[0]*C[2]-B[2]*C[0])
		  + A[2]*( B[0]*C[1]-B[1]*C[0]);
	}
	inline  TriangleAdjacent Adj(const TriangleAdjacent & a)
	  { return  a.Adj();}

	inline TriangleAdjacent Next(const TriangleAdjacent & ta) 
	  { return TriangleAdjacent(ta.t,NextEdge[ta.a]);}

	inline TriangleAdjacent Previous(const TriangleAdjacent & ta) 
	  { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);}
	inline void Adj(GeometricalEdge * & on,int &i) 
	  {int j=i;i=on->DirAdj[i];on=on->Adj[j];}
	inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc)
	  {
		Real4 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 {
			Real8 a = sqrt((Real8) (ac,ac)),
					b = sqrt((Real8) (bc,bc)),
					c = sqrt((Real8) (ab,ab)),
					p = a+b+c;
			Real8 h= Max(Max(a,b),c),ro=deta/p;
			ret = ro/h;}
			return ret;
	  }
	Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c){
		register  Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y; 
		register  Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y; 
		return  bax*cay - bay*cax;
	}
	inline  TriangleAdjacent FindTriangleAdjacent(Edge &E){
		Vertex * a = E.v[0];
		Vertex * b = E.v[1];

		Triangle * t = a->t;
		int i = a->vint;
		TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
		if (!t || i<0 || i>=3){
			throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
		}
		if ( a!=(*t)(i)){
			throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
		}
		int k=0;
		do { // turn around vertex in direct sens (trigo)
			k++;
			if (k>=20000){
				throw ErrorException(__FUNCT__,exprintf("k>=20000"));
			}
			//  in no crack => ta.EdgeVertex(1) == a otherwise ??? 
			if (ta.EdgeVertex(1) ==  a && ta.EdgeVertex(0) ==  b) return ta; // find 
			ta = ta.Adj();
			if (ta.EdgeVertex(0) ==  a && ta.EdgeVertex(1) ==  b) return ta; // find 
			--ta;
		} while (t != (Triangle *)ta);
		throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found"));
		return TriangleAdjacent(0,0);//for compiler
	}
	inline Vertex* TheVertex(Vertex * a){
		// give a unique vertex with smallest number
		// in case on crack in mesh 
		Vertex * r(a), *rr;
		Triangle * t = a->t;
		int i = a->vint;
		TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
		if (!t || i<0 || i>=3){
			throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
		}
		if ( a!=(*t)(i)){
			throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
		}
		int k=0;
		do { // turn around vertex in direct sens (trigo)
			k++;
			if (k>=20000){
				throw ErrorException(__FUNCT__,exprintf("k>=20000"));
			}
			//  in no crack => ta.EdgeVertex(1) == a
			if ((rr=ta.EdgeVertex(0)) < r) r = rr;
			ta = ta.Adj();
			if ((rr=ta.EdgeVertex(1)) < r) r =rr;
			--ta;
		} while (t != (Triangle*) ta);  
		return r;
	}

	/*}}}1*/

	/*INLINE functions of CLASS GeometricalVertex{{{1*/
	inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){ 
		*this  = rec;
	  }
	/*}}}1*/
	/*INLINE functions of CLASS VertexOnVertex{{{1*/
	inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) { 
		*this = Th.VertexOnBThVertex[i];  
		v = ThNew.vertices + Th.Number(v);

	  }
/*}}}1*/
	/*INLINE functions of CLASS Triangles{{{1*/
	inline  void  Triangles::ReMakeTriangleContainingTheVertex(){
		register Int4 i;
		for ( i=0;i<nbv;i++){
			vertices[i].vint=0;
			vertices[i].t=0;
		  }
		for ( i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
	  }

	inline  void  Triangles::UnMarkUnSwapTriangle()
	  {
		register Int4 i;
		for ( i=0;i<nbt;i++) 
		 for(int  j=0;j<3;j++)
		  triangles[i].SetUnMarkUnSwap(j);
	  }

	inline  void   Triangles::SetVertexFieldOn(){
		for (Int4 i=0;i<nbv;i++) 
		 vertices[i].onGeometry=0;
		for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ ) 
		 VerticesOnGeomVertex[j].SetOn();
		for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ ) 
		 VerticesOnGeomEdge[k].SetOn();
	}	       
	inline  void   Triangles::SetVertexFieldOnBTh(){
		for (Int4 i=0;i<nbv;i++) 
		 vertices[i].onGeometry=0;
		for (Int4 j=0;j<NbVertexOnBThVertex;j++ ) 
		 VertexOnBThVertex[j].SetOnBTh();
		for (Int4 k=0;k<NbVertexOnBThEdge;k++ ) 
		 VertexOnBThEdge[k].SetOnBTh();

	}
	/*}}}1*/
	/*INLINE functions of CLASS Triangle{{{1*/
	inline Triangle* Triangle::Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const
	  {
		// return the other triangle of the quad if a quad or 0 if not a quat
		Triangle * t =0;
		if (link) {
			int a=-1;
			if (TriaAdjSharedEdge[0] & 16 ) a=0;
			if (TriaAdjSharedEdge[1] & 16 ) a=1;
			if (TriaAdjSharedEdge[2] & 16 ) a=2;
			if (a>=0) {
				t = TriaAdjTriangles[a];
				//  if (t-this<0) return 0;
				v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];
				v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];
				v1 = TriaVertices[OppositeEdge[a]];
				v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
			}
		}
		return t;
	  }

	inline   double   Triangle::QualityQuad(int a,int option) const
	  { // first do the logique part 
		double q;
		if (!link || TriaAdjSharedEdge[a] &4)
		 q=  -1;
		else {
			Triangle * t = TriaAdjTriangles[a];
			if (t-this<0) q=  -1;// because we do 2 times 
			else if (!t->link ) q=  -1;
			else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16  || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 )
			 q= -1;
			else if(option) 
			  { 
				const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
				const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
				const Vertex & v1 = *TriaVertices[OppositeEdge[a]];
				const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
				q =  QuadQuality(v0,v1,v2,v3); // do the float part
			  }
			else q= 1;
		}
		return  q;
	  }
	inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew)
	  { 
		*this = rec;
		if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices +  Th.Number(TriaVertices[0]);
		if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices +  Th.Number(TriaVertices[1]);
		if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices +  Th.Number(TriaVertices[2]);
		if(TriaAdjTriangles[0]) TriaAdjTriangles[0] =  ThNew.triangles + Th.Number(TriaAdjTriangles[0]);
		if(TriaAdjTriangles[1]) TriaAdjTriangles[1] =  ThNew.triangles + Th.Number(TriaAdjTriangles[1]);
		if(TriaAdjTriangles[2]) TriaAdjTriangles[2] =  ThNew.triangles + Th.Number(TriaAdjTriangles[2]);
		if (link  >= Th.triangles && link  < Th.triangles + Th.nbt)
		 link = ThNew.triangles + Th.Number(link);
	  }
	inline  Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) {
		Vertex *v=Th->vertices;
		Int4 nbv = Th->nbv;
		if (i<0 || j<0 || k<0){
			throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0"));
		}
		if (i>=nbv || j>=nbv || k>=nbv){
			throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv"));
		}
		TriaVertices[0]=v+i;
		TriaVertices[1]=v+j;
		TriaVertices[2]=v+k;
		TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
		TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
		det=0;
	}
	inline  Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){
		TriaVertices[0]=v0;
		TriaVertices[1]=v1;
		TriaVertices[2]=v2;
		TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
		TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
		if (v0) det=0;
		else {
			det=-1;
			link=NULL;};  
	}
	inline    Real4 Triangle::qualite()
	  {
		return det < 0 ? -1 :  bamg::qualite(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
	  }

	/*}}}1*/
	 /*INLINE functions of CLASS Edge{{{1*/
	inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
	  { 
		*this = Th.edges[i];
		v[0] = ThNew.vertices + Th.Number(v[0]);    
		v[1] = ThNew.vertices + Th.Number(v[1]);
		if (onGeometry) 
		 onGeometry =  ThNew.Gh.edges+Th.Gh.Number(onGeometry);
		if (adj[0]) adj[0] =   ThNew.edges +   Th.Number(adj[0]);
		if (adj[1]) adj[1] =   ThNew.edges +   Th.Number(adj[1]);

	  }

	/*}}}1*/
	 /*INLINE functions of CLASS GeometricalEdge{{{1*/
	inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew)
	  { 
		*this = rec;
		v[0] = GhNew.vertices + Gh.Number(v[0]);    
		v[1] = GhNew.vertices + Gh.Number(v[1]); 
		if (Adj[0]) Adj[0] =  GhNew.edges + Gh.Number(Adj[0]);     
		if (Adj[1]) Adj[1] =  GhNew.edges + Gh.Number(Adj[1]);     
	  }
	/*}}}1*/
	 /*INLINE functions of CLASS Curve{{{1*/
	inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew)
	  {
		*this = rec;
		be = GhNew.edges + Gh.Number(be);    
		ee = GhNew.edges + Gh.Number(ee); 
		if(next) next= GhNew.curves + Gh.Number(next); 
	  }
	/*}}}1*/
	 /*INLINE functions of CLASS SubDomain{{{1*/
	inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
	  {
		*this = Th.subdomains[i];
		if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){
			throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt"));
		}
		head = ThNew.triangles + Th.Number(head) ; 
		if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{
			throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe"));
		}
		edge = ThNew.edges+ Th.Number(edge);
	  }
	/*}}}1*/
	 /*INLINE functions of CLASS GeometricalSubDomain{{{1*/
	inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew)
	  {
		*this = rec;
		edge = Gh.Number(edge) + GhNew.edges;
	  }
	/*}}}1*/
	/*INLINE functions of CLASS Vertex{{{1*/
	Int4 inline Vertex::Optim(int i,int koption){ 
		Int4 ret=0;
		if ( t && (vint >= 0 ) && (vint <3) ){
			ret = t->Optim(vint,koption);
			if(!i){
				t =0; // for no future optime 
				vint= 0;
			}
		}
		return ret;
	}
	inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){ 
		*this  = rec;
	}
	/*}}}1*/
	/*INLINE functions of CLASS VertexOnEdge{{{1*/
	inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
	  {
		*this = Th.VertexOnBThEdge[i];  
		v = ThNew.vertices + Th.Number(v);
	  }
	/*}}}1*/
	/*INLINE functions of CLASS VertexOnGeom{{{1*/
	inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew)
	  {
		*this = rec;  
		mv = ThNew.vertices + Th.Number(mv);
		if (gv)
		 if (abscisse < 0 )
		  gv = ThNew.Gh.vertices + Th.Gh.Number(gv);
		 else
		  ge = ThNew.Gh.edges + Th.Gh.Number(ge);

	  }
	/*}}}1*/
	 /*INLINE functions of CLASS TriangleAdjacent{{{1*/
	inline  void  TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l  ){
		//set Adjacent Triangle of a triangle
		if(t) {
			t->TriaAdjTriangles[a]=ta.t;
			t->TriaAdjSharedEdge[a]=ta.a|l;
		}
		if(ta.t) {
			ta.t->TriaAdjTriangles[ta.a] = t ;
			ta.t->TriaAdjSharedEdge[ta.a] = a| l ;
		}
	}
	inline int  TriangleAdjacent::Locked() const
	  { return t->TriaAdjSharedEdge[a] &4;}
	inline int  TriangleAdjacent::Cracked() const
	  { return t->TriaAdjSharedEdge[a] &32;}
	inline int  TriangleAdjacent::GetAllFlag_UnSwap() const
	  { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap
	inline int  TriangleAdjacent::MarkUnSwap() const
	  { return t->TriaAdjSharedEdge[a] &8;}
	inline void  TriangleAdjacent::SetLock(){ t->SetLocked(a);}
	inline void  TriangleAdjacent::SetCracked() { t->SetCracked(a);}
	inline  TriangleAdjacent TriangleAdjacent::Adj() const
	  { return  t->Adj(a);}
	inline Vertex  * TriangleAdjacent::EdgeVertex(const int & i) const
	  {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; }
	inline Vertex  * TriangleAdjacent::OppositeVertex() const
	  {return t->TriaVertices[bamg::OppositeVertex[a]]; }
	inline Icoor2 &  TriangleAdjacent::det() const
	  { return t->det;}
	int inline TriangleAdjacent::swap()
	  { return  t->swap(a);}
	/*}}}1*/

	/*Other prototypes {{{1*/
	TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ;
	TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b);
	Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside);
	void  swap(Triangle *t1,Int1 a1,
				Triangle *t2,Int1 a2,
				Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2);
	int SwapForForcingEdge(Vertex   *  & pva ,Vertex  * &   pvb ,
				TriangleAdjacent & tt1,Icoor2 & dets1,
				Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
	int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ;
	/*}}}1*/

}
#endif