Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2936)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2937)
@@ -71,4 +71,5 @@
 		throw ErrorException(__FUNCT__,exprintf("maxsubdiv %g is should be in ]1,%g]",maxsubdiv,boundmaxsubdiv));
 	}
+	// if iso, then anisomax = 1;
 	if (iso) anisomax=1;
 	// no metric -> no smoothing 
@@ -76,4 +77,6 @@
 		NbSmooth=0; 
 	}
+
+	/*If no mesh in input, generate one*/
 
 	if(bamgmesh->NumTriangles==0){
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2936)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2937)
@@ -137,4 +137,5 @@
 		//  ~Vertex(){}
 		Real8  Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,Real8 =1);
+		void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,BamgOpts* bamgopts);
 		int ref() const { return ReferenceNumber;}
 
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2936)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2937)
@@ -900,33 +900,13 @@
 		/*Options*/
 		const int dim = 2;
-		int AbsError;
 		double* s=NULL;
 		Int4 nbsol;
 		int* typsols=NULL;
-		Real8 hmin1;
-		Real8 hmax1;
-		Real8 coef;
-		Real8 anisomax;
-		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;
-		anisomax=bamgopts->anisomax;
-
-		/*process options*/
-		if (AbsError) CutOff=0.0;
-		coef=sqrt(bamgopts->err)*coef;
 
 		/*Get and process fields*/
@@ -955,8 +935,4 @@
 		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;
@@ -974,9 +950,4 @@
 		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");
 		}
 
@@ -1049,5 +1020,4 @@
 			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; 
@@ -1078,8 +1048,7 @@
 			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);
+			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
 
 			//skip constant field
@@ -1236,69 +1205,9 @@
 				}
 
-				//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(hmax);
-					Vp.Minh(hmin);
-
-					//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));
-				}
+				/*Compute Metric from Hessian*/
+				for ( iv=0;iv<nbv;iv++){
+					vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts);
+				}
+
 			} //  end of for all field
 		}// end for all solution 
@@ -1325,33 +1234,13 @@
 		/*Options*/
 		const int dim = 2;
-		int AbsError;
 		double* s=NULL;
 		Int4 nbsol;
 		int* typsols=NULL;
-		Real8 hmin1;
-		Real8 hmax1;
-		Real8 coef;
-		Real8 anisomax;
-		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;
-		anisomax=bamgopts->anisomax;
-
-		/*process options*/
-		if (AbsError) CutOff=0.0;
-		coef=sqrt(bamgopts->err)*coef;
 
 		/*Get and process fields*/
@@ -1380,8 +1269,4 @@
 		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;
@@ -1399,9 +1284,4 @@
 		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");
 		}
 
@@ -1476,5 +1356,4 @@
 			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; 
@@ -1505,8 +1384,7 @@
 			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);
+			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
 
 			//skip constant field
@@ -1657,9 +1535,4 @@
 						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;
@@ -1739,69 +1612,9 @@
 				}
 
-				//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(hmax);
-					Vp.Minh(hmin);
-
-					//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));
-				}
+				/*Compute Metric from Hessian*/
+				for ( iv=0;iv<nbv;iv++){
+					vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts);
+				}
+
 			} //  end of for all field
 		}// end for all solution 
@@ -1828,24 +1641,10 @@
 		int*   typsols=NULL;
 		int    verbosity;
-		Real8  hmin1;
-		Real8  hmax1;
-		double power;
-		double anisomax;
 
 		int   i,j,k,iA,iB,iC;
 		int   iv,nbfield;
-		double coef;
 
 		/*Recover options*/
 		verbosity=bamgopts->verbose;
-		power=bamgopts->power;
-		hmin1=bamgopts->hmin;
-		hmax1=bamgopts->hmax;
-
-		/*Recover options*/
-		verbosity=bamgopts->verbose;
-		hmin1=bamgopts->hmin;
-		hmax1=bamgopts->hmax;
-		anisomax=bamgopts->anisomax;
 
 		/*Get and process fields*/
@@ -1872,6 +1671,4 @@
 
 		//initialization of some variables
-		Real8   hmin = Max(hmin1,MinimalHmin());
-		Real8   hmax = Min(hmax1,MaximalHmax());
 		double* ss=(double*)s;
 		double  sA,sB,sC;
@@ -1953,4 +1750,40 @@
 		for (Int4 nusol=0;nusol<nbsol;nusol++) {
 			int   nbfield=typsols?sizeoftype[typsols[nusol]]:1; 
+			Real8 smin=ss[0],smax=ss[0];
+
+			//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 ){
+					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));
+
+			//display info
+			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,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;
+			}
+
 
 			//loop over all the fields of the solution
@@ -2017,35 +1850,8 @@
 
 				/*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);
-
-					//move eigen valuse to their absolute values
-					Vp.Abs();
-
-					//Apply a power if requested by user
-					if(power!=1.0) Vp.pow(power);
-
-					//modify eigen values according to hmin and hmax
-					Vp.Maxh(hmax);
-					Vp.Minh(hmin);
-
-					//Bound anisotropy by 1/(anisomax)^2
-					Vp.BoundAniso2(1/(anisomax*anisomax));
-
-					//rebuild Metric from Vp
-					Metric MVp(Vp);
-
-					//Apply Metric to vertex
-					vertices[iv].m.IntersectWith(MVp);
-				}
+					vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[n*iv],bamgopts);
+				}
+
 			}//for all fields
 		}//for all solutions
Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2936)
+++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2937)
@@ -127,4 +127,89 @@
 	}
 	/*}}}1*/
+	/*FUNCTION Vertex::MetricFromHessian{{{1*/
+	void Vertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,BamgOpts* bamgopts){
+		/*Compute Metric from Hessian*/
+
+		/*get options*/
+		double power=(bamgopts->power);
+		double anisomax=(bamgopts->anisomax);
+		double CutOff=bamgopts->cutoff;
+		double hmin=(bamgopts->hmin);
+		double hmax=(bamgopts->hmax);
+		double coef=bamgopts->coef;
+		int    Metrictype=(bamgopts->Metrictype);
+		int    Rescaling=(bamgopts->Rescaling);
+
+		/*Intermediary*/
+		double ci;
+
+		/*compute multiplicative coefficient depending on Metric Type (2/9 because it is 2d)*/
+
+		//Absolute Error
+		/*
+		 *            2         1       
+		 *Metric M = ---  ------------   Abs(Hessian)
+		 *            9   err * coeff^2  
+		 */
+		if (Metrictype==0){
+			if (Rescaling){
+				ci= 1/(bamgopts->err*coef*coef) * 1/(smax-smin);
+			}
+			else{
+				ci=1/(bamgopts->err*coef*coef);
+			}
+		}
+
+		//Relative Error
+		/*
+		 *            2         1            Abs(Hessian)
+		 *Metric M = ---  ------------  ----------------------
+		 *            9   err * coeff^2  max( |s| , cutoff*max(|s|) )
+		 *
+		 */
+		else if (Metrictype==1){
+			ci=1/(bamgopts->err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
+		}
+
+		//Interpolation Error
+		/*
+		 *            2         1       
+		 *Metric M = ---  ------------   Abs(Hessian)
+		 *            9   err * coeff^2  
+		 */
+		else if (Metrictype==2){
+			ci=(2.0/9.0)*(1.0/(bamgopts->err));
+		}
+		else{
+			throw ErrorException(__FUNCT__,exprintf("Metrictype %i not supported yet (use 0,1 or 2(default))",Metrictype));
+		}
+
+		//initialize metric Miv with ci*H
+		Metric Miv(Hxx*ci,Hyx*ci,Hyy*ci);
+
+		//Get eigen values and vectors of Miv
+		MatVVP2x2 Vp(Miv);
+
+		//move eigen valuse to their absolute values
+		Vp.Abs();
+
+		//Apply a power if requested by user
+		if(power!=1.0) Vp.pow(power);
+
+		//modify eigen values according to hmin and hmax
+		Vp.Maxh(hmax);
+		Vp.Minh(hmin);
+
+		//Bound anisotropy by 1/(anisomax)^2
+		Vp.BoundAniso2(1/(anisomax*anisomax));
+
+		//rebuild Metric from Vp
+		Metric MVp(Vp);
+
+		//Apply Metric to vertex
+		m.IntersectWith(MVp);
+
+	}
+	/*}}}1*/
 
 	/*Intermediary*/
Index: /issm/trunk/src/c/objects/BamgOpts.h
===================================================================
--- /issm/trunk/src/c/objects/BamgOpts.h	(revision 2936)
+++ /issm/trunk/src/c/objects/BamgOpts.h	(revision 2937)
@@ -16,4 +16,5 @@
 	int     NbSmooth;
 	int     Rescaling;
+	int     Metrictype;
 	int     nbjacobi;
 	int     AbsError;
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 2936)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 2937)
@@ -84,5 +84,6 @@
 bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1);
 bamg_options.AbsError=getfieldvalue(options,'AbsError',0);
-bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
+bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',2);
+bamg_options.Metrictype=getfieldvalue(options,'Metrictype',2);
 bamg_options.NbSmooth=getfieldvalue(options,'NbSmooth',3);
 bamg_options.omega=getfieldvalue(options,'omega',1.8);
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2936)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2937)
@@ -36,5 +36,5 @@
 	double err,errg,coef;
 	double power;
-	int    Hessiantype,NbSmooth;
+	int    Hessiantype,Metrictype,NbSmooth;
 	int    Rescaling,nbjacobi,AbsError;
 	double omega;
@@ -117,4 +117,6 @@
 	FetchData(&Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype"));
 	bamgopts.Hessiantype=Hessiantype;
+	FetchData(&Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype"));
+	bamgopts.Metrictype=Metrictype;
 	FetchData(&power,mxGetField(BAMGOPTIONS,0,"power"));
 	bamgopts.power=power;