Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2916)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2917)
@@ -108,15 +108,17 @@
 
 	inline void  MatVVP2x2::BoundAniso2(const Real8 coef){
-		if (coef<=1.00000000001) 
-		 if (lambda1 < lambda2)
-		  lambda1 = bamg::Max(lambda1,lambda2*coef);
-		 else
-		  lambda2 = bamg::Max(lambda2,lambda1*coef);
-		else  // a verifier 
-		 if (lambda1 > lambda2)
-		  lambda1 = bamg::Min(lambda1,lambda2*coef);
-		 else
-		  lambda2 = bamg::Min(lambda2,lambda1*coef);
-	  }
+		if (coef<=1.00000000001){
+			if (lambda1 < lambda2)
+			 lambda1 = bamg::Max(lambda1,lambda2*coef);
+			else
+			 lambda2 = bamg::Max(lambda2,lambda1*coef);
+		}
+		else{  //TO BE CHECKED
+			if (lambda1 > lambda2)
+			 lambda1 = bamg::Min(lambda1,lambda2*coef);
+			else
+			 lambda2 = bamg::Min(lambda2,lambda1*coef);
+		}
+	}
 
 	void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
Index: /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2916)
+++ /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2917)
@@ -24,24 +24,26 @@
 		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 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)
-		  { 
+
+		//if norm(M)<10^30 -> M=zeros(2,2)
+		if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0;
+
+		//if ???
+		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;
+		}
+
+		//general case: construction of 2 eigen vectors
+		else {
+			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/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2916)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2917)
@@ -854,33 +854,35 @@
 	/*FUNCTION Triangles::BoundAnisotropy{{{1*/
 	void  Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso) {
+
 		long int verbosity=0;
-
 		double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
+
+		//display info
 		if (verbosity > 1)  printf("   BoundAnisotropy by %g\n",anisomax);
-		Real8 h1=1.e30,h2=1e-30,rx=0;
+
+		Real8 h1=1.e30,h2=1e-30;
 		Real8 coef = 1./(anisomax*anisomax);
-		Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-		for (Int4 i=0;i<nbv;i++)
-		  {
-
+		Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0;  
+
+		//loop over all vertices
+		for (Int4 i=0;i<nbv;i++){
 			MatVVP2x2 Vp(vertices[i]);
 			double lmax=Vp.lmax();
-			h1=Min(h1,Vp.lmin());
-			h2=Max(h2,Vp.lmax());
-			rx = Max(rx,Vp.Aniso2());
-
-			Vp *= Min(lminaniso,lmax)/lmax;
-
+			h1 =Min(h1,Vp.lmin());
+			h2 =Max(h2,Vp.lmax());
+			Vp*=Min(lminaniso,lmax)/lmax;
 			Vp.BoundAniso2(coef);
-
-			hn1=Min(hn1,Vp.lmin());
-			hn2=Max(hn2,Vp.lmax());
-			rnx = Max(rnx,Vp.Aniso2());
-
-
 			vertices[i].m = Vp;
 
-		  }
-
+			//info to be displayed
+			if (verbosity>2){
+				hn1=Min(hn1,Vp.lmin());
+				hn2=Max(hn2,Vp.lmax());
+				rx =Max(rx,Vp.Aniso2());
+				rnx= Max(rnx,Vp.Aniso2());
+			}
+		}
+
+		//display info
 		if (verbosity>2){
 			printf("      input:  Hmin = %g, Hmax = %g, factor of anisotropy max  = %g\n",pow(h2,-0.5),pow(h1,-0.5),pow(rx,0.5));
@@ -890,6 +892,5 @@
 	/*}}}1*/
 	/*FUNCTION Triangles::ConsGeometry{{{1*/
-	void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges) // construct a geometry if no geo 
-	  {
+	void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges){
 		//  if equiedges existe taille nbe 
 		//   equiedges[i]/2 == i  original
@@ -2729,21 +2730,17 @@
 /*FUNCTION Triangles::IntersectConsMetric{{{1*/
 void Triangles::IntersectConsMetric(BamgOpts* bamgopts){
-	   /*const double* s,const Int4 nbsol,const int * typsols,
-		 *      const  Real8 hmin1,const Real8 hmax1,const Real8 coef,
-		 *           const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,
-		 *                const int Rescaling,const double power,const int Hessiantype)*/
 	//  the array of solution s is store    
 	// sol0,sol1,...,soln    on vertex 0
 	//  sol0,sol1,...,soln   on vertex 1
 	//  etc.
-	//  choise = 0 =>  H is computed with green formula
-	//  otherwise  =>  H is computed from P2 on 4T 
+	//  Hessiantype = 0 =>  H is computed with green formula
+	//  Hessiantype = 1 =>  H is computed from P2 on 4T 
 
 	/*Options*/
 	const int dim = 2;
 	int AbsError;
-	double* s;
+	double* s=NULL;
 	Int4 nbsol;
-	int * typsols;
+	int* typsols=NULL;
 	Real8 hmin1;
 	Real8 hmax1;
@@ -2759,7 +2756,4 @@
 	/*Recover options*/
 	verbosity=bamgopts->verbose;
-	s=bamgopts->field;
-	nbsol=1;   //for now, only one field
-	typsols=0; // only one dof per node
 	AbsError=bamgopts->AbsError;   
 	CutOff=bamgopts->cutoff;
@@ -2776,45 +2770,56 @@
 	coef=sqrt(bamgopts->err)*coef;
 
+	/*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 nb of field 
+	// 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
-
+	//n is the total number of fields
 	const Int4 n = ntmp;
 
-	Int4 i,k,iA,iB,iC,iv;
-	R2 O(0,0);
-	int RelativeMetric = CutOff>1e-30;
-	Real8 hmin = Max(hmin1,MinimalHmin());
-	Real8 hmax = Min(hmax1,MaximalHmax());
-	Real8 coef2 = 1/(coef*coef);
-
+	//initialization of some variables
+	Int4    i,k,iA,iB,iC,iv;
+	R2      O(0,0);
+	int     RelativeMetric = CutOff>1e-30;
+	Real8   hmin = Max(hmin1,MinimalHmin());
+	Real8   hmax = Min(hmax1,MaximalHmax());
+	Real8   coef2 = 1/(coef*coef);
+	double* ss=(double*)s;
+	double  sA,sB,sC;
+	Real8*  detT = new Real8[nbt];
+	Real8*  Mmass= new Real8[nbv];
+	Real8*  Mmassxx= new Real8[nbv];
+	Real8*  dxdx= new Real8[nbv];
+	Real8*  dxdy= new Real8[nbv];
+	Real8*  dydy= new Real8[nbv];
+	Real8*  workT= new Real8[nbt];
+	Real8*  workV= new Real8[nbv];
+	int*    OnBoundary = new int[nbv];
+
+	//display infos
 	if(verbosity>1) {
 		printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
 		printf("      coef = %g\n",coef); 
 		printf("      hmin = %g hmax = %g\n",hmin,hmax); 
-		printf("      anisomax = %g nb Jacobi = %i, power = %i\n",anisomax,NbJacobi,power); 
+		printf("      anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power); 
 		if (RelativeMetric) printf("      RelativeErr with CutOff= %g\n",CutOff);
 		else printf("      Absolute error\n");
 	}
-	double* ss=(double*)s;
-
-	double sA,sB,sC;
-
-	Real8 *detT = new Real8[nbt];
-	Real8 *Mmass= new Real8[nbv];
-	Real8 *Mmassxx= new Real8[nbv];
-	Real8 *dxdx= new Real8[nbv];
-	Real8 *dxdy= new Real8[nbv];
-	Real8 *dydy= new Real8[nbv];
-	Real8 *workT= new Real8[nbt];
-	Real8 *workV= new Real8[nbv];
-	int *OnBoundary = new int[nbv];
+
+	//initialize Mmass, OnBoundary and Massxx by zero
 	for (iv=0;iv<nbv;iv++){
 		Mmass[iv]=0;
@@ -2823,67 +2828,89 @@
 	}
 
-	for (i=0;i<nbt;i++) 
-	 if(triangles[i].link){ // the real triangles 
-		 const Triangle &t=triangles[i];
-		 // coor of 3 vertices 
-		 R2 A=t[0];
-		 R2 B=t[1];
-		 R2 C=t[2];
-
-
-		 // number of the 3 vertices
-		 iA = Number(t[0]);
-		 iB = Number(t[1]);
-		 iC = Number(t[2]);
-
-		 Real8 dett = bamg::Area2(A,B,C);
-		 detT[i]=dett;
-		 dett /= 6;
-
-		 // construction of on boundary 
-		 int nbb =0;
-		 for(int j=0;j<3;j++)
-			{
-			 Triangle *ta=t.Adj(j);
-			 if ( ! ta || !ta->link) // no adj triangle => edge on boundary
-			  OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1,
-				 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1,
-				 nbb++;
-			}
-
-		 workT[i] = nbb;
-		 Mmass[iA] += dett;
-		 Mmass[iB] += dett;
-		 Mmass[iC] += dett;
-
-		 if((nbb==0)|| !Hessiantype){
-			 Mmassxx[iA] += dett;
-			 Mmassxx[iB] += dett;
-			 Mmassxx[iC] += dett;
-			}
-		}
-	 else workT[i]=-1;
-
+	//Build detT Mmas Mmassxx workT and OnBoundary
+	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];
+
+			// number of the 3 vertices
+			iA = Number(t[0]);
+			iB = Number(t[1]);
+			iC = Number(t[2]);
+
+			//compute triangle determinant (2*Area)
+			Real8 dett = bamg::Area2(A,B,C);
+			detT[i]=dett;
+			dett /= 6;
+
+			// construction of OnBoundary (flag=1 if on boundary, else 0)
+			int nbb=0;
+			for(int j=0;j<3;j++){
+				//get adjacent triangle
+				Triangle *ta=t.Adj(j);
+				//if there is no adjacent triangle, the edge of the triangle t is on boundary
+				if ( !ta || !ta->link){
+					//mark the two vertices of the edge as OnBoundary
+					OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1;
+					OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1;
+					nbb++;
+				}
+			}
+
+			//number of vertices on boundary for current triangle t
+			workT[i] = nbb;
+
+			//Build Mmass Mmass[i] = Mmass[i] + Area/3
+			Mmass[iA] += dett;
+			Mmass[iB] += dett;
+			Mmass[iC] += dett;
+
+			//Build Massxx = Mmass
+			if((nbb==0) || (Hessiantype==0)){
+				Mmassxx[iA] += dett;
+				Mmassxx[iB] += dett;
+				Mmassxx[iC] += dett;
+			}
+		}
+
+		//else: the triangle is a boundary triangle -> workT=-1
+		else workT[i]=-1;
+	}
+
+	//for all Solution  
 	for (Int4 nusol=0;nusol<nbsol;nusol++) {
-		//for all Solution  
 
 		Real8 smin=ss[0],smax=ss[0];
-
 		Real8 h1=1.e30,h2=1e-30,rx=0;
 		Real8 coef = 1./(anisomax*anisomax);
 		Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-		int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 
-		if (nbfield == 1) 
-		 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
-			 dxdx[iv]=dxdy[iv]=dydy[iv]=0;
-			 smin=Min(smin,ss[k]);
-			 smax=Max(smax,ss[k]);
-		 }
+		int   nbfield=typsols?sizeoftype[typsols[nusol]]:1; 
+
+		//only one field
+		if (nbfield == 1) {
+			//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
+			for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
+				dxdx[iv]=dxdy[iv]=dydy[iv]=0;
+				smin=Min(smin,ss[k]);
+				smax=Max(smax,ss[k]);
+			}
+		}
+
+		//vectorial case
 		else{
-			//  cas vectoriel 
-			for ( iv=0,k=0; iv<nbv; iv++,k+=n ){	
+			for (iv=0,k=0;iv<nbv;iv++,k+=n ){
+				//compute v = √sum(s[i]^2)
 				double v=0;		     
-				for (int i=0;i<nbfield;i++) 
-				 v += ss[k+i]*ss[k+i];
+				for (int i=0;i<nbfield;i++){
+					v += ss[k+i]*ss[k+i];
+				}
 				v = sqrt(v);
 				smin=Min(smin,v);
@@ -2891,178 +2918,202 @@
 			}
 		}
-		Real8 sdelta = smax-smin;
+		Real8 sdelta=smax-smin;
 		Real8 absmax=Max(Abs(smin),Abs(smax));
-		Real8 cnorm = Rescaling ? coef2/sdelta : coef2;
-
+		Real8 cnorm =Rescaling ? coef2/sdelta : coef2;
+
+		//display info
 		if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield);
 
-		if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) {
+		//skip constant field
+		if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){
 			if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
 			continue;
 		}
 
-		double *sf  = ss; 
+		//pointer toward ss that is also a pointer toward s (solutions)
+		double* sf=ss; 
+
+		//loop over all the fields of the solution
 		for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){
-			for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-			 dxdx[iv]=dxdy[iv]=dydy[iv]=0;
-			for (i=0;i<nbt;i++) 
-			 if(triangles[i].link)
-				{// for real all triangles 
-				 // coor of 3 vertices 
-				 R2 A=triangles[i][0];
-				 R2 B=triangles[i][1];
-				 R2 C=triangles[i][2];
-
-
-				 // warning the normal is internal and the 
-				 //   size is the length of the edge
-				 R2 nAB = Orthogonal(B-A);
-				 R2 nBC = Orthogonal(C-B);
-				 R2 nCA = Orthogonal(A-C);
-				 // remark :  nAB + nBC + nCA == 0 
-
-				 // number of the 3 vertices
-				 iA = Number(triangles[i][0]);
-				 iB = Number(triangles[i][1]);
-				 iC = Number(triangles[i][2]);
-
-				 // for the test of  boundary edge
-				 // the 3 adj triangles 
-				 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
-				 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
-				 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);
-
-				 // value of the P1 fonction on 3 vertices 
-				 sA = ss[iA*n];
-				 sB = ss[iB*n];
-				 sC = ss[iC*n];
-
-				 R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ;
-				 if(Hessiantype){
-					 int nbb = 0;
-					 Real8 dd = detT[i];
-					 Real8 lla,llb,llc,llf;
-					 Real8  taa[3][3],bb[3];
-					 // construction of the trans of lin system
-					 for (int j=0;j<3;j++)
-						{
-						 int ie = OppositeEdge[j];
-						 TriangleAdjacent ta = triangles[i].Adj(ie);
-						 Triangle *tt = ta;
-						 if (tt && tt->link)
-							{
-							 Vertex &v = *ta.OppositeVertex();
-							 R2 V = v;
-							 Int4 iV = Number(v);
-							 Real8 lA  = bamg::Area2(V,B,C)/dd;
-							 Real8 lB  = bamg::Area2(A,V,C)/dd;
-							 Real8 lC  = bamg::Area2(A,B,V)/dd;
-							 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 ) ;
+			//ss++ so that for each iteration ss points toward the right field
+
+			//initialize the hessian matrix
+			for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
+
+			//loop over the triangles
+			for (i=0;i<nbt;i++){
+
+				//for real all triangles 
+				if(triangles[i].link){
+
+					// coor of 3 vertices 
+					R2 A=triangles[i][0];
+					R2 B=triangles[i][1];
+					R2 C=triangles[i][2];
+
+					//warning: the normal is internal and the size is the length of the edge
+					R2 nAB = Orthogonal(B-A);
+					R2 nBC = Orthogonal(C-B);
+					R2 nCA = Orthogonal(A-C);
+					//note that :  nAB + nBC + nCA == 0 
+
+					// number of the 3 vertices
+					iA = Number(triangles[i][0]);
+					iB = Number(triangles[i][1]);
+					iC = Number(triangles[i][2]);
+
+					// for the test of  boundary edge
+					// the 3 adj triangles 
+					Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
+					Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
+					Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);
+
+					// value of the P1 fonction on 3 vertices 
+					sA = ss[iA*n];
+					sB = ss[iB*n];
+					sC = ss[iC*n];
+
+					/*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))
+					  beta_A = (xC-xB)/(2*Area(T))
+					  and therefore:
+					  grad N_A = nA / detT
+					  for an interpolation of a solution s:
+					  grad(s) = s * sum_{i=A,B,C} grad(N_i) */
+
+					R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i];
+
+					//from P2 on 4T to compute Hessian
+					if(Hessiantype==1){
+						int   nbb=0;
+						Real8 dd = detT[i];
+						Real8 lla,llb,llc,llf;
+						Real8 taa[3][3],bb[3];
+
+						// construction of the transpose of lin system
+						for (int j=0;j<3;j++){
+							int              ie = OppositeEdge[j];
+							TriangleAdjacent ta = triangles[i].Adj(ie);
+							Triangle*        tt = ta;
+
+							//if the adjacent triangle is not a boundary triangle:
+							if (tt && tt->link){
+								Vertex &v = *ta.OppositeVertex();
+								R2     V = v;
+								Int4   iV = Number(v);
+								Real8  lA = bamg::Area2(V,B,C)/dd;
+								Real8  lB = bamg::Area2(A,V,C)/dd;
+								Real8  lC = bamg::Area2(A,B,V)/dd;
+								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) ;
 							}
-						 else
-							{
-							 nbb++;
-							 taa[0][j]=0;
-							 taa[1][j]=0;
-							 taa[2][j]=0;
-							 taa[j][j]=1;
-							 bb[j]=0;
+							else{
+								nbb++;
+								taa[0][j]=0;
+								taa[1][j]=0;
+								taa[2][j]=0;
+								taa[j][j]=1;
+								bb[j]=0;
 							}
 						}
 
-					 // resolution of 3x3 lineaire system transpose
-					 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"));
-					 }
-					 // computation of the gradient 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)
-					 Real8 Hxx = cAB * ( nBC.x*nCA.x) +  cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);
-					 Real8 Hyy = cAB * ( nBC.y*nCA.y) +  cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);
-					 Real8 Hxy = cAB * ( nBC.y*nCA.x) +  cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 
-						+ cAB * ( nBC.x*nCA.y) +  cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);
-					 Real8 coef = 1.0/(3*dd*det33);
-					 Real8 coef2 = 2*coef;
-					 Hxx *= coef2;
-					 Hyy *= coef2;
-					 Hxy *= coef2;
-					 if(nbb==0)
-						{
-						 dxdx[iA] += Hxx;
-						 dydy[iA] += Hyy;
-						 dxdy[iA] += Hxy;
-
-						 dxdx[iB] += Hxx;
-						 dydy[iB] += Hyy;
-						 dxdy[iB] += Hxy;
-
-						 dxdx[iC] += Hxx;
-						 dydy[iC] += Hyy;
-						 dxdy[iC] += Hxy;
+						// resolution of 3x3 linear system transpose
+						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"));
 						}
-
+						// 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)
+						Real8 Hxx = cAB * ( nBC.x*nCA.x) +  cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);
+						Real8 Hyy = cAB * ( nBC.y*nCA.y) +  cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);
+						Real8 Hxy = cAB * ( nBC.y*nCA.x) +  cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 
+						  + cAB * ( nBC.x*nCA.y) +  cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);
+						Real8 coef = 1.0/(3*dd*det33);
+						Real8 coef2 = 2*coef;
+						Hxx *= coef2;
+						Hyy *= coef2;
+						Hxy *= coef2;
+						if(nbb==0){
+							dxdx[iA] += Hxx;
+							dydy[iA] += Hyy;
+							dxdy[iA] += Hxy;
+
+							dxdx[iB] += Hxx;
+							dydy[iB] += Hyy;
+							dxdy[iB] += Hxy;
+
+							dxdx[iC] += Hxx;
+							dydy[iC] += Hyy;
+							dxdy[iC] += Hxy;
+						}
 					}
-				 else {
-
-					 // if edge on boundary no contribution  => normal = 0
-					 if ( ! tBC || ! tBC->link ) nBC = O;
-					 if ( ! tCA || ! tCA->link ) nCA = O;
-					 if ( ! tAB || ! tAB->link ) nAB = O;
-
-					 // remark we forgot a 1/2 because
-					 //       $\\int_{edge} w_i = 1/2 $ if $i$ is in edge 
-					 //                          0  if not
-					 // if we don't take the  boundary 
-					 // dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
-
-					 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
-					 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
-					 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
-
-					 // warning optimization (1) the divide by 2 is done on the metrix construction
-					 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
-					 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
-					 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 
-
-					 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
-					 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
-					 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
+
+					//Use Green to compute Hessian Matrix
+					else if (Hessiantype==0){
+
+						// if edge on boundary no contribution  => normal = 0
+						if ( !tBC || !tBC->link ) nBC=O;
+						if ( !tCA || !tCA->link ) nCA=O;
+						if ( !tAB || !tAB->link ) nAB=O;
+
+						// remark we forgot a 1/2 because
+						//       int_{edge} w_i = 1/2 if i is in edge 
+						//                         0  if not
+						// if we don't take the  boundary 
+						dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
+						dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
+						dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
+
+						//warning optimization (1) the division by 2 is done on the metric construction
+						dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
+						dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
+						dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 
+
+						dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
+						dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
+						dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
 					}
 
 				} // for real all triangles 
+			}
+			
 			Int4 kk=0;
-			for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )
-			 if(Mmassxx[iv]>0) 
-				{
-				 dxdx[iv] /= 2*Mmassxx[iv];
-				 // warning optimization (1) on term dxdy[iv]*ci/2 
-				 dxdy[iv] /= 4*Mmassxx[iv];
-				 dydy[iv] /= 2*Mmassxx[iv];
-				 // Compute the matrix with abs(eigen value)
-				 Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
-				 MatVVP2x2 Vp(M);
-				 Vp.Abs();
-				 M = Vp;
-				 dxdx[iv] = M.a11;
-				 dxdy[iv] = M.a21;
-				 dydy[iv] = M.a22;
+			for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
+				if(Mmassxx[iv]>0){
+					dxdx[iv] /= 2*Mmassxx[iv];
+					// warning optimization (1) on term dxdy[iv]*ci/2 
+					dxdy[iv] /= 4*Mmassxx[iv];
+					dydy[iv] /= 2*Mmassxx[iv];
+					// Compute the matrix with abs(eigen value)
+					Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
+					MatVVP2x2 Vp(M);
+					Vp.Abs();
+					M = Vp;
+					dxdx[iv] = M.a11;
+					dxdy[iv] = M.a21;
+					dydy[iv] = M.a22;
 				}
-			 else kk++;
+				else kk++;
+			}
 
 			// correction of second derivative
 			// by a laplacien
-			Real8 *d2[3] = { dxdx, dxdy, dydy};
-			Real8 *dd;
+			Real8* d2[3] = {dxdx, dxdy, dydy};
+			Real8* dd;
 			for (int xy = 0;xy<3;xy++) {
 				dd = d2[xy];
@@ -3070,6 +3121,5 @@
 				for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
 					for (i=0;i<nbt;i++) 
-					 if(triangles[i].link) // the real triangles 
-						{
+					 if(triangles[i].link){// the real triangles 
 						 // number of the 3 vertices
 						 iA = Number(triangles[i][0]);
@@ -3080,75 +3130,94 @@
 						  cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
 						 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
+					 }
+					for (iv=0;iv<nbv;iv++) workV[iv]=0;
+
+					for (i=0;i<nbt;i++){ 
+						if(triangles[i].link){ // the real triangles 
+							// number of the 3 vertices
+							iA = Number(triangles[i][0]);
+							iB = Number(triangles[i][1]);
+							iC = Number(triangles[i][2]);
+							Real8 cc =  workT[i]*detT[i];
+							workV[iA] += cc;
+							workV[iB] += cc;
+							workV[iC] += cc;
 						}
-					for (iv=0;iv<nbv;iv++)
-					 workV[iv]=0;
-
-					for (i=0;i<nbt;i++) 
-					 if(triangles[i].link) // the real triangles 
-						{
-						 // number of the 3 vertices
-						 iA = Number(triangles[i][0]);
-						 iB = Number(triangles[i][1]);
-						 iC = Number(triangles[i][2]);
-						 Real8 cc =  workT[i]*detT[i];
-						 workV[iA] += cc;
-						 workV[iB] += cc;
-						 workV[iC] += cc;
+					}
+
+					for (iv=0;iv<nbv;iv++){
+						if( ijacobi<NbJacobi || OnBoundary[iv]){
+							dd[iv] = workV[iv]/(Mmass[iv]*6);
 						}
-
-					for (iv=0;iv<nbv;iv++)
-					 if( ijacobi<NbJacobi || OnBoundary[iv])
-					  dd[iv] = workV[iv]/(Mmass[iv]*6);
+					}
 				}
 			}
 
-			// constuction  of the metrix from the Hessian dxdx. dxdy,dydy
+			//constuction of the metric from the Hessian dxdx. dxdy,dydy
 			Real8 rCutOff=CutOff*absmax;// relative cut off 
 
-			for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )
-			  { // for all vertices 
+			//loop over the nodes
+			for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
+
 				MetricIso Miso;
 				Real8 ci ;
-				if (RelativeMetric)
-				  { //   compute the norm of the solution
+
+				//   compute norm of the solution
+				if (RelativeMetric){
 					double xx =0,*sfk=sf+k; 
-					for (int ifield=0;ifield<nbfield;ifield++,sfk++)
-					 xx += *sfk* *sfk;	       
+					for (int ifield=0;ifield<nbfield;ifield++,sfk++){
+						xx += *sfk* *sfk;	       
+					}
 					xx=sqrt(xx);
-					ci = coef2/Max(xx,rCutOff);
-				  }
-				else ci = cnorm;
-
-				Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci,  dydy[iv]*ci);
+					ci=coef2/Max(xx,rCutOff);
+				}
+				else ci=cnorm;
+
+				//initialize metric Miv with ci*H
+				Metric    Miv(dxdx[iv]*ci, dxdy[iv]*ci,  dydy[iv]*ci);
+
+				//Get eigen values and vectors of Miv
 				MatVVP2x2 Vp(Miv);
 
+				//move eigen valuse to their absolute values
 				Vp.Abs();
+
+				//Allpy a power if requested by user
 				if(power!=1.0) Vp.pow(power);
 
+				//get minimum and maximum eigen values  
 				h1=Min(h1,Vp.lmin());
 				h2=Max(h2,Vp.lmax());
 
+				//modify eigen values according to hmin and hmax
 				Vp.Maxh(hmin);
 				Vp.Minh(hmax);
 
+				//multiply eigen values by coef
+				Vp.BoundAniso2(coef);
+
+				//rebuild Metric from Vp
+				Metric MVp(Vp);
+
+				//Apply Metric to vertex
+				vertices[iv].m.IntersectWith(MVp);
+
+				//info to be displayed
+				//rx = max(lmax/lmin) (anisotropy ratio)
 				rx = Max(rx,Vp.Aniso2());
-
-				Vp.BoundAniso2(coef);
-
 				hn1=Min(hn1,Vp.lmin());
 				hn2=Max(hn2,Vp.lmax());
 				rnx = Max(rnx,Vp.Aniso2());
-
-				Metric MVp(Vp);
-				vertices[iv].m.IntersectWith(MVp);
-				//vertices[iv].m=Vp; TEST: work like that
-			  }// for all vertices 
-			//vertices[nbv-1].m.Echo();
-			if (verbosity>-1) { 
+			}
+
+			//display info
+			if (verbosity>2) { 
+
+
 				printf("      Field %i of solution %i\n",nufield,nusol);
 				printf("         before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5));
 				printf("         after  bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5));
 			}
-		  } //  end of for all field
+		} //  end of for all field
 	}// end for all solution