Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2919)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2920)
@@ -889,4 +889,908 @@
 			printf("      output: Hmin = %g, Hmax = %g, factor of anisotropy max  = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5));
 		}
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::BuildMetric0 (Green formula){{{1*/
+	void Triangles::BuildMetric0(BamgOpts* bamgopts){
+		//  the array of solution s is store    
+		// sol0,sol1,...,soln    on vertex 0
+		//  sol0,sol1,...,soln   on vertex 1
+		//  etc.
+
+		/*Options*/
+		const int dim = 2;
+		int AbsError;
+		double* s=NULL;
+		Int4 nbsol;
+		int* typsols=NULL;
+		Real8 hmin1;
+		Real8 hmax1;
+		Real8 coef;
+		Real8 anisomax=1e30;
+		Real8 CutOff;
+		int NbJacobi;
+		int Rescaling;
+		double power;
+		int verbosity;
+
+		/*Recover options*/
+		verbosity=bamgopts->verbose;
+		AbsError=bamgopts->AbsError;   
+		CutOff=bamgopts->cutoff;
+		hmin1=bamgopts->hmin;
+		hmax1=bamgopts->hmax;
+		coef=bamgopts->coef;
+		NbJacobi=bamgopts->nbjacobi;
+		Rescaling=bamgopts->Rescaling; //do normalization
+		power=bamgopts->power;
+
+		/*process options*/
+		if (AbsError) CutOff=0.0;
+		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 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
+		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 = %g\n",anisomax,NbJacobi,power); 
+			if (RelativeMetric) printf("      RelativeErr with CutOff= %g\n",CutOff);
+			else printf("      Absolute error\n");
+		}
+
+		//initialize Mmass, OnBoundary and Massxx by zero
+		for (iv=0;iv<nbv;iv++){
+			Mmass[iv]=0;
+			OnBoundary[iv]=0;
+			Mmassxx[iv]=0;
+		}
+
+		//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
+				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++) {
+
+			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; 
+
+			//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{
+				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];
+					}
+					v = sqrt(v);
+					smin=Min(smin,v);
+					smax=Max(smax,v);
+				}
+			}
+			Real8 sdelta=smax-smin;
+			Real8 absmax=Max(Abs(smin),Abs(smax));
+			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);
+
+			//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;
+			}
+
+			//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++){
+				//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];
+
+						//Use Green to compute Hessian Matrix
+
+						// 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;
+					}
+					else kk++;
+				}
+
+				// correction of second derivative
+				// by a laplacien
+				Real8* d2[3] = {dxdx, dxdy, dydy};
+				Real8* dd;
+				for (int xy = 0;xy<3;xy++) {
+					dd = d2[xy];
+					// do leat 2 iteration for boundary problem
+					for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
+						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=3;
+							 if(ijacobi==0)
+							  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++){
+							if( ijacobi<NbJacobi || OnBoundary[iv]){
+								dd[iv] = workV[iv]/(Mmass[iv]*6);
+							}
+						}
+					}
+				}
+
+				//constuction of the metric from the Hessian dxdx. dxdy,dydy
+				Real8 rCutOff=CutOff*absmax;// relative cut off 
+
+				//loop over the nodes
+				for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
+
+					MetricIso Miso;
+					Real8 ci ;
+
+					//   compute norm of the solution
+					if (RelativeMetric){
+						double xx =0,*sfk=sf+k; 
+						for (int ifield=0;ifield<nbfield;ifield++,sfk++){
+							xx += *sfk* *sfk;	       
+						}
+						xx=sqrt(xx);
+						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());
+					hn1=Min(hn1,Vp.lmin());
+					hn2=Max(hn2,Vp.lmax());
+					rnx = Max(rnx,Vp.Aniso2());
+				}
+
+				//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 for all solution 
+
+		delete [] detT;
+		delete [] Mmass;
+		delete [] dxdx;
+		delete [] dxdy;
+		delete [] dydy;
+		delete []  workT;
+		delete [] workV;
+		delete [] Mmassxx;
+		delete []  OnBoundary;
+
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::BuildMetric1 (from P2 on 4T){{{1*/
+	void Triangles::BuildMetric1(BamgOpts* bamgopts){
+		//  the array of solution s is store    
+		// sol0,sol1,...,soln    on vertex 0
+		//  sol0,sol1,...,soln   on vertex 1
+		//  etc.
+
+		/*Options*/
+		const int dim = 2;
+		int AbsError;
+		double* s=NULL;
+		Int4 nbsol;
+		int* typsols=NULL;
+		Real8 hmin1;
+		Real8 hmax1;
+		Real8 coef;
+		Real8 anisomax=1e30;
+		Real8 CutOff;
+		int NbJacobi;
+		int Rescaling;
+		double power;
+		int verbosity;
+
+		/*Recover options*/
+		verbosity=bamgopts->verbose;
+		AbsError=bamgopts->AbsError;   
+		CutOff=bamgopts->cutoff;
+		hmin1=bamgopts->hmin;
+		hmax1=bamgopts->hmax;
+		coef=bamgopts->coef;
+		NbJacobi=bamgopts->nbjacobi;
+		Rescaling=bamgopts->Rescaling; //do normalization
+		power=bamgopts->power;
+
+		/*process options*/
+		if (AbsError) CutOff=0.0;
+		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 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
+		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 = %g\n",anisomax,NbJacobi,power); 
+			if (RelativeMetric) printf("      RelativeErr with CutOff= %g\n",CutOff);
+			else printf("      Absolute error\n");
+		}
+
+		//initialize Mmass, OnBoundary and Massxx by zero
+		for (iv=0;iv<nbv;iv++){
+			Mmass[iv]=0;
+			OnBoundary[iv]=0;
+			Mmassxx[iv]=0;
+		}
+
+		//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){
+					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++) {
+
+			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; 
+
+			//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{
+				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];
+					}
+					v = sqrt(v);
+					smin=Min(smin,v);
+					smax=Max(smax,v);
+				}
+			}
+			Real8 sdelta=smax-smin;
+			Real8 absmax=Max(Abs(smin),Abs(smax));
+			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);
+
+			//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;
+			}
+
+			//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++){
+				//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
+
+						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;
+							}
+						}
+
+						// 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;
+						}
+					} // 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;
+					}
+					else kk++;
+				}
+
+				// correction of second derivative
+				// by a laplacien
+				Real8* d2[3] = {dxdx, dxdy, dydy};
+				Real8* dd;
+				for (int xy = 0;xy<3;xy++) {
+					dd = d2[xy];
+					// do leat 2 iteration for boundary problem
+					for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
+						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=3;
+							 if(ijacobi==0)
+							  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++){
+							if( ijacobi<NbJacobi || OnBoundary[iv]){
+								dd[iv] = workV[iv]/(Mmass[iv]*6);
+							}
+						}
+					}
+				}
+
+				//constuction of the metric from the Hessian dxdx. dxdy,dydy
+				Real8 rCutOff=CutOff*absmax;// relative cut off 
+
+				//loop over the nodes
+				for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
+
+					MetricIso Miso;
+					Real8 ci ;
+
+					//   compute norm of the solution
+					if (RelativeMetric){
+						double xx =0,*sfk=sf+k; 
+						for (int ifield=0;ifield<nbfield;ifield++,sfk++){
+							xx += *sfk* *sfk;	       
+						}
+						xx=sqrt(xx);
+						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());
+					hn1=Min(hn1,Vp.lmin());
+					hn2=Max(hn2,Vp.lmax());
+					rnx = Max(rnx,Vp.Aniso2());
+				}
+
+				//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 for all solution 
+
+		delete [] detT;
+		delete [] Mmass;
+		delete [] dxdx;
+		delete [] dxdy;
+		delete [] dydy;
+		delete []  workT;
+		delete [] workV;
+		delete [] Mmassxx;
+		delete []  OnBoundary;
+
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::BuildMetric2{{{1*/
+	void Triangles::BuildMetric2(BamgOpts* bamgopts){
+
+		throw ErrorException(__FUNCT__,exprintf("not supported yet"));
 	}
 	/*}}}1*/
@@ -2738,498 +3642,15 @@
 
 	/*Options*/
-	const int dim = 2;
-	int AbsError;
-	double* s=NULL;
-	Int4 nbsol;
-	int* typsols=NULL;
-	Real8 hmin1;
-	Real8 hmax1;
-	Real8 coef;
-	Real8 anisomax=1e30;
-	Real8 CutOff;
-	int NbJacobi;
-	int Rescaling;
-	double power;
-	int Hessiantype;
-	int verbosity;
-
-	/*Recover options*/
-	verbosity=bamgopts->verbose;
-	AbsError=bamgopts->AbsError;   
-	CutOff=bamgopts->cutoff;
-	hmin1=bamgopts->hmin;
-	hmax1=bamgopts->hmax;
-	coef=bamgopts->coef;
-	NbJacobi=bamgopts->nbjacobi;
-	Rescaling=bamgopts->Rescaling; //do normalization
-	power=bamgopts->power;
-	Hessiantype=bamgopts->Hessiantype;
-
-	/*process options*/
-	if (AbsError) CutOff=0.0;
-	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 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]];
+	int Hessiantype=bamgopts->Hessiantype;
+
+	if (Hessiantype==0){
+		BuildMetric0(bamgopts);
 	}
-	else ntmp = nbsol;
-
-	//n is the total number of fields
-	const Int4 n = ntmp;
-
-	//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 = %g\n",anisomax,NbJacobi,power); 
-		if (RelativeMetric) printf("      RelativeErr with CutOff= %g\n",CutOff);
-		else printf("      Absolute error\n");
+	else if (Hessiantype==1){
+		BuildMetric1(bamgopts);
 	}
-
-	//initialize Mmass, OnBoundary and Massxx by zero
-	for (iv=0;iv<nbv;iv++){
-		Mmass[iv]=0;
-		OnBoundary[iv]=0;
-		Mmassxx[iv]=0;
+	else{
+		throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (0->use Green formula, 1-> from P2 on 4T, 2-> double P2 projection)",Hessiantype));
 	}
-
-	//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++) {
-
-		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; 
-
-		//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{
-			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];
-				}
-				v = sqrt(v);
-				smin=Min(smin,v);
-				smax=Max(smax,v);
-			}
-		}
-		Real8 sdelta=smax-smin;
-		Real8 absmax=Max(Abs(smin),Abs(smax));
-		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);
-
-		//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;
-		}
-
-		//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++){
-			//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;
-							}
-						}
-
-						// 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;
-						}
-					}
-
-					//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;
-				}
-				else kk++;
-			}
-
-			// correction of second derivative
-			// by a laplacien
-			Real8* d2[3] = {dxdx, dxdy, dydy};
-			Real8* dd;
-			for (int xy = 0;xy<3;xy++) {
-				dd = d2[xy];
-				// do leat 2 iteration for boundary problem
-				for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
-					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=3;
-						 if(ijacobi==0)
-						  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++){
-						if( ijacobi<NbJacobi || OnBoundary[iv]){
-							dd[iv] = workV[iv]/(Mmass[iv]*6);
-						}
-					}
-				}
-			}
-
-			//constuction of the metric from the Hessian dxdx. dxdy,dydy
-			Real8 rCutOff=CutOff*absmax;// relative cut off 
-
-			//loop over the nodes
-			for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
-
-				MetricIso Miso;
-				Real8 ci ;
-
-				//   compute norm of the solution
-				if (RelativeMetric){
-					double xx =0,*sfk=sf+k; 
-					for (int ifield=0;ifield<nbfield;ifield++,sfk++){
-						xx += *sfk* *sfk;	       
-					}
-					xx=sqrt(xx);
-					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());
-				hn1=Min(hn1,Vp.lmin());
-				hn2=Max(hn2,Vp.lmax());
-				rnx = Max(rnx,Vp.Aniso2());
-			}
-
-			//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 for all solution 
-
-	delete [] detT;
-	delete [] Mmass;
-	delete [] dxdx;
-	delete [] dxdy;
-	delete [] dydy;
-	delete []  workT;
-	delete [] workV;
-	delete [] Mmassxx;
-	delete []  OnBoundary;
-
 }
 /*}}}1*/