Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2923)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2924)
@@ -148,6 +148,6 @@
 
 		public:
-			Triangle * t; // le triangle 
-			int  a; // le numero de l arete
+			Triangle* t; // le triangle 
+			int  a;      // le numero de l arete
 
 			TriangleAdjacent(Triangle  * tt,int  aa): t(tt),a(aa &3) {};
Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2923)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2924)
@@ -91,6 +91,6 @@
 		void  Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;}
 
-		void Minh(double h) {Max(1.0/(h*h));}
-		void Maxh(double h) {Min(1.0/(h*h));}
+		void Minh(double h) {Min(1.0/(h*h));}
+		void Maxh(double h) {Max(1.0/(h*h));}
 		void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
 		MatVVP2x2(const MetricAnIso );
@@ -122,5 +122,5 @@
 	}
 
-	void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
+	void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
 	MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
 
Index: /issm/trunk/src/c/Bamgx/R2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/R2.h	(revision 2923)
+++ /issm/trunk/src/c/Bamgx/R2.h	(revision 2924)
@@ -7,4 +7,8 @@
 		  public:  
 			  R x,y;
+			  void Echo(){
+				  printf("   x: %g\n",x);
+				  printf("   y: %g\n",y);
+			  }
 			  P2 () :x(0),y(0) {};
 			  P2 (R a,R b)  :x(a),y(b)  {}
Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2923)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2924)
@@ -71,4 +71,9 @@
 	/*FUNCTION MetricAnIso::IntersectWith{{{1*/
 	int MetricAnIso::IntersectWith(const 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
+		 * */
+
 		int r=0;
 		MetricAnIso & M1 = *this;
@@ -76,5 +81,6 @@
 		double l1,l2;
 
-		ReductionSimultanee(*this,M2,l1,l2,M);
+		SimultaneousMatrixReduction(*this,M2,l1,l2,M);
+
 		R2 v1(M.x.x,M.y.x);
 		R2 v2(M.x.y,M.y.y);
@@ -97,6 +103,6 @@
 
 	/*Intermediary*/
-	/*FUNCTION ReductionSimultanee{{{1*/
-	void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
+	/*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;
@@ -113,6 +119,5 @@
 		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 ) )
-		  {
+		if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
 			// racine double;
 			if (Abs(a) < 1.e-30 )
@@ -153,5 +158,5 @@
 		D2xD2 M;
 		double l1,l2;
-		ReductionSimultanee(M1,M2,l1,l2,M);
+		SimultaneousMatrixReduction(M1,M2,l1,l2,M);
 		R2 v0(M.x.x,M.y.x);
 		R2 v1(M.x.y,M.y.y);
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2923)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2924)
@@ -700,6 +700,6 @@
 					MatVVP2x2 Vp(M/coef);
 
-					Vp.Maxh(hmin);
-					Vp.Minh(hmax);
+					Vp.Maxh(hmax);
+					Vp.Minh(hmin);
 					vertices[i].m = Vp;
 				}
@@ -1272,6 +1272,6 @@
 
 					//modify eigen values according to hmin and hmax
-					Vp.Maxh(hmin);
-					Vp.Minh(hmax);
+					Vp.Maxh(hmax);
+					Vp.Minh(hmin);
 
 					//multiply eigen values by coef
@@ -1577,5 +1577,4 @@
 						int   nbb=0;
 						Real8 dd = detT[i];
-						Real8 lla,llb,llc,llf;
 						Real8 taa[3][3],bb[3];
 
@@ -1586,23 +1585,53 @@
 							Triangle*        tt = ta;
 
+
+							/* V is the vertex opposed to the edge shared by the
+							 * current triangle i and its neighbor ta
+							 *         C       
+							 *       / | \     
+							 *      /  |  \    
+							 *     /   |   \   
+							 * V  / ta | i  \ B
+							 *    \    |    /  
+							 *     \   |   /   
+							 *      \  |  /
+							 *       \ | /     
+							 *         A
+							 * lA=area(VBC)/Area(ABC)
+							 * lB=area(AVC)/Area(ABC)
+							 * lC=area(ABV)/Area(ABC)
+							 */
+
+							//first, get V
+							Vertex &v = *ta.OppositeVertex();
+							Int4   iV = Number(v);
+
 							//if the adjacent triangle is not a boundary triangle:
 							if (tt && tt->link){
 								Vertex &v = *ta.OppositeVertex();
-								R2     V = v;
+								R2      V = v;
 								Int4   iV = Number(v);
+
+								//get lA lB and lC
 								Real8  lA = bamg::Area2(V,B,C)/dd;
 								Real8  lB = bamg::Area2(A,V,C)/dd;
 								Real8  lC = bamg::Area2(A,B,V)/dd;
+
+								//fill taa
 								taa[0][j] = lB*lC;
 								taa[1][j] = lC*lA;
 								taa[2][j] = lA*lB;
-								lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;
-								bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ;
+
+								//fill bb
+								bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ; //value of the solution on V minus...
 							}
 							else{
+								//there might be a problem here: if 2 nodes are inline lA, lB or lC
+								//is 0. This means that one column has only one term.
+								//if instead of 0.1 we use 0, and taa[j][j]=1, we might have det(taa)=0
 								nbb++;
-								taa[0][j]=0;
-								taa[1][j]=0;
-								taa[2][j]=0;
+								taa[0][j]=0.1;
+								taa[1][j]=0.1;
+								taa[2][j]=0.1;
 								taa[j][j]=1;
 								bb[j]=0;
@@ -1611,14 +1640,13 @@
 
 						// resolution of 3x3 linear system transpose
-						Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);		
+						Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);	
 						Real8 cBC   =  det3x3(bb,taa[1],taa[2]);
 						Real8 cCA   =  det3x3(taa[0],bb,taa[2]);
 						Real8 cAB   =  det3x3(taa[0],taa[1],bb);
-
 						if (!det33){
-							throw ErrorException(__FUNCT__,exprintf("!det33"));
+							throw ErrorException(__FUNCT__,exprintf("for triangle %i, det33==0! cannot compute theHessian matrix using Hessiantype=1",i+1));
 						}
+
 						// computation of the Hessian in the element 
-
 						// H( li*lj) = grad li grad lj + grad lj grad lj
 						// grad li = njk  / detT ; with i j k =(A,B,C)
@@ -1673,5 +1701,5 @@
 				for (int xy = 0;xy<3;xy++) {
 					dd = d2[xy];
-					// do leat 2 iteration for boundary problem
+					// do least 2 iteration for boundary problem
 					for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
 						for (i=0;i<nbt;i++) 
@@ -1746,6 +1774,6 @@
 
 					//modify eigen values according to hmin and hmax
-					Vp.Maxh(hmin);
-					Vp.Minh(hmax);
+					Vp.Maxh(hmax);
+					Vp.Minh(hmin);
 
 					//multiply eigen values by coef
@@ -1789,8 +1817,258 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::BuildMetric2{{{1*/
+	/*FUNCTION Triangles::BuildMetric2 (double P2 projection){{{1*/
 	void Triangles::BuildMetric2(BamgOpts* bamgopts){
 
-		throw ErrorException(__FUNCT__,exprintf("not supported yet"));
+		/*Options*/
+		const int dim = 2;
+		double* s=NULL;
+		Int4 nbsol;
+		int* typsols=NULL;
+		int verbosity;
+		Real8 hmin1;
+		Real8 hmax1;
+		int   i,j,k,iA,iB,iC;
+		int   iv,nbfield;
+
+		/*Recover options*/
+		verbosity=bamgopts->verbose;
+		hmin1=bamgopts->hmin;
+		hmax1=bamgopts->hmax;
+
+		/*Recover options*/
+		verbosity=bamgopts->verbose;
+		hmin1=bamgopts->hmin;
+		hmax1=bamgopts->hmax;
+
+		/*Get and process fields*/
+		s=bamgopts->field;
+		nbsol=1;       //for now, only one field
+		typsols=(int*)xmalloc(1*sizeof(int));
+		typsols[0]=0; // only one dof per node
+
+		//sizeoftype = {1,2,3,4}
+		int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 
+
+		// computation of the number of fields
+		Int4 ntmp = 0;
+		if (typsols){
+			//if there is more than one dof per vertex for one field
+			//increase ntmp to take into account all the fields
+			//if nbsol=1
+			for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
+		}
+		else ntmp = nbsol;
+
+		//n is the total number of fields
+		const Int4 n = ntmp;
+
+		//initialization of some variables
+		Real8   hmin = Max(hmin1,MinimalHmin());
+		Real8   hmax = Min(hmax1,MaximalHmax());
+		double* ss=(double*)s;
+		double  sA,sB,sC;
+		Real8*  detT = new Real8[nbt];
+		Real8*  sumareas = new Real8[nbv];
+		Real8*  alpha= new Real8[nbt*3];
+		Real8*  beta = new Real8[nbt*3];
+		Real8*  dx_elem    = new Real8[nbt];
+		Real8*  dy_elem    = new Real8[nbt];
+		Real8*  dx_vertex  = new Real8[nbv];
+		Real8*  dy_vertex  = new Real8[nbv];
+		Real8*  dxdx_elem  = new Real8[nbt];
+		Real8*  dxdy_elem  = new Real8[nbt];
+		Real8*  dydy_elem  = new Real8[nbt];
+		Real8*  dxdx_vertex= new Real8[nbv];
+		Real8*  dxdy_vertex= new Real8[nbv];
+		Real8*  dydy_vertex= new Real8[nbv];
+
+		//display infos
+		if(verbosity>1) {
+			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
+		}
+
+		//first, build the chains that will be used for the Hessian computation, as weel as the area of each element
+		int head_s[nbv];
+		int next_p[3*nbt];
+		int  p=0;
+		//initialization
+		for(i=0;i<nbv;i++){
+			sumareas[i]=0;
+			head_s[i]=-1;
+		}
+		for(i=0;i<nbt;i++){
+
+			//lopp over the real triangles (no boundary elements)
+			if(triangles[i].link){ 
+
+				//get current triangle t
+				const Triangle &t=triangles[i];
+
+				// coor of 3 vertices 
+				R2 A=t[0];
+				R2 B=t[1];
+				R2 C=t[2];
+
+				//compute triangle determinant (2*Area)
+				Real8 dett = bamg::Area2(A,B,C);
+				detT[i]=dett;
+
+				/*The nodal functions are such that for a vertex A:
+				 *    N_A(x,y)=alphaA x + beta_A y +gamma_A
+				 *    N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
+				 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
+				 * leads to:
+				 *    N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
+				 * and this gives:
+				 *    alpha_A = (yB-yC)/(2*Area(T))*/
+				alpha[i*3+0]=(B.y-C.y)/dett;
+				alpha[i*3+1]=(C.y-A.y)/dett;
+				alpha[i*3+2]=(A.y-B.y)/dett;
+				beta[ i*3+0]=(C.x-B.x)/dett;
+				beta[ i*3+1]=(A.x-C.x)/dett;
+				beta[ i*3+2]=(B.x-A.x)/dett;
+
+				//compute chains
+				for(j=0;j<3;j++){
+					k=Number(triangles[i][j]);
+					next_p[p]=head_s[k];
+					head_s[k]=p++;
+
+					//add area to sumareas
+					sumareas[k]+=dett;
+				}
+
+			}
+		}
+
+		//for all Solutions
+		for (Int4 nusol=0;nusol<nbsol;nusol++) {
+			int   nbfield=typsols?sizeoftype[typsols[nusol]]:1; 
+
+			//loop over all the fields of the solution
+			for (Int4 nufield=0;nufield<nbfield;nufield++,s++){
+				//ss++ so that for each iteration ss points toward the right field
+
+				//initialize the hessian and gradient matrices
+				for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
+
+				//1: Compute gradient for each element (exact)
+				for (i=0;i<nbt;i++){
+					if(triangles[i].link){
+						// number of the 3 vertices
+						iA = Number(triangles[i][0]);
+						iB = Number(triangles[i][1]);
+						iC = Number(triangles[i][2]);
+
+						// value of the P1 fonction on 3 vertices 
+						sA = ss[iA*n];
+						sB = ss[iB*n];
+						sC = ss[iC*n];
+
+						//gradient = (sum alpha_i s_i, sum_i beta_i s_i)
+						dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
+						dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
+					}
+				}
+
+				//2: then compute a gradient for each vertex using a P2 projection
+				for(i=0;i<nbv;i++){
+					for(p=head_s[i];p!=-1;p=next_p[p]){
+						//Get triangle number
+						k=(Int4)(p/3);
+						dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
+						dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
+					}
+				}
+
+				//3: compute Hessian matrix on each element
+				for (i=0;i<nbt;i++){
+					if(triangles[i].link){
+						// number of the 3 vertices
+						iA = Number(triangles[i][0]);
+						iB = Number(triangles[i][1]);
+						iC = Number(triangles[i][2]);
+
+						//Hessian
+						dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
+						dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
+						dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
+					}
+				}
+
+				//4: finaly compute Hessian on each vertex using the second P2 projection
+				for(i=0;i<nbv;i++){
+					for(p=head_s[i];p!=-1;p=next_p[p]){
+						//Get triangle number
+						k=(Int4)(p/3);
+						dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
+						dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
+						dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
+					}
+				}
+
+
+				/*Compute Metric from Hessian*/
+
+				//compute multiplicative coefficient (2/9 because it is 2d)
+				Real8 ci=(2.0/9.0)*(1.0/(bamgopts->err));
+
+				for ( iv=0;iv<nbv;iv++){
+
+					//initialize metric Miv with ci*H
+					Metric Miv(dxdx_vertex[iv]*ci,dxdy_vertex[iv]*ci,dydy_vertex[iv]*ci);
+
+					//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
+					Vp.Abs();
+
+					//modify eigen values according to hmin and hmax
+					Vp.Maxh(hmax);
+					Vp.Minh(hmin);
+
+					//rebuild Metric from Vp
+					Metric MVp(Vp);
+					if(iv==2){
+						MVp.Echo();
+						vertices[iv].m.Echo();
+					}
+
+					//Apply Metric to vertex
+					vertices[iv].m.IntersectWith(MVp);
+					if(iv==2){
+						vertices[iv].m.Echo();
+					}
+				}
+
+				for(i=0;i<3;i++){
+					vertices[i].m.Echo();
+				}
+			}//for all fields
+		}//for all solutions
+		for(i=0;i<5;i++){
+			vertices[i].m.Echo();
+		}
+
+		//clean up
+		delete [] detT;
+		delete [] alpha;
+		delete [] beta;
+		delete [] sumareas;
+		delete [] dx_elem;
+		delete [] dy_elem;
+		delete [] dx_vertex;
+		delete [] dy_vertex;
+		delete [] dxdx_elem;
+		delete [] dxdy_elem;
+		delete [] dydy_elem;
+		delete [] dxdx_vertex;
+		delete [] dxdy_vertex;
+		delete [] dydy_vertex;
 	}
 	/*}}}1*/
@@ -3650,4 +3928,7 @@
 		BuildMetric1(bamgopts);
 	}
+	else if (Hessiantype==2){
+		BuildMetric2(bamgopts);
+	}
 	else{
 		throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (0->use Green formula, 1-> from P2 on 4T, 2-> double P2 projection)",Hessiantype));