Index: /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2838)
+++ /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2838)
@@ -0,0 +1,44 @@
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <ctime>
+
+#include "../Mesh2.h"
+#include "../QuadTree.h"
+#include "../SetOfE4.h"
+
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+namespace bamg {
+
+	/*Constructor{{{1*/
+	MatVVP2x2::MatVVP2x2(const MetricAnIso M){
+		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);
+		if ( n2 < 1e-30) 
+		 lambda1=lambda2=0,v.x=1,v.y=0;
+		else if (delta < eps*n2)
+		  { 
+			lambda1=lambda2=-b/2, v.x=1,v.y=0;
+		  }
+		else 
+		  {  //    ---  construction  de 2 vecteur dans (Im ( A - D(i) Id) ortogonal 
+			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/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2838)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2838)
@@ -0,0 +1,245 @@
+#include <cstdio>
+#include <string.h>
+#include <cmath>
+#include <time.h>
+#include <iostream>
+
+#include "../Mesh2.h"
+#include "../QuadTree.h"
+#include "../SetOfE4.h"
+
+using namespace std;
+
+namespace bamg {
+
+	SaveMetricInterpole  LastMetricInterpole;
+
+	/*Constructor/Destructor*/
+	/*FUNCTION MetricAnIso::MetricAnIso(const Real8  a[3],const  MetricAnIso m0, const  MetricAnIso m1,const  MetricAnIso  m2 ){{{1*/ 
+	MetricAnIso::MetricAnIso(const Real8  a[3],const  MetricAnIso m0, const  MetricAnIso m1,const  MetricAnIso m2 ){
+		MetricAnIso mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11,
+					a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21,
+					a[0]*m0.a22 + a[1]*m1.a22 + a[2]*m2.a22);
+
+		MatVVP2x2 vab(mab);
+
+		R2 v1(vab.v.x,vab.v.y);
+		R2 v2(-v1.y,v1.x);
+
+		Real8 h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1);
+		Real8 h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2);
+
+		vab.lambda1 =  1 / (h1*h1);
+		vab.lambda2 =  1 / (h2*h2);
+		*this = vab;
+	}
+	/*}}}1*/
+	/*FUNCTION MetricAnIso::MetricAnIso( Real8  a,const  MetricAnIso ma, Real8  b,const  MetricAnIso mb){{{1*/
+	MetricAnIso::MetricAnIso( Real8  a,const  MetricAnIso ma, Real8  b,const  MetricAnIso mb) { 
+		MetricAnIso mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
+		MatVVP2x2 vab(mab);
+
+		R2 v1(vab.v.x,vab.v.y);
+		R2 v2(-v1.y,v1.x);
+
+
+		Real8 h1 = a / ma(v1) + b / mb(v1);
+		Real8 h2 = a / ma(v2) + b / mb(v2);
+		vab.lambda1 =  1 / (h1*h1);
+		vab.lambda2 =  1 / (h2*h2);
+		*this = vab;
+	}
+	/*}}}1*/
+
+	/*Others*/
+	/*FUNCTION ReductionSimultanee{{{1*/
+	void ReductionSimultanee( 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;
+		ReductionSimultanee(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);
+	}
+	/*}}}1*/
+	/*FUNCTION MetricAnIso::IntersectWith{{{1*/
+	int MetricAnIso::IntersectWith(const MetricAnIso M2) {
+		int r=0;
+		MetricAnIso & M1 = *this;
+		D2xD2 M;
+		double l1,l2;
+
+		ReductionSimultanee(*this,M2,l1,l2,M);
+		R2 v1(M.x.x,M.y.x);
+		R2 v2(M.x.y,M.y.y);
+		double l11=M1(v1,v1);
+		double l12=M1(v2,v2);
+		double l21=M2(v1,v1);
+		double l22=M2(v2,v2);
+		if ( l11 < l21 )  r=1,l11=l21;
+		if ( l12 < l22 )  r=1,l12=l22; 
+		if (r) { // change
+			D2xD2 M_1(M.inv());
+			D2xD2 D(l11,0,0,l12); 
+			D2xD2 Mi(M_1.t()*D*M_1);
+			a11=Mi.x.x;
+			a21=0.5*(Mi.x.y+Mi.y.x);
+			a22=Mi.y.y; }
+			return r;
+	}
+	/*}}}1*/
+	/*FUNCTION LengthInterpole{{{1*/
+	Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB) {
+		Real8 k=1./2.;
+		int level=0;
+		static int kkk=0;
+		static  Metric Ms1[32],Ms2[32];
+		static Real8 lMs1[32],lMs2[32];
+		static double K[32];
+		Real8 l=0,sss=0;
+		Ms1[level]=Ma;
+		Ms2[level]=Mb;
+		Real8 sa =  Ma(AB);
+		Real8 sb =  Mb(AB);
+		lMs1[level]=sa;
+		lMs2[level]=sb;
+		K[level]=k;
+		level++;
+		int i=0;
+		Real8 * L= LastMetricInterpole.L, *S = LastMetricInterpole.S;
+		Real8  sstop = 0.1; // Max(0.6,(sa+sb)/5000);
+		while (level) {
+			level--;
+			Metric M1=Ms1[level];
+			Metric M2=Ms2[level];
+			k=K[level];
+			Real8 s1=  lMs1[level];
+			Real8 s2=  lMs2[level];
+
+			Real8 s= (s1+s2)*k;
+			if( s > sstop   && level < 30 && i < 500-level ) {
+				Metric Mi(0.5,M1,0.5,M2);
+				Real8 si = Mi(AB);
+				if( Abs((s1+s2)-(si+si)) > s1*0.001) 
+				  {
+					k=k/2;
+					// we begin by the end to walk in the correct sens from a to b
+					// due to the stack 
+					Ms1[level]=Mi;
+					Ms2[level]=M2;
+					lMs1[level]=si;
+					lMs2[level]=s2;
+					K[level]=k;
+					level++;
+					Ms1[level]=M1;
+					Ms2[level]=Mi;
+					lMs1[level]=s1;
+					lMs2[level]=si;
+					K[level]=k;
+					level++;
+				  }
+				else
+				 L[i]= l += s,S[i]=sss+=k,i++;
+			}
+			else 
+			 L[i]= l += s,S[i]=sss+=k,i++;//cout << i << " l = " << l << " sss = " << sss << endl;
+		}
+		// warning for optimisation S is in [0:0.5] not in [0:1]
+		assert(i<512);
+		LastMetricInterpole.lab=l;
+		LastMetricInterpole.opt=i;
+		if (i>200 && kkk++<10)
+		 cout << "Warning LengthInterpole: ( i = " << i << " l = " << l << " sss " << sss << " ) " << sstop <<endl;
+		return l;
+	}
+	/*}}}1*/
+	/*FUNCTION abscisseInterpole{{{1*/
+	Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim) { 
+		if(!optim)  LengthInterpole(Ma,Mb,AB);
+		Real8 l  = s* LastMetricInterpole.lab,r;
+		int j=LastMetricInterpole.opt-1,i=0,k;
+
+		Real8 * L= LastMetricInterpole.L, *S = LastMetricInterpole.S;
+		// warning for optimisation S is the abcisse in [0:0.5]
+		// and L is le lenght 
+		if(l<=L[0])
+		 r=2*S[0]*l/L[0];
+		else if (l>=L[j])
+		 r=1;
+		else 
+		  {
+			while (j-i>1)
+			  {
+				k= (i+j)/2;
+				if(l<=L[k])
+				 j=k;// l<=L[j] 
+				else
+				 i=k; //  L[i]<l
+			  };
+			//   cout  << i << " " << j  <<" " << L[i] << " " << L[j] << " " << S[i] << " " <<  S[j]  << " l=" << l << endl;
+			if (i==j)
+			 r = 2*S[i];
+			else
+			 r =  2*(S[i]*(L[j]-l)+ S[j]*(l-L[i]))/(L[j]-L[i]);
+		  }
+		assert(r<=1 && r>=0);
+		return r ;
+	}
+	/*}}}1*/
+
+}