Index: issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3300)
+++ issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3301)
@@ -16,16 +16,12 @@
 
 	/*Bamg options*/
-	double anisomax;
 	int    maxnbv;
-	double hmin,hmax;
-	double coef,maxsubdiv;
+	double coef;
 	int    verbosity;
 	int    NbSmooth;
-	double omega;
 
 	/*intermediary*/
 	int noerr=1;
 	int i,j,num;
-	int allquad=0;
 	double costheta=2;
 	double hminaniso=1e-100; 
@@ -35,12 +31,7 @@
 
 	/*Bamg options*/
-	anisomax=bamgopts->anisomax;
 	NbSmooth=bamgopts->NbSmooth;
-	omega=bamgopts->omega;
 	coef=bamgopts->coef;
 	maxnbv=bamgopts->maxnbv;
-	maxsubdiv=bamgopts->maxsubdiv;
-	hmin=bamgopts->hmin;
-	hmax=bamgopts->hmax;
 	verbosity=bamgopts->verbose;
 
@@ -55,22 +46,20 @@
 		/*Mesh generation {{{1*/
 
+		//Step1: generate geometry Gh
 		if (verbosity>0) printf("Construction of a mesh from a given geometry\n");
 		if (verbosity>1) printf("   Processing geometry...\n");
 		Geometry Gh(bamggeom_in,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());
-
-		//build metric if not given in input
+		bamgopts->hmin = Max(bamgopts->hmin,Gh.MinimalHmin());
+		bamgopts->hmax = Min(bamgopts->hmax,Gh.MaximalHmax());
+
+		//build metric using geometry
 		if (verbosity>1) printf("   Generating Metric...\n");
 		for(i=0;i<Gh.nbv;i++){
 			Metric M=Gh[i];
 			MatVVP2x2 Vp(M/coef);
-			Vp.Maxh(hmax);
-			Vp.Minh(hmin);
+			Vp.Maxh(bamgopts->hmax);
+			Vp.Minh(bamgopts->hmin);
 			Gh.vertices[i].m = Vp;
 		}
@@ -79,5 +68,9 @@
 		if (verbosity>1) printf("   Generating Mesh...\n");
 		Triangles Th(maxnbv,Gh,bamgopts);
+
+		//Split corners if requested
 		if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
+
+		//Renumbering
 		Th.ReNumberingTheTriangleBySubDomain();
 
@@ -101,54 +94,95 @@
 		if (verbosity>1) printf("   Processing initial mesh and geometry...\n");
 		Triangles BTh(bamggeom_in,bamgmesh_in,bamgopts); 
-		hmin=Max(hmin,BTh.MinimalHmin());
-		hmax=Min(hmax,BTh.MaximalHmax());
 
 		//Make Quadtree from background mesh
 		BTh.MakeQuadTree();
 
-		//build metric if not given in input
+		//Bound hmin and hmax
+		bamgopts->hmin=Max(bamgopts->hmin,BTh.MinimalHmin());
+		bamgopts->hmax=Min(bamgopts->hmax,BTh.MaximalHmax());
+
+		//Generate initial metric
 		if (bamgopts->metric){
 			if (verbosity>1) printf("   Processing Metric...\n");
-			BTh.ReadMetric(bamgopts,hmin,hmax,coef);
+			BTh.ReadMetric(bamgopts);
 		}
 		else { 
-
-			// init with hmax 
-			Metric Mhmax(hmax);
+			if (verbosity>1) printf("   Generating initial Metric...\n");
+			Metric Mhmax(bamgopts->hmax);
 			for (int iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;
-			// change using hVertices if required
-			if (bamgmesh_in->hVertices){
-				for (i=0;i<BTh.nbv;i++){
-					if (!isnan(bamgmesh_in->hVertices[i])){
-						BTh[i].m=Metric((float)bamgmesh_in->hVertices[i]);
-					}
+		}
+
+		//use present fields to generate metric if present
+		if (bamgopts->field){
+			if (verbosity>1) printf("   Merge metric with field provided...\n");
+			BTh.AddMetric(bamgopts);
+		}
+
+		// change using hVertices if provided
+		if (bamgmesh_in->hVertices){
+			if (verbosity>1) printf("   Merging Metric with hVertices...\n");
+			for (i=0;i<BTh.nbv;i++){
+				if (!isnan(bamgmesh_in->hVertices[i])){
+					BTh[i].m=Metric((float)bamgmesh_in->hVertices[i]);
 				}
-
-			}
-
-			//use present fields to generate metric if present
-			if (bamgopts->field){
-				if (verbosity>1) printf("   Generating Metric from solution field...\n");
-				BTh.AddMetric(bamgopts);
-			}
-		}
-
-		BTh.AddGeometryMetric(bamgopts);
-		if (verbosity>1) printf("   Smoothing Metric...");
-		if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation);
-		if (verbosity>1) printf(" done\n");
-		BTh.MaxSubDivision(maxsubdiv);
-		BTh.BoundAnisotropy(anisomax,hminaniso);
-
-		//Build new mesh Thr
+			}
+		}
+
+		// change using hminVertices if provided
+		if (bamgopts->hminVertices){
+			if (verbosity>1) printf("   Merging Metric with hminVertices...\n");
+			for (i=0;i<BTh.nbv;i++){
+				if (!isnan(bamgopts->hminVertices[i])){
+					Metric M=BTh.vertices[i].m;
+					MatVVP2x2 Vp(M/coef);
+					Vp.Minh(bamgopts->hminVertices[i]);
+					BTh.vertices[i].m=Vp;
+				}
+			}
+		}
+
+		// change using hmaxVertices if provided
+		if (bamgopts->hmaxVertices){
+			if (verbosity>1) printf("   Merging Metric with hmaxVertices...\n");
+			for (i=0;i<BTh.nbv;i++){
+				if (!isnan(bamgopts->hmaxVertices[i])){
+					Metric M=BTh.vertices[i].m;
+					MatVVP2x2 Vp(M/coef);
+					Vp.Minh(bamgopts->hmaxVertices[i]);
+					BTh.vertices[i].m=Vp;
+				}
+			}
+		}
+
+		//Add geometry metric if provided
+		if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts);
+
+		//Smoothe metric
+		BTh.SmoothMetric(bamgopts->gradation);
+
+		//Control element subdivision
+		BTh.MaxSubDivision(bamgopts->maxsubdiv);
+
+		//Bound anisotropy
+		BTh.BoundAnisotropy(bamgopts->anisomax,hminaniso);
+
+		//Build new mesh
 		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,bamgopts,bamgopts->KeepVertices)));
 		if (Thr != &BTh) delete Thr;
+
+		//Make quadrangle if requested
 		if(costheta<=1.0) Th.MakeQuadrangles(costheta);
-		if (allquad) Th.SplitElement(allquad==2);
+		//if () Th.SplitElement(2);
+
+		//Split corners if requested
 		if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
+
+		//Renumber by subdomain
 		Th.ReNumberingTheTriangleBySubDomain();
-		if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
+
+		//Smooth vertices
+		if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,bamgopts->omega);
 
 		//display info
Index: issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
===================================================================
--- issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp	(revision 3300)
+++ issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp	(revision 3301)
@@ -25,71 +25,67 @@
 		double dca,dcb,dcta,dctb;
 		double ddca,ddcb,ddcta,ddctb;
-		// double t1 = 1 -theta;
-		// double t1t1 = t1*t1;
 		double tt = theta*theta;
-		if ( theta<0){
-			throw ErrorException(__FUNCT__,exprintf("theta<0"));
+
+		//check theta
+		assert(theta>=0 && theta<=1);
+
+		if (TgA()){ 
+			if (TgB()){
+				// Tangent A and B provided:
+				// interpolation d'hermite
+				dcb = 6*theta*(1-theta);
+				ddcb = 6*(1-2*theta);
+				dca = -dcb;
+				ddca = -ddcb;
+				dcta =  (3*theta - 4)*theta + 1;
+				ddcta=6*theta-4;
+				dctb = 3*tt - 2*theta;
+				ddctb = 6*theta-2;
+			}
+			else {
+				//Tangent A provided but tangent B not provided
+				// 1-t*t, t-t*t, t*t
+				double t = theta;
+				dcb = 2*t;
+				ddcb = 2;
+				dca = -dcb;
+				ddca = -2;
+				dcta = 1-dcb;
+				ddcta = -ddcb;
+				dctb=0;    
+				ddctb=0;    
+			}
 		}
-		if ( theta>1){
-			throw ErrorException(__FUNCT__,exprintf("theta>1"));
+		else{
+			if (TgB()){
+				//Tangent B provided but tangent A not provided
+				double t = 1-theta;
+				dca = -2*t;
+				ddca = 2;
+				dcb = -dca;
+				ddcb = -2;
+				dctb = 1+dca;
+				ddctb= ddca;
+				dcta =0;
+				ddcta =0;
+			}
+			else {
+				//Neither thangent A nor tangent B provided
+				// lagrange P1
+				t=B-A;
+				return 0;
+			} 
 		}
-		if (TgA()) 
-		 if (TgB()) // interpolation d'hermite
-			{ //cb =  theta*theta*(3-2*theta);
-			 dcb = 6*theta*(1-theta);
-			 ddcb = 6*(1-2*theta);
-			 //ca =  1-cb;     
-			 dca = -dcb;
-			 ddca = -ddcb;
-
-			 // cta = (1-theta)*(1-theta)*theta;
-			 dcta =  (3*theta - 4)*theta + 1;
-			 ddcta=6*theta-4;
-
-			 //ctb = (theta-1)*theta*theta ;
-			 dctb = 3*tt - 2*theta;
-			 ddctb = 6*theta-2;
-			}
-		 else { // 1-t*t, t-t*t, t*t
-			 double t = theta;
-			 // cb = t*t;
-			 dcb = 2*t;
-			 ddcb = 2;
-			 //ca = 1-cb;
-			 dca = -dcb;
-			 ddca = -2;
-			 // cta= t-cb;
-			 dcta = 1-dcb;
-			 ddcta = -ddcb;
-			 // ctb =0;
-			 dctb=0;    
-			 ddctb=0;    
-		 }    
-		else
-		 if (TgB()){
-			 double t = 1-theta;
-			 //ca = t*t;
-			 dca = -2*t;
-			 ddca = 2;
-			 //cb = 1-ca;
-			 dcb = -dca;
-			 ddcb = -2;
-			 //ctb= -t+ca;
-			 dctb = 1+dca;
-			 ddctb= ddca;
-			 //cta=0;    
-			 dcta =0;
-			 ddcta =0;
-		 }
-		 else {t=B-A;return 0;} // lagrange P1
-		R2 d =  A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb;
-
+		R2 d  =  A*dca  + B*dcb  + tg[0]* dcta  + tg[1] * dctb;
 		R2 dd =  A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;
 		double d2=(d,d);
 		double sd2 = sqrt(d2);
 		t=d;
-		if(d2>1.0e-20) {t/=sd2;return Abs(Det(d,dd))/(d2*sd2);}
+		if(d2>1.0e-20){
+			t/=sd2;
+			return Abs(Det(d,dd))/(d2*sd2);
+		}
 		else return 0;
-	  }
+	}
 	/*}}}1*/
 	/*FUNCTION GeometricalEdge::F{{{1*/
Index: issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3300)
+++ issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3301)
@@ -742,5 +742,4 @@
 		// generation of  all the tangents
 		for (i=0;i<nbe;i++) {
-
 			//Get AB = vertex1->vertex2
 			R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;        
@@ -748,10 +747,9 @@
 			double lAB=Norme2(AB);
 			//initialize tangent
-			double ltg2[2];
-			ltg2[0]=0;ltg2[1]=0;
+			double ltg2[2]={0.0};
 
 			//loop over the 2 vertices of the edge
 			for (jj=0;jj<2;jj++) {
-				R2    tg =edges[i].tg[jj];
+				R2     tg =edges[i].tg[jj];
 				double ltg=Norme2(tg);
 
@@ -760,15 +758,22 @@
 					//if the current vertex of the edge is not a corner
 					if(!edges[i].v[jj]->Corner()){
-						tg =  edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r;
-						ltg=  Norme2(tg);
-						tg =  tg *(lAB/ltg),ltg=lAB;
+						tg = edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; //previous point and next point on curve
+						ltg= Norme2(tg);
+						tg = tg *(lAB/ltg);
+						ltg= lAB;
 					}
 					//else:  a Corner with no tangent => nothing to do    
 				}
 				else{
-					tg = tg *(lAB/ltg),ltg=lAB;
-				}
+					//tangent has already been computed
+					tg = tg*(lAB/ltg),ltg=lAB;
+				}
+
 				ltg2[jj] = ltg;
-				if ((tg,AB)<0) tg = -tg;
+
+				if ((tg,AB)<0){
+					tg = -tg;
+				}
+
 				edges[i].tg[jj] = tg;
 			}
Index: issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3300)
+++ issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3301)
@@ -822,11 +822,13 @@
 	/*}}}1*/
 	/*FUNCTION Triangles::ReadMetric{{{1*/
-	void Triangles::ReadMetric(BamgOpts* bamgopts,const double hmin1=1.0e-30,const double hmax1=1.0e30,const double coef=1) {
+	void Triangles::ReadMetric(const BamgOpts* bamgopts) {
+
+		/*Intermediary*/
 		int  i,j;
 
 		if(bamgopts->verbose>3) printf("      processing metric\n");
-
-		double hmin = Max(hmin1,MinimalHmin());
-		double hmax = Min(hmax1,MaximalHmax());
+		double hmin = Max(bamgopts->hmin,MinimalHmin());
+		double hmax = Min(bamgopts->hmax,MaximalHmax());
+		double coef = bamgopts->coef;
 
 		//for now we only use j==3
@@ -904,5 +906,5 @@
 				Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
 
-				GeometricalEdge * eg = GV;
+				GeometricalEdge* eg = GV;
 				double s = GV;
 				R2 tg;
@@ -914,4 +916,5 @@
 				}
 				double hn=Min(hmax,ht*anisomax);
+
 				if (ht<=0 || hn<=0){
 					throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
@@ -4418,10 +4421,10 @@
 				}  
 			}
-		if (quadtree) {delete quadtree;
+		if (quadtree){
+			delete quadtree;
 			quadtree = new QuadTree(this);
 		}
 
 		for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
-		printf("end renumbervertex\n");
 	}
 	/*}}}1*/
Index: issm/trunk/src/c/Bamgx/objects/Triangles.h
===================================================================
--- issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3300)
+++ issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3301)
@@ -136,5 +136,5 @@
 			void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
 			void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
-			void ReadMetric(BamgOpts* bamgopts,const double hmin,const double hmax,const double coef);
+			void ReadMetric(const BamgOpts* bamgopts);
 			void WriteMetric(BamgOpts* bamgopts);
 			void AddMetric(BamgOpts* bamgopts);
Index: issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp
===================================================================
--- issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp	(revision 3300)
+++ issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp	(revision 3301)
@@ -11,9 +11,4 @@
 #include "../toolkits/toolkits.h"
 
-/*Bamg: */
-#include <unistd.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
 #include "../Bamgx/objects/BamgObjects.h"
 
Index: issm/trunk/src/c/objects/BamgOpts.cpp
===================================================================
--- issm/trunk/src/c/objects/BamgOpts.cpp	(revision 3300)
+++ issm/trunk/src/c/objects/BamgOpts.cpp	(revision 3301)
@@ -19,7 +19,10 @@
 	bamgopts->hmin=0;
 	bamgopts->hmax=0;
+	bamgopts->hminVertices=NULL;
+	bamgopts->hmaxVertices=NULL;
 	bamgopts->gradation=0;
 	bamgopts->cutoff=0;
 	bamgopts->splitcorners=0;
+	bamgopts->geometricalmetric=0;
 	bamgopts->verbose=0;
 	bamgopts->err=NULL;
@@ -43,4 +46,5 @@
 	if (bamgopts->errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0"));
 	if (bamgopts->nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0"));
+	if (bamgopts->geometricalmetric!=0  && bamgopts->geometricalmetric!=1) throw ErrorException(__FUNCT__,exprintf("'geometricalmetric' supported options are 0 and 1"));
 	if (bamgopts->NbSmooth<=0) throw ErrorException(__FUNCT__,exprintf("'NbSmooth' option should be >0"));
 	if (bamgopts->maxnbv<3) throw ErrorException(__FUNCT__,exprintf("'maxnbv' option should be >3"));
Index: issm/trunk/src/c/objects/BamgOpts.h
===================================================================
--- issm/trunk/src/c/objects/BamgOpts.h	(revision 3300)
+++ issm/trunk/src/c/objects/BamgOpts.h	(revision 3301)
@@ -22,10 +22,13 @@
 	double  hmin;
 	double  hmax;
+	double* hminVertices;
+	double* hmaxVertices;
 	double  gradation;
 	double  cutoff;
 	int     splitcorners;
+	int     geometricalmetric;
 	int     verbose;
-	double*  err;
-	double   errg;
+	double* err;
+	double  errg;
 	double  coef;
 	double* metric;
Index: issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- issm/trunk/src/m/classes/public/bamg.m	(revision 3300)
+++ issm/trunk/src/m/classes/public/bamg.m	(revision 3301)
@@ -4,8 +4,10 @@
 %   Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):
 %
-%   - domain : followed by an ARGUS file that prescribes the domain outline (rifts and holes)
-%   - hmin   : minimum edge length (default is 10^-100)
-%   - hmax   : maximum esge length (default is 10^100)
-%   - hVertices : imposed edge length for each vertex (geometry or mesh)
+%   - domain: followed by an ARGUS file that prescribes the domain outline (rifts and holes)
+%   - hmin  : minimum edge length (default is 10^-100)
+%   - hmax  : maximum esge length (default is 10^100)
+%   - hVertices   : imposed edge length for each vertex (geometry or mesh)
+%   - hminVertices: minimum edge length for each vertex (mesh)
+%   - hmaxVertices: maximum edge length for each vertex (mesh)
 %
 %   - anisomax    : maximum ration between the smallest and largest edges (default is 10^30)
@@ -31,5 +33,6 @@
 %   - omega       : relaxation parameter of the smoothing procedure (default is 1.8)
 %   - power       : power applied to the metric (default is 1)
-%   - splitcorner : split triangles whuch have 3 vertices on the outline (default is 1)
+%   - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1)
+%   - geometricalmetric : Take the geometry into account to generate the metric (default is 0)
 %   - verbose     : level of verbosity (default is 1)
 %
@@ -227,4 +230,6 @@
 bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
 bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
+bamg_options.hminVertices=getfieldvalue(options,'hminVertices',[]);
+bamg_options.hmaxVertices=getfieldvalue(options,'hmaxVertices',[]);
 bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1);
 bamg_options.MaximalAngleOfCorner=getfieldvalue(options,'MaximalAngleOfCorner',10);
@@ -238,17 +243,7 @@
 bamg_options.power=getfieldvalue(options,'power',1);
 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
+bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0);
 bamg_options.verbose=getfieldvalue(options,'verbose',1);
 %}}}
-
-%check bamg options
-fieldsize=size(bamg_options.field);
-if fieldsize(1),
-	if fieldsize(1)~=md.numberofgrids,
-		error(['bamg error message, ''field'' should have ' num2str(md.numberofgrids) ' lines (numberofgrids)']);
-	end
-	if (fieldsize(2)~=0 & length(bamg_options.err)~=fieldsize(2)),
-		error(['bamg error message, ''err'' should be of length ' num2str(fieldsize(2)) ' (As many as fields)']);
-	end
-end
 
 %call Bamg
Index: issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3300)
+++ issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3301)
@@ -82,5 +82,10 @@
 	FetchData(&bamgopts.verbose,mxGetField(BAMGOPTIONS,0,"verbose"));
 	FetchData(&bamgopts.splitcorners,mxGetField(BAMGOPTIONS,0,"splitcorners"));
+	FetchData(&bamgopts.geometricalmetric,mxGetField(BAMGOPTIONS,0,"geometricalmetric"));
 	FetchData(&bamgopts.MaximalAngleOfCorner,mxGetField(BAMGOPTIONS,0,"MaximalAngleOfCorner"));
+	FetchData(&bamgopts.hminVertices,&lines,&cols,mxGetField(BAMGOPTIONS,0,"hminVertices"));
+	if (bamgopts.hminVertices && (cols!=1 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hminVertices' should be [%i %i]",bamgmesh_in.NumVertices,1));}
+	FetchData(&bamgopts.hmaxVertices,&lines,&cols,mxGetField(BAMGOPTIONS,0,"hmaxVertices"));
+	if (bamgopts.hmaxVertices && (cols!=1 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hmaxVertices' should be [%i %i]",bamgmesh_in.NumVertices,1));}
 	FetchData(&bamgopts.metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric"));
 	if (bamgopts.metric && (cols!=3 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'metric' should be [%i %i]",bamgmesh_in.NumVertices,3));}
