Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2779)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2780)
@@ -20,9 +20,9 @@
 #include <new>
 #include <cassert>
+#include <iomanip>
+#include <fstream>
 #include "Meshio.h"
-#include <iomanip>
 #include "Mesh2.h"
 #include "QuadTree.h"
-#include <fstream>
 
 using namespace bamg;
@@ -43,5 +43,5 @@
 	hinterpole=1; // by def interpolation a h 
 	int fileout=0;
-	int nbvx = 50000;
+	int nbvx = 1000000;
 	int iso =0,AbsError=0,nbjacoby=1,allquad=0;
 	int NoMeshReconstruction=0;
@@ -50,5 +50,5 @@
 	Real8 cutoffradian=-1;
 	double anisomax = 1e6;
-	double err=0.01,errg=0.1,coef=1,hmin=1.e-100,hmax=1.e17,raison=0,CutOff=1e-5;
+	double err=0.01,errg=0.1,coef=1,hmin=1.e-100,hmax=1.e17,gradation=1.2,cutoff=1e-5;
 	int KeepBackVertices=1;
 	double hminaniso=1e-100; 
@@ -76,9 +76,11 @@
 	/*testing*/
 	int splitpbedge=1;
-	int nel=0;
-	int nods=0;
-	double* x=NULL;
-	double* y=NULL;
-	double* elements=NULL;
+
+	/*Bamg options*/
+	iso=bamgargs->iso;
+	hmin=bamgargs->hmin;
+	hmax=bamgargs->hmax;
+	gradation=bamgargs->gradation;
+	cutoff=bamgargs->cutoff;
 
 	/*Recover options from inputs: */
@@ -108,10 +110,11 @@
 		rBBeqMBB = !strcmp(rBB,fMBB);
 
-	if(1){
+	if(bamgmesh->numberofelements==0){
+		/*Mesh generation {{{1*/
 		if (verbosity>0) printf("Construction of a mesh from a given geometry\n");
 		Geometry Gh(bamggeom);
 		hmin = Max(hmin,Gh.MinimalHmin());
 		hmax = Min(hmax,Gh.MaximalHmax());
-		
+
 		//build metric if not given in input
 		for(i=0;i<Gh.nbv;i++)
@@ -127,5 +130,5 @@
 		Triangles Th(nbvx,Gh);
 
-		/*Build output {{{1*/
+		//Build output
 		nelout=Th.nbt-Th.NbOutT; //number of triangles - number of external triangles
 		nodsout=Th.nbv;
@@ -150,250 +153,265 @@
 			}
 		}
-
-		/*Assign output pointers*/
-		bamgmesh->numberofelements=nelout;
-		bamgmesh->numberofnodes=nodsout;
-		bamgmesh->index=elementsout;
-		bamgmesh->x=xout;
-		bamgmesh->y=yout;
-
-		return noerr;
-		/*}}}1*/
-
+		/*}}}*/
 	}
 
-	if (0){
-	/*Anisotropic mesh adaptation {{{1*/
-	/*Read background mesh from simple delaunay triangulation: */
-	Triangles BTh(elements,x,y,nel,nods,nbvx,cutoffradian); // read the background mesh 
-
-	hmin = Max(hmin,BTh.MinimalHmin());
-	hmax = Min(hmax,BTh.MaximalHmax());
-
-	throw ErrorException(__FUNCT__," done");
-
-	BTh.MakeQuadTree();
-	if (fmetrix) 
-		BTh.ReadMetric(fmetrix,hmin,hmax,coef);
-	else
-	{ // init with hmax 
-		Metric Mhmax(hmax);
-		for (Int4 iv=0;iv<BTh.nbv;iv++)
-			BTh[iv].m = Mhmax;
+	else{
+		/*Anisotropic mesh adaptation {{{1*/
+		if (verbosity>0) printf("Anisotropic mesh adaptation\n");
+
+		/*Read background mesh from simple delaunay triangulation: */
+		Triangles BTh(bamgmesh->index,bamgmesh->x,bamgmesh->y,bamgmesh->numberofelements,bamgmesh->numberofnodes,nbvx,cutoffradian); // read the background mesh 
+
+		hmin = Max(hmin,BTh.MinimalHmin());
+		hmax = Min(hmax,BTh.MaximalHmax());
+
+		throw ErrorException(__FUNCT__," done");
+
+		BTh.MakeQuadTree();
+		if (fmetrix){
+			BTh.ReadMetric(fmetrix,hmin,hmax,coef);
+		}
+		else { // init with hmax 
+			Metric Mhmax(hmax);
+			for (Int4 iv=0;iv<BTh.nbv;iv++)
+			 BTh[iv].m = Mhmax;
+		}
+		if (fMbb) {
+			solMbb = ReadbbFile(fMbb,nbsolbb,lsolbb,2,2);
+			if (lsolbb != BTh.nbv) 
+			  {
+				cerr << "fatal error  nbsol " << nbsolbb << " " << lsolbb<< " =! " << BTh.nbv << endl;
+				cerr << "  size of sol incorrect " << endl;
+				MeshError(99);
+			  }
+			assert(lsolbb==BTh.nbv);
+			BTh.IntersectConsMetric(solMbb,nbsolbb,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:cutoff,
+						nbjacoby,Rescaling,power,ChoiseHessien);
+			if(!rbbeqMbb)
+			 delete [] solMbb,solMbb =0;
+
+		}
+		if (fMBB) {
+			solMBB = ReadBBFile(fMBB,nbsolBB,lsolBB,typesolsBB,2,2);
+			if (lsolBB != BTh.nbv) 
+			  {
+				cerr << "fatal error  nbsol " << nbsolBB << " " << lsolBB << " =! " << BTh.nbv << endl;
+				cerr << "  size of sol incorrect " << endl;
+				MeshError(99);
+			  }
+			assert(lsolBB==BTh.nbv);
+			BTh.IntersectConsMetric(solMBB,nbsolBB,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:cutoff,
+						nbjacoby,Rescaling,ChoiseHessien);
+			if(!rBBeqMBB)
+			 delete [] solMBB,solMBB =0;
+
+		}
+
+		BTh.IntersectGeomMetric(errg,iso);
+
+		if(gradation) BTh.SmoothMetric(gradation);
+		BTh.MaxSubDivision(maxsubdiv);
+		BTh.BoundAnisotropy(anisomax,hminaniso);
+
+		if(foM) {
+			if(verbosity >2)
+			 cout << " -- write Metric  file " << foM <<endl;
+
+			ofstream f(foM);
+			if(f) BTh.WriteMetric(f,iso);
+		}
+
+		if ( fileout) {
+			if ( NoMeshReconstruction)
+			 if (( fmeshback == fmeshr) || (!strcmp(fmeshback,fmeshr)))
+			  Thr=&BTh,Thb=0; // back and r mesh are the same 
+			 else
+			  Thr= new Triangles(fmeshr,cutoffradian),Thb=&BTh; // read the new 
+
+
+
+			Triangles & Th( *(NoMeshReconstruction 
+							?  new Triangles(*Thr,&Thr->Gh,Thb,nbvx) // copy the mesh + free space to modification  
+							:  new Triangles(nbvx,BTh,KeepBackVertices)     // construct a new mesh 
+							));
+			if (Thr != &BTh) delete Thr;
+
+			if(costheta<=1.0)
+			 Th.MakeQuadrangles(costheta);
+			if (allquad)
+			 Th.SplitElement(allquad==2);
+			if(SplitEdgeWith2Boundary)
+			 Th.SplitInternalEdgeWithBorderVertices();
+			Th.ReNumberingTheTriangleBySubDomain();
+
+			if(verbosity>3) Th.ShowHistogram();
+			if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
+			if(verbosity>3 && NbSmooth>0) Th.ShowHistogram();
+
+			Th.BTh.ReMakeTriangleContainingTheVertex();
+
+			if (fmeshout) Th.Write(fmeshout  ,Triangles::BDMesh);
+			if (famfmt)   Th.Write(famfmt    ,Triangles::am_fmtMesh);
+			if (fam)      Th.Write(fam       ,Triangles::amMesh);
+			if (famdba)   Th.Write(famdba    ,Triangles::amdbaMesh);
+			if (fftq)     Th.Write(fftq      ,Triangles::ftqMesh);
+			if (fmsh)     Th.Write(fmsh      ,Triangles::mshMesh);
+			if (fnopo)    Th.Write(fnopo     ,Triangles::NOPOMesh);
+
+			if ( ( rbb && wbb)  ||( rBB && wBB)){  // the code for interpolation 
+
+				if(verbosity >1)
+				  {
+					if (rbb ) 
+					 cout << " -- interpolation P1  files : " << rbb << " -> " << wbb <<endl;
+					if (rBB ) 
+					 cout << " -- interpolation P1  files : " << rBB<< " -> " << wBB <<endl;
+				  }
+				const int dim=2;
+				// optimisation read only si rbb != fMbb
+
+				double *solbb=0;
+				double *solBB=0;
+
+				if (rbb) 
+				 solbb =  rbbeqMbb? solMbb : ReadbbFile(rbb,nbsolbb,lsolbb,2,2);
+				if (rBB) 
+				 solBB =  rBBeqMBB? solMBB : ReadBBFile(rBB,nbsolBB,lsolBB,typesolsBB,2,2);
+
+				// cout << " " << rbbeqMbb << " " <<  sol << " " << nbsol << " " << lsol << endl; 
+				if (!solBB && !solbb ) 
+				  {
+					cerr << " Fatal Error "  << "solBB =  " <<  solBB << " solbb= " << solbb << endl;
+					exit(2);}
+
+					ofstream *fbb = wbb ? new ofstream(wbb) :0;
+					ofstream *fBB = wBB ? new ofstream(wBB) :0;
+					Int4   nbfieldBB = 0, nbfieldbb = nbsolbb;
+					if (fbb)
+					 *fbb  << dim << " " << nbsolbb << " " << Th.nbv <<" " << 2 << endl; 
+					if (fBB)
+					  {
+						int i;
+						*fBB  << dim << " " << nbsolBB ;
+						for ( i=0;i<nbsolBB;i++)
+						 *fBB << " " << (typesolsBB ?typesolsBB[i]+1:1) ;
+						*fBB << " " << Th.nbv <<" " << 2 << endl; 
+						assert(dim==2);
+						for ( i=0;i<nbsolBB;i++)
+						 nbfieldBB += typesolsBB ? typesolsBB[i]+1 : 1;
+						// this code is good only if dim == 2 
+					  }
+					cout << "nb of field BB " << nbfieldBB << endl;
+					//		 if(fBB) fBB->precision(15);
+					//if(fbb) fbb->precision(15);
+					for(i=0;i<Th.nbv;i++)
+					  {
+						Int4 i0,i1,i2;
+						double a[3];
+						Icoor2 dete[3];
+						I2 I = Th.BTh.toI2(Th.vertices[i].r);
+						Triangle & tb = *Th.BTh.FindTriangleContening(I,dete);
+
+						if (tb.det>0) { // internal point 
+							a[0]= (Real8) dete[0]/ tb.det;
+							a[1]= (Real8) dete[1] / tb.det;
+							a[2] = (Real8) dete[2] / tb.det;
+							i0=Th.BTh.Number(tb[0]);
+							i1=Th.BTh.Number(tb[1]);
+							i2=Th.BTh.Number(tb[2]);}
+						else {
+							double aa,bb;
+
+							TriangleAdjacent  ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
+
+							int k = ta;
+							Triangle & tc = *(Triangle *)ta;
+							i0=Th.BTh.Number(tc[0]);
+							i1=Th.BTh.Number(tc[1]);
+							i2=Th.BTh.Number(tc[2]);
+							a[VerticesOfTriangularEdge[k][1]] =aa;
+							a[VerticesOfTriangularEdge[k][0]] = bb;
+							a[OppositeVertex[k]] = 1- aa -bb;}
+
+							Int4  ibb0 = nbfieldbb*i0;
+							Int4  ibb1 = nbfieldbb*i1;
+							Int4  ibb2 = nbfieldbb*i2;
+
+							Int4  iBB0 = nbfieldBB*i0;
+							Int4  iBB1 = nbfieldBB*i1;
+							Int4  iBB2 = nbfieldBB*i2;
+							Int4 j=0;
+
+							for ( j=0;j<nbfieldbb;j++) 
+							 *fbb  << " " << ( a[0] * solbb[ibb0++] + a[1] * solbb[ibb1++] + a[2]* solbb[ibb2++]) ;
+
+							for (j=0;j<nbfieldBB;j++) 
+							 *fBB  << " " << ( a[0] * solBB[iBB0++] + a[1] * solBB[iBB1++] + a[2]* solBB[iBB2++]) ;
+
+							if (fbb) *fbb  << endl;
+							if (fBB) *fBB  << endl;
+					  }
+					if (fbb)
+					 delete fbb  ; // close
+					if (fBB)
+					 delete fBB  ; // close
+					if (solbb) 
+					 delete [] solbb;
+					if (solBB) 
+					 delete [] solBB;
+
+			}
+
+			if(verbosity>0) 
+			  {
+				if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2) 
+				 cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
+				if (Th.NbOfQuad ) 
+				 cout  << " Nb Quadrilaterals = " << Th.NbOfQuad  ;
+				cout   << endl;
+			  }
+
+			//Build output
+			nelout=Th.nbt-Th.NbOutT; //number of triangles - number of external triangles
+			nodsout=Th.nbv;
+
+			xout=(double*)xmalloc(nodsout*sizeof(double));
+			yout=(double*)xmalloc(nodsout*sizeof(double));
+			for (i=0;i<nodsout;i++){
+				xout[i]=Th.vertices[i].r.x;
+				yout[i]=Th.vertices[i].r.y;
+			}
+
+			elementsout=(double*)xmalloc(3*nelout*sizeof(double));
+			num=0;
+			for(i=0;i<Th.nbt;i++){ 
+				Triangle &t=Th.triangles[i];
+				if (t.link){
+					//write the element only if it is part of the mesh (no boundary element)
+					elementsout[3*num+0]=Th.Number(t[0])+1;
+					elementsout[3*num+1]=Th.Number(t[1])+1;
+					elementsout[3*num+2]=Th.Number(t[2])+1;
+					num=num+1;
+				}
+			}
+
+			/*clean up*/
+			delete &Th;
+		}
+		/*}}}*/
 	}
-	if (fMbb)
-	{
-		solMbb = ReadbbFile(fMbb,nbsolbb,lsolbb,2,2);
-		if (lsolbb != BTh.nbv) 
-		{
-			cerr << "fatal error  nbsol " << nbsolbb << " " << lsolbb<< " =! " << BTh.nbv << endl;
-			cerr << "  size of sol incorrect " << endl;
-			MeshError(99);
-		}
-		assert(lsolbb==BTh.nbv);
-		BTh.IntersectConsMetric(solMbb,nbsolbb,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:CutOff,
-				nbjacoby,Rescaling,power,ChoiseHessien);
-		if(!rbbeqMbb)
-			delete [] solMbb,solMbb =0;
-
-	}
-	if (fMBB)
-	{
-		solMBB = ReadBBFile(fMBB,nbsolBB,lsolBB,typesolsBB,2,2);
-		if (lsolBB != BTh.nbv) 
-		{
-			cerr << "fatal error  nbsol " << nbsolBB << " " << lsolBB << " =! " << BTh.nbv << endl;
-			cerr << "  size of sol incorrect " << endl;
-			MeshError(99);
-		}
-		assert(lsolBB==BTh.nbv);
-		BTh.IntersectConsMetric(solMBB,nbsolBB,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:CutOff,
-				nbjacoby,Rescaling,ChoiseHessien);
-		if(!rBBeqMBB)
-			delete [] solMBB,solMBB =0;
-
-	}
-
-	BTh.IntersectGeomMetric(errg,iso);
-	if(raison) 
-		BTh.SmoothMetric(raison);
-	BTh.MaxSubDivision(maxsubdiv);
-	//
-	if (iso) 
-		anisomax=1;
-	BTh.BoundAnisotropy(anisomax,hminaniso);
-	if(foM)
-	{
-		if(verbosity >2)
-			cout << " -- write Metric  file " << foM <<endl;
-
-		ofstream f(foM);
-		if(f) BTh.WriteMetric(f,iso);
-	}
-
-	if ( fileout) 
-	{
-		if ( NoMeshReconstruction)
-			if (( fmeshback == fmeshr) || (!strcmp(fmeshback,fmeshr)))
-				Thr=&BTh,Thb=0; // back and r mesh are the same 
-			else
-				Thr= new Triangles(fmeshr,cutoffradian),Thb=&BTh; // read the new 
-
-
-
-		Triangles & Th( *(NoMeshReconstruction 
-					?  new Triangles(*Thr,&Thr->Gh,Thb,nbvx) // copy the mesh + free space to modification  
-					:  new Triangles(nbvx,BTh,KeepBackVertices)     // construct a new mesh 
-					));
-		if (Thr != &BTh) delete Thr;
-
-		if(costheta<=1.0)
-			Th.MakeQuadrangles(costheta);
-		if (allquad)
-			Th.SplitElement(allquad==2);
-		if(SplitEdgeWith2Boundary)
-			Th.SplitInternalEdgeWithBorderVertices();
-		Th.ReNumberingTheTriangleBySubDomain();
-
-		if(verbosity>3) Th.ShowHistogram();
-		if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
-		if(verbosity>3 && NbSmooth>0) Th.ShowHistogram();
-		
-		Th.BTh.ReMakeTriangleContainingTheVertex();
-
-		if (fmeshout) Th.Write(fmeshout  ,Triangles::BDMesh);
-		if (famfmt)   Th.Write(famfmt    ,Triangles::am_fmtMesh);
-		if (fam)      Th.Write(fam       ,Triangles::amMesh);
-		if (famdba)   Th.Write(famdba    ,Triangles::amdbaMesh);
-		if (fftq)     Th.Write(fftq      ,Triangles::ftqMesh);
-		if (fmsh)     Th.Write(fmsh      ,Triangles::mshMesh);
-		if (fnopo)    Th.Write(fnopo     ,Triangles::NOPOMesh);
-
-		if ( ( rbb && wbb)  ||( rBB && wBB)){  // the code for interpolation 
-
-			if(verbosity >1)
-			{
-				if (rbb ) 
-					cout << " -- interpolation P1  files : " << rbb << " -> " << wbb <<endl;
-				if (rBB ) 
-					cout << " -- interpolation P1  files : " << rBB<< " -> " << wBB <<endl;
-			}
-			const int dim=2;
-			// optimisation read only si rbb != fMbb
-
-			double *solbb=0;
-			double *solBB=0;
-
-			if (rbb) 
-				solbb =  rbbeqMbb? solMbb : ReadbbFile(rbb,nbsolbb,lsolbb,2,2);
-			if (rBB) 
-				solBB =  rBBeqMBB? solMBB : ReadBBFile(rBB,nbsolBB,lsolBB,typesolsBB,2,2);
-
-			// cout << " " << rbbeqMbb << " " <<  sol << " " << nbsol << " " << lsol << endl; 
-			if (!solBB && !solbb ) 
-			{
-				cerr << " Fatal Error "  << "solBB =  " <<  solBB << " solbb= " << solbb << endl;
-				exit(2);}
-
-				ofstream *fbb = wbb ? new ofstream(wbb) :0;
-				ofstream *fBB = wBB ? new ofstream(wBB) :0;
-				Int4   nbfieldBB = 0, nbfieldbb = nbsolbb;
-				if (fbb)
-					*fbb  << dim << " " << nbsolbb << " " << Th.nbv <<" " << 2 << endl; 
-				if (fBB)
-				{
-					int i;
-					*fBB  << dim << " " << nbsolBB ;
-					for ( i=0;i<nbsolBB;i++)
-						*fBB << " " << (typesolsBB ?typesolsBB[i]+1:1) ;
-					*fBB << " " << Th.nbv <<" " << 2 << endl; 
-					assert(dim==2);
-					for ( i=0;i<nbsolBB;i++)
-						nbfieldBB += typesolsBB ? typesolsBB[i]+1 : 1;
-					// this code is good only if dim == 2 
-				}
-				cout << "nb of field BB " << nbfieldBB << endl;
-				//		 if(fBB) fBB->precision(15);
-				//if(fbb) fbb->precision(15);
-				for(i=0;i<Th.nbv;i++)
-				{
-					Int4 i0,i1,i2;
-					double a[3];
-					Icoor2 dete[3];
-					I2 I = Th.BTh.toI2(Th.vertices[i].r);
-					Triangle & tb = *Th.BTh.FindTriangleContening(I,dete);
-
-					if (tb.det>0) { // internal point 
-						a[0]= (Real8) dete[0]/ tb.det;
-						a[1]= (Real8) dete[1] / tb.det;
-						a[2] = (Real8) dete[2] / tb.det;
-						i0=Th.BTh.Number(tb[0]);
-						i1=Th.BTh.Number(tb[1]);
-						i2=Th.BTh.Number(tb[2]);}
-					else {
-						double aa,bb;
-
-						TriangleAdjacent  ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
-
-						int k = ta;
-						Triangle & tc = *(Triangle *)ta;
-						i0=Th.BTh.Number(tc[0]);
-						i1=Th.BTh.Number(tc[1]);
-						i2=Th.BTh.Number(tc[2]);
-						a[VerticesOfTriangularEdge[k][1]] =aa;
-						a[VerticesOfTriangularEdge[k][0]] = bb;
-						a[OppositeVertex[k]] = 1- aa -bb;}
-
-						Int4  ibb0 = nbfieldbb*i0;
-						Int4  ibb1 = nbfieldbb*i1;
-						Int4  ibb2 = nbfieldbb*i2;
-
-						Int4  iBB0 = nbfieldBB*i0;
-						Int4  iBB1 = nbfieldBB*i1;
-						Int4  iBB2 = nbfieldBB*i2;
-						Int4 j=0;
-
-						for ( j=0;j<nbfieldbb;j++) 
-							*fbb  << " " << ( a[0] * solbb[ibb0++] + a[1] * solbb[ibb1++] + a[2]* solbb[ibb2++]) ;
-
-						for (j=0;j<nbfieldBB;j++) 
-							*fBB  << " " << ( a[0] * solBB[iBB0++] + a[1] * solBB[iBB1++] + a[2]* solBB[iBB2++]) ;
-
-						if (fbb) *fbb  << endl;
-						if (fBB) *fBB  << endl;
-				}
-				if (fbb)
-					delete fbb  ; // close
-				if (fBB)
-					delete fBB  ; // close
-				if (solbb) 
-					delete [] solbb;
-				if (solBB) 
-					delete [] solBB;
-
-		}
-
-		if(verbosity>0) 
-		{
-			if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2) 
-				cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
-			if (Th.NbOfQuad ) 
-				cout  << " Nb Quadrilaterals = " << Th.NbOfQuad  ;
-			cout   << endl;
-		}
-
-		// cout << "delete &Th " ;
-		delete &Th;
-		//cout << " end " << endl;
-	}
-	/*Assign output pointers:*/
+
+	/*Assign output pointers*/
 	bamgmesh->numberofelements=nelout;
 	bamgmesh->numberofnodes=nodsout;
+	xfree((void**)&bamgmesh->index);
 	bamgmesh->index=elementsout;
+	xfree((void**)&bamgmesh->x);
 	bamgmesh->x=xout;
+	xfree((void**)&bamgmesh->y);
 	bamgmesh->y=yout;
 
+	/*No error return*/
 	return noerr;
-	/*}}}*/
-	}
+
 }
Index: /issm/trunk/src/c/objects/BamgArgs.h
===================================================================
--- /issm/trunk/src/c/objects/BamgArgs.h	(revision 2779)
+++ /issm/trunk/src/c/objects/BamgArgs.h	(revision 2780)
@@ -8,5 +8,9 @@
 struct BamgArgs{
 
-	char* geomfile;
+	int    iso;
+	double hmin;
+	double hmax;
+	double gradation;
+	double cutoff;
 
 };
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 2779)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 2780)
@@ -13,35 +13,65 @@
 bamg_mesh=struct();
 
-%get domainoutline
-domainfile=getfieldvalueerr(options,'domain');
-resolution=getfieldvalue(options,'resolution',5000);
-if ~exist(domainfile,'file')
-	error(['bamg error message: file ' domainfile ' not found ']);
-else
-	domain=expread(domainfile);
+% Bamg Geometry parameters {{{1
+bamg_geometry.NumVertices=0;
+bamg_geometry.Vertices=zeros(0,3);
+bamg_geometry.NumEdges=0;
+bamg_geometry.Edges=zeros(0,3);
+bamg_geometry.hVertices=zeros(0,1);
+bamg_geometry.NumSubDomain=0;
+bamg_geometry.SubDomain=zeros(0,3);
+if exist(options,'domain'),
+	domainfile=getfieldvalueerr(options,'domain');
+	if ~exist(domainfile,'file')
+		error(['bamg error message: file ' domainfile ' not found ']);
+	else
+		domain=expread(domainfile);
 
-	%build geometry
-	count=0;
-	bamg_geometry.Vertices=zeros(0,3);
-	bamg_geometry.Edges=zeros(0,3);
-	bamg_geometry.hVertices=zeros(0,1);
-	bamg_geometry.SubDomain=zeros(0,3);
-	for i=1:length(domain),
-		nods=domain(i).nods-1; %the domain are closed 1=end;
-		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,
-			bamg_geometry.SubDomain=[2 count+1 1 1];
+		%build geometry
+		count=0;
+
+		for i=1:length(domain),
+			nods=domain(i).nods-1; %the domain are closed 1=end;
+			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;
+				bamg_geometry.SubDomain=[2 count+1 clockwise 1];
+			end
+			count=count+nods;
 		end
-		count=count+nods;
-	end
-	bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
-	bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
-	bamg_geometry.NumSubDomain=size(bamg_geometry.SubDomain,1);
+		bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
+		bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
+		bamg_geometry.NumSubDomain=size(bamg_geometry.SubDomain,1);
+	end 
 end
+%}}}
+
+% Bamg Mesh parameters {{{1
+bamg_mesh.numberofelements=0;
+bamg_mesh.numberofnodes=0;
+bamg_mesh.x=[];
+bamg_mesh.y=[];
+bamg_mesh.index=zeros(0,3);
+if (~exist(options,'domain') & md.numberofgrids~=0 & strcmpi(md.type,'2d')),
+	bamg_mesh.numberofelements=md.numberofelements;
+	bamg_mesh.numberofnodes=md.numberofgrids;
+	bamg_mesh.x=md.x;
+	bamg_mesh.y=md.y;
+	bamg_mesh.index=md.elements;
+end
+%}}}
+
+% Bamg Options {{{1
+bamg_options.iso=getfieldvalue(options,'iso',0);
+bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
+bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
+bamg_options.gradation=getfieldvalue(options,'gradation',1.2);
+bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
+%}}}
 
 %call Bamg
-[elements,x,y]=Bamg(bamg_geometry,bamg_options);
+[elements,x,y]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
 
 % plug results onto model
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2779)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2780)
@@ -13,6 +13,12 @@
 	BamgGeom bamggeom;
 
-	/*inputs: */
-	char*   geomfile=NULL;
+	/*Mesh inputs*/
+	int     numberofnodes;
+	int     numberofelements;
+	double* x;
+	double* y;
+	double* index;
+
+	/*Geom inputs: */
 	int     NumVertices;
 	double* Vertices=NULL;
@@ -23,4 +29,9 @@
 	double* SubDomain=NULL;
 
+	/*Options inputs*/
+	int    iso;
+	double hmin,hmax;
+	double gradation;
+	double cutoff;
 
 	/*Boot module: */
@@ -30,20 +41,14 @@
 	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&BamgUsage);
 
-	/*process inputs*/
-	//FetchData(&geomfile,mxGetField(BAMGOPTIONS,0,"fgeom"));
-
+	/*create bamg geometry input*/
 	FetchData(&NumVertices,mxGetField(BAMGGEOMETRY,0,"NumVertices"));
+	bamggeom.NumVertices=NumVertices;
 	FetchData(&Vertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Vertices"));
+	bamggeom.Vertices=Vertices;
 	FetchData(&NumEdges,mxGetField(BAMGGEOMETRY,0,"NumEdges"));
+	bamggeom.NumEdges=NumEdges;
 	FetchData(&Edges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges"));
+	bamggeom.Edges=Edges;
 	FetchData(&hVertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"hVertices"));
-	FetchData(&NumSubDomain,mxGetField(BAMGGEOMETRY,0,"NumSubDomain"));
-	FetchData(&SubDomain,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"SubDomain"));
-
-	/*create bamg geometry input*/
-	bamggeom.NumVertices=NumVertices;
-	bamggeom.Vertices=Vertices;
-	bamggeom.NumEdges=NumEdges;
-	bamggeom.Edges=Edges;
 	bamggeom.hVertices=hVertices;
 	bamggeom.Edges=Edges;
@@ -59,11 +64,32 @@
 	bamggeom.NumRequiredEdges=0;
 	bamggeom.RequiredEdges=NULL;
+	FetchData(&NumSubDomain,mxGetField(BAMGGEOMETRY,0,"NumSubDomain"));
 	bamggeom.NumSubDomain=NumSubDomain;
+	FetchData(&SubDomain,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"SubDomain"));
 	bamggeom.SubDomain=SubDomain;
 
 	/*create bamg mesh input*/
-	//bamgargs.geomfile=geomfile;
+	FetchData(&numberofnodes,mxGetField(BAMGMESH,0,"numberofnodes"));
+	bamgmesh.numberofnodes=numberofnodes;
+	FetchData(&numberofelements,mxGetField(BAMGMESH,0,"numberofelements"));
+	bamgmesh.numberofelements=numberofelements;
+	FetchData(&x,NULL,NULL,mxGetField(BAMGMESH,0,"x"));
+	bamgmesh.x=x;
+	FetchData(&y,NULL,NULL,mxGetField(BAMGMESH,0,"y"));
+	bamgmesh.y=y;
+	FetchData(&index,NULL,NULL,mxGetField(BAMGMESH,0,"index"));
+	bamgmesh.index=index;
 
-	/*create bamg mesh input*/
+	/*create bamg options input*/
+	FetchData(&iso,mxGetField(BAMGOPTIONS,0,"iso"));
+	bamgargs.iso=iso;
+	FetchData(&hmin,mxGetField(BAMGOPTIONS,0,"hmin"));
+	bamgargs.hmin=hmin;
+	FetchData(&hmax,mxGetField(BAMGOPTIONS,0,"hmax"));
+	bamgargs.hmax=hmax;
+	FetchData(&gradation,mxGetField(BAMGOPTIONS,0,"gradation"));
+	bamgargs.gradation=gradation;
+	FetchData(&cutoff,mxGetField(BAMGOPTIONS,0,"cutoff"));
+	bamgargs.cutoff=cutoff;
 
 	/*!Generate internal degree of freedom numbers: */
@@ -76,5 +102,4 @@
 
 	/*Free ressources: */
-	//xfree((void**)&geomfile);
 	xfree((void**)&Vertices);
 	xfree((void**)&Edges);
@@ -89,5 +114,5 @@
 {
 	_printf_("\n");
-	_printf_("   usage: [elements,x,y]=%s(elements,x,y,metric,gradation,splitpbedge,nbv);\n",__FUNCT__);
+	_printf_("   usage: [elements,x,y]=%s(bamgmesh,bamggeom,bamgoptions,nbv);\n",__FUNCT__);
 	_printf_("\n");
 }
Index: /issm/trunk/src/mex/Bamg/Bamg.h
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.h	(revision 2779)
+++ /issm/trunk/src/mex/Bamg/Bamg.h	(revision 2780)
@@ -18,6 +18,7 @@
     
 /* serial input macros: */
-#define BAMGGEOMETRY (mxArray*)prhs[0]
-#define BAMGOPTIONS  (mxArray*)prhs[1]
+#define BAMGMESH     (mxArray*)prhs[0]
+#define BAMGGEOMETRY (mxArray*)prhs[1]
+#define BAMGOPTIONS  (mxArray*)prhs[2]
 
 /* serial output macros: */
@@ -30,5 +31,5 @@
 #define NLHS  3
 #undef NRHS
-#define NRHS  2
+#define NRHS  3
 
 #endif  /* _BAMG_H */