Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2927)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2928)
@@ -96,6 +96,6 @@
 			MetricAnIso M=Gh[i];
 			MatVVP2x2 Vp(M/coef);
-			Vp.Maxh(hmin);
-			Vp.Minh(hmax);
+			Vp.Maxh(hmax);
+			Vp.Minh(hmin);
 			Gh.vertices[i].m = Vp;
 		}
Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2927)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2928)
@@ -63,5 +63,6 @@
 		MetricAnIso( Real8  a,const  MetricAnIso ma,
 					Real8  b,const  MetricAnIso mb);
-		int  IntersectWith(const MetricAnIso M);
+		int  IntersectWith(const MetricAnIso M2);
+		int  IntersectWith_new(const MetricAnIso M2);
 		MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return  MetricAnIso(a11*c2,a21*c2,a22*c2);} 
 		MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return  MetricAnIso(a11*c2,a21*c2,a22*c2);} 
@@ -122,6 +123,7 @@
 	}
 
-	void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
-	MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
+	void SimultaneousMatrixReduction(const MetricAnIso M1,const MetricAnIso M2, MetricAnIso* M);
+	void  SimultaneousMatrixReduction_bamg( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V);
+	void eigen2(double* M,double lambda[2],double vp[2][2]);;
 
 	inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) {
@@ -152,5 +154,5 @@
 
 	class SaveMetricInterpole {
-		friend  Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
+		friend Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
 		friend Real8 abscisseInterpole(const MetricAnIso ,const  MetricAnIso , R2 ,Real8 ,int );
 		int opt;
Index: /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2927)
+++ /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2928)
@@ -30,5 +30,5 @@
 		if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0;
 
-		//if ???
+		//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;
Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2927)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2928)
@@ -69,4 +69,34 @@
 	}
 	/*}}}*/
+	/*FUNCTION MetricAnIso::IntersectWith_new{{{1*/
+	int MetricAnIso::IntersectWith_new(MetricAnIso M2) {
+		/*Get a new metric from an existing metric (M1=this)
+		 * and a new metric given in input M2 using a 
+		 * Simultaneous Matrix Reduction
+		 * */
+
+		MetricAnIso M;
+
+		//Compute M intersection of M1 and M2
+		SimultaneousMatrixReduction(*this,M2,&M);
+
+		//check wether something has been done
+		if(M.a11!=a11 || M.a21!=a21 || M.a22!=a22){
+
+			//update this with M
+			a11=M.a11;
+			a21=M.a21;
+			a22=M.a22;
+
+			//return 1 -> something has been done
+			return 1;
+		}
+		else{
+
+			//nothing has changed, return 0
+			return 0;
+		}
+	}
+	/*}}}1*/
 	/*FUNCTION MetricAnIso::IntersectWith{{{1*/
 	int MetricAnIso::IntersectWith(const MetricAnIso M2) {
@@ -81,5 +111,5 @@
 		double l1,l2;
 
-		SimultaneousMatrixReduction(*this,M2,l1,l2,M);
+		SimultaneousMatrixReduction_bamg(*this,M2,l1,l2,M);
 
 		R2 v1(M.x.x,M.y.x);
@@ -103,66 +133,69 @@
 
 	/*Intermediary*/
-	/*FUNCTION SimultaneousMatrixReduction{{{1*/
-	void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
-		double a11=M1.a11,a21=M1.a21,a22=M1.a22;
-		double b11=M2.a11,b21=M2.a21,b22=M2.a22;
-		//  M1 v = l M2 v
-		// (M1 - l M2) v =0
-		// det (M1 - l M2) =0
-		// det (M1 - l M2) = a l^2 + b l + c;
-		// = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
-		const double /*c11 = a11*a11,*/ c21= a21*a21;
-		const double /*d11 = b11*b11,*/ d21= b21*b21;
-		const double a=b11*b22 - d21;
-		const double b=-a11*b22-a22*b11+2*a21*b21;
-		const double c=-c21+a11*a22;
-		const double bb = b*b,ac= a*c;
-		const double delta = bb - 4 * ac;
-		if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
-			// racine double;
-			if (Abs(a) < 1.e-30 )
-			 l1 = l2 = 0;
-			else 
-			 l1=l2=-b/(2*a); 
-			V= D2xD2(1,0,0,1);
-		  }
-		else {
-			const double delta2 = sqrt(delta);
-			l1= (-b - delta2)/(2*a);
-			l2= (-b + delta2)/(2*a);
-			// M1 v = l M2 v
-			//  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
-			// so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
-			double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
-			double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-			double vp1x,vp1y,vp2x,vp2y;
-
-			if(s1 < s0)
-			 s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
-			else
-			 s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
-
-			v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
-			s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-			if(s1 < s0)
-			 s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
-			else
-			 s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
-			V=D2xD2(vp1x,vp2x,vp1y,vp2y);
-		}
-		return;
-	}
-	/*}}}1*/
-	/*FUNCTION Intersection{{{1*/ 
-	MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) {
-		D2xD2 M;
-		double l1,l2;
-		SimultaneousMatrixReduction(M1,M2,l1,l2,M);
-		R2 v0(M.x.x,M.y.x);
-		R2 v1(M.x.y,M.y.y);
-		D2xD2 M_1(M.inv());
-		D2xD2 D(Max(M1(v0,v0),M2(v0,v0)),0,0,Max(M1(v1,v1),M2(v1,v1)));
-		D2xD2 Mi(M_1.t()*D*M_1);
-		return MetricAnIso(Mi.x.x,0.5*(Mi.x.y+Mi.y.x),Mi.y.y);
+	/*FUNCTION eigen2{{{1*/
+	void eigen2(double* M,double* lambda,double vp[2][2]) {
+		/*m=[a b]
+		 *  [c d] */
+
+		double normM,detM,trM;
+		double norm;
+		const double a=M[2*0+0];
+		const double b=M[2*0+1];
+		const double c=M[2*1+0];
+		const double d=M[2*1+1];
+
+		//compute norm(M)
+		normM=pow(pow(a,2.0)+pow(b,2.0)+pow(c,2.0)+pow(d,2.0),0.5);
+
+		//if normM<10-30, M=zeros(2,2)
+		if (normM<1.0e-30){
+			lambda[0]=0.0;
+			lambda[1]=0.0;
+			vp[0][0]=vp[0][1]=vp[1][0]=vp[1][1]=0.0;
+			return;
+		}
+
+		//check that matrix is already diagonal if det(M)=0
+		if(Abs(b)<(1.0e-12*normM) & Abs(c)<(1.0e-12*normM)){
+			lambda[0]=a;
+			lambda[1]=d;
+			vp[0][0]=1;
+			vp[0][1]=0;
+			vp[1][0]=0;
+			vp[1][1]=1;
+			return;
+		}
+
+		//compute det(M)
+		detM=a*d-c*b;
+
+		//not symetrical and still det<10^-30
+		if (Abs(detM)<1.0e-30){
+			throw ErrorException(__FUNCT__,exprintf("Cannot find eigen values of a non invertible matrix"));
+		}
+
+		//See http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html
+		trM=a+d;
+
+		//general case
+		lambda[0] = 0.5*(trM + pow( trM*trM - 4.0*detM ,0.5));
+		lambda[1] = 0.5*(trM - pow( trM*trM - 4.0*detM ,0.5));
+
+		if(Abs(b)>1.0e-30){ //if b!=0
+			norm=pow(b*b+pow(lambda[0]-a,2.0),0.5);
+			vp[0][0]=b/norm;
+			vp[1][0]=(lambda[0]-a)/norm;
+			norm=pow(b*b+pow(lambda[1]-a,2.0),0.5);
+			vp[0][1]=b/norm;
+			vp[1][1]=(lambda[1]-a)/norm;
+		}
+		else{
+			norm=pow(c*c+pow(lambda[0]-d,2.0),0.5);
+			vp[0][0]=(lambda[0]-d)/norm;
+			vp[1][0]=c/norm;
+			norm=pow(c*c+pow(lambda[1]-d,2.0),0.5);
+			vp[0][1]=(lambda[1]-d)/norm;
+			vp[1][1]=c/norm;
+		}
 	}
 	/*}}}1*/
@@ -233,4 +266,127 @@
 	}
 	/*}}}1*/
+	/*FUNCTION SimultaneousMatrixReduction{{{1*/
+	void SimultaneousMatrixReduction(MetricAnIso M1,MetricAnIso M2,MetricAnIso* pM) {
+
+		/*intermediary*/
+		double a11,a21,a22;
+		double b11,b21,b22;
+		double mu1,mu2;
+		double lambda1,lambda2;
+		double detM1,detP;
+		double N[2][2];
+		double lambda[2];
+		double vp[2][2];
+		double invP[2][2];
+		MetricAnIso M;
+
+		//get components of both metrics
+		a11=M1.a11; a21=M1.a21; a22=M1.a22;
+		b11=M2.a11; b21=M2.a21; b22=M2.a22;
+
+		//check diagnonal matrices (easy)
+		if (Abs(a21)< 1.0e-20 && Abs(b21)< 1.0e-20){
+			M.a11=Max(a11,b11);
+			M.a22=Max(a22,b22);
+			M.a21=0;
+			*pM=M;
+			return;
+		}
+
+		//Build N=M1^-1 M2 (Alauzet2003 p16)
+		detM1=a11*a22-a21*a21;
+		if (!detM1){
+			printf("a11 a21 a22 = %g %g %g\n",a11,a21,a22);
+			throw ErrorException(__FUNCT__,exprintf("One of the metric has a zero eigen value or cannot be inverted... (detP=%g)",detP));
+		}
+		N[0][0]= 1.0/detM1 * ( a22*b11-a21*b21);
+		N[0][1]= 1.0/detM1 * ( a22*b21-a21*b22);
+		N[1][0]= 1.0/detM1 * (-a21*b11+a11*b21);
+		N[1][1]= 1.0/detM1 * (-a21*b21+a11*b22);
+
+		//now, we must find the eigen values and vectors of N
+		eigen2(&N[0][0],lambda,vp);
+
+		//Now that we have P=vp, we can recompute M = M1 inter M2
+		//
+		//      -T [ max(lambda1, mu1)          0         ]  -1
+		// M = P   [                                      ] P
+		//         [        0            max(lambda2, mu2)]
+
+		//get det(P)
+		detP=vp[0][0]*vp[1][1]-vp[0][1]*vp[1][0];
+		if (!detP){
+			printf("vp= [%g %g ; %g %g]\n",vp[0][0],vp[0][1],vp[1][0],vp[1][1]);
+			throw ErrorException(__FUNCT__,exprintf("One of the metric has a zero eigen value or cannot be inverted... (detP=%g)",detM1));
+		}
+		invP[0][0]= 1.0/detP * ( vp[1][1]);
+		invP[0][1]= 1.0/detP * (-vp[0][1]);
+		invP[1][0]= 1.0/detP * (-vp[1][0]);
+		invP[1][1]= 1.0/detP * ( vp[0][0]);
+
+		//Compute lambdai and mui using Rayleigh formula: lambdai = ei' M1 ei
+		lambda1=vp[0][0]*(a11*vp[0][0]+a21*vp[1][0]) + vp[1][0]*(a21*vp[0][0]+a22*vp[1][0]);
+		lambda2=vp[0][1]*(a11*vp[0][1]+a21*vp[1][1]) + vp[1][1]*(a21*vp[0][1]+a22*vp[1][1]);
+		mu1    =vp[0][0]*(b11*vp[0][0]+b21*vp[1][0]) + vp[1][0]*(b21*vp[0][0]+b22*vp[1][0]);
+		mu2    =vp[0][1]*(b11*vp[0][1]+b21*vp[1][1]) + vp[1][1]*(b21*vp[0][1]+b22*vp[1][1]);
+
+		//finally compute M
+		M.a11=invP[0][0]*Max(lambda1,mu1)*invP[0][0] + invP[1][0]*Max(lambda2,mu2)*invP[1][0];
+		M.a21=invP[0][0]*Max(lambda1,mu1)*invP[0][1] + invP[1][0]*Max(lambda2,mu2)*invP[1][1];
+		M.a22=invP[0][1]*Max(lambda1,mu1)*invP[0][1] + invP[1][1]*Max(lambda2,mu2)*invP[1][1];
+		*pM=M;
+	}
+	/*}}}1*/
+	/*FUNCTION SimultaneousMatrixReduction_bamg{{{1*/
+	void SimultaneousMatrixReduction_bamg( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
+		double a11=M1.a11,a21=M1.a21,a22=M1.a22;
+		double b11=M2.a11,b21=M2.a21,b22=M2.a22;
+		//  M1 v = l M2 v
+		// (M1 - l M2) v =0
+		// det (M1 - l M2) =0
+		// det (M1 - l M2) = a l^2 + b l + c;
+		// = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
+		const double /*c11 = a11*a11,*/ c21= a21*a21;
+		const double /*d11 = b11*b11,*/ d21= b21*b21;
+		const double a=b11*b22 - d21;
+		const double b=-a11*b22-a22*b11+2*a21*b21;
+		const double c=-c21+a11*a22;
+		const double bb = b*b,ac= a*c;
+		const double delta = bb - 4 * ac;
+		if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
+			// racine double;
+			if (Abs(a) < 1.e-30 )
+			 l1 = l2 = 0;
+			else 
+			 l1=l2=-b/(2*a); 
+			V= D2xD2(1,0,0,1);
+		}
+		else {
+			const double delta2 = sqrt(delta);
+			l1= (-b - delta2)/(2*a);
+			l2= (-b + delta2)/(2*a);
+			// M1 v = l M2 v
+			//  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
+			// so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
+			double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
+			double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
+			double vp1x,vp1y,vp2x,vp2y;
+
+			if(s1 < s0)
+			 s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
+			else
+			 s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
+
+			v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
+			s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
+			if(s1 < s0)
+			 s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
+			else
+			 s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
+			V=D2xD2(vp1x,vp2x,vp1y,vp2y);
+		}
+		return;
+	}
+	/*}}}1*/
 	/*FUNCTION abscisseInterpole{{{1*/
 	Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim) { 
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2927)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2928)
@@ -2021,8 +2021,4 @@
 					//Get eigen values and vectors of Miv
 					MatVVP2x2 Vp(Miv);
-					if(iv==2){
-						printf("Hessien' = [%g %g %g]\n",dxdx_vertex[iv]*ci,dxdy_vertex[iv]*ci,dydy_vertex[iv]*ci);
-						Vp.Echo();
-					}
 
 					//move eigen valuse to their absolute values
@@ -2037,5 +2033,4 @@
 					if(iv==2){
 						MVp.Echo();
-						vertices[iv].m.Echo();
 					}
 
@@ -2044,9 +2039,6 @@
 					if(iv==2){
 						vertices[iv].m.Echo();
+						throw ErrorException(__FUNCT__,exprintf(""));
 					}
-				}
-
-				for(i=0;i<3;i++){
-					vertices[i].m.Echo();
 				}
 			}//for all fields