Index: /issm/trunk/src/c/modules/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Bamgx/Bamgx.cpp	(revision 5123)
+++ /issm/trunk/src/c/modules/Bamgx/Bamgx.cpp	(revision 5124)
@@ -67,5 +67,5 @@
 
 		//Renumbering
-		Th.ReNumberingTheTriangleBySubDomain();
+		Th.TrianglesRenumberBySubDomain();
 
 		//Crack mesh if requested
@@ -177,5 +177,5 @@
 
 		//Renumber by subdomain
-		Th.ReNumberingTheTriangleBySubDomain();
+		Th.TrianglesRenumberBySubDomain();
 
 		//Smooth vertices
Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp	(revision 5123)
+++ /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp	(revision 5124)
@@ -61,5 +61,5 @@
 	if (verbose) printf("Reading mesh\n");
 	Mesh Th(index_data,x_data,y_data,nods_data,nels_data); 
-	Th.ReMakeTriangleContainingTheVertex();
+	Th.CreateSingleVertexToTriangleConnectivity();
 
 	//Loop over output nodes
@@ -85,5 +85,5 @@
 
 			//Find triangle holding r/I
-			Triangle &tb=*Th.FindTriangleContaining(I,dete);
+			Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
 
 			// internal point 
Index: /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp	(revision 5124)
@@ -48,5 +48,5 @@
 		I2 IBTh  = BTh.toI2(PNew);
 
-		tstart=BTh.FindTriangleContaining(IBTh,deta,tstart);  
+		tstart=BTh.TriangleFindFromCoord(IBTh,deta,tstart);  
 
 		if (tstart->det <0){ // outside
Index: /issm/trunk/src/c/objects/Bamg/Edge.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Edge.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/Edge.h	(revision 5124)
@@ -28,5 +28,5 @@
 
 			//Methods
-			void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
+			void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
 				if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
 				if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
Index: /issm/trunk/src/c/objects/Bamg/GeometricalVertex.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/GeometricalVertex.cpp	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/GeometricalVertex.cpp	(revision 5124)
@@ -16,5 +16,5 @@
 	/*FUNCTION GeometricalVertex::Corner {{{1*/
 	int  GeometricalVertex::Corner() const {
-		return cas & 4;
+		return type & 4;
 	}
 	/*}}}*/
@@ -22,33 +22,17 @@
 	int  GeometricalVertex::Required()const {
 		// a corner is required
-		return cas & 6;
+		return type & 6;
 	}
-	/*}}}*/
-	/*FUNCTION GeometricalVertex::IsThe {{{1*/
-	int  GeometricalVertex::IsThe() const {
-		return link == this;
-	}  
 	/*}}}*/
 	/*FUNCTION GeometricalVertex::SetCorner {{{1*/
 	void GeometricalVertex::SetCorner(){
-		cas |= 4;
+		type |= 4;
 	}
 	/*}}}*/
 	/*FUNCTION GeometricalVertex::SetRequired {{{1*/
 	void GeometricalVertex::SetRequired(){
-		cas |= 2;
+		type |= 2;
 	}
 	/*}}}*/
-	/*FUNCTION GeometricalVertex::Set {{{1*/
-	void GeometricalVertex::Set(){
-		cas=0;
-	}
-	/*}}}*/
-	/*FUNCTION GeometricalVertex::The {{{1*/
-	GeometricalVertex* GeometricalVertex::The(){
-		// return a unique vertex
-		ISSMASSERT(link);
-		return link;
-	}/*}}}*/
 
 } 
Index: /issm/trunk/src/c/objects/Bamg/GeometricalVertex.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/GeometricalVertex.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/GeometricalVertex.h	(revision 5124)
@@ -14,23 +14,15 @@
 			friend class Geometry;
 
-			int cas;
-			GeometricalVertex* link; //  link all the same GeometricalVertex circular (Crack) 
+			int type;
 
 			//Constructors
-			GeometricalVertex() :cas(0), link(this) {};
+			GeometricalVertex():type(0){};
 
 			//Methods
 			int  Corner() const;
 			int  Required()const;
-			int  IsThe() const;
 			void SetCorner();
 			void SetRequired();
-			void Set();
-			GeometricalVertex* The();
 
-			//Inline methods
-			inline void Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){
-				*this=rec;
-			}
 	};
 
Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5124)
@@ -33,6 +33,4 @@
 		curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
 		subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;
-		for (i=0;i<nbv;i++)
-		 vertices[i].Set(Gh.vertices[i],Gh,*this);
 		for (i=0;i<nbe;i++)
 		 edges[i].Set(Gh.edges[i],Gh,*this);
@@ -89,5 +87,5 @@
 				vertices[i].DirOfSearch=NoDirOfSearch;
 				vertices[i].color =0;
-				vertices[i].Set();
+				vertices[i].type=0;
 			}
 			/*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
@@ -424,6 +422,4 @@
 
 		k=0;
-		//link all vertices to themselves by default
-		for (i=0;i<nbv;i++) vertices[i].link = vertices+i;
 
 		//build quadtree for this geometry
@@ -443,18 +439,9 @@
 			//if there is a vertex found that is to close to vertices[i] -> error
 			if( v && Norme1(v->r - vertices[i]) < eps ){
-				// mama's old trick to get j 
-				GeometricalVertex* vg=(GeometricalVertex*) (void*)v;
-				j=vg-v0g;
-				//check that the clostest vertex is not itself...
-				if ( v !=  &(BamgVertex &) vertices[j]){
-					ISSMERROR(" v !=  &(BamgVertex &) vertices[j]");
-				}
-				vertices[i].link = vertices + j;
-				k++;	      
-			}
-
-			//The nearest vertex was non existent or far enough from vertices[i]
+				printf("WARNING: two points of the geometry are very closed to each other\n");
+			}
+
 			//Add vertices[i] to the quadtree
-			else  quadtree.Add(vertices[i]);
+			quadtree.Add(vertices[i]);
 		}
 
@@ -463,7 +450,4 @@
 			printf("number of distinct vertices= %i, over %i\n",nbv - k,nbv);
 			printf("List of duplicate vertices:\n");
-			for (i=0;i<nbv;i++){
-				if (!vertices[i].IsThe()) printf("  %i and %i\n",i,Number(vertices[i].The()));
-			}
 			ISSMERROR("See above");
 		}
Index: /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 5124)
@@ -34,5 +34,5 @@
 		else {// not optimisation 
 			init();
-			t=tbegin = Bh.FindTriangleContaining(a,deta);
+			t=tbegin = Bh.TriangleFindFromCoord(a,deta);
 			if( t->det>=0)
 			 ilast=NewItem(t,double(deta[0])/t->det,double(deta[1])/t->det,double(deta[2])/t->det);
Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5124)
@@ -1188,7 +1188,7 @@
 		tt[i2]->SetAdj2(i1,tt[i1],i2);
 
-		tt[0]->SetTriangleContainingTheVertex();
-		tt[1]->SetTriangleContainingTheVertex();
-		tt[2]->SetTriangleContainingTheVertex();
+		tt[0]->SetSingleVertexToTriangleConnectivity();
+		tt[1]->SetSingleVertexToTriangleConnectivity();
+		tt[2]->SetSingleVertexToTriangleConnectivity();
 
 
@@ -2256,5 +2256,5 @@
 
 		//Now, find the triangles that hold a cracked edge
-		ReMakeTriangleContainingTheVertex();
+		CreateSingleVertexToTriangleConnectivity();
 
 		long* Edgeflags=new long[NbCrackedEdges];
@@ -2537,5 +2537,5 @@
 				NbSubDomains =Gh.NbSubDomains;
 				long err=0;
-				ReMakeTriangleContainingTheVertex();
+				CreateSingleVertexToTriangleConnectivity();
 				long * mark = new long[nbt];
 				Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge();
@@ -2623,100 +2623,4 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Mesh::FindTriangleContaining{{{1*/
-	Triangle * Mesh::FindTriangleContaining(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/
-
-		Triangle * t=0;	
-		int j,jp,jn,jj;
-		int counter;
-
-		/*Get starting triangle. Take tsart if provided*/
-		if (tstart){
-			t=tstart;
-		}
-		/*Else find the closest Triangle using the quadtree*/
-		else {
-
-			/*Check that the quadtree does exist*/
-			if (!quadtree) ISSMERROR("no starting triangle provided and no quadtree available");
-
-			/*Call NearestVertex*/
-			BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
-
-			/*Check output (Vertex a)*/
-			if (!a)    ISSMERROR("problem while trying to find nearest vertex from a given point. No output found");
-			if (!a->t) ISSMERROR("no triangle is associated to vertex number %i (another call to ReMakeTriangleContainingTheVertex is required)",Number(a)+1);
-			ISSMASSERT(a>=vertices && a<vertices+nbv);
-
-			/*Get starting triangle*/
-			t = a->t;
-			ISSMASSERT(t>=triangles && t<triangles+nbt);
-		}
-
-		Icoor2  detop ;
-
-		/*initialize number of test triangle*/
-		counter=0; 
-
-		/*The initial triangle might be outside*/
-		while (t->det < 0){ 
-
-			/*Get a real vertex from this triangle (k0)*/
-			int k0=(*t)(0)?(((*t)(1)?((*t)(2)?-1:2):1)):0;
-			ISSMASSERT(k0>=0);// k0 the NULL vertex
-			int k1=NextVertex[k0],k2=PreviousVertex[k0];
-			dete[k0]=det(B,(*t)[k1],(*t)[k2]);
-			dete[k1]=dete[k2]=-1;     
-			if (dete[k0] > 0) // outside B 
-			 return t; 
-			t = t->TriangleAdj(OppositeEdge[k0]);
-			counter++;
-			ISSMASSERT(counter<2);
-		}
-
-		jj=0;
-		detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
-
-		while(t->det>0) { 
-
-			/*Increase counter*/
-			if (++counter>=10000) ISSMERROR("Maximum number of iteration reached (threshold = %i).",counter);
-
-			j= OppositeVertex[jj];
-			dete[j] = detop;  //det(*b,*s1,*s2);
-			jn = NextVertex[j];
-			jp = PreviousVertex[j];
-			dete[jp]= det(*(*t)(j),*(*t)(jn),B);
-			dete[jn] = t->det-dete[j] -dete[jp];
-
-			// count the number k of  dete <0
-			int k=0,ii[3];
-			if (dete[0] < 0 ) ii[k++]=0; 
-			if (dete[1] < 0 ) ii[k++]=1;
-			if (dete[2] < 0 ) ii[k++]=2;
-			// 0 => ok
-			// 1 => go in way 1
-			// 2 => two way go in way 1 or 2 randomly
-
-			if (k==0) break;
-			if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]);
-			if ( k>=3){
-				ISSMERROR("k>=3");
-			}
-			TriangleAdjacent t1 = t->Adj(jj=ii[0]);
-			if ((t1.det() < 0 ) && (k == 2))
-			 t1 = t->Adj(jj=ii[1]);
-			t=t1;
-			j=t1;// for optimisation we now the -det[OppositeVertex[j]];
-			detop = -dete[OppositeVertex[jj]];
-			jj = j;
-		}
-
-		if (t->det<0) // outside triangle 
-		 dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
-		//  NbOfTriangleSearchFind += counter;  
-		return t;
-	}
-	/*}}}1*/
 	/*FUNCTION Mesh::Init{{{1*/
 	void Mesh::Init(long maxnbv_in) {
@@ -2861,6 +2765,6 @@
 		triangles[1].det = -1;  //boundary triangle: det = -1
 
-		triangles[0].SetTriangleContainingTheVertex();
-		triangles[1].SetTriangleContainingTheVertex();
+		triangles[0].SetSingleVertexToTriangleConnectivity();
+		triangles[1].SetSingleVertexToTriangleConnectivity();
 
 		triangles[0].link=&triangles[1];
@@ -2884,5 +2788,5 @@
 			//Find the triangle in which newvertex is located
 			Icoor2 dete[3];
-			Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
+			Triangle* tcvi = TriangleFindFromCoord(newvertex->i,dete); //(newvertex->i = integer coordinates)
 
 			//Add newvertex to the quadtree
@@ -2918,5 +2822,5 @@
 			if(!NbSwap) break;
 		}
-		ReMakeTriangleContainingTheVertex(); 
+		CreateSingleVertexToTriangleConnectivity(); 
 		// because we break the TriangleContainingTheVertex
 #endif
@@ -2975,5 +2879,5 @@
 				}
 				vj.ReferenceNumber=0; 
-				Triangle *tcvj=FindTriangleContaining(vj.i,dete);
+				Triangle *tcvj=TriangleFindFromCoord(vj.i,dete);
 				if (tcvj && !tcvj->link){
 					tcvj->Echo();
@@ -3155,5 +3059,5 @@
 		I2 a = toI2(A);
 		Icoor2 deta[3];
-		Triangle * t =FindTriangleContaining(a,deta);
+		Triangle * t =TriangleFindFromCoord(a,deta);
 		if (t->det <0) { // outside
 			double ba,bb;
@@ -3202,8 +3106,8 @@
 				}
 			}
-			Bh.ReMakeTriangleContainingTheVertex();     
+			Bh.CreateSingleVertexToTriangleConnectivity();     
 			InsertNewPoints(nbvold,NbTSwap)   ;            
 		}  
-		else Bh.ReMakeTriangleContainingTheVertex();     
+		else Bh.CreateSingleVertexToTriangleConnectivity();     
 
 		// generation of the list of next Triangle 
@@ -3370,5 +3274,5 @@
 		double abscisse = -1;
 
-		for (int cas=0;cas<2;cas++){
+		for (int step=0;step<2;step++){
 			// 2 times algo:
 			//    1 for computing the length l
@@ -3385,19 +3289,13 @@
 
 				kkk=kkk+1;
-				if (kkk>=100){
-					ISSMERROR("kkk>=100");
-				}
-				if (!eee){
-					ISSMERROR("!eee");
-				}
+				ISSMASSERT(kkk<100);
+				ISSMASSERT(eee);
 				double lg0 = lg;
 				double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
 				lg += dp;
-				if (cas && abscisse <= lg) { // ok we find the geom edge 
+				if (step && abscisse <= lg) { // ok we find the geom edge 
 					double sss  =   (abscisse-lg0)/dp;
 					double thetab = te0*(1-sss)+ sss*iii;
-					if (thetab<0 || thetab>1){
-						ISSMERROR("thetab<0 || thetab>1");
-					}
+					ISSMASSERT(thetab>=0 && thetab<=1);
 					BR = VertexOnEdge(&R,eee,thetab);
 					return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
@@ -3410,7 +3308,5 @@
 
 				double lg0 = lg;
-				if (!eee){
-					ISSMERROR("!eee");
-				}
+				ISSMASSERT(eee);
 				v1 = pvB;
 				double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
@@ -3421,7 +3317,5 @@
 					double sss  =   (abscisse-lg0)/dp;
 					double thetab = te0*(1-sss)+ sss*tB;
-					if (thetab<0 || thetab>1){
-						ISSMERROR("thetab<0 || thetab>1");
-					}
+					ISSMASSERT(thetab>=0 && thetab<=1);
 					BR = VertexOnEdge(&R,eee,thetab);
 					return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
@@ -3594,6 +3488,6 @@
 	triangles[1].det = -1;  // boundary triangles
 
-	triangles[0].SetTriangleContainingTheVertex();
-	triangles[1].SetTriangleContainingTheVertex();
+	triangles[0].SetSingleVertexToTriangleConnectivity();
+	triangles[1].SetSingleVertexToTriangleConnectivity();
 
 	triangles[0].link=&triangles[1];
@@ -3610,5 +3504,5 @@
 		BamgVertex *vi  = ordre[icount];
 		Icoor2 dete[3];
-		Triangle *tcvi = FindTriangleContaining(vi->i,dete);
+		Triangle *tcvi = TriangleFindFromCoord(vi->i,dete);
 		quadtree->Add(*vi); 
 		AddVertex(*vi,tcvi,dete);
@@ -3742,6 +3636,6 @@
 }
 /*}}}1*/
-	/*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/
-	void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){
+	/*FUNCTION Mesh::TrianglesRenumberBySubDomain{{{1*/
+	void Mesh::TrianglesRenumberBySubDomain(bool justcompress){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
 
@@ -3786,5 +3680,5 @@
 		// do the change on all the pointeur 
 		for ( it=0;it<nbt;it++)
-		 triangles[it].ReNumbering(triangles,te,renu);
+		 triangles[it].Renumbering(triangles,te,renu);
 
 		for ( i=0;i<NbSubDomains;i++)
@@ -3812,6 +3706,6 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Mesh::ReNumberingVertex{{{1*/
-	void Mesh::ReNumberingVertex(long * renu) {
+	/*FUNCTION Mesh::VerticesRenumber{{{1*/
+	void Mesh::VerticesRenumber(long * renu) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
 
@@ -3824,9 +3718,9 @@
 		printf("renumbering triangles\n");
 		for ( it=0;it<nbt;it++) 
-		 triangles[it].ReNumbering(vertices,ve,renu);
+		 triangles[it].Renumbering(vertices,ve,renu);
 
 		printf("renumbering edges\n");
 		for ( ie=0;ie<nbe;ie++) 
-		 edges[ie].ReNumbering(vertices,ve,renu);
+		 edges[ie].Renumbering(vertices,ve,renu);
 
 		printf("renumbering vertices on geom\n");
@@ -4067,5 +3961,5 @@
 	if (quadtree) delete quadtree;
 	quadtree=0;
-	ReMakeTriangleContainingTheVertex();
+	CreateSingleVertexToTriangleConnectivity();
 	Triangle vide; // a triangle to mark the boundary vertex
 	Triangle   ** tstart= new Triangle* [nbv];
@@ -4109,5 +4003,5 @@
 	if(raisonmax<1.1) return;
 	if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
-	ReMakeTriangleContainingTheVertex();
+	CreateSingleVertexToTriangleConnectivity();
 	long i,j,kch,kk,ip;
 	long *first_np_or_next_t0 = new long[nbv];
@@ -4188,5 +4082,5 @@
 		const  int withBackground = &BTh != this && &BTh;
 
-		ReNumberingTheTriangleBySubDomain();
+		TrianglesRenumberBySubDomain();
 		//int nswap =0;
 		const long nfortria( choice ? 4 : 6);
@@ -4638,7 +4532,7 @@
 		SetIntCoor("In SplitElement"); 
 
-		ReMakeTriangleContainingTheVertex();
+		CreateSingleVertexToTriangleConnectivity();
 		if(withBackground)
-		 BTh.ReMakeTriangleContainingTheVertex();
+		 BTh.CreateSingleVertexToTriangleConnectivity();
 
 		delete [] edges;
@@ -4662,5 +4556,5 @@
 		NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
 		NbVertexOnBThEdge=newNbVertexOnBThEdge;
-		//  ReMakeTriangleContainingTheVertex();
+		//  CreateSingleVertexToTriangleConnectivity();
 
 		ReconstructExistingMesh();
@@ -4721,5 +4615,5 @@
 		  }
 	}
-	ReMakeTriangleContainingTheVertex();    
+	CreateSingleVertexToTriangleConnectivity();    
 	if (nbvold!=nbv){
 		long  iv = nbvold;
@@ -4734,5 +4628,5 @@
 			vi.ReferenceNumber=0; 
 			vi.DirOfSearch =NoDirOfSearch;
-			Triangle *tcvi = FindTriangleContaining(vi.i,dete);
+			Triangle *tcvi = TriangleFindFromCoord(vi.i,dete);
 			if (tcvi && !tcvi->link) {
 				printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
@@ -4758,4 +4652,150 @@
 
 return  NbSplitEdge;
+}
+/*}}}1*/
+/*FUNCTION Mesh::TriangleFindFromCoord{{{1*/
+Triangle * Mesh::TriangleFindFromCoord(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/
+
+	Triangle * t=0;	
+	int j,jp,jn,jj;
+	int counter;
+
+	/*Get starting triangle. Take tsart if provided*/
+	if (tstart) t=tstart;
+
+	/*Else find the closest Triangle using the quadtree*/
+	else {
+
+		/*Check that the quadtree does exist*/
+		if (!quadtree) ISSMERROR("no starting triangle provided and no quadtree available");
+
+		/*Call NearestVertex*/
+		BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
+
+		/*Check output (Vertex a)*/
+		if (!a)    ISSMERROR("problem while trying to find nearest vertex from a given point. No output found");
+		if (!a->t) ISSMERROR("no triangle is associated to vertex number %i (orphan?)",Number(a)+1);
+		ISSMASSERT(a>=vertices && a<vertices+nbv);
+
+		/*Get starting triangle*/
+		t = a->t;
+		ISSMASSERT(t>=triangles && t<triangles+nbt);
+	}
+
+	Icoor2  detop ;
+
+	/*initialize number of test triangle*/
+	counter=0; 
+
+	/*The initial triangle might be outside*/
+	while (t->det < 0){ 
+
+		/*Get a real vertex from this triangle (k0)*/
+		int k0=(*t)(0)?(((*t)(1)?((*t)(2)?-1:2):1)):0;
+		ISSMASSERT(k0>=0);// k0 the NULL vertex
+		int k1=NextVertex[k0],k2=PreviousVertex[k0];
+		dete[k0]=det(B,(*t)[k1],(*t)[k2]);
+		dete[k1]=dete[k2]=-1;     
+		if (dete[k0] > 0) // outside B 
+		 return t; 
+		t = t->TriangleAdj(OppositeEdge[k0]);
+		counter++;
+		ISSMASSERT(counter<2);
+	}
+
+	jj=0;
+	detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
+
+	while(t->det>0) { 
+
+		/*Increase counter*/
+		if (++counter>=10000) ISSMERROR("Maximum number of iteration reached (threshold = %i).",counter);
+
+		j= OppositeVertex[jj];
+		dete[j] = detop;  //det(*b,*s1,*s2);
+		jn = NextVertex[j];
+		jp = PreviousVertex[j];
+		dete[jp]= det(*(*t)(j),*(*t)(jn),B);
+		dete[jn] = t->det-dete[j] -dete[jp];
+
+		// count the number k of  dete <0
+		int k=0,ii[3];
+		if (dete[0] < 0 ) ii[k++]=0; 
+		if (dete[1] < 0 ) ii[k++]=1;
+		if (dete[2] < 0 ) ii[k++]=2;
+		// 0 => ok
+		// 1 => go in way 1
+		// 2 => two way go in way 1 or 2 randomly
+
+		if (k==0) break;
+		if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]);
+		if ( k>=3){
+			ISSMERROR("k>=3");
+		}
+		TriangleAdjacent t1 = t->Adj(jj=ii[0]);
+		if ((t1.det() < 0 ) && (k == 2))
+		 t1 = t->Adj(jj=ii[1]);
+		t=t1;
+		j=t1;// for optimisation we now the -det[OppositeVertex[j]];
+		detop = -dete[OppositeVertex[jj]];
+		jj = j;
+	}
+
+	if (t->det<0) // outside triangle 
+	 dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
+	//  NbOfTriangleSearchFind += counter;  
+	return t;
+}
+/*}}}1*/
+/*FUNCTION Mesh::TriangleIntNumbering{{{1*/
+void Mesh::TriangleIntNumbering(long* renumbering){
+
+	long num=0;
+	for (int i=0;i<nbt;i++){
+		if (triangles[i].det>0) renumbering[i]=num++;
+		else renumbering[i]=-1;
+	}
+	return;   
+}
+/*}}}1*/
+/*FUNCTION Mesh::TriangleReferenceList{{{1*/
+long  Mesh::TriangleReferenceList(long* reft) const {
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
+
+	long int verbose=0;
+	register Triangle *t0,*t;
+	register long k=0, num;   
+
+	//initialize all triangles as -1 (outside)
+	for (int it=0;it<nbt;it++) reft[it]=-1;
+
+	//loop over all subdomains
+	for (int i=0;i<NbSubDomains;i++){ 
+
+		//first triangle of the subdomain i
+		t=t0=subdomains[i].head;
+
+		//check that the subdomain is not empty
+		if (!t0){ ISSMERROR("At least one subdomain is empty");}
+
+		//loop
+		do{
+			k++;
+
+			//get current triangle number
+			num = Number(t);
+
+			//check that num is in [0 nbt[
+			if (num<0 || num>=nbt){ ISSMERROR("num<0 || num>=nbt");}
+
+			//reft of this triangle is the subdomain number
+			reft[num]=i;
+
+		} while (t0 != (t=t->link));
+		//stop when all triangles of subdomains have been tagged
+
+	}
+	return k;   
 }
 /*}}}1*/
@@ -4793,5 +4833,5 @@
 		//Compute number of vertices on geometrical vertex
 		for (i=0;i<Gh.nbv;i++){
-			if (Gh[i].Required() && Gh[i].IsThe()) NbVerticesOnGeomVertex++;
+			if (Gh[i].Required()) NbVerticesOnGeomVertex++;
 		}
 
@@ -4805,9 +4845,6 @@
 		nbv=0;
 		for (i=0;i<Gh.nbv;i++){
-			/* Add vertex only if required and not duplicate 
-			 * (IsThe is a method of GeometricalVertex that checks
-			 * that we are taking into acount only one vertex is
-			 * 2 vertices are superimposed) */
-			if (Gh[i].Required() && Gh[i].IsThe()) {//Gh  vertices Required
+			/* Add vertex only if required*/
+			if (Gh[i].Required()) {//Gh  vertices Required
 
 				//Add the vertex (provided that nbv<maxnbv)
@@ -4873,6 +4910,6 @@
 								else{ 
 									e=&ei;
-									a=ei(0)->The();
-									b=ei(1)->The();
+									a=ei(0);
+									b=ei(1);
 
 									//check that edges has been allocated
@@ -4906,5 +4943,5 @@
 								k=j;            // k = vertex number in edge (0 or 1)
 								e=&ei;          // e = reference of current edge
-								a=ei(k)->The(); // a = pointer toward the kth vertex of the current edge
+								a=ei(k);        // a = pointer toward the kth vertex of the current edge
 								va = a->to;     // va = pointer toward vertex associated
 								e->SetMark();   // Mark edge
@@ -4914,5 +4951,5 @@
 								for(;;){ 
 									k = 1-k;
-									b = (*e)(k)->The();// b = pointer toward the other vertex of the current edge
+									b = (*e)(k);// b = pointer toward the other vertex of the current edge
 									AB= b->r - a->r;   // AB = vector of the current edge
 									Metric MA = background ? BTh.MetricAt(a->r) :a->m ;  //Get metric associated to A
@@ -5403,55 +5440,4 @@
 	}
 	/*}}}1*/
-/*FUNCTION Mesh::TriangleReferenceList{{{1*/
-long  Mesh::TriangleReferenceList(long* reft) const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
-
-	long int verbose=0;
-	register Triangle *t0,*t;
-	register long k=0, num;   
-
-	//initialize all triangles as -1 (outside)
-	for (int it=0;it<nbt;it++) reft[it]=-1;
-
-	//loop over all subdomains
-	for (int i=0;i<NbSubDomains;i++){ 
-
-		//first triangle of the subdomain i
-		t=t0=subdomains[i].head;
-
-		//check that the subdomain is not empty
-		if (!t0){ ISSMERROR("At least one subdomain is empty");}
-
-		//loop
-		do{
-			k++;
-
-			//get current triangle number
-			num = Number(t);
-
-			//check that num is in [0 nbt[
-			if (num<0 || num>=nbt){ ISSMERROR("num<0 || num>=nbt");}
-
-			//reft of this triangle is the subdomain number
-			reft[num]=i;
-
-		} while (t0 != (t=t->link));
-		//stop when all triangles of subdomains have been tagged
-
-	}
-	return k;   
-}
-/*}}}1*/
-/*FUNCTION Mesh::TriangleIntNumbering{{{1*/
-void Mesh::TriangleIntNumbering(long* renumbering){
-
-	long num=0;
-	for (int i=0;i<nbt;i++){
-		if (triangles[i].det>0) renumbering[i]=num++;
-		else renumbering[i]=-1;
-	}
-	return;   
-}
-/*}}}1*/
 
 	/*Intermediary*/
@@ -5762,6 +5748,6 @@
 	t2->det = det2;
 
-	t1->SetTriangleContainingTheVertex();
-	t2->SetTriangleContainingTheVertex();
+	t1->SetSingleVertexToTriangleConnectivity();
+	t2->SetSingleVertexToTriangleConnectivity();
 } // end swap 
 /*}}}1*/
Index: /issm/trunk/src/c/objects/Bamg/Mesh.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.h	(revision 5124)
@@ -107,6 +107,6 @@
 			void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1);
 			long InsertNewPoints(long nbvold,long & NbTSwap) ; 
-			void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
-			void ReNumberingVertex(long * renu);
+			void TrianglesRenumberBySubDomain(bool justcompress=false);
+			void VerticesRenumber(long * renu);
 			void SmoothingVertex(int =3,double=0.3);
 			Metric MetricAt (const R2 &) const;
@@ -120,5 +120,5 @@
 			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;
+			Triangle* TriangleFindFromCoord(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);
@@ -135,7 +135,7 @@
 
 			//Inline methods
-			inline  void ReMakeTriangleContainingTheVertex(){
+			inline  void CreateSingleVertexToTriangleConnectivity(){
 				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();
+				for (int i=0;i<nbt;i++) triangles[i].SetSingleVertexToTriangleConnectivity();
 			}
 			inline  void  UnMarkUnSwapTriangle(){
Index: /issm/trunk/src/c/objects/Bamg/Metric.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Metric.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/Metric.h	(revision 5124)
@@ -13,5 +13,4 @@
 	class Metric;
 	class MatVVP2x2;
-
 
 	class Metric{
Index: /issm/trunk/src/c/objects/Bamg/R2.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/R2.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/R2.h	(revision 5124)
@@ -7,5 +7,5 @@
 namespace bamg {
 
-	template <class R,class RR> class P2 {
+	template <class R,class RR> class P2{
 
 		  public:  
@@ -37,5 +37,5 @@
 	  };
 
-	template <class R,class RR> class P2xP2 {
+	template <class R,class RR> class P2xP2{
 
 		  private:
Index: /issm/trunk/src/c/objects/Bamg/Triangle.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Triangle.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/Triangle.h	(revision 5124)
@@ -53,5 +53,5 @@
 			Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
 			Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ;
-			void  ReNumbering(Triangle *tb,Triangle *te, long *renu){
+			void  Renumbering(Triangle *tb,Triangle *te, long *renu){
 				if (link  >=tb && link  <te) link  = tb + renu[link -tb];
 				if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
@@ -59,5 +59,5 @@
 				if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];    
 			}
-			void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
+			void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
 				if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
 				if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
@@ -81,5 +81,5 @@
 				}
 			}
-			void SetTriangleContainingTheVertex() { 
+			void SetSingleVertexToTriangleConnectivity() { 
 				if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
 				if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);
Index: /issm/trunk/src/c/objects/Bamg/typedefs.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/typedefs.h	(revision 5123)
+++ /issm/trunk/src/c/objects/Bamg/typedefs.h	(revision 5124)
@@ -16,7 +16,5 @@
 	/*I2 and R2*/
 	typedef P2<Icoor1,Icoor2>  I2;
-	typedef P2xP2<short,long>  I2xI2;
 	typedef P2<double,double>  R2;
-	typedef P2<double,double>  R2xR2;
 }