Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5338)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5339)
@@ -478,11 +478,11 @@
 		 * a vertex v?
 		 * To do so, we use the same chaining algorithm but there is a difficulty
-		 * coming from the fact that F we have a couple of vertices and not one 
-		 * vertices.
-		 * To overcome this difficulty, we use a global indices exactly like in 
+		 * 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
+		 * j=p%2      column number of p
 		 *
 		 * Algorithm:
@@ -506,7 +506,5 @@
 			double lv10=Norme2(v10);
 			//check that its length is not 0
-			if(lv10==0) {
-				ISSMERROR("Length of edge %i is 0",i);
-			}
+			if(lv10==0)ISSMERROR("Length of edge %i is 0",i);
 			//compute angle in [-Pi Pi]
 			eangle[i] = atan2(v10.y,v10.x);
@@ -521,22 +519,22 @@
 		//sort head_v by order of increasing edges angle
 		for (i=0;i<nbv;i++) {
-			int exch=1, ord=0;      
+			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;               // pointer toward head_v[vertex i]
-				long *po=p;                     // copy of pointer p
-				long  n=*p;                     // next value of edge holding i
+				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
+				// loop over the edges that contain the vertex i (till -1)
 				while (n >=0){
 					ord++;
-					register long  i1=n/2;       // i1 = floor (n/2)
-					register long  j1=n%2;       // j1 = 1 if n is odd
-					register long* pn=next_p+n;  // pointer to next_p[n]
-
-					//  n = next_p[n] = position in edge of next vertex i
+					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;                       
 
@@ -562,6 +560,6 @@
 			}
 
-			// angular test on current vertex to guess whether it is a corner (ord = number of edges horlding i)
-			if(ord == 2) { 
+			// 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];
@@ -574,5 +572,5 @@
 					vertices[i].SetCorner() ; 
 				}
-				// if the ReferenceNumber a changing then is     SetRequired();
+				// if the edge type/referencenumber a changing then is SetRequired();
 				if (edges[i1].type != edges[i2].type || edges[i1].Required()){
 					vertices[i].SetRequired();
@@ -586,19 +584,21 @@
 			}
 
-			// close the list around the vertex 
+			/*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];
-			// now the list around the vertex is circular
-		}
-
+		}
+
+		/*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]) {
-					ISSMERROR("Bug Adj edge");
-				}
+
+				if( edges[i1].v[j1] != edges[i].v[j]) ISSMERROR("Problem while processing edges: check the edge list");
+
 				edges[i1].Adj[j1] = edges + i;
 				edges[i1].DirAdj[j1] = j;
@@ -608,26 +608,26 @@
 		// generation of  all the tangents
 		for (i=0;i<nbe;i++) {
-			//Get AB = vertex1->vertex2
-			R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;        
-			//Get length of AB
-			double lAB=Norme2(AB);
-			//initialize tangent
-			double ltg2[2]={0.0};
+			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 (jj=0;jj<2;jj++) {
-				R2     tg =edges[i].tg[jj];
+			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(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; //previous point and next point on curve
+					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].DirAdj[j]]->r;
 						ltg= Norme2(tg);
 						tg = tg *(lAB/ltg);
 						ltg= lAB;
 					}
-					//else:  a Corner with no tangent => nothing to do    
+					//else:  a Corner no tangent => nothing to do    
 				}
 				else{
@@ -636,11 +636,9 @@
 				}
 
-				ltg2[jj] = ltg;
-
-				if ((tg,AB)<0){
-					tg = -tg;
-				}
-
-				edges[i].tg[jj] = tg;
+				ltg2[j] = ltg;
+
+				if ((tg,AB)<0) tg = -tg;
+
+				edges[i].tg[j]=tg;
 			}
 			if (ltg2[0]!=0) edges[i].SetTgA();
@@ -649,42 +647,46 @@
 
 		for (int step=0;step<2;step++){
+
 			for (i=0;i<nbe;i++) edges[i].SetUnMark();
+
 			nbcurves = 0;
 			long  nbgem=0;
-			for (int level=0;level < 2 && nbgem != nbe;level++)
-			 for (i=0;i<nbe;i++) {
-				 GeometricalEdge & ei = edges[i];   
-				 for(jj=0;jj<2;jj++){
-					 if (!ei.Mark() && (level || ei[jj].Required())) { 
-						 // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++) 
-						 int k0=jj,k1;
-						 GeometricalEdge *e = & ei;
-						 GeometricalVertex *a=(*e)(k0); // begin 
-						 if(curves){
-							 curves[nbcurves].be=e;
-							 curves[nbcurves].kb=k0;
-						 }
-						 int nee=0;
-						 for(;;) { 
-							 nee++;
-							 k1 = 1-k0; // next vertex of the edge 
-							 e->SetMark();
-							 nbgem++;
-							 e->CurveNumber=nbcurves;
-							 if(curves) {
-								 curves[nbcurves].ee=e;
-								 curves[nbcurves].ke=k1;
-							 }
-
-							 GeometricalVertex *b=(*e)(k1);
-							 if (a == b ||  b->Required() ) break;
-							 k0 = e->DirAdj[k1];//  vertex in next edge
-							 e = e->Adj[k1]; // next edge
-						 }
-						 nbcurves++;
-						 if(level) a->SetRequired();
-					 }
-				 }
-			 } 
+
+			for (int level=0;level < 2 && nbgem != nbe;level++){
+				for (i=0;i<nbe;i++){
+					GeometricalEdge & ei = edges[i];   
+					for(j=0;j<2;j++){
+						if (!ei.Mark() && (level || ei[j].Required())) { 
+							// warning ei.Mark() can be change in loop for(j=0;j<2;j++) 
+							int k0=j,k1;
+							GeometricalEdge *e = & ei;
+							GeometricalVertex *a=(*e)(k0); // begin 
+							if(curves){
+								curves[nbcurves].be=e;
+								curves[nbcurves].kb=k0;
+							}
+							int nee=0;
+							for(;;) { 
+								nee++;
+								k1 = 1-k0; // next vertex of the edge 
+								e->SetMark();
+								nbgem++;
+								e->CurveNumber=nbcurves;
+								if(curves) {
+									curves[nbcurves].ee=e;
+									curves[nbcurves].ke=k1;
+								}
+
+								GeometricalVertex *b=(*e)(k1);
+								if (a == b ||  b->Required() ) break;
+								k0 = e->DirAdj[k1];//  vertex in next edge
+								e = e->Adj[k1]; // next edge
+							}
+							nbcurves++;
+							if(level) a->SetRequired();
+						}
+					}
+				} 
+			}
 			ISSMASSERT(nbgem && nbe);
 			if(step==0){