Index: issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 21821)
+++ issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 21837)
@@ -980,55 +980,4 @@
 
 	/*Methods*/
-	void Mesh::AddGeometryMetric(BamgOpts* bamgopts){/*{{{*/
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectGeomMetric)*/
-
-		/*Get options*/
-		double anisomax = bamgopts->anisomax;
-		double errg     = bamgopts->errg;
-
-		double ss[2]={0.00001,0.99999};
-		double errC = 2*sqrt(2*errg);
-		double hmax = Gh.MaximalHmax();
-		double hmin = Gh.MinimalHmin();
-
-		//check that hmax is positive
-		_assert_(hmax>0); 
-
-		//errC cannot be higher than 1
-		if(errC>1) errC=1;
-
-		//Set all vertices to "on"
-		SetVertexFieldOn();
-
-		//loop over all the vertices on edges
-		for(int  i=0;i<nbe;i++){
-			for (int j=0;j<2;j++){
-
-				BamgVertex V;
-				VertexOnGeom GV;
-				Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
-
-				GeomEdge* eg = GV;
-				double s = GV;
-				R2 tg;
-				double  R1= eg->R1tg(s,tg);
-				double  ht=hmax;
-				// err relative to the length of the edge
-				if (R1>1.0e-20) {  
-					ht = Min(Max(errC/R1,hmin),hmax);
-				}
-				double hn=Min(hmax,ht*anisomax);
-
-				if (ht<=0 || hn<=0){
-					_error_("ht<=0 || hn<=0");
-				}
-				EigenMetric Vp(1/(ht*ht),1/(hn*hn),tg);
-				Metric MVp(Vp);
-				edges[i][j].m.IntersectWith(MVp);
-			}
-		}
-		// the problem is for the vertex on vertex 
-	}
-	/*}}}*/
 	void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/
 		//  Hessiantype = 0 =>  H is computed using double L2 projection
@@ -1231,10 +1180,8 @@
 		int i,j,k,kk,it,jt;
 		int    verbose=0;
-		double cutoffradian=10*Pi/180;
 
 		/*Recover options*/
-		if (bamgopts){
+		if(bamgopts){
 			verbose=bamgopts->verbose;
-			cutoffradian=bamgopts->MaxCornerAngle*Pi/180;
 		}
 
@@ -1243,10 +1190,5 @@
 
 		//check that the mesh is not empty
-		if (nbt<=0 || nbv <=0 ) {
-			_error_("nbt or nbv is negative (Mesh empty?)");
-		}
-
-		//Gh is the geometry of the mesh (this), initialize MaxCornerAngle
-		if (cutoffradian>=0) Gh.MaxCornerAngle = cutoffradian;
+		if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)");
 
 		/*Construction of the edges*/
@@ -1447,7 +1389,5 @@
 
 		//check that nbsubdomains is empty
-		if (nbsubdomains){
-			_error_("nbsubdomains should be 0");
-		}
+		if(nbsubdomains) _error_("nbsubdomains should be 0");
 		nbsubdomains=0;
 
@@ -2649,49 +2589,49 @@
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/
 
-		/* initialize random seed: */
-		srand(19999999);
-
 		/*Initialize fields*/
-		NbRef=0;
-		quadtree=NULL;
-		nbv=0;
-		nbt=0;
-		nbe=0;
-		edges=NULL;
-		nbsubdomains=0;
-		subdomains=NULL;
-		maxnbv=maxnbv_in;
-		maxnbt=2 *maxnbv_in-2;
-		NbVertexOnBThVertex=0;
-		VertexOnBThVertex=NULL;
-		NbVertexOnBThEdge=0;
-		VertexOnBThEdge=NULL;
-		NbCrackedVertices=0;
-		CrackedVertices =NULL;
-		NbCrackedEdges =0;
-		CrackedEdges =NULL;
-		NbVerticesOnGeomVertex=0;
-		VerticesOnGeomVertex=NULL;
-		NbVerticesOnGeomEdge=0;
-		VerticesOnGeomEdge=NULL;
+		this->NbRef                  = 0;
+		this->quadtree               = NULL;
+		this->nbv                    = 0;
+		this->nbt                    = 0;
+		this->nbe                    = 0;
+		this->edges                  = NULL;
+		this->nbsubdomains           = 0;
+		this->subdomains             = NULL;
+		this->maxnbv                 = maxnbv_in;
+		this->maxnbt                 = 2 *maxnbv_in-2;
+		this->NbVertexOnBThVertex    = 0;
+		this->VertexOnBThVertex      = NULL;
+		this->NbVertexOnBThEdge      = 0;
+		this->VertexOnBThEdge        = NULL;
+		this->NbCrackedVertices      = 0;
+		this->CrackedVertices        = NULL;
+		this->NbCrackedEdges         = 0;
+		this->CrackedEdges           = NULL;
+		this->NbVerticesOnGeomVertex = 0;
+		this->VerticesOnGeomVertex   = NULL;
+		this->NbVerticesOnGeomEdge   = 0;
+		this->VerticesOnGeomEdge     = NULL;
+
+		/*Initialize random seed*/
+		this->randomseed = 1;
 
 		/*Allocate if maxnbv_in>0*/
 		if(maxnbv_in){
-			vertices=new BamgVertex[maxnbv];
-			_assert_(vertices);
-			orderedvertices=new BamgVertex* [maxnbv];
-			_assert_(orderedvertices);
-			triangles=new Triangle[maxnbt];
-			_assert_(triangles);
+			this->vertices=new BamgVertex[this->maxnbv];
+			this->orderedvertices=new BamgVertex* [this->maxnbv];
+			this->triangles=new Triangle[this->maxnbt];
+         _assert_(this->vertices);
+         _assert_(this->orderedvertices);
+			_assert_(this->triangles);
 		}
 		else{
-			vertices=NULL;
-			orderedvertices=NULL;
-			triangles=NULL;
-			maxnbt=0;
+			this->vertices        = NULL;
+			this->orderedvertices = NULL;
+			this->triangles       = NULL;
+			this->maxnbt          = 0;
 		} 
 	}
 	/*}}}*/
-	void Mesh::Insert(bool random) {/*{{{*/
+	void Mesh::Insert(){/*{{{*/
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
 
@@ -2740,6 +2680,5 @@
 		//Get Prime number
 		const long PrimeNumber= BigPrimeNumber(nbv);
-		int temp = rand(); if(!random) temp = 756804110;
-		int  k0=temp%nbv; 
+		int  k0=this->RandomNumber(nbv);
 
 		//Build orderedvertices
@@ -2825,5 +2764,5 @@
 	}
 	/*}}}*/
-	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap,bool random) {/*{{{*/
+	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
 
@@ -2845,7 +2784,5 @@
 		/*construction of a random order*/
 		const long PrimeNumber= BigPrimeNumber(nbv)  ;
-		//remainder of the division of rand() by nbvnew
-		int  temp = rand(); if(!random) temp = 1098566905;
-		long k3 = temp%nbvnew;
+		long k3 = long(this->RandomNumber(nbvnew));
 		//loop over the new points
 		for (int is3=0; is3<nbvnew; is3++){
@@ -3088,5 +3025,5 @@
 			//if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");
 			Bh.CreateSingleVertexToTriangleConnectivity();     
-			InsertNewPoints(nbvold,NbTSwap,bamgopts->random);
+			InsertNewPoints(nbvold,NbTSwap);
 		}
 		else Bh.CreateSingleVertexToTriangleConnectivity();     
@@ -3146,5 +3083,5 @@
 			}// for triangle   
 
-			if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break;
+			if (!InsertNewPoints(nbvold,NbTSwap)) break;
 			for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
 			Headt = nbt; // empty list 
@@ -4095,5 +4032,5 @@
 
 		/*Insert Vertices*/
-		Insert(true);
+		Insert();
 	}
 	/*}}}*/
@@ -4396,5 +4333,5 @@
 		if (verbose>3) _printf_("      Creating initial Constrained Delaunay Triangulation...\n");
 		if (verbose>3) _printf_("         Inserting boundary points\n");
-		Insert(bamgopts->random);
+		Insert();
 
 		//Force the boundary
@@ -4722,5 +4659,5 @@
 		if (verbose>3) _printf_("      Creating initial Constrained Delaunay Triangulation...\n");
 		if (verbose>3) _printf_("         Inserting boundary points\n");
-		Insert(bamgopts->random);
+		Insert();
 
 		//Force the boundary
@@ -4735,4 +4672,173 @@
 		NewPoints(BTh,bamgopts,KeepVertices) ;
 		if (verbose>4) _printf_("      -- current number of vertices = " << nbv << "\n");
+	}
+	/*}}}*/
+	int  Mesh::RandomNumber(int max){/*{{{*/
+		/*  Generate a random number from 0 to max-1 using linear congruential
+		 *  random number generator*/
+
+		this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l;
+		return this->randomseed / (714025l / max + 1);
+	} /*}}}*/
+	int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret)  { /*{{{*/
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/
+
+		int NbSwap =0;
+		if (!a.t || !b.t){ // the 2 vertex is in a mesh
+			_error_("!a.t || !b.t");
+		}
+		int k=0;
+		taret=AdjacentTriangle(0,0); // erreur 
+
+		AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]);
+		BamgVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
+		// we turn around a in the  direct direction  
+
+		long long det2 = v2 ? det(*v2,a,b): -1 , det1;
+		if(v2) // normal case 
+		 det2 = det(*v2,a,b);
+		else { // no chance infini vertex try the next
+			tta= Previous(Adj(tta));
+			v2 = tta.EdgeVertex(0);
+			vbegin =v2;
+			if (!v2){
+				_error_("!v2");
+			}
+			det2 = det(*v2,a,b);
+		}
+
+		while (v2 != &b) {
+			AdjacentTriangle tc = Previous(Adj(tta));    
+			v1 = v2; 
+			v2 = tc.EdgeVertex(0);
+			det1 = det2;
+			det2 =  v2 ? det(*v2,a,b): det2; 
+
+			if((det1 < 0) && (det2 >0)) { 
+				// try to force the edge 
+				BamgVertex * va = &a, *vb = &b;
+				tc = Previous(tc);
+				if (!v1 || !v2){
+					_error_("!v1 || !v2");
+				}
+				long long detss = 0,l=0;
+				while ((this->SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
+				 if(l++ > 10000000) {
+					 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
+				 }
+				BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
+				if (((aa == &a ) && (bb == &b)) ||((bb ==  &a ) && (aa == &b))){
+					tc.SetLock();
+					a.Optim(1,0);
+					b.Optim(1,0);
+					taret = tc;
+					return NbSwap;
+				}
+				else 
+				  {
+					taret = tc;
+					return -2; // error  boundary is crossing
+				  }
+			}
+			tta = tc;
+			k++;
+			if (k>=2000){
+				_error_("k>=2000");
+			}
+			if ( vbegin == v2 ) return -1;// error 
+		}
+
+		tta.SetLock();
+		taret=tta;
+		a.Optim(1,0);
+		b.Optim(1,0);
+		return NbSwap; 
+	}
+	/*}}}*/
+	int Mesh::SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/
+		// l'arete ta coupe l'arete pva pvb
+		// de cas apres le swap sa coupe toujours
+		// on cherche l'arete suivante 
+		// on suppose que detsa >0 et detsb <0
+		// attention la routine echange pva et pvb 
+
+		if(tt1.Locked()) return 0; // frontiere croise 
+
+		AdjacentTriangle tt2 = Adj(tt1);
+		Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
+		short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
+		if ( a1<0 || a1>=3 ){
+			_error_("a1<0 || a1>=3");
+		}
+
+		BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
+		BamgVertex & s1= (*t1)[OppositeVertex[a1]];
+		BamgVertex & s2= (*t2)[OppositeVertex[a2]];
+
+		long long dets2 = det(*pva,*pvb,s2);
+		long long det1=t1->det , det2=t2->det ;
+		long long detT = det1+det2;
+		if ((det1<=0 ) || (det2<=0)){
+			_error_("(det1<=0 ) || (det2<=0)");
+		}
+		if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb
+			_error_("(detsa>=0) || (detsb<=0)");
+		}
+		long long ndet1 = bamg::det(s1,sa,s2);
+		long long ndet2 = detT - ndet1;
+
+		int ToSwap =0; //pas de swap
+		if ((ndet1 >0) && (ndet2 >0)) 
+		  { // on peut swaper  
+			if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
+			 ToSwap =1; 
+			else // swap alleatoire 
+			 if (BinaryRand()) 
+			  ToSwap =2; 
+		  }
+		if (ToSwap) NbSwap++,
+		 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
+
+		int ret=1;
+
+		if (dets2 < 0) {// haut
+			dets1 = ToSwap ? dets1 : detsa ;
+			detsa = dets2; 
+			tt1 =  Previous(tt2) ;}
+		else if (dets2 > 0){// bas 
+			dets1 = ToSwap ? dets1 : detsb ;
+			detsb = dets2;
+			//xxxx tt1 = ToSwap ? tt1 : Next(tt2);
+			if(!ToSwap) tt1 =  Next(tt2);
+		}
+		else { // changement de direction 
+			ret = -1;
+			Exchange(pva,pvb);
+			Exchange(detsa,detsb);
+			Exchange(dets1,dets2);
+			Exchange(tt1,tt2);
+			dets1=-dets1;
+			dets2=-dets2;
+			detsa=-detsa;
+			detsb=-detsb;
+
+			if(ToSwap){
+				if (dets2 < 0) {// haut
+					dets1 = (ToSwap ? dets1 : detsa) ;
+					detsa = dets2; 
+					tt1 =  Previous(tt2) ;}
+				else if(dets2 > 0){// bas 
+					dets1 = (ToSwap ? dets1 : detsb) ;
+					detsb =  dets2;
+					if(!ToSwap) tt1 =  Next(tt2);
+				}
+				else {// on a fin ???
+					tt1 = Next(tt2);
+					ret =0;}
+			}
+
+		}
+		return ret;
 	}
 	/*}}}*/
@@ -4776,79 +4882,4 @@
 							return edge;
 		} 
-	}
-	/*}}}*/
-	int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret)  { /*{{{*/
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/
-
-		int NbSwap =0;
-		if (!a.t || !b.t){ // the 2 vertex is in a mesh
-			_error_("!a.t || !b.t");
-		}
-		int k=0;
-		taret=AdjacentTriangle(0,0); // erreur 
-
-		AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]);
-		BamgVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
-		// we turn around a in the  direct direction  
-
-		long long det2 = v2 ? det(*v2,a,b): -1 , det1;
-		if(v2) // normal case 
-		 det2 = det(*v2,a,b);
-		else { // no chance infini vertex try the next
-			tta= Previous(Adj(tta));
-			v2 = tta.EdgeVertex(0);
-			vbegin =v2;
-			if (!v2){
-				_error_("!v2");
-			}
-			det2 = det(*v2,a,b);
-		}
-
-		while (v2 != &b) {
-			AdjacentTriangle tc = Previous(Adj(tta));    
-			v1 = v2; 
-			v2 = tc.EdgeVertex(0);
-			det1 = det2;
-			det2 =  v2 ? det(*v2,a,b): det2; 
-
-			if((det1 < 0) && (det2 >0)) { 
-				// try to force the edge 
-				BamgVertex * va = &a, *vb = &b;
-				tc = Previous(tc);
-				if (!v1 || !v2){
-					_error_("!v1 || !v2");
-				}
-				long long detss = 0,l=0;
-				while ((SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
-				 if(l++ > 10000000) {
-					 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
-				 }
-				BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
-				if (((aa == &a ) && (bb == &b)) ||((bb ==  &a ) && (aa == &b))){
-					tc.SetLock();
-					a.Optim(1,0);
-					b.Optim(1,0);
-					taret = tc;
-					return NbSwap;
-				}
-				else 
-				  {
-					taret = tc;
-					return -2; // error  boundary is crossing
-				  }
-			}
-			tta = tc;
-			k++;
-			if (k>=2000){
-				_error_("k>=2000");
-			}
-			if ( vbegin == v2 ) return -1;// error 
-		}
-
-		tta.SetLock();
-		taret=tta;
-		a.Optim(1,0);
-		b.Optim(1,0);
-		return NbSwap; 
 	}
 	/*}}}*/
@@ -4897,90 +4928,3 @@
 	} // end swap 
 	/*}}}*/
-	int SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/
-		// l'arete ta coupe l'arete pva pvb
-		// de cas apres le swap sa coupe toujours
-		// on cherche l'arete suivante 
-		// on suppose que detsa >0 et detsb <0
-		// attention la routine echange pva et pvb 
-
-		if(tt1.Locked()) return 0; // frontiere croise 
-
-		AdjacentTriangle tt2 = Adj(tt1);
-		Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
-		short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
-		if ( a1<0 || a1>=3 ){
-			_error_("a1<0 || a1>=3");
-		}
-
-		BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
-		BamgVertex & s1= (*t1)[OppositeVertex[a1]];
-		BamgVertex & s2= (*t2)[OppositeVertex[a2]];
-
-		long long dets2 = det(*pva,*pvb,s2);
-		long long det1=t1->det , det2=t2->det ;
-		long long detT = det1+det2;
-		if ((det1<=0 ) || (det2<=0)){
-			_error_("(det1<=0 ) || (det2<=0)");
-		}
-		if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb
-			_error_("(detsa>=0) || (detsb<=0)");
-		}
-		long long ndet1 = bamg::det(s1,sa,s2);
-		long long ndet2 = detT - ndet1;
-
-		int ToSwap =0; //pas de swap
-		if ((ndet1 >0) && (ndet2 >0)) 
-		  { // on peut swaper  
-			if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
-			 ToSwap =1; 
-			else // swap alleatoire 
-			 if (BinaryRand()) 
-			  ToSwap =2; 
-		  }
-		if (ToSwap) NbSwap++,
-		 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
-
-		int ret=1;
-
-		if (dets2 < 0) {// haut
-			dets1 = ToSwap ? dets1 : detsa ;
-			detsa = dets2; 
-			tt1 =  Previous(tt2) ;}
-		else if (dets2 > 0){// bas 
-			dets1 = ToSwap ? dets1 : detsb ;
-			detsb = dets2;
-			//xxxx tt1 = ToSwap ? tt1 : Next(tt2);
-			if(!ToSwap) tt1 =  Next(tt2);
-		}
-		else { // changement de direction 
-			ret = -1;
-			Exchange(pva,pvb);
-			Exchange(detsa,detsb);
-			Exchange(dets1,dets2);
-			Exchange(tt1,tt2);
-			dets1=-dets1;
-			dets2=-dets2;
-			detsa=-detsa;
-			detsb=-detsb;
-
-			if(ToSwap){
-				if (dets2 < 0) {// haut
-					dets1 = (ToSwap ? dets1 : detsa) ;
-					detsa = dets2; 
-					tt1 =  Previous(tt2) ;}
-				else if(dets2 > 0){// bas 
-					dets1 = (ToSwap ? dets1 : detsb) ;
-					detsb =  dets2;
-					if(!ToSwap) tt1 =  Next(tt2);
-				}
-				else {// on a fin ???
-					tt1 = Next(tt2);
-					ret =0;}
-			}
-
-		}
-		return ret;
-	}
-	/*}}}*/
 }
