Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2964)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2965)
@@ -114,5 +114,5 @@
 			Direction DirOfSearch;
 			union {
-				Triangle * t; // one triangle which contained  the vertex
+				Triangle* t; // one triangle which contained  the vertex
 				Int4 color;
 				Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed 
@@ -366,10 +366,9 @@
 		}
 
-		void SetTriangleContainingTheVertex()
-		  { 
+		void SetTriangleContainingTheVertex() { 
 			if (ns[0]) (ns[0]->t=this,ns[0]->vint=0);
 			if (ns[1]) (ns[1]->t=this,ns[1]->vint=1);
 			if (ns[2]) (ns[2]->t=this,ns[2]->vint=2);
-		  }
+		}
 
 		int swap(Int2 a1,int=0);
@@ -782,5 +781,5 @@
 
 			Vertex * NearestVertex(Icoor1 i,Icoor1 j) ;
-			Triangle * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
+			Triangle * FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
 
 			void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2964)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2965)
@@ -234,5 +234,5 @@
 		Real8 Hmin = HUGE_VAL;// the infinie value 
 		Int4 hvertices=0;
-		int i,j,k,n;
+		int i,j,k,n,i1,i2;
 
 		//initialize some variables
@@ -290,25 +290,31 @@
 		//Edges
 		if (bamggeom->Edges){
-			int i1,i2;
-			R2 zero2(0,0);
-			Real4 *len =0;
+			R2     zero2(0,0);
+			Real4* len=NULL;
 
 			if(verbose>3) printf("      processing Edges\n");
 			edges = new GeometricalEdge[nbe];
 
+			//if hvertices==0, initialize len (length of each edge)
 			if (!hvertices) {
 				len = new Real4[nbv];
-				for(i=0;i<nbv;i++)
-				 len[i]=0;
+				for(i=0;i<nbv;i++) len[i]=0;
 			}
 
 			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].ref=(Int4)bamggeom->Edges[i*3+2];
-				edges[i].v[0]=  vertices + i1;
-				edges[i].v[1]=  vertices + i2;
-				R2 x12 = vertices[i2].r-vertices[i1].r;
+
+				edges[i].v[0]= vertices + i1;
+				edges[i].v[1]= vertices + i2;
+
+				//get length of edge
+				R2    x12=vertices[i2].r-vertices[i1].r;
 				Real8 l12=sqrt((x12,x12));
+				Hmin=Min(Hmin,l12);
+
+				//initialize other fields
 				edges[i].tg[0]=zero2;
 				edges[i].tg[1]=zero2;
@@ -322,6 +328,4 @@
 					len[i2] += l12;
 				}
-
-				Hmin = Min(Hmin,l12);
 			}
 
@@ -330,7 +334,7 @@
 				for (i=0;i<nbv;i++) 
 				 if (vertices[i].color > 0) 
-				  vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
+				  vertices[i].m=Metric(len[i] /(Real4) vertices[i].color);
 				 else 
-				  vertices[i].m=  Metric(Hmin);
+				  vertices[i].m=Metric(Hmin);
 				delete [] len;
 			}
@@ -460,7 +464,8 @@
 				j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
 				if (j>nbe-1 || j<0){
-					throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbv));
-				}
-				edges[j].SetRequired();  }
+					throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbe));
+				}
+				edges[j].SetRequired();  
+			}
 		}
 		else{
@@ -741,18 +746,26 @@
 		}
 
-		// generation of  all the tangente 
+		// generation of  all the tangents
 		for (i=0;i<nbe;i++) {
-			R2 AB = edges[i].v[1]->r -edges[i].v[0]->r;        
-			Real8 lAB = Norme2(AB); // length of current edge AB
+
+			//Get AB = vertex1->vertex2
+			R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;        
+			//Get length of AB
+			Real8 lAB=Norme2(AB);
+			//initialize tangent
 			Real8 ltg2[2];
 			ltg2[0]=0;ltg2[1]=0;
+
+			//loop over the 2 vertices of the edge
 			for (jj=0;jj<2;jj++) {
-				R2 tg =  edges[i].tg[jj];
-				Real8 ltg = Norme2(tg); // length of tg
-				if(ltg == 0) {// no tg
-					if( ! edges[i].v[jj]->Corner())   { // not a Corner       
-						tg =  edges[i].v[1-jj]->r 
-						  - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r;
-						ltg =  Norme2(tg);
+				R2    tg =edges[i].tg[jj];
+				Real8 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[jj]->Corner()){
+						tg =  edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r;
+						ltg=  Norme2(tg);
 						tg =  tg *(lAB/ltg),ltg=lAB;
 					}
Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 2964)
+++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 2965)
@@ -41,5 +41,5 @@
 		else {// not optimisation 
 			init();
-			t=tbegin = Bh.FindTriangleContening(a,deta);
+			t=tbegin = Bh.FindTriangleContaining(a,deta);
 			if( t->det>=0)
 			 ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det);
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2964)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2965)
@@ -730,5 +730,5 @@
 	/*Methods*/
 	/*FUNCTION Triangles::Add{{{1*/
-	void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3) {
+	void Triangles::Add( Vertex &s,Triangle* t, Icoor2* det3) {
 		// -------------------------------------------
 		//             s2
@@ -746,21 +746,33 @@
 		//-------------------------------------------- 
 
-		Triangle * tt[3]; // the 3 new Triangles
-		Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2];
-		Icoor2  det3local[3];
-		int infv = &s0 ?  ((  &s1 ? ( &s2  ? -1 : 2) : 1  )) : 0;
-		// infv = ordre of the infini vertex (null)
-		register int nbd0 =0; // number of zero det3
-		register int izerodet=-1,iedge; // izerodet = egde contening the vertex s
+		//the three triangles tt
+		Triangle* tt[3];
+		//three vertices of t
+		Vertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];
+		//three determinants
+		Icoor2 det3local[3];
+		// number of zero det3
+		register int nbzerodet =0; 
+		// izerodet = egde contening the vertex s
+		register int izerodet=-1,iedge; 
+		//determinant of t
 		Icoor2 detOld = t->det;
 
-		if (( infv <0 ) && (detOld <0) ||  ( infv >=0  ) && (detOld >0) ){
-			throw ErrorException(__FUNCT__,exprintf("infv=%g det=%g"));
-		}
-
-		// if det3 do not exist then constuct det3
+		// infinitevertexpos = order of the infinite vertex (NULL)
+		// if no infinite vertex (NULL) infinitevertexpos=-1
+		// else if v_i is infinite, infinitevertexpos=i
+		int infinitevertexpos = &s0 ?  ((  &s1 ? ( &s2  ? -1 : 2) : 1  )) : 0;
+
+		//some checks
+		if (( infinitevertexpos <0 ) && (detOld <0) ||  ( infinitevertexpos >=0  ) && (detOld >0) ){
+			throw ErrorException(__FUNCT__,exprintf("bug in Triangles::Add, bad configuration"));
+		}
+
+		// if det3 does not exist, build it 
 		if (!det3) { 
-			det3 = det3local; // alloc 
-			if ( infv<0 ) {
+			//allocate
+			det3 = det3local;
+			//if no infinite vertex
+			if (infinitevertexpos<0 ) {
 				det3[0]=bamg::det(s ,s1,s2);
 				det3[1]=bamg::det(s0,s ,s2);
@@ -768,84 +780,91 @@
 			else { 
 				// one of &s1  &s2  &s0 is NULL so (&si || &sj) <=> !&sk
-				det3[0]=  &s0 ? -1  : bamg::det(s ,s1,s2) ;
-				det3[1]=  &s1 ? -1 : bamg::det(s0,s ,s2) ;
-				det3[2]=  &s2 ? -1 : bamg::det(s0,s1,s ) ;}}
-
-
-				if (!det3[0]) izerodet=0,nbd0++;
-				if (!det3[1]) izerodet=1,nbd0++;
-				if (!det3[2]) izerodet=2,nbd0++;
-
-				if  (nbd0 >0 ) // point s on a egde or on a vertex 
-				 if (nbd0 == 1) {
-					 iedge = OppositeEdge[izerodet];
-					 TriangleAdjacent ta = t->Adj(iedge);
-
-					 // the point is on the edge 
-					 // if the point is one the boundary 
-					 // add the point in outside part 
-					 if ( t->det >=0) { // inside triangle
-						 if ((( Triangle *) ta)->det < 0 ) {
-							 // add in outside triangle 
-							 Add(s,( Triangle *) ta);
-							 return;}
-					 }}
-				 else {
-					 printf("bug (%i): Bug double points in\n",nbd0);
-					 throw ErrorException(__FUNCT__,exprintf("See above"));
-				 }
-
-				// remove de MarkUnSwap edge
-				t->SetUnMarkUnSwap(0);     
-				t->SetUnMarkUnSwap(1);     
-				t->SetUnMarkUnSwap(2);
-
-				tt[0]= t;
-				tt[1]= &triangles[nbt++];
-				tt[2]= &triangles[nbt++];
-
-				if (nbt>nbtx) {
-					throw ErrorException(__FUNCT__,exprintf("Not ebough triangles"));
-				}
-
-				*tt[1]=   *tt[2]= *t;
-				// gestion of the link
-				tt[0]->link=tt[1];
-				tt[1]->link=tt[2]; 
-
-				(* tt[0])(OppositeVertex[0])=&s;
-				(* tt[1])(OppositeVertex[1])=&s;
-				(* tt[2])(OppositeVertex[2])=&s;
-
-				tt[0]->det=det3[0];
-				tt[1]->det=det3[1];
-				tt[2]->det=det3[2];         
-
-				//  update adj des triangles externe 
-				tt[0]->SetAdjAdj(0);
-				tt[1]->SetAdjAdj(1);
-				tt[2]->SetAdjAdj(2);
-				//  update des adj des 3 triangle interne
-				const int i0 = 0;
-				const int i1= NextEdge[i0];
-				const int i2 = PreviousEdge[i0];
-
-				tt[i0]->SetAdj2(i2,tt[i2],i0);
-				tt[i1]->SetAdj2(i0,tt[i0],i1);
-				tt[i2]->SetAdj2(i1,tt[i1],i2);
-
-				tt[0]->SetTriangleContainingTheVertex();
-				tt[1]->SetTriangleContainingTheVertex();
-				tt[2]->SetTriangleContainingTheVertex();
-
-
-				// swap if the point s is on a edge
-				if(izerodet>=0) {
-					int rswap =tt[izerodet]->swap(iedge);
-
-					if (!rswap) {
-						throw ErrorException(__FUNCT__,exprintf("swap the point s is on a edge"));
+				det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ;
+				det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ;
+				det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;
+			}
+		}
+
+		if (!det3[0]) izerodet=0,nbzerodet++;
+		if (!det3[1]) izerodet=1,nbzerodet++;
+		if (!det3[2]) izerodet=2,nbzerodet++;
+
+		//if nbzerodet>0, point s is on an egde or on a vertex 
+		if  (nbzerodet >0 ){ 
+			if (nbzerodet == 1) {
+				iedge = OppositeEdge[izerodet];
+				TriangleAdjacent ta = t->Adj(iedge);
+
+				// the point is on the edge 
+				// if the point is one the boundary 
+				// add the point in outside part 
+				if ( t->det >=0) { // inside triangle
+					if ((( Triangle *) ta)->det < 0 ) {
+						// add in outside triangle 
+						Add(s,( Triangle *) ta);
+						return;
 					}
 				}
+			}
+			else {
+				t->Echo();
+				printf("\nproblem while trying to add:\n");
+				s.Echo();
+				throw ErrorException(__FUNCT__,exprintf("Bug in Triangles::Add points duplicated %i times",nbzerodet));
+			}
+		}
+
+		// remove de MarkUnSwap edge
+		t->SetUnMarkUnSwap(0);     
+		t->SetUnMarkUnSwap(1);     
+		t->SetUnMarkUnSwap(2);
+
+		tt[0]= t;
+		tt[1]= &triangles[nbt++];
+		tt[2]= &triangles[nbt++];
+
+		if (nbt>nbtx) {
+			throw ErrorException(__FUNCT__,exprintf("Not ebough triangles"));
+		}
+
+		*tt[1]=   *tt[2]= *t;
+		// gestion of the link
+		tt[0]->link=tt[1];
+		tt[1]->link=tt[2]; 
+
+		(* tt[0])(OppositeVertex[0])=&s;
+		(* tt[1])(OppositeVertex[1])=&s;
+		(* tt[2])(OppositeVertex[2])=&s;
+
+		tt[0]->det=det3[0];
+		tt[1]->det=det3[1];
+		tt[2]->det=det3[2];         
+
+		//  update adj des triangles externe 
+		tt[0]->SetAdjAdj(0);
+		tt[1]->SetAdjAdj(1);
+		tt[2]->SetAdjAdj(2);
+		//  update des adj des 3 triangle interne
+		const int i0 = 0;
+		const int i1= NextEdge[i0];
+		const int i2 = PreviousEdge[i0];
+
+		tt[i0]->SetAdj2(i2,tt[i2],i0);
+		tt[i1]->SetAdj2(i0,tt[i0],i1);
+		tt[i2]->SetAdj2(i1,tt[i1],i2);
+
+		tt[0]->SetTriangleContainingTheVertex();
+		tt[1]->SetTriangleContainingTheVertex();
+		tt[2]->SetTriangleContainingTheVertex();
+
+
+		// swap if the point s is on a edge
+		if(izerodet>=0) {
+			int rswap =tt[izerodet]->swap(iedge);
+
+			if (!rswap) {
+				throw ErrorException(__FUNCT__,exprintf("swap the point s is on a edge"));
+			}
+		}
 	}
 	/*}}}1*/
@@ -2011,19 +2030,23 @@
 		/*}}}1*/
 	/*FUNCTION Triangles::ForceBoundary{{{1*/
-	void Triangles::ForceBoundary() {
-		long int verbosity=2;
-		if (verbosity > 2) printf("   ForceBoundary  nb of edge: %i\n",nbe);
-		int k=0;
-		Int4  nbfe=0,nbswp=0,Nbswap=0;
-		for (Int4 t = 0; t < nbt; t++){
-			if (!triangles[t].det) k++;
-		}
-		if (k!=0) {
-			throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k));
-		}
-
+		void Triangles::ForceBoundary() {
+			long int verbosity=2;
+			int k=0;
+			int nbfe=0,nbswp=0,Nbswap=0;
+
+			//display
+			if (verbosity > 2) printf("   ForceBoundary  nb of edge: %i\n",nbe);
+
+			//check that there is no triangle with 0 determinant
+			for (Int4 t = 0; t < nbt; t++){
+				if (!triangles[t].det) k++;
+			}
+			if (k!=0) {
+				throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k));
+			}
+
+			//Force Edges
 			TriangleAdjacent ta(0,0);
-			for (Int4 i = 0; i < nbe; i++) 
-			  {
+			for (Int4 i = 0; i < nbe; i++){
 				nbswp =  ForceEdge(edges[i][0],edges[i][1],ta);
 
@@ -2033,17 +2056,17 @@
 				if ( nbswp < 0 && k < 5){
 					throw ErrorException(__FUNCT__,exprintf("Missing  Edge %i, v0=%i,v1=%i",i ,Number(edges[i][0]),Number(edges[i][1])));
-				  }
+				}
 				if ( nbswp >=0 && edges[i].onGeometry->Cracked())
 				 ta.SetCracked();
-			  }
-
+			}
 
 			if (k!=0) {
 				throw ErrorException(__FUNCT__,exprintf("There are %i lost edges, the boundary might be crossing",k));
 			}
-			for (Int4 j=0;j<nbv;j++)
-			 Nbswap +=  vertices[j].Optim(1,0);
+			for (Int4 j=0;j<nbv;j++){
+				Nbswap +=  vertices[j].Optim(1,0);
+			}
 			if (verbosity > 3) printf("      number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 
-	}
+		}
 	/*}}}1*/
 	/*FUNCTION Triangles::FillHoleInMesh{{{1*/
@@ -2209,5 +2232,5 @@
 				Vertex *vi  = ordre[icount];
 				Icoor2 dete[3];
-				Triangle *tcvi = FindTriangleContening(vi->i,dete);
+				Triangle *tcvi = FindTriangleContaining(vi->i,dete);
 				quadtree->Add(*vi); 
 				Add(*vi,tcvi,dete);
@@ -2584,6 +2607,6 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::FindTriangleContening{{{1*/
-	Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
+	/*FUNCTION Triangles::FindTriangleContaining{{{1*/
+	Triangle * Triangles::FindTriangleContaining(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
 		Triangle * t=0;	
 		int j,jp,jn,jj;
@@ -2599,5 +2622,5 @@
 					printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a));
 				}
-				throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContening"));
+				throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContaining"));
 			}
 			if (a<vertices || a>=vertices+nbv){
@@ -2682,6 +2705,8 @@
 		Gh.NbRef++;// add a ref to GH
 
-		Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve;
+		int i,j,k;
+		int NbOfCurves=0,NbNewPoints,NbEdgeCurve;
 		Real8 lcurve,lstep,s;
+		const int MaxSubEdge = 10;
 
 		R2 AB;
@@ -2710,5 +2735,5 @@
 		}
 
-		//compute nbv (number of vertices)
+		//Add all the geometrical vertices to the mesh
 		nbv=0;
 		for (i=0;i<Gh.nbv;i++){
@@ -2730,5 +2755,4 @@
 				Gh[i].to=vertices+nbv;
 
-
 				//Build VerticesOnGeomVertex for current point
 				VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);
@@ -2753,6 +2777,4 @@
 			nbe=0;
 			Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
-
-
 			Gh.UnMarkEdges();	
 			NbOfCurves = 0;
@@ -2761,4 +2783,5 @@
 			for (i=0;i<Gh.nbe;i++) {
 
+				//ei = current Geometrical edge
 				GeometricalEdge &ei=Gh.edges[i];   
 
@@ -2773,5 +2796,6 @@
 							Edge* PreviousNewEdge=NULL;
 
-							lstep = -1;//do not create points
+							//do not create points if required
+							lstep = -1;
 							if(ei.Required()){
 								if (j==0){
@@ -2781,6 +2805,8 @@
 										a=ei(0)->The();
 										b=ei(1)->The();
+
+										//check that edges has been allocated
 										if (!edges){
-											throw ErrorException(__FUNCT__,exprintf("!edges"));
+											throw ErrorException(__FUNCT__,exprintf("edges has not been allocated..."));
 										}
 										edges[nbe].v[0]=a->to;
@@ -2794,7 +2820,11 @@
 								}
 							}
-							else { // on curve ------
+
+							//if not required: on curve
+							else {
 								for ( int kstep=0;kstep<= step;kstep++){
-									// if 2nd step where 2 step
+									//step=0, nohing
+									//step=1,
+									//step=2
 									// 1 compute the length of the curve
 									// 2 create the points and edge
@@ -2807,31 +2837,48 @@
 									lcurve =0;
 									s = lstep;
-									int k=j;
-									e = & ei;
-									a=ei(k)->The();
-									va = a->to;
-									e->SetMark();
-
-									// if SameGeo  We have go in the background geometry 
-									// to find the discretisation of the curve
+
+									// i = edge number, j=[0;1] vertex number in edge
+									
+									k=j;            // k = vertex number in edge
+									e=&ei;          // e = address of current edge
+									a=ei(k)->The(); // a = pointer toward the kth vertex of the current edge
+									va = a->to;     // va= pointer toward vertex of ???
+									e->SetMark();   // Mark e
+
+									//if SameGeo We have go to the background geometry 
+									//to find the discretisation of the curve
 									for(;;){ 
 										k = 1-k;
-										b= (*e)(k)->The();
-										AB = b->r - a->r;
-										Metric MA = background ? BTh.MetricAt(a->r) :a->m ;
-										Metric MB =  background ? BTh.MetricAt(b->r) :b->m ;
-										Real8 ledge = (MA(AB) + MB(AB))/2;
-
-										const int MaxSubEdge = 10;
+										b = (*e)(k)->The(); // 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
+										Metric MB =  background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
+										Real8 ledge = (MA(AB) + MB(AB))/2;                   //Get edge length in metric
+
+										/* We are now creating the edges of the mesh from the
+										 * geometrical edge selected above.
+										 * The edge will be divided according to the metric
+										 * previously computed and cannot be divided more
+										 * than 10 times (MaxSubEdge). */
+
+										//By default, there is only one subedge that is the geometrical edge itself
 										int NbSubEdge = 1;
+
+										//initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
 										Real8 lSubEdge[MaxSubEdge];
-										R2 A,B;
+
+										//Build Subedges according to the edge length
+										//if ledge < 1.5 (between one and 2), take the edge as is
 										if (ledge < 1.5) lSubEdge[0] = ledge;
+										//else, divide the edge
 										else {
+											//compute number of subedges (division of the edge)
 											NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
-											A= a->r;
-											Metric MAs =MA,MBs;
-											ledge = 0; 
-											Real8 x =0, xstep= 1. /  NbSubEdge;
+											//A and B are the position of points on the edge
+											R2 A,B;
+											A=a->r;
+											Metric MAs=MA,MBs;
+											ledge=0; 
+											Real8 x =0, xstep= 1./NbSubEdge;
 											for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
 												x += xstep;
@@ -2844,48 +2891,48 @@
 
 										Real8 lcurveb = lcurve+ ledge ;
+
+										//Now, create corresponding points
+										//s = lstep
 										while (lcurve<=s && s <= lcurveb && nbv < nbvend){
-											// New points
-
-											// Real8 aa=(lcurveb-s)/ledge;
-											// Real8 bb=(s-lcurve)/ledge;
 
 											Real8 ss = s-lcurve;
-											// 1) find the SubEdge containing ss by dichotomie
+
+											/*Find the SubEdge containing ss using Dichotomy*/
+
 											int kk0=-1,kk1=NbSubEdge-1,kkk;
 											Real8 ll0=0,ll1=ledge,llk;
-											while (kk1-kk0>1)
-											  {
+											while (kk1-kk0>1){
 												if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
 												 kk1=kkk,ll1=llk;
 												else
-												 kk0=kkk,ll0=llk;}
-												if (kk1 == kk0){
-													throw ErrorException(__FUNCT__,exprintf("kk1 == kk0"));
-												}
-
-												Real8 sbb = (ss-ll0  )/(ll1-ll0);
-												Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb;
-
-												// new vertex on edge
-												vb = &vertices[nbv++];
-												vb->m = Metric(aa,a->m,bb,b->m);
-												vb->ReferenceNumber = e->ref;
-												vb->DirOfSearch =NoDirOfSearch;
-												Real8 abcisse = k ? bb : aa;
-												vb->r =  e->F( abcisse );
-												VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);        
-
-												// to take into account the direction of the edge
-												s += lstep;
-												edges[nbe].v[0]=va;
-												edges[nbe].v[1]=vb;
-												edges[nbe].ref = e->ref;
-												edges[nbe].onGeometry = e;
-												edges[nbe].adj[0] = PreviousNewEdge;
-												if(PreviousNewEdge)
-												 PreviousNewEdge->adj[1] = &edges[nbe];
-												PreviousNewEdge = edges + nbe;
-												nbe++;
-												va = vb;
+												 kk0=kkk,ll0=llk;
+											}
+											if (kk1 == kk0){
+												throw ErrorException(__FUNCT__,exprintf("kk1 == kk0"));
+											}
+
+											Real8 sbb = (ss-ll0  )/(ll1-ll0);
+											Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb;
+
+											// new vertex on edge
+											vb = &vertices[nbv++];
+											vb->m = Metric(aa,a->m,bb,b->m);
+											vb->ReferenceNumber = e->ref;
+											vb->DirOfSearch =NoDirOfSearch;
+											Real8 abcisse = k ? bb : aa;
+											vb->r =  e->F( abcisse );
+											VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);        
+
+											// to take into account the direction of the edge
+											s += lstep;
+											edges[nbe].v[0]=va;
+											edges[nbe].v[1]=vb;
+											edges[nbe].ref = e->ref;
+											edges[nbe].onGeometry = e;
+											edges[nbe].adj[0] = PreviousNewEdge;
+											if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
+											PreviousNewEdge=edges+nbe;
+											nbe++;
+											va = vb;
 										}
 										lcurve = lcurveb;
@@ -2947,6 +2994,11 @@
 		}
 
+		//Insert points inside existing triangles
 		Insert();
+
+		//Force the boundary
 		ForceBoundary();
+
+		//Extract SubDomains
 		FindSubDomain();
 
@@ -3369,19 +3421,58 @@
 	/*FUNCTION Triangles::Insert{{{1*/
 	void Triangles::Insert() {
+		/*Insert points in the existing Geometry*/
+
+		//Intermediary
+		int i;
+
+		/*Get options*/
 		long int verbosity=2;
+
+		//Display info
 		if (verbosity>2) printf("   Insert initial %i vertices\n",nbv);
+
+		//Compute integer coordinates and determinants for the existing vertices (from Geometry)
 		SetIntCoor();
-		Int4 i;
-
-		//fill out ordre with the adresses of each vertex
+
+		/*Now we want to build a list (ordre) of the vertices in a random
+		 * order. To do so, we use the following method:
+		 *
+		 * From an initial k0 in [0 nbv[ random (vertex number)
+		 * the next k (vertex number) is computed using a big
+		 * prime number (PN>>nbv) following:
+		 *
+		 * k_{i+1} = k_i + PN  [nbv]
+		 *
+		 * let's show that:
+		 *
+		 *   for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j
+		 *
+		 * Let's assume that there are 2 distinct j1 and j2 such that
+		 * k_j1 = k_j2
+		 *
+		 * This means that
+		 *  
+		 *  k0+j1*PN = k0+j2*PN [nbv]
+		 *  (j1-j2)*PN =0       [nbv]
+		 * since PN is a prime number larger than nbv, and nbv!=1
+		 *  j1-j2=0             [nbv]
+		 * BUT
+		 *  j1 and j2 are in [0 nbv[ which is impossible.
+		 *
+		 *  We hence have built a random list of nbv elements of
+		 *  [0 nbv[ all distincts*/
 		for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ;
-
-		//Build a random order
 		const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
-		Int4 k3 = rand()%nbv ; 
-		for (int is3=0; is3<nbv; is3++) 
-		 ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
-
+		int   k0=rand()%nbv; 
+		for (int is3=0; is3<nbv; is3++){
+			ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv];
+		}
+
+		/*Modify ordre such that the first 3 vertices form a triangle*/
+
+		//get first vertex i such that [0,1,i] are not aligned
 		for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
+			//if i is higher than nbv, it means that all the determinants are 0,
+			//all vertices are aligned!
 			if  ( ++i >= nbv) {
 				throw ErrorException(__FUNCT__,exprintf("all the vertices are aligned"));
@@ -3389,23 +3480,21 @@
 		}
 		// exchange i et 2 in "ordre" so that
-		// the first 3 vertices are not aligned
+		// the first 3 vertices are not aligned (real triangle)
 		Exchange(ordre[2], ordre[i]);
 
-		// We add a vertex at the infinity to build the mesh
-		// so that we have a simple definition of the boundary edges
+		/*Take the first edge formed by the first two vertices and build
+		 * the initial simple mesh from this edge and 2 boundary triangles*/
+
+		Vertex *  v0=ordre[0], *v1=ordre[1];
+
 		nbt = 2;
-
-		// We build a trivial mesh from one edge
-		// and 2 triangles built with the 2 oriented
-		// et
-		Vertex *  v0=ordre[0], *v1=ordre[1];
-
 		triangles[0](0) = 0; //infinite vertex
 		triangles[0](1) = v0;
 		triangles[0](2) = v1;
-
 		triangles[1](0) = 0;//infinite vertex
 		triangles[1](2) = v0;
 		triangles[1](1) = v1;
+
+		//Build adjacence
 		const int e0 = OppositeEdge[0];
 		const int e1 = NextEdge[e0];
@@ -3429,18 +3518,29 @@
 		quadtree->Add(*v1);
 
-		//the vertices are added one by one
+		/*Now, add the vertices One by One*/
+
 		Int4 NbSwap=0;
-
-
 		if (verbosity>3) printf("   Begining of insertion process...\n");
 
 		for (Int4 icount=2; icount<nbv; icount++) {
-			Vertex *vi=ordre[icount];
+
+			//Get new vertex
+			Vertex *newvertex=ordre[icount];
+
+			//Find the triangle in which newvertex is located
 			Icoor2 dete[3];
-			Triangle *tcvi = FindTriangleContening(vi->i,dete);
-			quadtree->Add(*vi); 
-			Add(*vi,tcvi,dete);
-			NbSwap += vi->Optim(1,0);
-		}
+			Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
+
+			//Add newvertex to the quadtree
+			quadtree->Add(*newvertex); 
+
+			//Add newvertex to the existing mesh
+			Add(*newvertex,tcvi,dete);
+
+			//Make the mesh Delaunay around newvertex by swaping the edges
+			NbSwap += newvertex->Optim(1,0);
+		}
+
+		//Display info
 		if (verbosity>3) {
 			printf("      NbSwap of insertion: %i\n",NbSwap);
@@ -3450,14 +3550,12 @@
 		}
 		NbUnSwap = 0;
-		// construction d'un ordre aleatoire 
-		//  const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ;
+
 #ifdef NBLOOPOPTIM
 
-		k3 = rand()%nbv ; 
+		k0 = rand()%nbv ; 
 		for (int is4=0; is4<nbv; is4++) 
-		 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
-
-		for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++) 
-		  {
+		 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];
+
+		for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){
 			Int4  NbSwap = 0;
 			for (int is1=0; is1<nbv; is1++) 
@@ -3471,5 +3569,5 @@
 			NbUnSwap = 0;
 			if(!NbSwap) break;
-		  }
+		}
 		ReMakeTriangleContainingTheVertex(); 
 		// because we break the TriangleContainingTheVertex
@@ -3528,5 +3626,5 @@
 				}
 				vj.ReferenceNumber=0; 
-				Triangle *tcvj=FindTriangleContening(vj.i,dete);
+				Triangle *tcvj=FindTriangleContaining(vj.i,dete);
 				if (tcvj && !tcvj->link){
 					tcvj->Echo();
@@ -3703,5 +3801,5 @@
 		I2 a = toI2(A);
 		Icoor2 deta[3];
-		Triangle * t =FindTriangleContening(a,deta);
+		Triangle * t =FindTriangleContaining(a,deta);
 		if (t->det <0) { // outside
 			double ba,bb;
@@ -4221,8 +4319,9 @@
 /*FUNCTION Triangles::SetIntCoor{{{1*/
 void Triangles::SetIntCoor(const char * strfrom) {
+	/*Set integer coordinate for existing vertices*/
+
+	//Get extrema coordinates of the existing vertices
 	pmin =  vertices[0].r;
 	pmax =  vertices[0].r;
-
-	// recherche des extrema des vertices pmin,pmax
 	Int4 i;
 	for (i=0;i<nbv;i++) {
@@ -4235,7 +4334,9 @@
 	pmin = pmin-DD;
 	pmax = pmax+DD; 
+
+	//Compute coefIcoor
 	coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
 	if (coefIcoor<=0){
-		throw ErrorException(__FUNCT__,exprintf("coefIcoor<=0"));
+		throw ErrorException(__FUNCT__,exprintf("coefIcoor should be positive, a problem in the geometry is likely"));
 	}
 
@@ -4251,6 +4352,7 @@
 		Vertex & v1 = triangles[i][1];
 		Vertex & v2 = triangles[i][2];
-		if ( &v0 && &v1 &&  &v2 ) // a good triangles;
-		  {
+
+		//If this is not a boundary triangle
+		if ( &v0 && &v1 &&  &v2 ){
 			triangles[i].det= det(v0,v1,v2);
 			if (triangles[i].det <=0 && Nberr++ <10){
@@ -4263,6 +4365,8 @@
 				 }
 			}
-		  }
-		else triangles[i].det= -1; // Null triangle; 
+		}
+
+		//else, set as -1
+		else triangles[i].det=-1;
 	}
 }
@@ -5043,5 +5147,5 @@
 			vi.ReferenceNumber=0; 
 			vi.DirOfSearch =NoDirOfSearch;
-			Triangle *tcvi = FindTriangleContening(vi.i,dete);
+			Triangle *tcvi = FindTriangleContaining(vi.i,dete);
 			if (tcvi && !tcvi->link) {
 				printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
@@ -5244,5 +5348,5 @@
 		I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); 
 		Icoor2 dete[3];
-		Triangle & tb = *Th.FindTriangleContening(I,dete);
+		Triangle & tb = *Th.FindTriangleContaining(I,dete);
 
 		if  (tb.link) 
Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2964)
+++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2965)
@@ -55,5 +55,5 @@
 		I2 IBTh  = BTh.toI2(PNew);
 
-		tstart=BTh.FindTriangleContening(IBTh,deta,tstart);  
+		tstart=BTh.FindTriangleContaining(IBTh,deta,tstart);  
 
 		if (tstart->det <0){ // outside