Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3278)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3279)
@@ -53,9 +53,4 @@
 	hmax=bamgopts->hmax;
 	verbosity=bamgopts->verbose;
-
-	// some verification
-	if ( maxsubdiv > 10 || maxsubdiv <= 1.0){
-		throw ErrorException(__FUNCT__,exprintf("maxsubdiv (%g) should be in ]1,10]",maxsubdiv));
-	}
 
 	// no metric -> no smoothing 
@@ -115,6 +110,6 @@
 		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());
+		hmin=Max(hmin,BTh.MinimalHmin());
+		hmax=Min(hmax,BTh.MaximalHmax());
 
 		//Make Quadtree from background mesh
@@ -138,4 +133,5 @@
 					}
 				}
+
 			}
 
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3278)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3279)
@@ -3,4 +3,5 @@
 #include <cmath>
 #include <ctime>
+#include <assert.h>
 
 #include "BamgObjects.h"
@@ -227,5 +228,4 @@
 
 		double Hmin = HUGE_VAL;// the infinie value 
-		long hvertices=0;
 		int i,j,k,n,i1,i2;
 
@@ -290,9 +290,7 @@
 			edges = new GeometricalEdge[nbe];
 
-			//if hvertices==0, initialize len (length of each edge)
-			if (!hvertices) {
-				len = new double[nbv];
-				for(i=0;i<nbv;i++) len[i]=0;
-			}
+			//initialize len (length of each edge)
+			len = new double[nbv];
+			for(i=0;i<nbv;i++) len[i]=0;
 
 			for (i=0;i<nbe;i++){
@@ -316,21 +314,20 @@
 				edges[i].Adj[0] = edges[i].Adj[1] = 0;
 				edges[i].flag = 0;
-				if (!hvertices) {
-					vertices[i1].color++;
-					vertices[i2].color++;
-					len[i1] += l12;
-					len[i2] += l12;
-				}
+
+				//prepare metric
+				vertices[i1].color++;
+				vertices[i2].color++;
+				len[i1] += l12;
+				len[i2] += l12;
 			}
 
 			// definition  the default of the given mesh size 
-			if (!hvertices){
-				for (i=0;i<nbv;i++) 
-				 if (vertices[i].color > 0) 
-				  vertices[i].m=Metric(len[i] /(double) vertices[i].color);
-				 else 
-				  vertices[i].m=Metric(Hmin);
-				delete [] len;
-			}
+			for (i=0;i<nbv;i++) 
+			 if (vertices[i].color > 0) 
+			  vertices[i].m=Metric(len[i] /(double) vertices[i].color);
+			 else 
+			  vertices[i].m=Metric(Hmin);
+			delete [] len;
+			
 		}
 		else{
@@ -353,5 +350,7 @@
 			if(verbose>5) printf("      processing hVertices\n");
 			for (i=0;i< nbv;i++){
-				vertices[i].m=Metric((double)bamggeom->hVertices[i]);
+				if (!isnan(bamggeom->hVertices[i])){
+					vertices[i].m=Metric((double)bamggeom->hVertices[i]);
+				}
 			}
 		}
@@ -363,5 +362,4 @@
 		if(bamggeom->MetricVertices){
 			if(verbose>5) printf("      processing MetricVertices\n");
-			hvertices=1;
 			for (i=0;i< nbv;i++) {
 				vertices[i].m = Metric((double)bamggeom->MetricVertices[i*3+0],(double)bamggeom->MetricVertices[i*3+1],(double)bamggeom->MetricVertices[i*3+2]);
@@ -376,5 +374,4 @@
 			if(verbose>5) printf("      processing h1h2VpVertices\n");
 			double h1,h2,v1,v2;
-			hvertices =1;
 			for (i=0;i< nbv;i++) {
 				h1=(double)bamggeom->MetricVertices[i*4+0];
@@ -438,5 +435,5 @@
 		//RequiredVertices
 		if(bamggeom->RequiredVertices){
-			if(verbose>5) printf("      processing RequiredVertices");
+			if(verbose>5) printf("      processing RequiredVertices\n");
 			n=bamggeom->NumRequiredVertices;
 			for (i=0;i<n;i++) {     
@@ -453,5 +450,5 @@
 		//RequiredEdges
 		if(bamggeom->RequiredEdges){
-			if(verbose>5) printf("      processing RequiredEdges");
+			if(verbose>5) printf("      processing RequiredEdges\n");
 			n=bamggeom->NumRequiredEdges;
 			for (i=0;i<n;i++) {     
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3278)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3279)
@@ -261,7 +261,5 @@
 		long i1,i2,i3,iref;
 		long i,j;
-		long hvertices =0;
 		Metric M1(1);
-		int Version=1,dim=2;
 		int verbose=0;
 
@@ -311,8 +309,6 @@
 		long i1,i2,i3,iref;
 		long i,j;
-		long hvertices =0;
 		long ifgeom=0;
 		Metric M1(1);
-		int Version=1,dim=2;
 
 		verbose=bamgopts->verbose;
@@ -387,18 +383,4 @@
 		}
 
-		//hVertices
-		if(bamgmesh->hVertices){
-			if(verbose>5) printf("      processing hVertices\n");
-			hvertices=1;
-			for (i=0;i< nbv;i++){
-				if (!isnan(bamgmesh->hVertices[i])){
-					vertices[i].m=Metric((double)bamgmesh->hVertices[i]);
-				}
-			}
-		}
-		else{
-			if(verbose>5) printf("      no hVertices found\n");
-		}
-
 		//VerticesOnGeometricEdge
 		if(bamgmesh->VerticesOnGeometricEdge){
@@ -443,14 +425,12 @@
 			nbe=bamgmesh->NumEdges;
 			edges= new Edge[nbe];
-
-			if (!hvertices) {
-				len= new double[nbv];
-				for(i=0;i<nbv;i++)
-				 len[i]=0;
-			}
+			//initialize length of each edge (used to provided metric)
+			len= new double[nbv];
+			for(i=0;i<nbv;i++) len[i]=0;
 
 			for (i=0;i<nbe;i++){
 				i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing
 				i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing
+				edges[i].ref=(long)bamgmesh->Edges[i*3+2];
 				edges[i].v[0]= vertices +i1;
 				edges[i].v[1]= vertices +i2;
@@ -460,26 +440,24 @@
 				double l12=sqrt((x12,x12));
 
-				if (!hvertices){
-					vertices[i1].color++;
-					vertices[i2].color++;
-					len[i1]+=l12;
-					len[i2]+=l12;
-				}
+				//prepare metric
+				vertices[i1].color++;
+				vertices[i2].color++;
+				len[i1]+=l12;
+				len[i2]+=l12;
 				Hmin = Min(Hmin,l12);
 			}
 
 			// definition  the default of the given mesh size 
-			if (!hvertices){
-				for (i=0;i<nbv;i++) 
-				 if (vertices[i].color > 0) 
-				  vertices[i].m=Metric(len[i] /(double) vertices[i].color);
-				 else 
-				  vertices[i].m=Metric(Hmin);
-				delete [] len;
-			}
+			for (i=0;i<nbv;i++){
+				if (vertices[i].color>0) 
+				 vertices[i].m=Metric(len[i]/(double)vertices[i].color);
+				else 
+				 vertices[i].m=Metric(Hmin);
+			}
+			delete [] len;
 
 			// construction of edges[].adj 
 			for (i=0;i<nbv;i++){ 
-				vertices[i].color = (vertices[i].color ==2) ? -1 : -2;
+				vertices[i].color=(vertices[i].color ==2) ?-1:-2;
 			}
 			for (i=0;i<nbe;i++){
@@ -537,4 +515,17 @@
 		else{
 			if(verbose>5) printf("      no EdgesOnGeometricEdge found\n");
+		}
+
+		//hVertices
+		if(bamgmesh->hVertices){
+			if(verbose>5) printf("      processing hVertices\n");
+			for (i=0;i< nbv;i++){
+				if (!isnan(bamgmesh->hVertices[i])){
+					vertices[i].m=Metric((double)bamgmesh->hVertices[i]);
+				}
+			}
+		}
+		else{
+			if(verbose>5) printf("      no hVertices found\n");
 		}
 
@@ -1586,5 +1577,4 @@
 
 		/*Options*/
-		const int dim = 2;
 		double* s=NULL;
 		long    nbsol;
@@ -1790,5 +1780,4 @@
 
 		/*Options*/
-		const int dim = 2;
 		double* s=NULL;
 		long nbsol;
@@ -2622,310 +2611,4 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::ReconstructExistingMesh{{{1*/
-	void Triangles::ReconstructExistingMesh() {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
-
-		/*This routine reconstruct an existing mesh to make it CONVEX:
-		 * -all the holes are filled
-		 * -concave boundaries are filled
-		 * A convex mesh is required for a lot of operations. This is why every mesh
-		 * goes through this process.
-		 * This routine also generates mesh properties such as adjencies,...
-		 */
-
-		/*Intermediary*/
-		int verbosity=0;
-
-		// generation of the integer coordinate
-
-		// find extrema coordinates of vertices pmin,pmax
-		long i;
-		if(verbosity>2) printf("      Reconstruct mesh of %i vertices\n",nbv); 
-
-		//initialize ordre
-		assert(ordre);
-		for (i=0;i<nbv;i++) ordre[i]=0;
-
-		//Initialize NbSubDomains
-		NbSubDomains =0;
-
-		/* generation of triangles adjacency*/
-
-		//First add existing edges
-		long kk =0;
-		SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
-		for (i=0;i<nbe;i++){
-			kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
-		}
-		if (kk != nbe){ 
-			throw ErrorException(__FUNCT__,exprintf("There are %i double edges in the mesh",kk-nbe));
-		}
-
-		//Add edges of all triangles in existing mesh
-		long* st = new long[nbt*3];
-		for (i=0;i<nbt*3;i++) st[i]=-1;
-		for (i=0;i<nbt;i++){
-			for (int j=0;j<3;j++){
-
-				//Add current triangle edge to edge4
-				long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
-
-				long invisible=triangles[i].Hidden(j);
-
-				//If the edge has not been added to st, add it
-				if(st[k]==-1) st[k]=3*i+j;
-				
-				//If the edge already exists, add adjacency
-				else if(st[k]>=0) {
-					assert(!triangles[i].TriangleAdj(j));
-					assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
-
-					triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
-					if (invisible) triangles[i].SetHidden(j);
-					if (k<nbe)     triangles[i].SetLocked(j);
-
-					//Make st[k] negative so that it will throw an error message if it is found again
-					st[k]=-2-st[k]; 
-				}
-
-				//An edge belongs to 2 triangles
-				else {
-					throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",
-									Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));
-				}
-			}
-		}
-
-		//Display info if required
-		if(verbosity>5) {
-			printf("         info of Mesh:\n");
-			printf("            - number of vertices    = %i \n",nbv); 
-			printf("            - number of triangles   = %i \n",nbt); 
-			printf("            - number of given edges = %i \n",nbe); 
-			printf("            - number of all edges   = %i \n"  ,edge4->nb()); 
-			printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv); 
-		}
-
-		//check the consistency of edge[].adj and the geometrical required vertex
-		long k=0;
-		for (i=0;i<edge4->nb();i++){
-			if (st[i]>=0){ // edge alone 
-				if (i<nbe){
-					long i0=edge4->i(i);
-					ordre[i0] = vertices+i0;
-					long i1=edge4->j(i);
-					ordre[i1] = vertices+i1;
-				}
-				else {
-					k=k+1;
-					if (k<10) {
-						//print only 10 edges
-						printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
-					}
-					else if (k==10){
-						printf("Other lost boundary edges not shown...\n");
-					}
-				}
-			}
-		}
-		if(k) {
-			throw ErrorException(__FUNCT__,exprintf("%i boundary edges (from the geometry) are not defined as mesh edges",k));
-		}
-
-		/* mesh generation with boundary points*/
-		long nbvb=0;
-		for (i=0;i<nbv;i++){ 
-			vertices[i].t=0;
-			vertices[i].vint=0;
-			if (ordre[i]) ordre[nbvb++]=ordre[i];
-		}
-
-		Triangle* savetriangles=triangles;
-		long savenbt=nbt;
-		long savenbtx=nbtx;
-		SubDomain* savesubdomains=subdomains;
-		subdomains=0;
-
-		long  Nbtriafillhole=2*nbvb;
-		Triangle* triafillhole=new Triangle[Nbtriafillhole];
-		triangles = triafillhole;
-
-		nbt=2;
-		nbtx= Nbtriafillhole;
-
-		//Find a vertex that is not aligned with vertices 0 and 1
-		for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;) 
-		 if  (++i>=nbvb) {
-			 throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned"));
-		 }
-		//Move this vertex (i) to the 2d position in ordre
-		Exchange(ordre[2], ordre[i]);
-
-		/*Reconstruct mesh beginning with 2 triangles*/
-		Vertex *  v0=ordre[0], *v1=ordre[1];
-
-		triangles[0](0) = NULL; // Infinite vertex
-		triangles[0](1) = v0;
-		triangles[0](2) = v1;
-
-		triangles[1](0) = NULL;// Infinite vertex
-		triangles[1](2) = v0;
-		triangles[1](1) = v1;
-		const int e0 = OppositeEdge[0];
-		const int e1 = NextEdge[e0];
-		const int e2 = PreviousEdge[e0];
-		triangles[0].SetAdj2(e0, &triangles[1] ,e0);
-		triangles[0].SetAdj2(e1, &triangles[1] ,e2);
-		triangles[0].SetAdj2(e2, &triangles[1] ,e1);
-
-		triangles[0].det = -1;  // boundary triangles
-		triangles[1].det = -1;  // boundary triangles
-
-		triangles[0].SetTriangleContainingTheVertex();
-		triangles[1].SetTriangleContainingTheVertex();
-
-		triangles[0].link=&triangles[1];
-		triangles[1].link=&triangles[0];
-
-		if (!quadtree) delete quadtree; //ReInitialise;
-		quadtree = new QuadTree(this,0);
-		quadtree->Add(*v0);
-		quadtree->Add(*v1);
-
-		// vertices are added one by one
-		long NbSwap=0;
-		for (int icount=2; icount<nbvb; icount++) {
-			Vertex *vi  = ordre[icount];
-			Icoor2 dete[3];
-			Triangle *tcvi = FindTriangleContaining(vi->i,dete);
-			quadtree->Add(*vi); 
-			AddVertex(*vi,tcvi,dete);
-			NbSwap += vi->Optim(1,1);
-		}
-
-		//enforce the boundary 
-		TriangleAdjacent ta(0,0);
-		long nbloss = 0,knbe=0;
-		for ( i = 0; i < nbe; i++){
-			if (st[i] >=0){ //edge alone => on border
-				Vertex &a=edges[i][0], &b=edges[i][1];
-				if (a.t && b.t){
-					knbe++;
-					if (ForceEdge(a,b,ta)<0) nbloss++;
-				}
-			}
-		}
-		if(nbloss) {
-			throw ErrorException(__FUNCT__,exprintf("we lost %i existing edges other %i",nbloss,knbe));
-		}
-
-		FindSubDomain(1);
-		// remove all the hole 
-		// remove all the good sub domain
-		long krm =0;
-		for (i=0;i<nbt;i++){
-			if (triangles[i].link){ // remove triangles
-				krm++;
-				for (int j=0;j<3;j++){
-					TriangleAdjacent ta =  triangles[i].Adj(j);
-					Triangle &tta = *(Triangle*)ta;
-					//if edge between remove and not remove 
-					if(! tta.link){ 
-						// change the link of ta;
-						int ja = ta;
-						Vertex *v0= ta.EdgeVertex(0);
-						Vertex *v1= ta.EdgeVertex(1);
-						long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
-
-						assert(st[k]>=0);
-						tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
-						ta.SetLock();
-						st[k]=-2-st[k]; 
-					}
-				}
-			}
-		}
-		long NbTfillHoll =0;
-		for (i=0;i<nbt;i++){
-			if (triangles[i].link) {
-				triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
-				triangles[i].color=-1;
-			}
-			else{
-				triangles[i].color= savenbt+ NbTfillHoll++;
-			}
-		}
-		assert(savenbt+NbTfillHoll<=savenbtx);
-
-		// copy of the outside triangles in saveTriangles 
-		for (i=0;i<nbt;i++){
-			if(triangles[i].color>=0) {
-				savetriangles[savenbt]=triangles[i];
-				savetriangles[savenbt].link=0;
-				savenbt++;
-			}
-		}
-		// gestion of the adj
-		k =0;
-		Triangle * tmax = triangles + nbt;
-		for (i=0;i<savenbt;i++) { 
-			Triangle & ti = savetriangles[i];
-			for (int j=0;j<3;j++){
-				Triangle * ta = ti.TriangleAdj(j);
-				int aa = ti.NuEdgeTriangleAdj(j);
-				int lck = ti.Locked(j);
-				if (!ta) k++; // bug 
-				else if ( ta >= triangles && ta < tmax){
-					ta= savetriangles + ta->color;
-					ti.SetAdj2(j,ta,aa);
-					if(lck) ti.SetLocked(j);
-				}
-			}
-		}
-
-		// restore triangles;
-		nbt=savenbt;
-		nbtx=savenbtx;
-		delete [] triangles;
-		delete [] subdomains;
-		triangles = savetriangles;
-		subdomains = savesubdomains;
-		if (k) {
-			throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));
-		}
-		FindSubDomain();
-
-		delete edge4;
-		delete [] st;
-		for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
-
-		SetVertexFieldOn();
-
-		/*Check requirements consistency*/
-		for (i=0;i<nbe;i++){
-			/*If the current mesh edge is on Geometry*/
-			if(edges[i].onGeometry){
-				for(int j=0;j<2;j++){
-					/*Go through the edges adjacent to current edge (if on the same curve)*/
-					if (!edges[i].adj[j]){
-						/*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
-						/*Check that the 2 vertices are on geometry AND required*/
-						if(!edges[i][j].onGeometry->IsRequiredVertex()){
-							printf("ReconstructExistingMesh error message: problem with the edge %i: [%i][%i]\n",Number(edges[i][j])+1,i,j);
-							printf("   This edge is on geometry (Edge of geometry %i)",Gh.Number(edges[i].onGeometry));
-							if (edges[i][j].onGeometry->OnGeomVertex())
-							 printf(" Vertex %i\n",Gh.Number(edges[i][j].onGeometry->gv));
-							else if (edges[i][j].onGeometry->OnGeomEdge())
-							 printf(" Edges %i\n",Gh.Number(edges[i][j].onGeometry->ge));
-							else
-							 printf(" = %p\n",edges[i][j].onGeometry);
-							throw ErrorException(__FUNCT__,exprintf("See above (might be cryptic...)"));
-						}
-					}
-				}
-			}
-		}
-	}
-	/*}}}1*/
 	/*FUNCTION Triangles::GeomToTriangles0{{{1*/
 	void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){
@@ -3893,12 +3576,9 @@
 		long it,nbchange=0;    
 		double lmax=0;
-		for (it=0;it<nbt;it++)
-		  {
+		for (it=0;it<nbt;it++){
 			Triangle &t=triangles[it];
-			for (int j=0;j<3;j++)
-			  {
+			for (int j=0;j<3;j++){
 				Triangle &tt = *t.TriangleAdj(j);
-				if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)) 
-				  {
+				if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)){
 					Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];
 					Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];
@@ -3907,6 +3587,6 @@
 					double l = M(AB,AB);
 					lmax = Max(lmax,l);
-					if(l> maxsubdiv2)
-					  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
+					if(l> maxsubdiv2){
+					  R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
 						double lc = M(AC,AC);
 						D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
@@ -3916,10 +3596,10 @@
 						v0.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
 						nbchange++;
-					  }
+					}
 					M = v1;
 					l = M(AB,AB);
 					lmax = Max(lmax,l);
-					if(l> maxsubdiv2)
-					  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
+					if(l> maxsubdiv2){
+					  R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
 						double lc = M(AC,AC);
 						D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
@@ -3929,10 +3609,8 @@
 						v1.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
 						nbchange++;
-					  }
-
-
-				  }
-			  }
-		  }
+					}
+				}
+			}
+		}
 		if(verbosity>3){
 			printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
@@ -4285,6 +3963,314 @@
 	}                  
 	/*}}}1*/
+/*FUNCTION Triangles::ReconstructExistingMesh{{{1*/
+void Triangles::ReconstructExistingMesh(){
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
+
+	/*This routine reconstruct an existing mesh to make it CONVEX:
+	 * -all the holes are filled
+	 * -concave boundaries are filled
+	 * A convex mesh is required for a lot of operations. This is why every mesh
+	 * goes through this process.
+	 * This routine also generates mesh properties such as adjencies,...
+	 */
+
+	/*Intermediary*/
+	int verbosity=0;
+
+	// generation of the integer coordinate
+
+	// find extrema coordinates of vertices pmin,pmax
+	long i;
+	if(verbosity>2) printf("      Reconstruct mesh of %i vertices\n",nbv); 
+
+	//initialize ordre
+	assert(ordre);
+	for (i=0;i<nbv;i++) ordre[i]=0;
+
+	//Initialize NbSubDomains
+	NbSubDomains =0;
+
+	/* generation of triangles adjacency*/
+
+	//First add existing edges
+	long kk =0;
+	SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
+	for (i=0;i<nbe;i++){
+		kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
+	}
+	if (kk != nbe){ 
+		throw ErrorException(__FUNCT__,exprintf("There are %i double edges in the mesh",kk-nbe));
+	}
+
+	//Add edges of all triangles in existing mesh
+	long* st = new long[nbt*3];
+	for (i=0;i<nbt*3;i++) st[i]=-1;
+	for (i=0;i<nbt;i++){
+		for (int j=0;j<3;j++){
+
+			//Add current triangle edge to edge4
+			long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
+
+			long invisible=triangles[i].Hidden(j);
+
+			//If the edge has not been added to st, add it
+			if(st[k]==-1) st[k]=3*i+j;
+
+			//If the edge already exists, add adjacency
+			else if(st[k]>=0) {
+				assert(!triangles[i].TriangleAdj(j));
+				assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
+
+				triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
+				if (invisible) triangles[i].SetHidden(j);
+				if (k<nbe)     triangles[i].SetLocked(j);
+
+				//Make st[k] negative so that it will throw an error message if it is found again
+				st[k]=-2-st[k]; 
+			}
+
+			//An edge belongs to 2 triangles
+			else {
+				throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",
+								Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));
+			}
+		}
+	}
+
+	//Display info if required
+	if(verbosity>5) {
+		printf("         info of Mesh:\n");
+		printf("            - number of vertices    = %i \n",nbv); 
+		printf("            - number of triangles   = %i \n",nbt); 
+		printf("            - number of given edges = %i \n",nbe); 
+		printf("            - number of all edges   = %i \n"  ,edge4->nb()); 
+		printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv); 
+	}
+
+	//check the consistency of edge[].adj and the geometrical required vertex
+	long k=0;
+	for (i=0;i<edge4->nb();i++){
+		if (st[i]>=0){ // edge alone 
+			if (i<nbe){
+				long i0=edge4->i(i);
+				ordre[i0] = vertices+i0;
+				long i1=edge4->j(i);
+				ordre[i1] = vertices+i1;
+			}
+			else {
+				k=k+1;
+				if (k<10) {
+					//print only 10 edges
+					printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
+				}
+				else if (k==10){
+					printf("Other lost boundary edges not shown...\n");
+				}
+			}
+		}
+	}
+	if(k) {
+		throw ErrorException(__FUNCT__,exprintf("%i boundary edges (from the geometry) are not defined as mesh edges",k));
+	}
+
+	/* mesh generation with boundary points*/
+	long nbvb=0;
+	for (i=0;i<nbv;i++){ 
+		vertices[i].t=0;
+		vertices[i].vint=0;
+		if (ordre[i]) ordre[nbvb++]=ordre[i];
+	}
+
+	Triangle* savetriangles=triangles;
+	long savenbt=nbt;
+	long savenbtx=nbtx;
+	SubDomain* savesubdomains=subdomains;
+	subdomains=0;
+
+	long  Nbtriafillhole=2*nbvb;
+	Triangle* triafillhole=new Triangle[Nbtriafillhole];
+	triangles = triafillhole;
+
+	nbt=2;
+	nbtx= Nbtriafillhole;
+
+	//Find a vertex that is not aligned with vertices 0 and 1
+	for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;) 
+	 if  (++i>=nbvb) {
+		 throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned"));
+	 }
+	//Move this vertex (i) to the 2d position in ordre
+	Exchange(ordre[2], ordre[i]);
+
+	/*Reconstruct mesh beginning with 2 triangles*/
+	Vertex *  v0=ordre[0], *v1=ordre[1];
+
+	triangles[0](0) = NULL; // Infinite vertex
+	triangles[0](1) = v0;
+	triangles[0](2) = v1;
+
+	triangles[1](0) = NULL;// Infinite vertex
+	triangles[1](2) = v0;
+	triangles[1](1) = v1;
+	const int e0 = OppositeEdge[0];
+	const int e1 = NextEdge[e0];
+	const int e2 = PreviousEdge[e0];
+	triangles[0].SetAdj2(e0, &triangles[1] ,e0);
+	triangles[0].SetAdj2(e1, &triangles[1] ,e2);
+	triangles[0].SetAdj2(e2, &triangles[1] ,e1);
+
+	triangles[0].det = -1;  // boundary triangles
+	triangles[1].det = -1;  // boundary triangles
+
+	triangles[0].SetTriangleContainingTheVertex();
+	triangles[1].SetTriangleContainingTheVertex();
+
+	triangles[0].link=&triangles[1];
+	triangles[1].link=&triangles[0];
+
+	if (!quadtree) delete quadtree; //ReInitialise;
+	quadtree = new QuadTree(this,0);
+	quadtree->Add(*v0);
+	quadtree->Add(*v1);
+
+	// vertices are added one by one
+	long NbSwap=0;
+	for (int icount=2; icount<nbvb; icount++) {
+		Vertex *vi  = ordre[icount];
+		Icoor2 dete[3];
+		Triangle *tcvi = FindTriangleContaining(vi->i,dete);
+		quadtree->Add(*vi); 
+		AddVertex(*vi,tcvi,dete);
+		NbSwap += vi->Optim(1,1);
+	}
+
+	//enforce the boundary 
+	TriangleAdjacent ta(0,0);
+	long nbloss = 0,knbe=0;
+	for ( i = 0; i < nbe; i++){
+		if (st[i] >=0){ //edge alone => on border
+			Vertex &a=edges[i][0], &b=edges[i][1];
+			if (a.t && b.t){
+				knbe++;
+				if (ForceEdge(a,b,ta)<0) nbloss++;
+			}
+		}
+	}
+	if(nbloss) {
+		throw ErrorException(__FUNCT__,exprintf("we lost %i existing edges other %i",nbloss,knbe));
+	}
+
+	FindSubDomain(1);
+	// remove all the hole 
+	// remove all the good sub domain
+	long krm =0;
+	for (i=0;i<nbt;i++){
+		if (triangles[i].link){ // remove triangles
+			krm++;
+			for (int j=0;j<3;j++){
+				TriangleAdjacent ta =  triangles[i].Adj(j);
+				Triangle &tta = *(Triangle*)ta;
+				//if edge between remove and not remove 
+				if(! tta.link){ 
+					// change the link of ta;
+					int ja = ta;
+					Vertex *v0= ta.EdgeVertex(0);
+					Vertex *v1= ta.EdgeVertex(1);
+					long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
+
+					assert(st[k]>=0);
+					tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
+					ta.SetLock();
+					st[k]=-2-st[k]; 
+				}
+			}
+		}
+	}
+	long NbTfillHoll =0;
+	for (i=0;i<nbt;i++){
+		if (triangles[i].link) {
+			triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
+			triangles[i].color=-1;
+		}
+		else{
+			triangles[i].color= savenbt+ NbTfillHoll++;
+		}
+	}
+	assert(savenbt+NbTfillHoll<=savenbtx);
+
+	// copy of the outside triangles in saveTriangles 
+	for (i=0;i<nbt;i++){
+		if(triangles[i].color>=0) {
+			savetriangles[savenbt]=triangles[i];
+			savetriangles[savenbt].link=0;
+			savenbt++;
+		}
+	}
+	// gestion of the adj
+	k =0;
+	Triangle * tmax = triangles + nbt;
+	for (i=0;i<savenbt;i++) { 
+		Triangle & ti = savetriangles[i];
+		for (int j=0;j<3;j++){
+			Triangle * ta = ti.TriangleAdj(j);
+			int aa = ti.NuEdgeTriangleAdj(j);
+			int lck = ti.Locked(j);
+			if (!ta) k++; // bug 
+			else if ( ta >= triangles && ta < tmax){
+				ta= savetriangles + ta->color;
+				ti.SetAdj2(j,ta,aa);
+				if(lck) ti.SetLocked(j);
+			}
+		}
+	}
+
+	// restore triangles;
+	nbt=savenbt;
+	nbtx=savenbtx;
+	delete [] triangles;
+	delete [] subdomains;
+	triangles = savetriangles;
+	subdomains = savesubdomains;
+	if (k) {
+		throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));
+	}
+	FindSubDomain();
+
+	delete edge4;
+	delete [] st;
+	for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
+
+	SetVertexFieldOn();
+
+	/*Check requirements consistency*/
+	for (i=0;i<nbe;i++){
+ 	/*If the current mesh edge is on Geometry*/
+		if(edges[i].onGeometry){
+			for(int j=0;j<2;j++){
+				/*Go through the edges adjacent to current edge (if on the same curve)*/
+				if (!edges[i].adj[j]){
+					/*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
+					/*Check that the 2 vertices are on geometry AND required*/
+					if(!edges[i][j].onGeometry->IsRequiredVertex()){
+						printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1);
+						printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1);
+						if (edges[i][j].onGeometry->OnGeomVertex())
+						 printf("the vertex number %i of this edge is a geometric Vertex number %i\n",Gh.Number(edges[i][j].onGeometry->gv)+1);
+						else if (edges[i][j].onGeometry->OnGeomEdge())
+						 printf("the vertex number %i of this edge is a geometric Edge number %i\n",Gh.Number(edges[i][j].onGeometry->ge)+1);
+						else
+						 printf("Its pointer is %p\n",edges[i][j].onGeometry);
+
+						printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
+						throw ErrorException(__FUNCT__,exprintf("See above (might be cryptic...)"));
+					}
+				}
+			}
+		}
+	}
+}
+/*}}}1*/
 	/*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/
-	void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) {
+	void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
 
Index: /issm/trunk/src/c/objects/BamgOpts.cpp
===================================================================
--- /issm/trunk/src/c/objects/BamgOpts.cpp	(revision 3278)
+++ /issm/trunk/src/c/objects/BamgOpts.cpp	(revision 3279)
@@ -37,5 +37,5 @@
 
 	if (bamgopts->coef==0) throw ErrorException(__FUNCT__,exprintf("'coef' should be positive"));
-	if (bamgopts->maxsubdiv<1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >=1"));
+	if (bamgopts->maxsubdiv<=1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >1"));
 	if (bamgopts->Hessiantype!=0  && bamgopts->Hessiantype!=1) throw ErrorException(__FUNCT__,exprintf("'Hessiantype' supported options are 0 and 1"));
 	if (bamgopts->Metrictype!=0   && bamgopts->Metrictype!=1 && bamgopts->Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2"));
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3278)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3279)
@@ -30,4 +30,10 @@
 	FetchData(&bamggeom_in.NumEdges,mxGetField(BAMGGEOMETRY,0,"NumEdges"));
 	FetchData(&bamggeom_in.Edges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges"));
+	FetchData(&bamggeom_in.NumCorners,mxGetField(BAMGGEOMETRY,0,"NumCorners"));
+	FetchData(&bamggeom_in.Corners,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Corners"));
+	FetchData(&bamggeom_in.NumRequiredVertices,mxGetField(BAMGGEOMETRY,0,"NumRequiredVertices"));
+	FetchData(&bamggeom_in.RequiredVertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"RequiredVertices"));
+	FetchData(&bamggeom_in.NumRequiredEdges,mxGetField(BAMGGEOMETRY,0,"NumRequiredEdges"));
+	FetchData(&bamggeom_in.RequiredEdges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"RequiredEdges"));
 	FetchData(&bamggeom_in.hVertices,&lines,&cols,mxGetField(BAMGGEOMETRY,0,"hVertices"));
 	FetchData(&bamggeom_in.NumCrackedEdges,mxGetField(BAMGGEOMETRY,0,"NumCrackedEdges"));