Index: /issm/trunk/src/c/objects/Bamg/Curve.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Curve.cpp	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/Curve.cpp	(revision 5397)
@@ -11,5 +11,5 @@
 	/*Constructors/Destructors*/
 	/*FUNCTION Curve::Curve(){{{1*/
-	Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0),master(true){
+	Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0){
 
 	} 
Index: /issm/trunk/src/c/objects/Bamg/Curve.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Curve.h	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/Curve.h	(revision 5397)
@@ -16,5 +16,4 @@
 			int kb,ke;  //  begin vetex and end vertex
 			Curve* next; // next curve equi to this
-			bool master; // true => of equi curve point on this curve  
 
 			//Methods
Index: /issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp	(revision 5397)
@@ -152,25 +152,25 @@
 	/*FUNCTION GeometricalEdge::SetCracked{{{1*/
 	void   GeometricalEdge::SetCracked()     { 
-		type |= 1;
+		type |= 1;/*=>1st digit to 1*/
 	}/*}}}*/
 	/*FUNCTION GeometricalEdge::SetTgA{{{1*/
 	void   GeometricalEdge::SetTgA()         { 
-		type |=4; 
+		type |=4; /*=>2d digit to 1*/
 	}/*}}}*/
 	/*FUNCTION GeometricalEdge::SetTgB{{{1*/
 	void   GeometricalEdge::SetTgB()         { 
-		type |=8; 
+		type |=8; /*=> 3d digit to 1*/
 	}/*}}}*/
 	/*FUNCTION GeometricalEdge::SetMark{{{1*/
 	void   GeometricalEdge::SetMark()        { 
-		type |=16;
+		type |=16;/*=> 4th digiy to 1*/
 	}/*}}}*/
 	/*FUNCTION GeometricalEdge::SetUnMark{{{1*/
 	void   GeometricalEdge::SetUnMark()      { 
-		type &= 1007 /* 1023-16*/;
+		type &= 1007 /* 1023-16 = 000111110111 => 4th digit to 0*/;
 	}/*}}}*/
 	/*FUNCTION GeometricalEdge::SetRequired{{{1*/
 	void   GeometricalEdge::SetRequired()    { 
-		type |= 64; 
+		type |= 64;/*=>=>6th digit to 1*/ 
 	}/*}}}*/
 	  /*FUNCTION GeometricalEdge::Tg{{{1*/
Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5397)
@@ -22,5 +22,5 @@
 		Init();
 		ReadGeometry(bamggeom,bamgopts);
-		AfterRead();
+		PostRead();
 	}
 	/*}}}*/
@@ -165,5 +165,5 @@
 			}
 
-			// definition  the default of the given mesh size 
+			// definition the default of the given mesh size 
 			for (i=0;i<nbv;i++) {
 				if (vertices[i].color > 0) 
@@ -400,6 +400,6 @@
 
 	/*Methods*/
-	/*FUNCTION Geometry::AfterRead{{{1*/
-	void Geometry::AfterRead(){
+	/*FUNCTION Geometry::PostRead{{{1*/
+	void Geometry::PostRead(){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
 
@@ -606,5 +606,5 @@
 		}
 
-		// generation of  all the tangents
+		/* 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
@@ -646,19 +646,25 @@
 		} 
 
+		/* 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;
-			long  nbgem=0;
-
-			for (int level=0;level < 2 && nbgem != nbe;level++){
+
+			for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
 				for (i=0;i<nbe;i++){
-					GeometricalEdge & ei = edges[i];   
+
+					GeometricalEdge & 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())) { 
-							// warning ei.Mark() can be change in loop for(j=0;j<2;j++) 
 							int k0=j,k1;
-							GeometricalEdge *e = & ei;
+							GeometricalEdge   *e=&ei;
 							GeometricalVertex *a=(*e)(k0); // begin 
 							if(curves){
@@ -667,11 +673,11 @@
 							}
 							int nee=0;
-							for(;;) { 
+							for(;;){ 
 								nee++;
 								k1 = 1-k0; // next vertex of the edge 
 								e->SetMark();
-								nbgem++;
+								nb_marked_edges++;
 								e->CurveNumber=nbcurves;
-								if(curves) {
+								if(curves){
 									curves[nbcurves].ee=e;
 									curves[nbcurves].ke=k1;
@@ -679,7 +685,14 @@
 
 								GeometricalVertex *b=(*e)(k1);
-								if (a == b ||  b->Required() ) break;
-								k0 = e->AdjVertexNumber[k1];//  vertex in next edge
-								e = e->Adj[k1]; // next edge
+
+								//break if we have reached the other end of the curve
+								if (a==b || b->Required()){
+									break;
+								}
+								//else: go to next edge (adjacent)
+								else{
+									k0 = e->AdjVertexNumber[k1];//  vertex in next edge
+									e  = e->Adj[k1]; // next edge
+								}
 							}
 							nbcurves++;
@@ -689,13 +702,11 @@
 				} 
 			}
-			ISSMASSERT(nbgem && nbe);
-			if(step==0){
-				curves = new Curve[nbcurves];
-			}
+			ISSMASSERT(nb_marked_edges && nbe);
+			//allocate if first step
+			if(step==0) curves=new Curve[nbcurves];
 		} 
 		for(int i=0;i<nbcurves ;i++){
 			GeometricalEdge * be=curves[i].be, *eqbe=be;
 			//GeometricalEdge * ee=curves[i].ee, *eqee=be;
-			curves[i].master=true;
 		}
 
@@ -707,32 +718,5 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Geometry::Containing{{{1*/
-	GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
-
-		GeometricalEdge* on =start,* pon=0;
-		// walk with the cos on geometry
-		int counter=0;
-		while(pon != on){  
-			counter++;
-			ISSMASSERT(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;
-	}
-	/*}}}1*/
 	/*FUNCTION Geometry::Echo {{{1*/
-
 	void Geometry::Echo(void){
 
@@ -774,4 +758,7 @@
 	/*FUNCTION Geometry::MinimalHmin{{{1*/
 	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;
 	}/*}}}*/
@@ -800,4 +787,30 @@
 		return c - curves;
 	}/*}}}*/
+	/*FUNCTION Geometry::Containing{{{1*/
+	GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
+
+		GeometricalEdge* on =start,* pon=0;
+		// walk with the cos on geometry
+		int counter=0;
+		while(pon != on){  
+			counter++;
+			ISSMASSERT(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;
+	}
+	/*}}}1*/
 	/*FUNCTION Geometry::ProjectOnCurve {{{1*/
 	GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
Index: /issm/trunk/src/c/objects/Bamg/Geometry.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.h	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.h	(revision 5397)
@@ -53,5 +53,5 @@
 			void             ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
 			void             Init(void);
-			void             AfterRead();
+			void             PostRead();
 			long             GetId(const GeometricalVertex &t) const;
 			long             GetId(const GeometricalVertex *t) const;
Index: /issm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp	(revision 5397)
@@ -11,32 +11,80 @@
 	/*FUNCTION MatVVP2x2::MatVVP2x2(const Metric M){{{1*/
 	MatVVP2x2::MatVVP2x2(const Metric M){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MatVVP2x2)*/
+		/*From a metric (a11,a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
 
+		/*Intermediaries*/
 		double a11=M.a11,a21=M.a21,a22=M.a22;
-		const double eps = 1.e-5;
-		double c11 = a11*a11, c22 = a22*a22, c21= a21*a21;
-		double b=-a11-a22,c=-c21+a11*a22;
-		double delta = b*b - 4 * c ;
-		double n2=(c11+c22+c21);
+		double norm1,norm2,normM;
+		double delta,b;
 
-		//if norm(M)<10^30 -> M=zeros(2,2)
-		if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0;
+		/*To get the eigen values, we must solve the following equation:
+		 *     | a11 - lambda    a21        |
+		 * det |                            | = 0
+		 *     | a21             a22-lambda |
+		 *
+		 * We have to solve the following polynom:
+		 *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
 
-		//if matrix is already diagonal and has a 0 eigen value
-		else if (delta < eps*n2){ 
-			lambda1=lambda2=-b/2, v.x=1,v.y=0;
+		/*Compute polynom determinant*/
+		b=-a11-a22;
+		delta=b*b - 4*(a11*a22-a21*a21);
+
+
+		/*Compute norm of M to avoid round off errors*/
+		normM=a11*a11 + a22*a22 + a21*a21;
+
+		/*1: normM too small: eigen values = 0*/
+		if(normM<1.e-30){
+			lambda1=0;
+			lambda2=0;
+			v.x=1;
+			v.y=0;
+		}
+		/*2: delta is small -> double root*/
+		else if (delta < 1.e-5*normM){
+			lambda1=-b/2;
+			lambda2=-b/2;
+			v.x=1;
+			v.y=0;
+		}
+		/*3: general case -> two roots*/
+		else{
+			delta     = sqrt(delta);
+			lambda1   = (-b-delta)/2.0;
+			lambda2   = (-b+delta)/2.0;
+
+			/*Now, one must find the eigen vectors. For that we use the following property of the inner product
+			 *    <Ax,y> = <x,tAy>
+			 * Here, M'(M-lambda*Id) is symmetrical, which gives:
+			 *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
+			 * And we have the following:
+			 *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
+			 * We have shown that
+			 *    Im(M') ⊥ Ker(M')
+			 *
+			 * To find the eigen vectors of M, we only have to find two vectors
+			 * of the image of M' and take their perpendicular as long as they are
+			 * not 0.
+			 * To do that, we take the images (1,0) and (0,1):
+			 *  x1 = (a11 - lambda)      x2 = a21
+			 *  y1 = a21                 y2 = (a22-lambda)
+			 *
+			 * We take the vector that has the larger norm and take its perpendicular.*/
+
+			double norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 
+			double norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
+
+			if (norm2<norm1){
+				norm1=sqrt(norm1);
+				v.x = - a21/norm1;
+				v.y = (a11-lambda1)/norm1;
+			}
+			else{
+				norm2=sqrt(norm2);
+				v.x = - (a22-lambda1)/norm2;
+				v.y = a21/norm2;
+			}
 		}
 
-		//general case: construction of 2 eigen vectors
-		else {
-			delta     = sqrt(delta);
-			lambda1   = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;
-			double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;
-			double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-			if(s1 < s0)
-			 s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;
-			else
-			 s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;
-		};
 	}
 	/*}}}1*/
Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5397)
@@ -21,5 +21,5 @@
 		if(bamggeom->Edges) {
 			Gh.ReadGeometry(bamggeom,bamgopts);
-			Gh.AfterRead();
+			Gh.PostRead();
 		}
 
@@ -32,5 +32,5 @@
 			printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
 			BuildGeometryFromMesh(bamgopts);
-			Gh.AfterRead();
+			Gh.PostRead();
 		}
 
@@ -141,5 +141,5 @@
 		  //double cutoffradian = 10.0/180.0*Pi;
 		  BuildGeometryFromMesh(bamgopts);
-		  Gh.AfterRead(); 
+		  Gh.PostRead(); 
 		  SetIntCoor();
 		  ReconstructExistingMesh();
@@ -300,5 +300,5 @@
 		BuildGeometryFromMesh();
 		if (verbose) printf("Completing geometry\n");
-		Gh.AfterRead();
+		Gh.PostRead();
 	}
 	/*}}}1*/
@@ -5319,7 +5319,4 @@
 				int jedge=bcurve[icurve]%2;
 
-				/*Skip if we are on a equi curve (duplicate)*/
-				if(!Gh.curves[icurve].master) continue; 
-
 				/*Get edge of Bth with index iedge*/
 				Edge &ei = BTh.edges[iedge];
Index: /issm/trunk/src/c/objects/Bamg/Metric.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Metric.cpp	(revision 5396)
+++ /issm/trunk/src/c/objects/Bamg/Metric.cpp	(revision 5397)
@@ -269,5 +269,5 @@
 			for(i=0;i<2;i++){
 				/*Now, one must find the eigen vectors. For that we use the 
-				 * following property of the scalare product
+				 * following property of the inner product
 				 *    (Ax,b) = transp(b) Ax = transp(x) transp(A) b
 				 *           = (transp(A) b ,x)