Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2931)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2932)
@@ -149,5 +149,7 @@
 
 		BTh.IntersectGeomMetric(errg,iso);
+		if (verbosity>1) printf("   Smoothing Metric...");
 		if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation);
+		if (verbosity>1) printf(" done\n");
 		BTh.MaxSubDivision(maxsubdiv);
 		BTh.BoundAnisotropy(anisomax,hminaniso);
Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2931)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2932)
@@ -64,5 +64,4 @@
 					Real8  b,const  MetricAnIso mb);
 		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);} 
@@ -123,7 +122,5 @@
 	}
 
-	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]);;
+	void  SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V);
 
 	inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) {
Index: /issm/trunk/src/c/Bamgx/R2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/R2.h	(revision 2931)
+++ /issm/trunk/src/c/Bamgx/R2.h	(revision 2932)
@@ -8,4 +8,5 @@
 			  R x,y;
 			  void Echo(){
+				  printf("Member of P2:\n");
 				  printf("   x: %g\n",x);
 				  printf("   y: %g\n",y);
@@ -35,4 +36,9 @@
 		  public:
 		  P2<R,RR> x,y; 
+		  void Echo(){
+			  printf("Member of P2xP2:\n");
+			  printf("   x.x: %g   x.y: %g\n",x.x,x.y);
+			  printf("   y.x: %g   y.x: %g\n",y.x,y.y);
+		  }
 		  P2xP2 (): x(),y()  {}
 		  P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {}
Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2931)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2932)
@@ -69,34 +69,4 @@
 	}
 	/*}}}*/
-	/*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) {
@@ -111,5 +81,5 @@
 		double l1,l2;
 
-		SimultaneousMatrixReduction_bamg(*this,M2,l1,l2,M);
+		SimultaneousMatrixReduction(*this,M2,l1,l2,M);
 
 		R2 v1(M.x.x,M.y.x);
@@ -127,77 +97,11 @@
 			a11=Mi.x.x;
 			a21=0.5*(Mi.x.y+Mi.y.x);
-			a22=Mi.y.y; }
-			return r;
+			a22=Mi.y.y;
+		}
+		return r;
 	}
 	/*}}}1*/
 
 	/*Intermediary*/
-	/*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*/
 	/*FUNCTION LengthInterpole{{{1*/
 	Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB) {
@@ -267,77 +171,5 @@
 	/*}}}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) {
+	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;
@@ -353,6 +185,7 @@
 		const double c=-c21+a11*a22;
 		const double bb = b*b,ac= a*c;
+		const double b2a = b/(2*a);
 		const double delta = bb - 4 * ac;
-		if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
+		if ( (bb + Abs(ac) < 1.0e-34 ) ||  (delta < 1.0e-6*bb) ){
 			// racine double;
 			if (Abs(a) < 1.e-30 )
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2931)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2932)
@@ -2008,5 +2008,4 @@
 				}
 
-
 				/*Compute Metric from Hessian*/
 
@@ -2033,5 +2032,5 @@
 
 					//Apply Metric to vertex
-					vertices[iv].m.IntersectWith_new(MVp);
+					vertices[iv].m.IntersectWith(MVp);
 				}
 			}//for all fields