#include <cstdio>
#include <cstring>
#include <cmath>
#include <ctime>

#include "../objects.h"
#include "../../include/include.h"
#include "../../shared/Exceptions/exceptions.h"

namespace bamg {

	static const  Direction NoDirOfSearch=Direction();

	/*Constructors/Destructors*/
	/*FUNCTION Geometry::Geometry(){{{*/
	Geometry::Geometry(){
		Init();
	}
	/*}}}*/
	/*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{*/
	Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
		Init();
		ReadGeometry(bamggeom,bamgopts);
		PostRead();
	}
	/*}}}*/
	/*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{*/
	Geometry::Geometry(const Geometry & Gh) {
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/

		long i;
		*this = Gh;
		NbRef =0;
		quadtree=0;
		vertices = nbv ? new GeomVertex[nbv] : NULL;
		edges = nbe ? new GeomEdge[nbe]:NULL;
		curves= nbcurves ? new Curve[nbcurves]:NULL;
		subdomains = nbsubdomains ? new GeomSubDomain[nbsubdomains]:NULL;
		for (i=0;i<nbe;i++)
		 edges[i].Set(Gh.edges[i],Gh,*this);
		for (i=0;i<nbcurves;i++)
		 curves[i].Set(Gh.curves[i],Gh,*this);
		for (i=0;i<nbsubdomains;i++)
		 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
	}
	/*}}}*/
	/*FUNCTION Geometry::~Geometry(){{{*/
	Geometry::~Geometry() {
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
		if(NbRef>0){   _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
		if(vertices)   delete [] vertices;  vertices=0;
		if(edges)      delete [] edges;     edges=0;
		if(quadtree)   delete    quadtree;  quadtree=0;
		if(curves)     delete [] curves;    curves=0;nbcurves=0;
		if(subdomains) delete [] subdomains;subdomains=0;
		Init();
	}
	/*}}}*/

	/*IO*/
	/*FUNCTION Geometry::ReadGeometry{{{*/
	void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){

		int verbose;
		nbv=0;
		nbe=0;
		nbcurves=0;

		double Hmin = HUGE_VAL;// the infinie value 
		int i,j,k,n,i0,i1,i2,i3;

		/*initialize some variables*/
		verbose= bamgopts->verbose;
		nbv    = bamggeom->VerticesSize[0];
		nbe    = bamggeom->EdgesSize[0];

		//some checks
		if (bamggeom->Vertices==NULL) _error_("the domain provided does not contain any vertex");
		if (bamggeom->Edges==NULL) _error_("the domain provided does not contain any edge");

		//Vertices
		if (bamggeom->Vertices){
			if(verbose>5) _printLine_("      processing Vertices");
			if (bamggeom->VerticesSize[1]!=3) _error_("Vertices should have 3 columns");
			vertices = new GeomVertex[nbv];
			for (i=0;i<nbv;i++) {
				vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
				vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
				vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2];
				vertices[i].DirOfSearch=NoDirOfSearch;
				vertices[i].color =0;
				vertices[i].type=0;
			}
			/*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
			pmin =  vertices[0].r;
			pmax =  vertices[0].r;
			for (i=0;i<nbv;i++) {
				pmin.x = Min(pmin.x,vertices[i].r.x);
				pmin.y = Min(pmin.y,vertices[i].r.y);
				pmax.x = Max(pmax.x,vertices[i].r.x);
				pmax.y = Max(pmax.y,vertices[i].r.y);
			}
			/*Offset pmin and pmax to avoid round-off errors*/
			R2 offset = (pmax-pmin)*0.05;
			pmin -= offset;
			pmax += offset;
			/*coefIcoor is the coefficient used for integer coordinates:
			 *                       (x-pmin.x)
			 * Icoor x = (2^30 -1) ------------ 
			 *                          D
			 * where D is the longest side of the domain (direction x or y)
			 * so that (x-pmin.x)/D is in ]0 1[
			 *
			 * coefIcoor = (2^30 -1)/D
			 */
			coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
			if(coefIcoor<=0) _error_("coefIcoor should be positive");
		}
		else{
			_error_("No BamgVertex provided");
		}

		//Edges
		if (bamggeom->Edges){
			R2      zerovector(0,0);
			double* verticeslength=NULL;

			if(verbose>5) _printLine_("      processing Edges");
			if (bamggeom->EdgesSize[1]!=3) _error_("Edges should have 3 columns");
			edges = new GeomEdge[nbe];

			//initialize verticeslength (sum of the lengths of the edges holding vertex)
			verticeslength = new double[nbv];
			for(i=0;i<nbv;i++) verticeslength[i]=0;

			/*Loop over the edges*/
			for (i=0;i<nbe;i++){

				i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
				i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
				edges[i].v[0]= vertices + i1;     //pointer toward vertex i1 (=&vertices[i1])
				edges[i].v[1]= vertices + i2;     //pointer toward vertex i2
				edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2];

				//get length of edge
				R2     x12=vertices[i2].r-vertices[i1].r;
				double l12=sqrt((x12,x12));
				Hmin=Min(Hmin,l12);

				//initialize other fields
				edges[i].tg[0]=zerovector;
				edges[i].tg[1]=zerovector;
				edges[i].AdjVertexIndex[0] = edges[i].AdjVertexIndex[1] = -1;
				edges[i].Adj[0] = edges[i].Adj[1] = NULL;
				edges[i].type = 0;

				//Cracked?
				if (edges[i].ReferenceNumber!=1) edges[i].SetCracked();

				//prepare metric
				vertices[i1].color++;
				vertices[i2].color++;
				verticeslength[i1] += l12;
				verticeslength[i2] += l12;
			}

			// definition the default of the given mesh size 
			for (i=0;i<nbv;i++) {
				if (vertices[i].color > 0) 
				 vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color);
				else 
				 vertices[i].m=Metric(Hmin);
			}
			delete [] verticeslength;
			
		}
		else{
			_error_("No edges provided");
		}

		//hVertices
		if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){
			if(verbose>5) _printLine_("      processing hVertices");
			for (i=0;i< nbv;i++){
				if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
					vertices[i].m=Metric((double)bamgopts->hVertices[i]);
				}
			}
		}

		//MetricVertices
		if(bamgopts->metric && bamgopts->metric[0]==nbv){
			if(verbose>5) _printLine_("      processing MetricVertices");
			for (i=0;i< nbv;i++) {
				vertices[i].m = Metric((double)bamgopts->metric[i*3+0],(double)bamgopts->metric[i*3+1],(double)bamgopts->metric[i*3+2]);
			}
		}

		//MaxCornerAngle
		if (bamgopts->MaxCornerAngle){
			if(verbose>5) _printLine_("      processing MaxCornerAngle");
			MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;
		}

		//TangentAtEdges
		if (bamggeom->TangentAtEdges){
			if(verbose>5) _printString_("      processing TangentAtEdges");
			if (bamggeom->TangentAtEdgesSize[1]!=4) _error_("TangentAtEdges should have 4 columns");
			int n,i,j,k;
			R2 tg;

			n=bamggeom->TangentAtEdgesSize[0];
			for (k=0;k<n;k++) {
				i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
				j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
				tg.x=bamggeom->TangentAtEdges[k*4+2];
				tg.y=bamggeom->TangentAtEdges[k*4+3];
				if (i<0 || i>=nbe) _error_("TangentAtEdges first index exceeds matrix dimension");
				if (j!=0 && j!=1)  _error_("TangentAtEdges second index should be 1 or 2 only");
				edges[i].tg[j] = tg;
			}
		}

		//Corners
		if(bamggeom->Corners){
			if(verbose>5) _printString_("      processing Corners");
			if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column");
			n=bamggeom->CornersSize[0];
			for (i=0;i<n;i++) {     
				j=(int)bamggeom->Corners[i]-1; //for C indexing
				if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]");
				/*Required => at the same time SetRequired and SetCorner*/
				vertices[j].SetCorner();
				vertices[j].SetRequired();
			}
		}

		//RequiredVertices
		if(bamggeom->RequiredVertices){
			if(verbose>5) _printLine_("      processing RequiredVertices");
			if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column");
			n=bamggeom->RequiredVerticesSize[0];
			for (i=0;i<n;i++) {     
				j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
				if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess  definition: should in [0 " << nbv << "]");
				vertices[j].SetRequired();
			}
		}

		//RequiredEdges
		if(bamggeom->RequiredEdges){
			if(verbose>5) _printLine_("      processing RequiredEdges");
			if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column");
			n=bamggeom->RequiredEdgesSize[0];
			for (i=0;i<n;i++) {     
				j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
				if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]");
				edges[j].SetRequired();  
			}
		}

		//SubDomain
		if(bamggeom->SubDomains){
			if(verbose>5) _printLine_("      processing SubDomains");
			if (bamggeom->SubDomainsSize[1]!=4) _error_("SubDomains should have 4 columns");
			nbsubdomains=bamggeom->SubDomainsSize[0];
			subdomains = new GeomSubDomain[nbsubdomains];
			for (i=0;i<nbsubdomains;i++){
				i0=(int)bamggeom->SubDomains[i*4+0];
				i1=(int)bamggeom->SubDomains[i*4+1];
				i2=(int)bamggeom->SubDomains[i*4+2];
				i3=(int)bamggeom->SubDomains[i*4+3];
				if (i0!=2) _error_("Bad Subdomain definition: first number should be 2 (for Edges)");
				if (i1>nbe || i1<=0) _error_("Bad Subdomain definition: second number should in [1 " << nbe << "] (edge number)");
				subdomains[i].edge=edges + (i1-1);
				subdomains[i].direction = (int) i2;
				subdomains[i].ReferenceNumber = i3;
			}
		}
	}
	/*}}}*/
	/*FUNCTION Geometry::WriteGeometry{{{*/
	void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){

		int verbose;
		int nbreq=0;
		int nbreqv=0;
		int nbtan=0;
		int nbcracked=0;
		int i,count;

		/*Get options*/
		verbose=bamgopts->verbose;

		/*Vertices*/
		if(verbose>5) _printLine_("      writing Vertices");
		bamggeom->VerticesSize[0]=nbv;
		bamggeom->VerticesSize[1]=3;
		if (nbv){
			bamggeom->Vertices=xNew<double>(3*nbv);
			for (i=0;i<nbv;i++){
				bamggeom->Vertices[i*3+0]=vertices[i].r.x;
				bamggeom->Vertices[i*3+1]=vertices[i].r.y;
				bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber();

				//update counters
				if (vertices[i].Required()) nbreqv++;
			}
		}

		/*Edges*/
		if(verbose>5) _printLine_("      writing Edges");
		bamggeom->EdgesSize[0]=nbe;
		bamggeom->EdgesSize[1]=3;
		if (nbe){
			bamggeom->Edges=xNew<double>(3*nbe);
			for (i=0;i<nbe;i++){
				bamggeom->Edges[i*3+0]=GetId(edges[i][0])+1; //back to Matlab indexing
				bamggeom->Edges[i*3+1]=GetId(edges[i][1])+1; //back to Matlab indexing
				bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber;

				//update counters
				if (edges[i].Required()) nbreq++;
				if (edges[i].TgA() && edges[i][0].Corner()) nbtan++;
				if (edges[i].TgB() && edges[i][1].Corner()) nbtan++;
			}
		}

		/*RequiredEdges*/
		if(verbose>5) _printLine_("      writing " << nbreq << " RequiredEdges");
		bamggeom->RequiredEdgesSize[0]=nbreq;
		bamggeom->RequiredEdgesSize[1]=1;
		if (nbreq){
			bamggeom->RequiredEdges=xNew<double>(1*nbreq);
			count=0;
			for (i=0;i<nbe;i++){
				if (edges[i].Required()){
					bamggeom->RequiredEdges[count]=i+1; //back to Matlab indexing
					count=count+1;
				}
			}
		}

		//No corners

		/*RequiredVertices*/
		if(verbose>5) _printLine_("      writing " << nbreqv << " RequiredVertices");
		bamggeom->RequiredVerticesSize[0]=nbreqv;
		bamggeom->RequiredVerticesSize[1]=1;
		if (nbreqv){
			bamggeom->RequiredVertices=xNew<double>(1*nbreqv);
			count=0;
			for (i=0;i<nbv;i++){
				if (vertices[i].Required()){
					bamggeom->RequiredVertices[count]=i+1; //back to Matlab indexing
					count=count+1;
				}
			}
		}

		/*SubDomains*/
		if(verbose>5) _printLine_("      writing SubDomains");
		bamggeom->SubDomainsSize[0]=nbsubdomains;
		bamggeom->SubDomainsSize[1]=4;
		if (nbsubdomains){
			bamggeom->SubDomains=xNew<double>(4*nbsubdomains);
			for (i=0;i<nbsubdomains;i++){
				bamggeom->SubDomains[4*i+0]=2;
				bamggeom->SubDomains[4*i+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
				bamggeom->SubDomains[4*i+2]=subdomains[i].direction;
				bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber;
			}
		}

		/*TangentAtEdges*/
		if(verbose>5) _printLine_("      writing TangentAtEdges");
		bamggeom->TangentAtEdgesSize[0]=nbtan;
		bamggeom->TangentAtEdgesSize[1]=4;
		if (nbtan){
			bamggeom->TangentAtEdges=xNew<double>(4*nbtan);
			count=0;
			for (i=0;i<nbe;i++){
				if (edges[i].TgA() && edges[i][0].Corner()){
					bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
					bamggeom->TangentAtEdges[4*i+1]=1;
					bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x;
					bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y;
				}
				if (edges[i].TgB() && edges[i][1].Corner()){
					bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
					bamggeom->TangentAtEdges[4*i+1]=2;
					bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x;
					bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[1].y;
				}
				count=count+1;
			}
		}
	}
	/*}}}*/

	/*Methods*/
	/*FUNCTION Geometry::Echo {{{*/
	void Geometry::Echo(void){

		_printLine_("Geometry:");
		_printLine_("   nbv  (number of vertices) : " << nbv);
		_printLine_("   nbe  (number of edges)    : " << nbe);
		_printLine_("   nbsubdomains: " << nbsubdomains);
		_printLine_("   nbcurves: " << nbcurves);
		_printLine_("   vertices: " << vertices);
		_printLine_("   edges: " << edges);
		_printLine_("   quadtree: " << quadtree);
		_printLine_("   subdomains: " << subdomains);
		_printLine_("   curves: " << curves);
		_printLine_("   pmin (x,y): (" << pmin.x << " " << pmin.y << ")");
		_printLine_("   pmax (x,y): (" << pmax.x << " " << pmax.y << ")");
		_printLine_("   coefIcoor: " << coefIcoor);
		_printLine_("   MaxCornerAngle: " << MaxCornerAngle);

		return;
	}
	/*}}}*/
	/*FUNCTION Geometry::Init{{{*/
	void Geometry::Init(void){
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/

		NbRef=0;
		nbv=0;
		nbe=0;
		quadtree=NULL;
		curves=NULL;
		edges=NULL;
		vertices=NULL;
		nbsubdomains=0;
		nbcurves=0;
		subdomains=NULL;
		MaxCornerAngle = 10*Pi/180; //default is 10 degres
	}
	/*}}}*/
	/*FUNCTION Geometry::MinimalHmin{{{*/
	double Geometry::MinimalHmin() {
		/* coeffIcoor = (2^30-1)/D
		 * We cannot go beyond hmin = D/2^30 because of the quadtree
		 * hmin is therefore approximately 2/coeffIcoor */
		return 2.0/coefIcoor;
	}/*}}}*/
	/*FUNCTION Geometry::MaximalHmax{{{*/
	double Geometry::MaximalHmax() {
		return Max(pmax.x-pmin.x,pmax.y-pmin.y);
	}/*}}}*/
	/*FUNCTION Geometry::GetId(const GeomVertex &t){{{*/
	long Geometry::GetId(const GeomVertex & t) const  {
		return &t - vertices;
	}/*}}}*/
	/*FUNCTION Geometry::GetId(const GeomVertex * t){{{*/
	long Geometry::GetId(const GeomVertex * t) const  {
		return t - vertices;
	}/*}}}*/
	/*FUNCTION Geometry::GetId(const GeomEdge & t){{{*/
	long Geometry::GetId(const GeomEdge & t) const  {
		return &t - edges;
	}/*}}}*/
	/*FUNCTION Geometry::GetId(const GeomEdge * t){{{*/
	long Geometry::GetId(const GeomEdge * t) const  {
		return t - edges;
	}/*}}}*/
	/*FUNCTION Geometry::GetId(const Curve * c){{{*/
	long Geometry::GetId(const Curve * c) const  {
		return c - curves;
	}/*}}}*/
	/*FUNCTION Geometry::Containing{{{*/
	GeomEdge* Geometry::Containing(const R2 P,  GeomEdge * start) const {
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/

		GeomEdge* on =start,* pon=0;
		// walk with the cos on geometry
		int counter=0;
		while(pon != on){  
			counter++;
			_assert_(counter<100);
			pon = on;
			R2 A= (*on)[0];
			R2 B= (*on)[1];
			R2 AB = B-A;
			R2 AP = P-A;
			R2 BP = P-B;
			if ( (AB,AP) < 0) 
			 on = on->Adj[0];
			else if ( (AB,BP)  > 0) 
			 on = on->Adj[1];
			else
			 return on;
		}
		return on;
	}
	/*}}}*/
	/*FUNCTION Geometry::PostRead{{{*/
	void Geometry::PostRead(){
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/

		long               i,j,k;
		int                jj;
		long              *head_v   = new long[nbv];
		long              *next_p   = new long[2*nbe];
		float             *eangle   = new float[nbe];
		double             eps      = 1e-20;
		BamgQuadtree           quadtree; // build quadtree to find duplicates
		BamgVertex        *v0       = vertices;
		GeomVertex *v0g      = (GeomVertex*) (void*)v0;

		k=0;

		//build quadtree for this geometry
		for (i=0;i<nbv;i++){

			/*build integer coordinates (non unique)
			  these coordinates are used by the quadtree to group
			  the vertices by groups of 5:
			  All the coordinates are transformed to ]0,1[^2
			  then, the integer coordinates are computed using 
			  the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
			vertices[i].i=R2ToI2(vertices[i].r); 

			/*find nearest vertex already present in the quadtree (NULL if empty)*/
			BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 

			/*if there is a vertex found that is to close to vertices[i] -> error*/
			if( v && Norme1(v->r - vertices[i].r) < eps ){
				_printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
				_printLine_("Id: " << i+1);
				_printLine_("Coords: ["<<v->r.x<<" "<<v->r.y<<"] ["<<vertices[i].r.x<<" "<<vertices[i].r.y<<"]");

				delete [] next_p;
				delete [] head_v;
				delete [] eangle;
				_error_("two points of the geometry are very closed to each other (see reference numbers above)");
			}

			/*Add vertices[i] to the quadtree*/
			quadtree.Add(vertices[i]);
		}

		/* Here we use a powerful chaining algorithm
		 *
		 * 1. What is a chaining algorithm?
		 *
		 * If F is a function that goes from i in [0 n] to j in [0 m]
		 * and we want to compute the reciprocal function F-1 of F
		 * (what are the antecedents of a given j in Im(F) )
		 * We use 2 lists:
		 *    head_F[j] that holds the head of lists
		 *    next_F[i] that holds the list of elements that have the same image
		 *
		 * Example:
		 *    i1, i2, ..., ip in [0,n] are all antecedents of a given j in [0 m]
		 *    head_F[j] = ip
		 *    next_F[ip]= ip-1
		 *    ....
		 *    next_F[i2]= i1
		 *    next_F[i1]= -1  //end of the list
		 *
		 * Algorithm:
		 *    for(j=0;j<m;j++)  head_F[j] = -1 //initialization
		 *    for(i=0;i<n;i++){
		 *       j=F[i];
		 *       next_F[i]= head_F[j];
		 *       head_F[j]=i;
		 *    }
		 * 
		 *    Then, we can go through all the elements that have for image j:
		 *    for(i=head_F[j]; i!=-1; i=next_F[i])
		 *    initialization of i by i=head_F[j]
		 *    stop the loop when i=-1 (end of the chain)
		 *    iterate using i=next_F[i] (next element that have for image j)
		 * 
		 * 2. How to use this algorithm here?
		 * 
		 * Here F is a function that associates two vertices v0 and v1 for a given edge E
		 * We want to build the reciprocal function: what are the edges that contains
		 * a vertex v?
		 * To do so, we use the same chaining algorithm but there is a difficulty
		 * coming from the fact that for F we have a couple of vertices and not one 
		 * vertex.
		 * To overcome this difficulty, we use a global indexing exactly like in 
		 * C/C++ so that
		 * a member of a 2-column-table can be described by one index p=i*2+j
		 * i=(int)p/2 line number of p
		 * j=p%2      column number of p
		 *
		 * Algorithm:
		 *    for(i=0;i<nbv;i++)  head_v[i] = -1 //initialization
		 *    for(i=0;i<nbe;i++){
		 *       for(j=0;j<2;j++){
		 *          p=2*i+j;
		 *          v=edges(i,j);
		 *          next_p[p]= head_v[v];
		 *          head_v[v]=p;
		 *       }
		 *    }
		 */

		//initialize head_v as -1
		for (i=0;i<nbv;i++) head_v[i]=-1;
		k=0;
		for (i=0;i<nbe;i++) {
			//compute vector of edge i that goes from vertex 0 to vertex 1
			R2 v10=edges[i].v[1]->r - edges[i].v[0]->r;
			double lv10=Norme2(v10);
			//check that its length is not 0
			if(lv10==0){
				delete [] next_p;
				delete [] head_v;
				delete [] eangle;
				_error_("Length of edge " << i << " is 0");
			}
			//compute angle in [-Pi Pi]
			eangle[i] = atan2(v10.y,v10.x);
			//build chains head_v and next_p
			for (j=0;j<2;j++){
				long v=GetId(edges[i].v[j]);
				next_p[k]=head_v[v];
				head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1;
			}
		}

		//sort head_v by order of increasing edges angle
		for (i=0;i<nbv;i++) {
			int exch=1,ord=0;      

			//exchange vertices position in head_v and next_p till tey are sorted
			while (exch){
				long *p=head_v+i;               
				long *po=p;                     
				long  n=*p;                     
				register float angleold=-1000 ; // angle = - infinity
				ord=0; exch=0;

				// loop over the edges that contain the vertex i (till -1)
				while (n >=0){
					ord++;
					long  i1=n/2;       // i1 = floor (n/2)      -> Edge number
					long  j1=n%2;       // j1 = 1 if n is odd    -> Vertex index for this edge (0 or 1)
					long* pn=next_p+n;

					//Next vertex index
					n=*pn;                       

					//compute angle between horizontal axis and v0->v1
					float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1]; 

					//exchange if the current edge angle is smaller than the previous one
					if (angleold > angle){
						exch=1;
						*pn=*po;  // next_p[n] = n + 1
						*po=*p;   // 
						*p=n;     // next_p[n+1] = n
						po=pn;    // po now point toward pn (invert next and current)
					}

					//else, continue going to the next edge position
					else{                        //  to have : po -> p -> pn
						angleold=angle; // update maximum angle
						po=p;           // po now point toward p  (current position)
						p=pn;           // p  now point toward pn (next position)
					}
				}
			}

			// angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
			if(ord==2) { 
				long  n1 = head_v[i];
				long  n2 = next_p[n1];
				long  i1 = n1/2, i2 = n2/2; // edge number
				long  j1 = n1%2, j2 = n2%2; // vertex in the edge 
				float angle1=  j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
				float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
				float da12 = Abs(angle2-angle1);
				if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
					vertices[i].SetCorner() ; 
				}
				// if the edge type/referencenumber a changing then is SetRequired();
				if (edges[i1].type != edges[i2].type || edges[i1].Required()){
					vertices[i].SetRequired();
				}
				if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
					vertices[i].SetRequired();
				}
			}
			if(ord != 2) {
				vertices[i].SetCorner();
			}

			/*close the list around the vertex to have a circular loop*/
			long no=-1, ne = head_v[i];
			while (ne >=0) ne = next_p[no=ne];        
			if(no>=0) next_p[no] = head_v[i];
		}

		/*Check that the list makes sense (we have all the time the same vertex)
		 * and build adjacent edges*/
		k =0;
		for (i=0;i<nbe;i++){
			for (j=0;j<2;j++){

				long n1 = next_p[k++]; 
				long i1 = n1/2 ,j1=n1%2;

				if( edges[i1].v[j1] != edges[i].v[j]) _error_("Problem while processing edges: check the edge list");

				edges[i1].Adj[j1] = edges + i;
				edges[i1].AdjVertexIndex[j1] = j;
			}
		}

		/* generation of  all the tangents*/
		for (i=0;i<nbe;i++) {
			R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
			double lAB=Norme2(AB);                       // Get length of AB
			double ltg2[2]={0.0};                        // initialize tangent

			//loop over the 2 vertices of the edge
			for (j=0;j<2;j++) {
				R2     tg =edges[i].tg[j];
				double ltg=Norme2(tg);

				//by default, tangent=[0 0]
				if(ltg==0){
					//if the current vertex of the edge is not a corner
					if(!edges[i].v[j]->Corner()){
						/*The tangent is set as the vector between the
						 * previous and next vertices connected to current vertex
						 * normed by the edge length*/
						tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r;
						ltg= Norme2(tg);
						tg = tg *(lAB/ltg);
						ltg= lAB;
					}
					//else:  a Corner no tangent => nothing to do    
				}
				else{
					//tangent has already been computed
					tg = tg*(lAB/ltg),ltg=lAB;
				}

				ltg2[j] = ltg;

				if ((tg,AB)<0) tg = -tg;

				edges[i].tg[j]=tg;
			}
			if (ltg2[0]!=0) edges[i].SetTgA();
			if (ltg2[1]!=0) edges[i].SetTgB();
		} 

		/* generation of  all curves (from corner to corner)*/
		/*We proceed in 2 steps: first allocate, second build*/
		for (int step=0;step<2;step++){

			//unmark all edges
			for (i=0;i<nbe;i++) edges[i].SetUnMark();
			long  nb_marked_edges=0;

			//initialize number of curves
			nbcurves = 0;

			for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
				for (i=0;i<nbe;i++){

					GeomEdge & ei=edges[i];   
					for(j=0;j<2;j++){
						/*If current edge ei is unmarked and (level=1 or vertex i is required (corner)):
						 * we do have the first edge of a new curve*/
						if (!ei.Mark() && (level || ei[j].Required())) { 
							int k0=j,k1;
							GeomEdge   *e=&ei;
							GeomVertex *a=(*e)(k0); // begin 
							if(curves){
								curves[nbcurves].FirstEdge=e;
								curves[nbcurves].FirstVertexIndex=k0;
							}
							int nee=0;
							for(;;){ 
								nee++;
								k1 = 1-k0; // next vertex of the edge 
								e->SetMark();
								nb_marked_edges++;
								e->CurveNumber=nbcurves;
								GeomVertex *b=(*e)(k1);

								//break if we have reached the other end of the curve
								if (a==b || b->Required()){
									if(curves){
										curves[nbcurves].LastEdge=e;
										curves[nbcurves].LastVertexIndex=k1;
									}
									break;
								}
								//else: go to next edge (adjacent)
								else{
									k0 = e->AdjVertexIndex[k1];//  vertex in next edge
									e  = e->Adj[k1]; // next edge
								}
							}
							nbcurves++;
							if(level) a->SetRequired();
						}
					}
				} 
			}
			_assert_(nb_marked_edges && nbe);
			//allocate if first step
			if(step==0) curves=new Curve[nbcurves];
		} 

		/*clean up*/
		delete [] next_p;
		delete [] head_v;
		delete [] eangle;

	}
	/*}}}*/
	/*FUNCTION Geometry::ProjectOnCurve {{{*/
	GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
		/*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/

		double save_s=s;
		int NbTry=0;

retry:

		s=save_s;
		GeomEdge* on=e.GeomEdgeHook;
		if (!on){
			_error_("ProjectOnCurve error message: edge provided should be on geometry");
		}
		if (!e[0].GeomEdgeHook ||  !e[1].GeomEdgeHook){
			_error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
		}

		//Get the two vertices of the edge
		const BamgVertex &v0=e[0];
		const BamgVertex &v1=e[1];

		//Get position of V0, V1 and vector v0->v1
		R2 V0=v0,V1=v1,V01=V1-V0;

		//Get geometrical vertices corresponding to v0 and v1
		VertexOnGeom  vg0=*v0.GeomEdgeHook,  vg1=*v1.GeomEdgeHook;

		//build two pointers towrad current geometrical edge
		GeomEdge *eg0=on, *eg1=on;

		//Get edge direction and swap v0 and v1 if necessary
		R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; 
		int OppositeSens = (V01,AB)<0;
		int direction0=0,direction1=1;
		if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);

		//Compute metric of new vertex as a linear interpolation of the two others
		V.m=Metric(1.0-s,v0,s,v1);

		const int mxe=100;
		GeomEdge* ge[mxe+1];
		int     directionge[mxe+1];
		double  lge[mxe+1];
		int bge=mxe/2,tge=bge;
		ge[bge] = e.GeomEdgeHook;
		directionge[bge]=1;

		while (eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){ 
			if (bge<=0) {
				if(NbTry) {
					_printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
					_printLine_("That bug might come from:");
					_printLine_(" 1)  a mesh edge  containing more than " << mxe/2 << " geometrical edges");
					_printLine_(" 2)  code bug : be sure that we call   Mesh::SetVertexFieldOn() before");
					_printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
					_error_("see above");
				}
				NbTry++;
				goto retry;
			}
			GeomEdge* tmpge = eg0;
			ge[--bge] =eg0 = eg0->Adj[direction0];
			_assert_(bge>=0 && bge<=mxe);
			direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]);
		}
		while (eg1 != (GeomEdge*) vg1  &&  (*eg1)(direction1) != (GeomVertex*) vg1) { 
			if(tge>=mxe ) { 
				_printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
				NbTry++;
				if (NbTry<2) goto retry;
				_printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
				_printLine_("That bug might come from:");
				_printLine_(" 1)  a mesh edge  contening more than " << mxe/2 << " geometrical edges");
				_printLine_(" 2)  code bug : be sure that we call   Mesh::SetVertexFieldOn() before");
				_printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
				_error_("see above");
			}
			GeomEdge* tmpge = eg1;
			ge[++tge] =eg1 = eg1->Adj[direction1];
			directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1];
			_assert_(tge>=0 && tge<=mxe);
		}


		if ((*eg0)(direction0)==(GeomVertex*)vg0)
		 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce

		if ((*eg1)(direction1)==(GeomVertex*)vg1)
		 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1);

		double sg;
		if (eg0 == eg1) { 
			register double s0=vg0,s1=vg1;
			sg =  s0*(1.0-s) +  s*s1;
			on=eg0;
		}
		else {
			R2 AA=V0,BB;
			double s0,s1;
			int i=bge;
			double ll=0;
			for(i=bge;i<tge;i++){
				_assert_(i>=0 && i<=mxe);
				BB =  (*ge[i])[directionge[i]];
				lge[i]=ll += Norme2(AA-BB);
				AA=BB ;}
				lge[tge]=ll+=Norme2(AA-V1); 
				// search the geometrical edge
				_assert_(s<=1.0);
				double ls= s*ll;
				on =0;
				s0 = vg0;
				s1= directionge[bge];
				double l0=0,l1;
				i=bge;
				while (  (l1=lge[i]) < ls ) {
					_assert_(i>=0 && i<=mxe);
					i++,s0=1-(s1=directionge[i]),l0=l1;
				}
				on=ge[i];
				if (i==tge) 
				 s1=vg1;

				s  =(ls-l0)/(l1-l0);
				sg =s0*(1.0-s)+s*s1;    
		} 
		_assert_(on);
		V.r= on->F(sg);
		GV=VertexOnGeom(V,*on,sg);
		return on;
	}
	/*}}}*/
	/*FUNCTION Geometry::R2ToI2{{{*/
	I2 Geometry::R2ToI2(const R2 & P) const {
		/*coefIcoor is the coefficient used for integer coordinates:
		 *                       (x-pmin.x)
		 * Icoor x = (2^30 -1) ------------ 
		 *                          D
		 * where D is the longest side of the domain (direction x or y)
		 * so that (x-pmin.x)/D is in ]0 1[
		 *
		 * coefIcoor = (2^30 -1)/D
		 */
		return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
	}/*}}}*/
	/*FUNCTION Geometry::UnMarkEdges{{{*/
	void Geometry::UnMarkEdges() {
		for (int i=0;i<nbe;i++) edges[i].SetUnMark();
	}/*}}}*/
}