Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3021)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3022)
@@ -92,5 +92,5 @@
 		//generate mesh
 		if (verbosity>1) printf("   Generating Mesh...\n");
-		Triangles Th(maxnbv,Gh);
+		Triangles Th(maxnbv,Gh,bamgopts);
 		if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
 		Th.ReNumberingTheTriangleBySubDomain();
@@ -103,4 +103,8 @@
 		if (verbosity>1) printf("   Write Geometry...\n");
 		Gh.WriteGeometry(bamggeom,bamgopts);
+
+		//clean up
+	//	delete &Th;
+	//	delete &Gh;
 		/*}}}*/
 	}
@@ -183,5 +187,6 @@
 
 		/*clean up*/
-		delete &Th; //TEST crash
+		delete &Th;
+		//delete &BTh;
 		/*}}}*/
 	}
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3021)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3022)
@@ -655,6 +655,6 @@
 
 			//Genetare mesh from geometry
-			Triangles(Int4 nbvx,Geometry & G) :Gh(G),BTh(*this){
-				try { GeomToTriangles0(nbvx);}
+			Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
+				try { GeomToTriangles0(nbvx,bamgopts);}
 				catch(...) { this->~Triangles(); throw; }
 			}
@@ -680,5 +680,6 @@
 			void Insert();
 			void ForceBoundary();
-			void Heap();
+			//void Heap();
+			void Renumber(BamgOpts* bamgopts);
 			void FindSubDomain(int OutSide=0);
 			Int4 ConsRefTriangle(Int4 *) const;
@@ -742,5 +743,5 @@
 		private:
 			void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption
-			void GeomToTriangles0(Int4 nbvx);                   // the real constructor mesh generator
+			void GeomToTriangles0(Int4 nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
 			void PreInit(Int4);
 
@@ -755,6 +756,4 @@
 			Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux,
 			Int4 NbSubDomains; // 
-			Int4 NbEquiEdges;
-			Int4 NbCrackedEdges; 
 			//  Int4 nbtf;//  de triangle frontiere
 			Int4 NbOfCurves;
@@ -810,24 +809,63 @@
 	// to sort in descending order
 	template<class T>  inline void  HeapSort(T *c,long n){
-		long l,j,r,i;
-		T crit;
-		c--; // on decale de 1 pour que le tableau commence a 1
-		if( n <= 1) return;
-		l = n/2 + 1;
-		r = n;
-		while (1) { // label 2
-			if(l <= 1 ) { // label 20
-				crit = c[r];
-				c[r--] = c[1];
-				if ( r == 1 ) { c[1]=crit; return;}
-			} else  crit = c[--l]; 
+		/*Intermediary*/
+		int l,j,r,i;
+		T   crit;
+		c--;                    //the array must starts at 1 and not 0 
+		if(n<=1) return;        //return if size <=1
+		l=n/2+1;                //initialize l and r
+		r=n;
+		for(;;){
+			if(l<=1){
+				crit  =c[r];
+				c[r--]=c[1];
+				if (r==1){c[1]=crit; return;}
+			}
+			else  crit = c[--l]; 
 			j=l;
-			while (1) {// label 4
+			for(;;){
 				i=j;
 				j=2*j;
-				if  (j>r) {c[i]=crit;break;} // L8 -> G2
-				if ((j<r) && (c[j] < c[j+1])) j++; // L5
-				if (crit < c[j]) c[i]=c[j]; // L6+1 G4
-				else {c[i]=crit;break;} //L8 -> G2
+				if  (j>r) {c[i]=crit;break;}
+				if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
+				if (crit < c[j]) c[i]=c[j];       //c[j]  > crit -> stock this large value in i(<j)
+				else{c[i]=crit;break;}            //c[j]  < crit -> stock crit in i (<j)
+			}
+		}
+	}
+	// to sort in descending order and return ordering
+	template<class T>  inline void  HeapSort(int** porder,T* c,int n){
+		/*Intermediary*/
+		int  l,j,r,i;
+		T    crit;
+		int  pos;
+		int* order=NULL;
+		order=(int*)xmalloc(n*sizeof(int));
+		for(i=0;i<n;i++) order[i]=i+1;
+		c--;                    //the array must starts at 1 and not 0 
+		order--;
+		if(n<=1) return;        //return if size <=1
+		l=n/2+1;                //initialize l and r
+		r=n;
+		for(;;){
+			if(l<=1){
+				crit  =c[r]; pos=order[r];
+				c[r--]=c[1]; order[r+1]=order[1];
+				if (r==1){
+					c[1]=crit; order[1]=pos;
+					order++;
+					*porder=order;
+					return;
+				}
+			}
+			else  {crit=c[--l]; pos=order[l];}
+			j=l;
+			for(;;){
+				i=j;
+				j=2*j;
+				if  (j>r) {c[i]=crit;order[i]=pos;break;}
+				if ((j<r) && (c[j] < c[j+1]))j++;
+				if (crit < c[j]) {c[i]=c[j];order[i]=order[j];} 
+				else{c[i]=crit;order[i]=pos;break;}
 			}
 		}
@@ -958,14 +996,11 @@
 	/*INLINE functions of CLASS Triangles{{{1*/
 	inline Triangles::Triangles(Int4 i) :Gh(*new Geometry()),BTh(*this){PreInit(i);}
-	inline  void  Triangles::ReMakeTriangleContainingTheVertex()
-	  {
+	inline  void  Triangles::ReMakeTriangleContainingTheVertex(){
 		register Int4 i;
-		for ( i=0;i<nbv;i++) 
-		  {
-			vertices[i].vint = 0;
+		for ( i=0;i<nbv;i++){
+			vertices[i].vint=0;
 			vertices[i].t=0;
 		  }
-		for ( i=0;i<nbt;i++) 
-		 triangles[i].SetTriangleContainingTheVertex();
+		for ( i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
 	  }
 
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3021)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3022)
@@ -893,6 +893,4 @@
 		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);
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3021)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3022)
@@ -479,10 +479,15 @@
 
 		int i,j,k,num;
-		int verbose;
-
-		verbose=bamgopts->verbose;
+		int verbose=bamgopts->verbose;
+
+		//renumber
+		if (bamgopts->renumber){
+			if(verbose>3) printf("      Renumbering...");
+			Renumber(bamgopts);
+			if(verbose>3) printf(" done\n"); 
+		}
 
 		//Build reft that holds the number the subdomain number of each triangle
-		Int4 *reft = new Int4[nbt];
+		Int4* reft = new Int4[nbt];
 		Int4 nbInT = ConsRefTriangle(reft);
 
@@ -524,8 +529,13 @@
 
 		//Segments
+		bamgmesh->NumSegments=0;
 		if(verbose>3) printf("      writing Segments\n");
+
 		//chaining algorithm
-		int head_v[nbv];
-		int next_p[3*nbt];
+		int* head_v=NULL;
+		head_v=(int*)xmalloc(nbv*sizeof(int));
+		int* next_p=NULL;
+		next_p=(int*)xmalloc(3*nbt*sizeof(int));
+
 		for (i=0;i<nbv;i++) head_v[i]=-1;
 		k=0;
@@ -535,5 +545,7 @@
 				for (j=0;j<3;j++){
 					int v=Number(triangles[i][j]); //jth vertex of the ith triangle
+					if (k>3*nbt-1 || k<0) throw ErrorException(__FUNCT__,exprintf("k = %i, nbt = %i",k,nbt));
 					next_p[k]=head_v[v];
+					if (v>nbv-1 || v<0) throw ErrorException(__FUNCT__,exprintf("v = %i, nbv = %i",v,nbv));
 					head_v[v]=k++;
 				}
@@ -582,4 +594,7 @@
 			}
 		}
+		//clean up
+		xfree((void**)&head_v);
+		xfree((void**)&next_p);
 
 		//CrackedEdges
@@ -615,9 +630,4 @@
 					bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
 					num=num+1;
-					/*if (reft[i]==-1){
-					  if (!t(0)) printf("%i %i    %i\n",Number(t[1])+1,Number(t[2]),i);
-					  if (!t(1)) printf("%i %i    %i\n",Number(t[2])+1,Number(t[0]),i);
-					  if (!t(2)) printf("%i %i    %i\n",Number(t[0])+1,Number(t[1]),i);
-					  }*/
 				}
 			}
@@ -746,4 +756,8 @@
 			bamgmesh->EdgesOnGeometricEdge=NULL;
 		}
+
+
+		//clean up
+		delete [] reft;
 	}
 	/*}}}1*/
@@ -1462,6 +1476,8 @@
 
 		//first, build the chains that will be used for the Hessian computation, as weel as the area of each element
-		int head_s[nbv];
-		int next_p[3*nbt];
+		int* head_s=NULL;
+		head_s=(int*)xmalloc(nbv*sizeof(int));
+		int* next_p=NULL;
+		next_p=(int*)xmalloc(3*nbt*sizeof(int));
 		int  p=0;
 		//initialization
@@ -1602,4 +1618,6 @@
 
 		//clean up
+		xfree((void**)&head_s);
+		xfree((void**)&next_p);
 		delete [] detT;
 		delete [] alpha;
@@ -2769,5 +2787,5 @@
 	/*}}}1*/
 	/*FUNCTION Triangles::GeomToTriangles0{{{1*/
-	void Triangles::GeomToTriangles0(Int4 inbvx) {
+	void Triangles::GeomToTriangles0(Int4 inbvx,BamgOpts* bamgopts) {
 		/*Generate mesh from geometry*/
 
@@ -2784,5 +2802,9 @@
 		GeometricalEdge * e;
 
+		/*Get options*/
+		int verbose=bamgopts->verbose;
+
 		//initialize this
+		if (verbose>3) printf("      Generating vertices\n");
 		PreInit(inbvx);
 		nbv=0;
@@ -3064,13 +3086,17 @@
 
 		//Insert points inside existing triangles
+		if (verbose>3) printf("      Inserting vertices in mesh\n");
 		Insert();
 
 		//Force the boundary
+		if (verbose>3) printf("      Forcing boundaries\n");
 		ForceBoundary();
 
 		//Extract SubDomains
+		if (verbose>3) printf("      Extracting subdomains\n");
 		FindSubDomain();
 
 		// NewPointsOld(*this) ;
+		if (verbose>3) printf("      Generating mesh properties\n");
 		NewPoints(*this,0) ;
 	}
@@ -3381,8 +3407,4 @@
 					}
 				  }//for(phase;;)
-				/*  modif FH add Curve class  		  
-					 }}//for (iedge=0;iedge<BTh.nbe;iedge++) 
-					 */
-				// new code Curve class  	
 		  } //  end of curve loop 
 		// end new code	    
@@ -3886,29 +3908,33 @@
 	/*FUNCTION Triangles::NewPoints{{{1*/
 	void  Triangles::NewPoints(Triangles & Bh,int KeepVertices) {
+
 		long int verbosity=2;
 		Int4 nbtold(nbt),nbvold(nbv);
-		if (verbosity>2)  printf("   Triangles::NewPoints\n");
 		Int4 i,k;
 		int j;
 		Int4 *first_np_or_next_t = new Int4[nbtx];
 		Int4 NbTSwap =0;
-		//    insert old point
+
+		//display info
+		if (verbosity>2)  printf("   Triangles::NewPoints\n");
+
+		/*First, insert old points*/
+
 		nbtold = nbt;
 
 		if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){
 			//   Bh.SetVertexFieldOn();
-			for (i=0;i<Bh.nbv;i++)
-			  { 
-				Vertex & bv = Bh[i];
-				if (!bv.onGeometry) {
+			for (i=0;i<Bh.nbv;i++){ 
+				Vertex &bv=Bh[i];
+				if (!bv.onGeometry){
 					vertices[nbv].r = bv.r;
-					vertices[nbv++].m = bv.m;}
-			  }
+					vertices[nbv++].m = bv.m;
+				}
+			}
 			int nbv1=nbv;
 			Bh.ReMakeTriangleContainingTheVertex();     
 			InsertNewPoints(nbvold,NbTSwap)   ;            
 		}  
-		else 
-		 Bh.ReMakeTriangleContainingTheVertex();     
+		else Bh.ReMakeTriangleContainingTheVertex();     
 
 		Triangle *t;
@@ -3922,5 +3948,6 @@
 		// the next traingle on i is  -first_np_or_next_t[i]
 		int iter=0;
-		// Big loop 
+
+		// Big loop (most time consuming)
 		do {
 			iter++;
@@ -3929,5 +3956,4 @@
 
 			// default size of  IntersectionTriangle
-
 			i=Headt;
 			next_t=-first_np_or_next_t[i];
@@ -3999,20 +4025,7 @@
 		Int4 NbSwapf =0,NbSwp;
 
-		// bofbof 
-
-
 		NbSwp = NbSwapf;
 		for (i=0;i<nbv;i++)
 		 NbSwapf += vertices[i].Optim(0);
-		/*
-			for (i=0;i<nbv;i++)
-			NbSwapf += vertices[i].Optim(0);
-			for (i=0;i<nbv;i++)
-			NbSwapf += vertices[i].Optim(0);
-			for (i=0;i<nbv;i++)
-			NbSwapf += vertices[i].Optim(0);
-			for (i=0;i<nbv;i++)
-			NbSwapf += vertices[i].Optim(0);
-			*/
 		NbTSwap +=  NbSwapf ;
 	}
@@ -4212,4 +4225,84 @@
 	}                  
 	/*}}}1*/
+/*FUNCTION Triangles::Renumber {{{1*/
+void Triangles::Renumber(BamgOpts* bamgopts){
+
+	/*Intermediary*/
+	int i,j,k;
+	int*    r=NULL;
+	int*    rprime=NULL;
+	rprime  =(int*)   xmalloc(nbv*sizeof(int));
+
+	//renumbering with respect to lower left corner
+	if (bamgopts->renumber==1){
+
+		//compute distance to pmin
+		double* distance;
+		distance=(double*)xmalloc(nbv*sizeof(double));
+		for (i=0;i<nbv;i++) distance[i]=pow(vertices[i].r.x - pmin.x,2)+pow(vertices[i].r.y - pmin.y,2);
+
+		//get ordering
+		HeapSort(&r,distance,nbv);
+		xfree((void**)&distance);
+	}
+	//renumbering randomly
+	else if (bamgopts->renumber==2){
+
+		//allocate memory to r
+		r=(int*)xmalloc(nbv*sizeof(int));
+		int PrimeNumber= AGoodNumberPrimeWith(nbv) ;
+		int k0=rand()%nbv; 
+		for (int i=0; i<nbv; i++){
+			r[i]=(k0=(k0+PrimeNumber)%nbv)+1;
+		}
+	}
+	//else: error
+	else{
+		throw ErrorException(__FUNCT__,exprintf("renumbering option %i not supported yet",bamgopts->renumber)); 
+	}
+
+	//Keep copy of old vertices
+	Vertex* oldvertices=NULL;
+	oldvertices=(Vertex*)xmalloc(nbv*sizeof(Vertex));
+	for (i=0;i<nbv;i++){
+		oldvertices[i]=vertices[i];
+	}
+
+	//renumber vertices
+	for (i=0;i<nbv;i++){
+		//check r[i]
+		if (r[i]>nbv){
+			throw ErrorException(__FUNCT__,exprintf("r[i]>nbv (r[i]=%i and nbv=%i)",r[i],nbv));
+		}
+		if (r[i]<=0){
+			throw ErrorException(__FUNCT__,exprintf("r[i]<=0 (r[i]=%i)",r[i]));
+		}
+		vertices[i]=oldvertices[r[i]-1];
+		rprime[r[i]-1]=i+1;
+	}
+	//clean up
+	xfree((void**)&oldvertices);
+
+	//renumber triangles
+	for (i=0; i<nt;i++){
+		for (j=0;j<3;j++){
+			if (triangles[i](j)){
+				triangles[i](j)=&vertices[rprime[Number(triangles[i](j))] -1 ];
+			}
+		}
+	}
+
+	//renumber edges
+	for (i=0; i<nbe;i++){
+		for (j=0;j<2;j++){
+			edges[i].v[j]=&vertices[rprime[Number(edges[i].v[j])] -1 ];
+		}
+	}
+
+	//clean up
+	xfree((void**)&r);
+	xfree((void**)&rprime);
+}
+/*}}}1*/
 	/*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/
 	void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) {
Index: /issm/trunk/src/c/objects/BamgOpts.h
===================================================================
--- /issm/trunk/src/c/objects/BamgOpts.h	(revision 3021)
+++ /issm/trunk/src/c/objects/BamgOpts.h	(revision 3022)
@@ -14,4 +14,5 @@
 	int     Metrictype;
 	int     KeepVertices;
+	int     renumber;
 	double  maxsubdiv;
 	double  power;
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 3021)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 3022)
@@ -2,6 +2,43 @@
 %BAMG - mesh generation
 %
-%   Usage:
-%      md=bamg(md,varargin);
+%   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)
+%
+%   - anisomax    : maximum ration between the smallest and largest edges (default is 10^30)
+%   - coef        : coefficient applied to the metric (2-> twice as many elements, default is 1)
+%   - cutoff      : scalar used to compute the metric when metric type 2 or 3 are applied
+%   - err         : error used to generate the metric from a field
+%   - errg        : geometrical error (default is 0.1)
+%   - field       : field of the model that will be used to compute the metric
+%                   to apply several fields, use one column per field
+%   - gradation   : maximum ration between two adjacent edges
+%   - Hessiantype : 0 -> use double P2 projection (default)
+%                   1 -> use Green formula
+%   - KeepVertices: try to keep initial vertices when adaptation is done on an existing mesh (default 1)
+%   - MaximalAngleOfCorner: maximal angle of corners in degree (default is 10)
+%   - maxnbv      : maximum number of vertices used to allocate memory (default is 10^6)
+%   - maxsubdiv   : maximum subdivision of exisiting elements (default is 10)
+%   - metric      : matrix (numberofgrids x 3) used as a metric
+%   - Metrictype  : 1 -> absolute error          c/(err coeff^2) * Abs(H)        (default)
+%                   2 -> relative error          c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s))
+%                   3 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin)
+%   - nbjacoby    : correction used by Hessiantype=1 (default is 1)
+%   - NbSmooth    : number of metric smoothing procedure (default is 3)
+%   - omega       : relaxation parameter of the smoothing procedure (default is 1.8)
+%   - power       : power applied to the metric (default is 1)
+%   - renumber    : 0 -> no renumbering
+%                   1 -> renumber by distance to lower left corner (default)
+%                   2 -> random renumbering
+%   - splitcorner : split triangles whuch have 3 vertices on the outline (default is 1)
+%   - verbose     : level of verbosity (default is 1)
+%
+%   Examples:
+%      md=bamg(md,'domain','DomainOutline.exp','hmax',3000);
+%      md=bamg(md,'field',[md.vel_obs md.thickness],'hmax',20000,'hmin',1000);
+%      md=bamg(md,'metric',A,'hmin',1000,'hmax',20000,'gradation',3,'anisomax',1);
 
 %process options
@@ -104,26 +141,27 @@
 
 % Bamg Options {{{1
+bamg_options.anisomax=getfieldvalue(options,'anisomax',10^30);
+bamg_options.coef=getfieldvalue(options,'coef',1);
+bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
 bamg_options.err=getfieldvalue(options,'err',0.01);
 bamg_options.errg=getfieldvalue(options,'errg',0.1);
-bamg_options.coef=getfieldvalue(options,'coef',1);
+bamg_options.field=getfieldvalue(options,'field',[]);
+bamg_options.gradation=getfieldvalue(options,'gradation',1.5);
+bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
+bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
+bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
+bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1);
+bamg_options.MaximalAngleOfCorner=getfieldvalue(options,'MaximalAngleOfCorner',10);
+bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6);
 bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10);
-bamg_options.power=getfieldvalue(options,'power',1);
+bamg_options.metric=getfieldvalue(options,'metric',[]);
+bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
 bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1);
-bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
-bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
-bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1);
 bamg_options.NbSmooth=getfieldvalue(options,'NbSmooth',3);
 bamg_options.omega=getfieldvalue(options,'omega',1.8);
-bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6);
-bamg_options.MaximalAngleOfCorner=getfieldvalue(options,'MaximalAngleOfCorner',10);
-bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
-bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
-bamg_options.anisomax=getfieldvalue(options,'anisomax',10^30);
-bamg_options.gradation=getfieldvalue(options,'gradation',1.5);
-bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
+bamg_options.power=getfieldvalue(options,'power',1);
+bamg_options.renumber=getfieldvalue(options,'renumber',1);
+bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
 bamg_options.verbose=getfieldvalue(options,'verbose',1);
-bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
-bamg_options.metric=getfieldvalue(options,'metric',[]);
-bamg_options.field=getfieldvalue(options,'field',[]);
 %}}}
 
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3021)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3022)
@@ -10,4 +10,5 @@
 	int   i;
 	int   nods=0; //to be removed
+	int   lines,cols;
 	BamgOpts bamgopts;
 	BamgMesh bamgmesh;
@@ -41,5 +42,5 @@
 	double power;
 	int    Hessiantype,Metrictype,NbSmooth;
-	int    nbjacobi,KeepVertices;
+	int    nbjacobi,KeepVertices,renumber;
 	double omega;
 	double gradation,errg;
@@ -65,5 +66,6 @@
 	FetchData(&EdgesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges"));
 	bamggeom.Edges=EdgesGeom;
-	FetchData(&hVerticesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"hVertices"));
+	FetchData(&hVerticesGeom,&lines,&cols,mxGetField(BAMGGEOMETRY,0,"hVertices"));
+	if (hVerticesGeom && (cols!=1 || lines!=NumVerticesGeom)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",NumVerticesGeom,1));}
 	bamggeom.hVertices=hVerticesGeom;
 	bamggeom.MetricVertices=NULL;
@@ -95,5 +97,6 @@
 	FetchData(&TrianglesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles"));
 	bamgmesh.Triangles=TrianglesMesh;
-	FetchData(&hVerticesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"hVertices"));
+	FetchData(&hVerticesMesh,&lines,&cols,mxGetField(BAMGMESH,0,"hVertices"));
+	if (hVerticesMesh && (cols!=1 || lines!=NumVerticesMesh)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",NumVerticesMesh,1));}
 	bamgmesh.hVertices=hVerticesMesh;
 	bamgmesh.NumQuadrilaterals=0;
@@ -121,32 +124,48 @@
 	/*create bamg options input*/
 	FetchData(&coef,mxGetField(BAMGOPTIONS,0,"coef"));
+	if (coef==0) throw ErrorException(__FUNCT__,exprintf("'coef' should be positive"));
 	bamgopts.coef=coef;
 	FetchData(&maxsubdiv,mxGetField(BAMGOPTIONS,0,"maxsubdiv"));
+	if (maxsubdiv<1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >=1"));
 	bamgopts.maxsubdiv=maxsubdiv;
 	FetchData(&Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype"));
+	if (Hessiantype!=0 && Hessiantype!=1) throw ErrorException(__FUNCT__,exprintf("'Hessiantype' supported options are 0 and 1"));
 	bamgopts.Hessiantype=Hessiantype;
 	FetchData(&Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype"));
+	if (Metrictype!=0 && Metrictype!=1 && Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2"));
 	bamgopts.Metrictype=Metrictype;
 	FetchData(&KeepVertices,mxGetField(BAMGOPTIONS,0,"KeepVertices"));
+	if (KeepVertices!=0 && KeepVertices!=1) throw ErrorException(__FUNCT__,exprintf("'KeepVertices' supported options are 0 and 1"));
 	bamgopts.KeepVertices=KeepVertices;
+	FetchData(&renumber,mxGetField(BAMGOPTIONS,0,"renumber"));
+	if (renumber!=0 && renumber!=1 && renumber!=2) throw ErrorException(__FUNCT__,exprintf("'renumber' supported options are 0, 1 and 2"));
+	bamgopts.renumber=renumber;
 	FetchData(&power,mxGetField(BAMGOPTIONS,0,"power"));
 	bamgopts.power=power;
 	FetchData(&errg,mxGetField(BAMGOPTIONS,0,"errg"));
+	if (errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0"));
 	bamgopts.errg=errg;
 	FetchData(&nbjacobi,mxGetField(BAMGOPTIONS,0,"nbjacobi"));
+	if (nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0"));
 	bamgopts.nbjacobi=nbjacobi;
 	FetchData(&NbSmooth,mxGetField(BAMGOPTIONS,0,"NbSmooth"));
+	if (NbSmooth<=0) throw ErrorException(__FUNCT__,exprintf("'NbSmooth' option should be >0"));
 	bamgopts.NbSmooth=NbSmooth;
 	FetchData(&omega,mxGetField(BAMGOPTIONS,0,"omega"));
 	bamgopts.omega=omega;
 	FetchData(&maxnbv,mxGetField(BAMGOPTIONS,0,"maxnbv"));
+	if (maxnbv<3) throw ErrorException(__FUNCT__,exprintf("'maxnbv' option should be >3"));
 	bamgopts.maxnbv=maxnbv;
 	FetchData(&hmin,mxGetField(BAMGOPTIONS,0,"hmin"));
+	if (hmin<=0) throw ErrorException(__FUNCT__,exprintf("'hmin' option should be >0"));
 	bamgopts.hmin=hmin;
 	FetchData(&hmax,mxGetField(BAMGOPTIONS,0,"hmax"));
+	if (hmax<=0 || hmax<hmin) throw ErrorException(__FUNCT__,exprintf("'hmax' option should be between 0 and hmin=%g",hmin));
 	bamgopts.hmax=hmax;
 	FetchData(&anisomax,mxGetField(BAMGOPTIONS,0,"anisomax"));
+	if (anisomax<1) throw ErrorException(__FUNCT__,exprintf("'anisomax' option should be >=1"));
 	bamgopts.anisomax=anisomax;
 	FetchData(&gradation,mxGetField(BAMGOPTIONS,0,"gradation"));
+	if (anisomax<1) throw ErrorException(__FUNCT__,exprintf("'anisomax' option should be >=1"));
 	bamgopts.gradation=gradation;
 	FetchData(&cutoff,mxGetField(BAMGOPTIONS,0,"cutoff"));
@@ -158,16 +177,15 @@
 	FetchData(&MaximalAngleOfCorner,mxGetField(BAMGOPTIONS,0,"MaximalAngleOfCorner"));
 	bamgopts.MaximalAngleOfCorner=MaximalAngleOfCorner;
-	FetchData(&metric,NULL,NULL,mxGetField(BAMGOPTIONS,0,"metric"));
+	FetchData(&metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric"));
+	if (metric && (cols!=3 || lines!=NumVerticesMesh)){throw ErrorException(__FUNCT__,exprintf("the size of 'metric' should be [%i %i]",NumVerticesMesh,3));}
 	bamgopts.metric=metric;
-	FetchData(&field,NULL,&numfields,mxGetField(BAMGOPTIONS,0,"field"));
+	FetchData(&field,&lines,&numfields,mxGetField(BAMGOPTIONS,0,"field"));
+	if (field && lines!=NumVerticesMesh){throw ErrorException(__FUNCT__,exprintf("the size of 'field' should be [%i %i]",NumVerticesMesh,numfields));}
 	bamgopts.numfields=numfields;
 	bamgopts.field=field;
-	FetchData(&err,NULL,&numerr,mxGetField(BAMGOPTIONS,0,"err"));
+	FetchData(&err,NULL,&cols,mxGetField(BAMGOPTIONS,0,"err"));
+	if (numfields!=0 && cols!=numfields){throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));}
+	for (i=0;i<numfields;i++) {if (err[i]<=0) throw ErrorException(__FUNCT__,exprintf("'err' option should be >0"));};
 	bamgopts.err=err;
-
-	/*Some checks*/
-	if (numfields!=0 && numerr!=numfields){
-		throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));
-	}
 
 	/*!Generate internal degree of freedom numbers: */
@@ -192,5 +210,5 @@
 {
 	_printf_("\n");
-	_printf_("   usage: [elements,x,y]=%s(bamgmesh,bamggeom,bamgoptions,nbv);\n",__FUNCT__);
+	_printf_("   usage: [elements,vertices,segments,segmentmarkers,metric]=%s(bamgmesh,bamggeom,bamgoptions);\n",__FUNCT__);
 	_printf_("\n");
 }