Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2899)
@@ -52,4 +52,8 @@
 	double omega=1.8;
 	int NbSmooth=3;
+	int ChoiseHessien = 0;
+	double power=1;
+	int Rescaling=1;
+
 	Triangles* Thr=NULL;
 	Triangles* Thb=NULL;
@@ -79,4 +83,9 @@
 		if (verbosity>1) printf("   Processing geometry...\n");
 		Geometry Gh(bamggeom,bamgopts);
+		if (verbosity>10){
+			Gh.Echo();
+		}
+
+		//get hmin and hmax from geometry to generate the metric
 		hmin = Max(hmin,Gh.MinimalHmin());
 		hmax = Min(hmax,Gh.MaximalHmax());
@@ -118,16 +127,19 @@
 		hmax = Min(hmax,BTh.MaximalHmax());
 
-		//?????TEST
+		//Make Quadtree from background mesh
 		BTh.MakeQuadTree();
 
 		//build metric if not given in input
-		if (verbosity>1) printf("   Processing Metric...\n");
 		if (bamgopts->metric){
+			if (verbosity>1) printf("   Processing Metric...\n");
 			BTh.ReadMetric(bamgopts,hmin,hmax,coef);
+		}
+		else if (bamgopts->field){
+			if (verbosity>1) printf("   Generating Metric from solution field...\n");
+			BTh.IntersectConsMetric(bamgopts->field,1,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:cutoff,nbjacoby,Rescaling,power,ChoiseHessien);
 		}
 		else { // init with hmax 
 			Metric Mhmax(hmax);
-			for (Int4 iv=0;iv<BTh.nbv;iv++)
-			 BTh[iv].m = Mhmax;
+			for (Int4 iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;
 		}
 
@@ -138,4 +150,5 @@
 
 		//Build new mesh Thr
+		if (verbosity>1) printf("   Generating Mesh...\n");
 		Thr=&BTh,Thb=0;
 		Triangles & Th( *(0 ?  new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) :  new Triangles(maxnbv,BTh,KeepBackVertices)));
@@ -163,4 +176,6 @@
 		if (verbosity>1) printf("   Write Geometry...\n");
 		Th.Gh.WriteGeometry(bamggeom,bamgopts);
+		if (verbosity>1) printf("   Write Metric...\n");
+		BTh.WriteMetric(bamgopts);
 
 		/*clean up*/
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2899)
@@ -845,6 +845,5 @@
 			Real8 MaximalAngleOfCorner;
 
-			//  end of data 
-
+			void Echo();
 
 			I2 toI2(const R2 & P) const {
Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2899)
@@ -73,4 +73,5 @@
 		inline void  Box(Real4 &hx,Real4 &hy) const ;  
 		MetricAnIso(const MatVVP2x2);
+		void  Echo();
 	};
 
@@ -82,4 +83,5 @@
 		D2 v;
 
+		void  Echo();
 		MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
 
@@ -103,5 +105,5 @@
 		void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
 		void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
-	  };
+	};
 
 	inline void  MatVVP2x2::BoundAniso2(const Real8 coef){
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2899)
@@ -580,35 +580,56 @@
 			}
 		}
-
-		// bulle sort on the angle of edge  
+		// in our case, for 4 points and 4 edges (in parenthesis from now on):
+		// hv =  7 2 4 6 
+		// ev =  -1 -1 1 -1 3 -1 5 0 
+
+		//compute angle of edges
 		for (i=0;i<nbv;i++) {
+
 			int exch=1, ord=0;      
 			while (exch){
-				exch  =0;
-				Int4* p=hv+i;
-				Int4* po=p;
-				Int4  n=*p;
+				Int4 *p=hv+i;                   // pointer to hv[vertex i]  (p=hv)
+				Int4 *po=p;                     // copy of pointer p        (po=p)
+				Int4  n=*p;                     // value hv[vertex i]       (n=7)
 				register float angleold=-1000 ; // angle = - infinity
-				ord = 0;
-				while (n >=0){
+				ord=0; exch=0;
+				while (n >=0){// loop as long as ...
 					ord++;
-					register Int4 i1=n/2;
-					register Int4 j1=n%2;
-					register Int4 *pn=ev+n;
+					register Int4  i1=n/2;       // i1 = floor (n/2)          (i1 = 7/2 = 3)
+					register Int4  j1=n%2;       // j1 = 1 if n is odd        (j1 = 1)
+					register Int4* pn=ev+n;      // pointer to ev[n]          (pn = &ev[7])
+					//compute angle from eangle:
+					// if n odd: angle = eangle - Pi [2*Pi]
+					// else    : angle = eangle
 					float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1];
-					n = *pn;
-					if (angleold > angle) // exch to have : po -> pn -> p 
-					 exch=1,*pn = *po,*po=*p,*p=n,po = pn;
-					else //  to have : po -> p -> pn 
-					 angleold =  angle, po = p,p  = pn;
-				}
-			}
-
-			// angulare test to find a corner 
+					n=*pn;                       //  n = ev[n]                (n = ev[7] = 0)
+					if (angleold > angle){       // exch to have : po -> pn -> p 
+						exch=1, *pn=*po, *po=*p, *p=n, po=pn;
+					}
+					else{                        //  to have : po -> p -> pn
+						angleold=angle, po=p, p=pn;
+					}
+				}
+			}
+			//first iteration: i=0
+			//hv 0 2 4 6 
+			//ev 7 -1 1 -1 3 -1 5 -1 
+
+			// angular test on current vertex to guess whether it is a corner
 			if(ord == 2) { 
-				Int4 n1 = hv[i];
-				Int4 n2 = ev[n1];
-				Int4 i1 = n1 /2, i2 = n2/2; // edge number
-				Int4  j1 = n1 %2, j2 = n2%2; // vertex in the edge 
+				Int4  n1 = hv[i];
+				Int4  n2 = ev[n1];
+				Int4  i1 = n1/2, i2 = n2/2; // edge number
+				Int4  j1 = n1%2, j2 = n2%2; // vertex in the edge 
+				/*//first iteration i=0
+				//n1=hv[i]=0
+				//n2=ev[i]=7
+				//i1 = 0  i2 = 3
+				//j1 = 0  j2 = 1
+				printf("i=%i\n",i);
+				printf("n1=hv[i]=%i\n",n1);
+				printf("n2=ev[i]=%i\n",n2);
+				printf("i1 = %i  i2 = %i\n",i1,i2);
+				printf("j1 = %i  j2 = %i\n",j1,j2);*/
 				float angle1=  j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
 				float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
@@ -638,12 +659,12 @@
 		k =0;
 		for (i=0;i<nbe;i++){
-			for (jj=0;jj<2;jj++){
+			for (j=0;j<2;j++){
 				Int4 n1 = ev[k++]; 
 				Int4 i1 = n1/2 ,j1=n1%2;
-				if( edges[i1].v[j1] != edges[i].v[jj]) {
+				if( edges[i1].v[j1] != edges[i].v[j]) {
 					throw ErrorException(__FUNCT__,exprintf("Bug Adj edge"));
 				}
 				edges[i1].Adj[j1] = edges + i;
-				edges[i1].DirAdj[j1] = jj;
+				edges[i1].DirAdj[j1] = j;
 			}
 		}
@@ -776,4 +797,34 @@
 	}
 	/*}}}1*/
+	/*FUNCTION  Geometry::Echo {{{1*/
+
+	void Geometry::Echo(void){
+
+		printf("Geometry:\n");
+		printf("   name: %s\n",name);
+		printf("   nbvx (maximum number of vertices) : %i\n",nbvx);
+		printf("   nbtx (maximum number of triangles): %i\n",nbtx);
+		printf("   nbv  (number of vertices) : %i\n",nbv);
+		printf("   nbt  (number of triangles): %i\n",nbt);
+		printf("   nbe  (number of edges)    : %i\n",nbe);
+		printf("   nbv  (number of initial vertices) : %i\n",nbiv);
+		printf("   NbSubDomains: %i\n",NbSubDomains);
+		printf("   NbEquiEdges: %i\n",NbEquiEdges);
+		printf("   NbCrackedEdges: %i\n",NbCrackedEdges);
+		printf("   NbOfCurves: %i\n",NbOfCurves);
+		printf("   vertices: %p\n",vertices);
+		printf("   triangles: %p\n",triangles);
+		printf("   edges: %p\n",edges);
+		printf("   quadtree: %p\n",quadtree);
+		printf("   subdomains: %p\n",subdomains);
+		printf("   curves: %p\n",curves);
+		printf("   pmin (x,y): (%g %g)\n",pmin.x,pmin.y);
+		printf("   pmax (x,y): (%g %g)\n",pmax.x,pmax.y);
+		printf("   coefIcoor: %g\n",coefIcoor);
+		printf("   MaximalAngleOfCorner: %g\n",MaximalAngleOfCorner);
+
+		return;
+	}
+	/*}}}*/
 	/*FUNCTION  Geometry::EmptyGeometry(){{{1*/
 	void Geometry::EmptyGeometry() {
Index: /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 2899)
@@ -47,3 +47,18 @@
 	/*}}}1*/
 
+	/*Methods*/
+	/*FUNCTION MatVVP2x2::Echo {{{1*/
+
+	void MatVVP2x2::Echo(void){
+
+		printf("MatVVP2x2:\n");
+		printf("   lambda1: %g\n",lambda1);
+		printf("   lambda2: %g\n",lambda2);
+		printf("   v.x: %g\n",v.x);
+		printf("   v.y: %g\n",v.y);
+
+		return;
+	}
+	/*}}}*/
+
 } 
Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2899)
@@ -59,4 +59,14 @@
 
 	/*Methods*/
+	/*FUNCTION MetricAnIso::Echo {{{1*/
+
+	void MetricAnIso::Echo(void){
+
+		printf("MetricAnIso:\n");
+		printf("   [a11 a21 a22]: [%g %g %g]\n",a11,a21,a22);
+
+		return;
+	}
+	/*}}}*/
 	/*FUNCTION MetricAnIso::IntersectWith{{{1*/
 	int MetricAnIso::IntersectWith(const MetricAnIso M2) {
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2898)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2899)
@@ -325,5 +325,7 @@
 			hvertices=1;
 			for (i=0;i< nbv;i++){
-				vertices[i].m=Metric((Real4)bamgmesh->hVertices[i]);
+				if (!isnan(bamgmesh->hVertices[i])){
+					vertices[i].m=Metric((Real4)bamgmesh->hVertices[i]);
+				}
 			}
 		}
@@ -689,14 +691,17 @@
 			}
 			else if (j==3){
-				Real8 a,b,c;	     
-				a=bamgopts->metric[i*3+0];
-				b=bamgopts->metric[i*3+1];
-				c=bamgopts->metric[i*3+2];
-				MetricAnIso M(a,b,c);
-				MatVVP2x2 Vp(M/coef);
-
-				Vp.Maxh(hmin);
-				Vp.Minh(hmax);
-				vertices[i].m = Vp;
+				//do not erase metric computed by hVertices
+				if (vertices[i].m.a11==1 && vertices[i].m.a21==0 && vertices[i].m.a22==1){
+					Real8 a,b,c;	     
+					a=bamgopts->metric[i*3+0];
+					b=bamgopts->metric[i*3+1];
+					c=bamgopts->metric[i*3+2];
+					MetricAnIso M(a,b,c);
+					MatVVP2x2 Vp(M/coef);
+
+					Vp.Maxh(hmin);
+					Vp.Minh(hmax);
+					vertices[i].m = Vp;
+				}
 			}
 		}
@@ -2723,16 +2728,16 @@
 	/*}}}1*/
 /*FUNCTION Triangles::IntersectConsMetric{{{1*/
-void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols,
+void Triangles::IntersectConsMetric(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 DoNormalisation,const double power,const int choice)
-{ //  the array of solution s is store    
+			const int DoNormalisation,const double power,const int choice){
+	//  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 formule
-	//   otherwise  => H is computed from P2 on 4T 
+	//  choise = 0 =>  H is computed with green formula
+	//  otherwise  =>  H is computed from P2 on 4T 
+
 	const int dim = 2;
-
 	long int verbosity=0;
 
@@ -2741,11 +2746,8 @@
 	// computation of the nb of field 
 	Int4 ntmp = 0;
-	if (typsols)
-	  {
-		for (Int4 i=0;i<nbsol;i++)
-		 ntmp += sizeoftype[typsols[i]];
+	if (typsols){
+		for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
 	  }
-	else
-	 ntmp = nbsol;
+	else ntmp = nbsol;
 
 	// n is the total number of fields
@@ -2768,5 +2770,5 @@
 		else printf("      Absolute error\n");
 	}
-	double *ss=(double*)s;
+	double* ss=(double*)s;
 
 	double sA,sB,sC;
@@ -2781,14 +2783,12 @@
 	Real8 *workV= new Real8[nbv];
 	int *OnBoundary = new int[nbv];
-	for (iv=0;iv<nbv;iv++)
-	  {
+	for (iv=0;iv<nbv;iv++){
 		Mmass[iv]=0;
 		OnBoundary[iv]=0;
 		Mmassxx[iv]=0;
-	  }
+	}
 
 	for (i=0;i<nbt;i++) 
-	 if(triangles[i].link) // the real triangles 
-		{
+	 if(triangles[i].link){ // the real triangles 
 		 const Triangle &t=triangles[i];
 		 // coor of 3 vertices 
@@ -2823,6 +2823,5 @@
 		 Mmass[iC] += dett;
 
-		 if((nbb==0)|| !choice)
-			{
+		 if((nbb==0)|| !choice){
 			 Mmassxx[iA] += dett;
 			 Mmassxx[iB] += dett;
@@ -2830,9 +2829,8 @@
 			}
 		}
-	 else
-	  workT[i]=-1;
-
-	for (Int4 nusol=0;nusol<nbsol;nusol++)
-	  { //for all Solution  
+	 else workT[i]=-1;
+
+	for (Int4 nusol=0;nusol<nbsol;nusol++) {
+		//for all Solution  
 
 		Real8 smin=ss[0],smax=ss[0];
@@ -2843,15 +2841,12 @@
 		int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 
 		if (nbfield == 1) 
-		 for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-			{
+		 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]);
-			}
-		else
-		  {
+		 }
+		else{
 			//  cas vectoriel 
-			for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-			  {	
+			for ( iv=0,k=0; iv<nbv; iv++,k+=n ){	
 				double v=0;		     
 				for (int i=0;i<nbfield;i++) 
@@ -2860,6 +2855,6 @@
 				smin=Min(smin,v);
 				smax=Max(smax,v);
-			  }
-		  }
+			}
+		}
 		Real8 sdelta = smax-smin;
 		Real8 absmax=Max(Abs(smin),Abs(smax));
@@ -2874,6 +2869,5 @@
 
 		double *sf  = ss; 
-		for (Int4 nufield=0;nufield<nbfield;nufield++,ss++) 
-		  {
+		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;
@@ -3035,16 +3029,12 @@
 			 else kk++;
 
-
-			// correction of second derivate
+			// correction of second derivative
 			// by a laplacien
-
 			Real8 *d2[3] = { dxdx, dxdy, dydy};
 			Real8 *dd;
-			for (int xy = 0;xy<3;xy++)
-			  {
+			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 (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
 					for (i=0;i<nbt;i++) 
 					 if(triangles[i].link) // the real triangles 
@@ -3078,13 +3068,8 @@
 					 if( ijacobi<NbJacobi || OnBoundary[iv])
 					  dd[iv] = workV[iv]/(Mmass[iv]*6);
-
-
-				  }
-
-
-			  }
+				}
+			}
 
 			// constuction  of the metrix from the Hessian dxdx. dxdy,dydy
-
 			Real8 rCutOff=CutOff*absmax;// relative cut off 
 
@@ -3107,8 +3092,5 @@
 
 				Vp.Abs();
-				if(power!=1.0) 
-				 Vp.pow(power);
-
-
+				if(power!=1.0) Vp.pow(power);
 
 				h1=Min(h1,Vp.lmin());
@@ -3128,6 +3110,8 @@
 				Metric MVp(Vp);
 				vertices[iv].m.IntersectWith(MVp);
+				//vertices[iv].m=Vp; TEST: work like that
 			  }// for all vertices 
-			if (verbosity>2) { 
+			//vertices[nbv-1].m.Echo();
+			if (verbosity>-1) { 
 				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));
@@ -3135,5 +3119,5 @@
 			}
 		  } //  end of for all field
-	  }// end for all solution 
+	}// end for all solution 
 
 	delete [] detT;
Index: /issm/trunk/src/c/objects/BamgOpts.h
===================================================================
--- /issm/trunk/src/c/objects/BamgOpts.h	(revision 2898)
+++ /issm/trunk/src/c/objects/BamgOpts.h	(revision 2899)
@@ -18,4 +18,5 @@
 	int     verbose;
 	double* metric;
+	double* field;
 
 };
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 2898)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 2899)
@@ -35,5 +35,4 @@
 			bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
 			bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  ones(nods,1)]];
-			bamg_geometry.hVertices=[];%[bamg_geometry.hVertices; resolution*ones(nods,1)];
 			if i>1,
 				clockwise=-1;
@@ -41,4 +40,10 @@
 			end
 			count=count+nods;
+		end
+		if exist(options,'hVertices'),
+			bamg_geometry.hVertices=getfieldvalueerr(options,'hVertices');
+			if length(bamg_geometry.hVertices)~=nods,
+				error(['hVertices option should be of size [' num2str(nods) ',1]']);
+			end
 		end
 		bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
@@ -55,4 +60,5 @@
 bamg_mesh.Triangles=[];
 bamg_mesh.index=zeros(0,3);
+bamg_mesh.hVertices=[];
 if (~exist(options,'domain') & md.numberofgrids~=0 & strcmpi(md.type,'2d')),
 	bamg_mesh.NumVertices=md.numberofgrids;
@@ -60,4 +66,10 @@
 	bamg_mesh.NumTriangles=md.numberofelements;
 	bamg_mesh.Triangles=[md.elements ones(md.numberofelements,1)];
+	if exist(options,'hVertices'),
+		bamg_mesh.hVertices=getfieldvalueerr(options,'hVertices');
+		if length(bamg_mesh.hVertices)~=md.numberofgrids,
+			error(['hVertices option should be of size [' num2str(md.numberofgrids) ',1]']);
+		end
+	end
 end
 %}}}
@@ -74,8 +86,9 @@
 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
 bamg_options.metric=getfieldvalue(options,'metric',zeros(0,3));
+bamg_options.field=getfieldvalue(options,'field',zeros(0,1));
 %}}}
 
 %call Bamg
-[triangles vertices]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
+[triangles vertices metric]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
 
 % plug results onto model
@@ -83,4 +96,5 @@
 md.y=vertices(:,2);
 md.elements=triangles(:,1:3);
+md.dummy=metric;
 
 %Fill in rest of fields:
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2898)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2899)
@@ -9,4 +9,5 @@
 	int   noerr=1;
 	int   i;
+	int   nods=0; //to be removed
 	BamgOpts bamgopts;
 	BamgMesh bamgmesh;
@@ -18,4 +19,5 @@
 	int     NumTrianglesMesh;
 	double* TrianglesMesh=NULL;
+	double* hVerticesMesh=NULL;
 
 	/*Geom inputs: */
@@ -34,5 +36,6 @@
 	double gradation;
 	double cutoff;
-	double* metric;
+	double* metric=NULL;
+	double* field=NULL;
 
 	/*Boot module: */
@@ -79,5 +82,6 @@
 	FetchData(&TrianglesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles"));
 	bamgmesh.Triangles=TrianglesMesh;
-	bamgmesh.hVertices=NULL;
+	FetchData(&hVerticesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"hVertices"));
+	bamgmesh.hVertices=hVerticesMesh;
 	bamgmesh.NumQuadrilaterals=0;
 	bamgmesh.Quadrilaterals=NULL;
@@ -118,6 +122,9 @@
 	FetchData(&metric,NULL,NULL,mxGetField(BAMGOPTIONS,0,"metric"));
 	bamgopts.metric=metric;
+	FetchData(&field,NULL,NULL,mxGetField(BAMGOPTIONS,0,"field"));
+	bamgopts.field=field;
 
 	/*!Generate internal degree of freedom numbers: */
+	nods=bamgmesh.NumVertices;
 	Bamgx(&bamgmesh,&bamggeom,&bamgopts);
 
@@ -125,4 +132,5 @@
 	WriteData(TRIANGLESOUT,bamgmesh.Triangles,bamgmesh.NumTriangles,4);
 	WriteData(VERTICESOUT,bamgmesh.Vertices,bamgmesh.NumVertices,3);
+	WriteData(METRICOUT,bamgopts.metric,nods,3);
 
 	/*Free ressources: */
Index: /issm/trunk/src/mex/Bamg/Bamg.h
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.h	(revision 2898)
+++ /issm/trunk/src/mex/Bamg/Bamg.h	(revision 2899)
@@ -25,8 +25,9 @@
 #define TRIANGLESOUT (mxArray**)&plhs[0]
 #define VERTICESOUT  (mxArray**)&plhs[1]
+#define METRICOUT    (mxArray**)&plhs[2]
 
 /* serial arg counts: */
 #undef NLHS
-#define NLHS  2
+#define NLHS  3
 #undef NRHS
 #define NRHS  3