Index: /issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp	(revision 5488)
+++ /issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp	(revision 5489)
@@ -31,30 +31,34 @@
 		ISSMASSERT(theta>-1e-12 && theta<1+1e-12);
 
-		if (TgA()) 
-		 if (TgB()) // interpolation d'hermite
-			{ cb =  theta*theta*(3-2*theta);
-			 ca =  1-cb;     
-			 cta = (1-theta)*(1-theta)*theta;
-			 ctb = (theta-1)*theta*theta ;
+		if (TgA()){ 
+			if (TgB()){ //Hermite interpolation
+				cb =  theta*theta*(3-2*theta);
+				ca =  1-cb;     
+				cta = (1-theta)*(1-theta)*theta;
+				ctb = (theta-1)*theta*theta ;
 			}
-		 else { // 1-t*t, t-t*t, t*t
-			 double t = theta;
-			 cb = t*t;
-			 ca = 1-cb;
-			 cta= t-cb;
-			 ctb=0;    
-		 }    
-		else
-		 if (TgB()){
-			 double t = 1-theta;
-			 ca = t*t;
-			 cb = 1-ca;
-			 ctb= -t+ca;
-			 cta=0;    
-		 }
-		 else {
-			 ca =(1-theta),cb = theta,cta=ctb=0; // lagrange P1
-		 }
-		return A*ca + B*cb + tg[0]* cta + tg[1] * ctb;
+			else {
+				double t = theta;
+				cb = t*t;
+				ca = 1-cb;
+				cta= t-cb;
+				ctb=0;    
+			}
+		}
+		else{
+			if (TgB()){
+				double t = 1-theta;
+				ca = t*t;
+				cb = 1-ca;
+				ctb= -t+ca;
+				cta=0;    
+			}
+			else { // lagrange P1
+				ca =(1-theta);
+				cb = theta;
+				cta=ctb=0;
+			}
+		}
+		return A*ca + B*cb + tg[0]*cta + tg[1]*ctb;
 	  }
 	/*}}}1*/
Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5488)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5489)
@@ -5047,19 +5047,18 @@
 											x += xstep;
 											B =  e->F(k ? x : 1-x);
-											MBs= background ? BTh.MetricAt(B) : Metric(1-x, MA, x ,MB);
+											MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
 											AB = A-B;
-											lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);
+											lSubEdge[kk]= (ledge+=(MAs(AB)+MBs(AB))/2);
 										}
 									}
 
-									double lcurveb = lcurve+ ledge ;
-
-									//Now, create corresponding points
-									while (lcurve<=s && s <= lcurveb && nbv < nbvend){
+									double lcurveb = lcurve+ledge;
+
+									/*Now, create corresponding points*/
+									while (lcurve<=s && s<=lcurveb && nbv<nbvend){
 
 										double ss = s-lcurve;
 
 										/*Find the SubEdge containing ss using Dichotomy*/
-
 										int kk0=-1,kk1=NbSubEdge-1,kkk;
 										double ll0=0,ll1=ledge,llk;
Index: /issm/trunk/src/c/objects/Bamg/Metric.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Metric.cpp	(revision 5488)
+++ /issm/trunk/src/c/objects/Bamg/Metric.cpp	(revision 5489)
@@ -51,16 +51,18 @@
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/EigenMetric)*/
 
+		/*Compute metric (linear combination of ma and mb)*/
 		Metric mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
+
+		/*Get Eigen values and vectors*/
 		EigenMetric vab(mab);
-
 		R2 v1(vab.v.x,vab.v.y);
 		R2 v2(-v1.y,v1.x);
 
-
-		double h1 = a / ma(v1) + b / mb(v1);
-		double h2 = a / ma(v2) + b / mb(v2);
-		vab.lambda1 =  1 / (h1*h1);
-		vab.lambda2 =  1 / (h2*h2);
-		*this = vab;
+		/*Modify eigen values (a+b=1)*/
+		double h1 = a/ma(v1) + b/mb(v1);
+		double h2 = a/ma(v2) + b/mb(v2);
+		vab.lambda1 =  1/(h1*h1);
+		vab.lambda2 =  1/(h2*h2);
+		*this=vab;
 	}
 	/*}}}1*/
Index: /issm/trunk/src/c/objects/Bamg/Metric.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Metric.h	(revision 5488)
+++ /issm/trunk/src/c/objects/Bamg/Metric.h	(revision 5489)
@@ -46,5 +46,5 @@
 			Metric operator/(double c) const {double c2=1/(c*c);return  Metric(a11*c2,a21*c2,a22*c2);} 
 			operator D2xD2(){ return D2xD2(a11,a21,a21,a22);}
-			double  operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
+			double  operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};        // length of x in metric sqrt(<Mx,x>)
 			double  operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};