Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 5119)
+++ /issm/trunk/src/c/Makefile.am	(revision 5120)
@@ -69,6 +69,6 @@
 					./objects/Bamg/Mesh.cpp\
 					./objects/Bamg/Mesh.h\
-					./objects/Bamg/MeshVertex.cpp\
-					./objects/Bamg/MeshVertex.h\
+					./objects/Bamg/BamgVertex.cpp\
+					./objects/Bamg/BamgVertex.h\
 					./objects/Bamg/VertexOnEdge.h\
 					./objects/Bamg/VertexOnEdge.cpp\
@@ -605,6 +605,6 @@
 					./objects/Bamg/Mesh.h\
 					./objects/Bamg/Mesh.cpp\
-					./objects/Bamg/MeshVertex.cpp\
-					./objects/Bamg/MeshVertex.h\
+					./objects/Bamg/BamgVertex.cpp\
+					./objects/Bamg/BamgVertex.h\
 					./objects/Bamg/VertexOnEdge.h\
 					./objects/Bamg/VertexOnEdge.cpp\
Index: /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp	(revision 5120)
+++ /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp	(revision 5120)
@@ -0,0 +1,251 @@
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <ctime>
+
+#include "../objects.h"
+
+namespace bamg {
+
+	/*Methods*/
+	/*FUNCTION BamgVertex::Smoothing{{{1*/
+	double  BamgVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/
+
+		register BamgVertex* s=this;
+		BamgVertex &vP = *s,vPsave=vP;
+
+		register Triangle* tbegin= t , *tria = t , *ttc;
+
+		register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc;
+		R2 P(s->r),PNew(0,0);
+		do {
+			k++; 
+
+			if (!tria->Hidden(j)){
+				BamgVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]]; 
+
+				R2 Q = vQ,QP(P-Q);
+				double lQP = LengthInterpole(vP,vQ,QP);
+				PNew += Q+QP/Max(lQP,1e-20);
+				kk ++;
+			}
+			ttc =  tria->TriangleAdj(j);
+			jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
+			tria = ttc;
+			j = NextEdge[jc];
+			if (k>=2000){
+				ISSMERROR("k>=2000 (Maximum number of iterations reached)");
+			}
+		} while ( tbegin != tria); 
+		if (kk<4) return 0;
+		PNew = PNew/(double)kk;
+		R2 Xmove((PNew-P)*omega);
+		PNew = P+Xmove;
+		double delta=Norme2_2(Xmove); 
+
+		Icoor2 deta[3];
+		I2 IBTh  = BTh.toI2(PNew);
+
+		tstart=BTh.FindTriangleContaining(IBTh,deta,tstart);  
+
+		if (tstart->det <0){ // outside
+			double ba,bb;
+			TriangleAdjacent edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ;
+			tstart = edge;
+			vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));
+		}
+		else { // inside
+			double   aa[3];
+			double s = deta[0]+deta[1]+deta[2];
+			aa[0]=deta[0]/s;
+			aa[1]=deta[1]/s;
+			aa[2]=deta[2]/s;
+			vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]);
+		}
+
+		// recompute the det of the triangle
+		vP.r = PNew;
+
+		vP.i = Th.toI2(PNew);
+
+		BamgVertex vPnew = vP;
+
+		int ok=1;
+		int loop=1;
+		k=0;
+		while (ok){
+			ok =0;
+			do {
+				k++; 
+				double detold = tria->det;
+				tria->det =  bamg::det( (*tria)[0],(*tria)[1]  ,(*tria)[2]);
+				if (loop) {
+					BamgVertex *v0,*v1,*v2,*v3;
+					if (tria->det<0) ok =1;			       
+					else if (tria->Quadrangle(v0,v1,v2,v3))
+					  {
+						vP = vPsave;
+						double qold =QuadQuality(*v0,*v1,*v2,*v3);
+						vP = vPnew;
+						double qnew =QuadQuality(*v0,*v1,*v2,*v3);
+						if (qnew<qold) ok = 1;
+					  }
+					else if ( (double)tria->det < detold/2 ) ok=1;
+
+				}
+				tria->SetUnMarkUnSwap(0);
+				tria->SetUnMarkUnSwap(1);
+				tria->SetUnMarkUnSwap(2);
+				ttc =  tria->TriangleAdj(j);
+				jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
+				tria = ttc;
+				j = NextEdge[jc];
+				if (k>=2000){
+					ISSMERROR("k>=2000");
+				}
+			}while ( tbegin != tria); 
+
+			if (ok && loop) vP=vPsave; // no move 
+			loop=0;
+		}
+		return delta;
+	}
+	/*}}}1*/
+	/*FUNCTION BamgVertex::MetricFromHessian{{{1*/
+	void BamgVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){
+		/*Compute Metric from Hessian*/
+
+		/*get options*/
+		double power=(bamgopts->power);
+		double anisomax=(bamgopts->anisomax);
+		double CutOff=bamgopts->cutoff;
+		double hmin=(bamgopts->hmin);
+		double hmax=(bamgopts->hmax);
+		double coef=bamgopts->coeff;
+		int    Metrictype=(bamgopts->Metrictype);
+
+		/*Intermediary*/
+		double ci;
+
+		/*compute multiplicative coefficient depending on Metric Type (2/9 because it is 2d)*/
+
+		//Absolute Error
+		/*
+		 *            2         1       
+		 *Metric M = ---  ------------   Abs(Hessian)
+		 *            9   err * coeff^2  
+		 */
+		if (Metrictype==0){
+			ci= 2.0/9.0 * 1/(err*coef*coef);
+		}
+
+		//Relative Error
+		/*
+		 *            2         1            Abs(Hessian)
+		 *Metric M = ---  ------------  ----------------------
+		 *            9   err * coeff^2  max( |s| , cutoff*max(|s|) )
+		 *
+		 */
+		else if (Metrictype==1){
+			ci= 2.0/9.0 * 1/(err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
+		}
+
+		//Rescaled absolute error
+		/*
+		 *            2         1            Abs(Hessian)
+		 *Metric M = ---  ------------  ---------------------- 
+		 *            9   err * coeff^2       (smax-smin)
+		 */
+		else if (Metrictype==2){
+			ci= 2.0/9.0 * 1/(err*coef*coef) * 1/(smax-smin);
+		}
+		else{
+			ISSMERROR("Metrictype %i not supported yet (use 0,1 or 2(default))",Metrictype);
+		}
+
+		//initialize metric Miv with ci*H
+		Metric Miv(Hxx*ci,Hyx*ci,Hyy*ci);
+
+		//Get eigen values and vectors of Miv
+		MatVVP2x2 Vp(Miv);
+
+		//move eigen valuse to their absolute values
+		Vp.Abs();
+
+		//Apply a power if requested by user
+		if(power!=1.0) Vp.pow(power);
+
+		//modify eigen values according to hmin and hmax
+		Vp.Maxh(hmax);
+		Vp.Minh(hmin);
+
+		//Bound anisotropy by 1/(anisomax)^2
+		Vp.BoundAniso2(1/(anisomax*anisomax));
+
+		//rebuild Metric from Vp
+		Metric MVp(Vp);
+
+		//Apply Metric to vertex
+		m.IntersectWith(MVp);
+
+	}
+	/*}}}1*/
+	/*FUNCTION BamgVertex::Echo {{{1*/
+
+	void BamgVertex::Echo(void){
+
+		printf("Vertex:\n");
+		printf("  integer   coordinates i.x: %i, i.y: %i\n",i.x,i.y);
+		printf("  Euclidean coordinates r.x: %g, r.y: %g\n",r.x,r.y);
+		printf("  ReferenceNumber = %i\n",ReferenceNumber);
+		m.Echo();
+
+		return;
+	}
+	/*}}}*/
+	/*FUNCTION BamgVertex::Optim {{{1*/
+	long BamgVertex::Optim(int i,int koption){ 
+		long ret=0;
+		if ( t && (vint >= 0 ) && (vint <3) ){
+			ret = t->Optim(vint,koption);
+			if(!i){
+				t =0; // for no future optime 
+				vint= 0;
+			}
+		}
+		return ret;
+	}
+	/*}}}*/
+
+	/*Intermediary*/
+	/*FUNCTION QuadQuality{{{1*/
+	double QuadQuality(const BamgVertex & a,const BamgVertex &b,const BamgVertex &c,const BamgVertex &d) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/
+
+		// calcul de 4 angles --
+		R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d);
+		R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D);
+		//  Move(A),Line(B),Line(C),Line(D),Line(A);
+		const Metric & Ma  = a;
+		const Metric & Mb  = b;
+		const Metric & Mc  = c;
+		const Metric & Md  = d;
+
+		double lAB=Norme2(AB);
+		double lBC=Norme2(BC);
+		double lCD=Norme2(CD);
+		double lDA=Norme2(DA);
+		AB /= lAB;  BC /= lBC;  CD /= lCD;  DA /= lDA;
+		// version aniso 
+		double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB);
+		double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC);
+		double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD);
+		double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA);
+		double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA));
+		if (sinmin<=0) return sinmin;
+		return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA)));
+	}
+	/*}}}1*/
+
+} 
Index: /issm/trunk/src/c/objects/Bamg/BamgVertex.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/BamgVertex.h	(revision 5120)
+++ /issm/trunk/src/c/objects/Bamg/BamgVertex.h	(revision 5120)
@@ -0,0 +1,57 @@
+#ifndef _MESHVERTEX_H_
+#define _MESHVERTEX_H_
+
+#include "./include.h"
+#include "./Metric.h"
+#include "./Direction.h"
+#include "./BamgOpts.h"
+
+namespace bamg {
+
+	//classes
+	class Triangle;
+	class Mesh;
+	class VertexOnGeom;
+	class VertexOnEdge;
+
+	class BamgVertex {
+
+		public:
+
+			I2        i;                 // integer coordinates
+			R2        r;                 // real coordinates
+			Metric    m;
+			long      ReferenceNumber;
+			Direction DirOfSearch;
+			short     vint;              // the vertex number in triangle; varies between 0 and 2 in t
+
+			union {
+				Triangle     *t;                    // one triangle which is containing the vertex
+				long          color;
+				BamgVertex   *to;                   // used in geometry BamgVertex to know the Mesh BamgVertex associated
+				VertexOnGeom *onGeometry;           // if vint == 8; // set with Mesh::SetVertexFieldOn()
+				BamgVertex   *onBackgroundVertex;   // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()
+				VertexOnEdge *onBackgroundEdge;     // if vint == 32 on Background edge
+			};
+
+			//Operators
+			operator I2() const {return i;}             // Cast operator
+			operator const R2 & () const {return r;}    // Cast operator
+			operator Metric () const {return m;}        // Cast operator
+			double operator()(R2 x) const { return m(x);} // Get x in the metric m
+
+			//methods (No constructor and no destructors...)
+			double Smoothing(Mesh & ,const Mesh & ,Triangle  * & ,double =1);
+			void   MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
+			void   Echo();
+			int    ref() const { return ReferenceNumber;}
+			long   Optim(int =1,int =0); 
+
+			//inline functions
+			inline void Set(const BamgVertex &rec,const Mesh & ,Mesh & ){*this=rec;}
+	};
+
+	//Intermediary
+	double QuadQuality(const BamgVertex &,const BamgVertex &,const BamgVertex &,const BamgVertex &);
+}
+#endif
Index: /issm/trunk/src/c/objects/Bamg/Edge.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Edge.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Edge.h	(revision 5120)
@@ -3,5 +3,5 @@
 
 #include "./include.h"
-#include "./MeshVertex.h"
+#include "./BamgVertex.h"
 #include "../../include/include.h"
 #include "../../shared/Exceptions/exceptions.h"
@@ -16,5 +16,5 @@
 
 		public:
-			MeshVertex* v[2];
+			BamgVertex* v[2];
 			long ref;
 			GeometricalEdge* onGeometry;
@@ -22,11 +22,11 @@
 
 			//Operators
-			MeshVertex & operator[](int i){return *v[i];};
-			MeshVertex * operator()(int i){return v[i];};
+			BamgVertex & operator[](int i){return *v[i];};
+			BamgVertex * operator()(int i){return v[i];};
 			R2       operator()(double t) const; // return the point 
-			const MeshVertex & operator[](int i) const { return *v[i];};
+			const BamgVertex & operator[](int i) const { return *v[i];};
 
 			//Methods
-			void ReNumbering(MeshVertex *vb,MeshVertex *ve, long *renu){
+			void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
 				if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
 				if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
Index: /issm/trunk/src/c/objects/Bamg/GeometricalVertex.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/GeometricalVertex.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/GeometricalVertex.h	(revision 5120)
@@ -3,5 +3,5 @@
 
 #include "./include.h"
-#include "MeshVertex.h"
+#include "BamgVertex.h"
 
 namespace bamg {
@@ -9,5 +9,5 @@
 	class Geometry;
 
-	class GeometricalVertex : public MeshVertex { 
+	class GeometricalVertex : public BamgVertex { 
 
 		public:
Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5120)
@@ -13,5 +13,13 @@
 
 	/*Constructors/Destructors*/
-	/*FUNCTION  Geometry::Geometry(const Geometry & Gh){{{1*/
+	Geometry::Geometry(){
+		Init();
+	}
+	Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
+		Init();
+		ReadGeometry(bamggeom,bamgopts);
+		AfterRead();
+	}
+	/*FUNCTION  Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/
 	Geometry::Geometry(const Geometry & Gh) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/
@@ -22,5 +30,4 @@
 		quadtree=0;
 		vertices = nbv ? new GeometricalVertex[nbv] : NULL;
-		triangles = nbt ? new  Triangle[nbt]:NULL;
 		edges = nbe ? new GeometricalEdge[nbe]:NULL;
 		curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
@@ -34,7 +41,4 @@
 		for (i=0;i<NbSubDomains;i++)
 		 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
-		if (nbt);  {
-			ISSMERROR("nbt");
-		}
 	}
 	/*}}}1*/
@@ -42,19 +46,12 @@
 	Geometry::~Geometry() {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
-
-		long int verbose=0;
-
-		if (NbRef>0){
-			ISSMERROR("NbRef>0");
-		}
-		if(vertices)  delete [] vertices;vertices=0;
-		if(edges)     delete [] edges;edges=0;
-		// if(edgescomponante) delete [] edgescomponante; edgescomponante=0;
-		if(triangles) delete [] triangles;triangles=0;
-		if(quadtree)  delete  quadtree;quadtree=0;
-		if(curves)  delete  []curves;curves=0;NbOfCurves=0;
+		ISSMASSERT(NbRef<=0);
+		if(vertices)  delete [] vertices;   vertices=0;
+		if(edges)     delete [] edges;      edges=0;
+		if(triangles) delete [] triangles;  triangles=0;
+		if(quadtree)  delete  quadtree;     quadtree=0;
+		if(curves)    delete  []curves;     curves=0;NbOfCurves=0;
 		if(subdomains) delete [] subdomains;subdomains=0;
-		//  if(ordre)     delete [] ordre;
-		EmptyGeometry();
+		Init();
 	}
 	/*}}}1*/
@@ -65,6 +62,6 @@
 
 		int verbose;
-		nbiv=nbv=nbvx=0;
-		nbe=nbt=nbtx=0;
+		nbv=0;
+		nbe=0;
 		NbOfCurves=0;
 
@@ -72,19 +69,12 @@
 		int i,j,k,n,i1,i2;
 
-		//initialize some variables
-		int Version=1,dim=2;
-		nbv=bamggeom->VerticesSize[0];
-		nbe=bamggeom->EdgesSize[0];
-		nbvx = nbv;
-		nbiv = nbv;
+		/*initialize some variables*/
 		verbose=bamgopts->verbose;
+		nbv  = bamggeom->VerticesSize[0];
+		nbe  = bamggeom->EdgesSize[0];
 
 		//some checks
-		if (bamggeom->Vertices==NULL){
-			ISSMERROR("the domain provided does not contain any vertex");
-		}
-		if (bamggeom->Edges==NULL){
-			ISSMERROR("the domain provided does not contain any edge");
-		}
+		if (bamggeom->Vertices==NULL) ISSMERROR("the domain provided does not contain any vertex");
+		if (bamggeom->Edges==NULL) ISSMERROR("the domain provided does not contain any edge");
 
 		//Vertices
@@ -92,5 +82,5 @@
 			if(verbose>5) printf("      processing Vertices\n");
 			if (bamggeom->VerticesSize[1]!=3) ISSMERROR("Vertices should have 3 columns");
-			vertices = new GeometricalVertex[nbvx];
+			vertices = new GeometricalVertex[nbv];
 			for (i=0;i<nbv;i++) {
 				vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
@@ -101,6 +91,5 @@
 				vertices[i].Set();
 			}
-			//find domain extrema for pmin,pmax that will define the square box
-			//used for by the quadtree
+			/*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
 			pmin =  vertices[0].r;
 			pmax =  vertices[0].r;
@@ -111,16 +100,20 @@
 				pmax.y = Max(pmax.y,vertices[i].r.y);
 			}
-			R2 DD05 = (pmax-pmin)*0.05;
-			pmin -=  DD05;
-			pmax +=  DD05;
-			//coefIcoor is the coefficient used to have coordinates
-			//in ]0 1[^2:  x'=coefIcoor*(x-pmin.x)
-			coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
-			if(coefIcoor <=0){
-				ISSMERROR("coefIcoor should be positive");
-			}
+			/*Offset pmin and pmax to avoid round-off errors*/
+			R2 offset = (pmax-pmin)*0.05;
+			pmin -= offset;
+			pmax += offset;
+			/*coefIcoor is the coefficient used for integer coordinates:
+			 *                       (x-pmin.x)
+			 * Icoor x = (2^30 -1) ------------ 
+			 *                          D
+			 * where D is the longest side of the domain (direction x or y)
+			 * so that (x-pmin.x)/D is in ]0 1[
+			 */
+			coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
+			if(coefIcoor<=0) ISSMERROR("coefIcoor should be positive");
 		}
 		else{
-			ISSMERROR("No MeshVertex provided");
+			ISSMERROR("No BamgVertex provided");
 		}
 
@@ -427,5 +420,5 @@
 		double eps=1e-20;
 		QuadTree quadtree; // build quadtree to find duplicates
-		MeshVertex* v0=vertices; 
+		BamgVertex* v0=vertices; 
 		GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0;   
 
@@ -446,5 +439,5 @@
 
 			//find nearest vertex already present in the quadtree (NULL if empty)
-			MeshVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 
+			BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 
 
 			//if there is a vertex found that is to close to vertices[i] -> error
@@ -454,6 +447,6 @@
 				j=vg-v0g;
 				//check that the clostest vertex is not itself...
-				if ( v !=  &(MeshVertex &) vertices[j]){
-					ISSMERROR(" v !=  &(MeshVertex &) vertices[j]");
+				if ( v !=  &(BamgVertex &) vertices[j]){
+					ISSMERROR(" v !=  &(BamgVertex &) vertices[j]");
 				}
 				vertices[i].link = vertices + j;
@@ -748,11 +741,9 @@
 		GeometricalEdge* on =start,* pon=0;
 		// walk with the cos on geometry
-		int k=0;
+		int counter=0;
 		while(pon != on){  
-			k++;
+			counter++;
+			ISSMASSERT(counter<100);
 			pon = on;
-			if (k>=100){
-				ISSMERROR("k>=100");
-			}
 			R2 A= (*on)[0];
 			R2 B= (*on)[1];
@@ -775,10 +766,6 @@
 
 		printf("Geometry:\n");
-		printf("   nbvx (maximum number of vertices) : %i\n",nbvx);
-		printf("   nbtx (maximum number of triangles): %i\n",nbtx);
 		printf("   nbv  (number of vertices) : %i\n",nbv);
-		printf("   nbt  (number of triangles): %i\n",nbt);
 		printf("   nbe  (number of edges)    : %i\n",nbe);
-		printf("   nbv  (number of initial vertices) : %i\n",nbiv);
 		printf("   NbSubDomains: %i\n",NbSubDomains);
 		printf("   NbOfCurves: %i\n",NbOfCurves);
@@ -797,28 +784,24 @@
 	}
 	/*}}}*/
-	/*FUNCTION  Geometry::EmptyGeometry(){{{1*/
-	void Geometry::EmptyGeometry() {
+	/*FUNCTION  Geometry::Init{{{1*/
+	void Geometry::Init(void){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
 
 		NbRef=0;
+		nbv=0;
+		nbe=0;
 		quadtree=0;
 		curves=0;
-		// edgescomponante=0;
 		triangles=0;
 		edges=0;
 		vertices=0;
 		NbSubDomains=0;
-		//  nbtf=0;
-		//  BeginOfCurve=0;  
-		nbiv=nbv=nbvx=0;
-		nbe=nbt=nbtx=0;
 		NbOfCurves=0;
-		//  BeginOfCurve=0;
 		subdomains=0;
-		MaxCornerAngle = 10*Pi/180;
+		MaxCornerAngle = 10*Pi/180; //default is 10 degres
 	}
 	/*}}}1*/
 	/*FUNCTION  Geometry::ProjectOnCurve {{{1*/
-	GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,MeshVertex &V,VertexOnGeom &GV) const {
+	GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
 		/*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/
@@ -839,6 +822,6 @@
 
 		//Get the two vertices of the edge
-		const MeshVertex &v0=e[0];
-		const MeshVertex &v1=e[1];
+		const BamgVertex &v0=e[0];
+		const BamgVertex &v1=e[1];
 
 		//Get position of V0, V1 and vector v0->v1
@@ -906,8 +889,8 @@
 
 		if ((*eg0)(sens0)==(GeometricalVertex*)vg0)
-		 vg0=VertexOnGeom(*(MeshVertex*) vg0,*eg0,sens0); //vg0 = absisce
+		 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,sens0); //vg0 = absisce
 
 		if ((*eg1)(sens1)==(GeometricalVertex*)vg1)
-		 vg1=VertexOnGeom(*(MeshVertex*) vg1,*eg1,sens1);
+		 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,sens1);
 
 		double sg;
Index: /issm/trunk/src/c/objects/Bamg/Geometry.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.h	(revision 5120)
@@ -19,23 +19,25 @@
 
 		public:
-			long NbRef;     // counter of ref on the this class if 0 we can delete
-			long nbvx,nbtx; // maximum number of vertices
-			long nbv,nbt,nbiv,nbe; // number of vertices
-			long NbSubDomains; // 
-			long NbOfCurves;
-			GeometricalVertex* vertices;
-			Triangle* triangles; 
-			GeometricalEdge* edges;
-			QuadTree* quadtree;
-			GeometricalSubDomain* subdomains;
-			Curve* curves;
-			R2 pmin,pmax; // extrema
-			double coefIcoor;  // coef to integer Icoor1;
-			double MaxCornerAngle;
+
+			long                  NbRef;                         // counter of ref on the this class if 0 we can delete
+			long                  nbv;                           // number of vertices
+			long                  nbe;                           // number of edges
+			long                  NbSubDomains;
+			long                  NbOfCurves;
+			GeometricalVertex    *vertices;
+			Triangle             *triangles;
+			GeometricalEdge      *edges;
+			QuadTree             *quadtree;
+			GeometricalSubDomain *subdomains;
+			Curve                *curves;
+			R2                    pmin,pmax;                     // domain extrema coordinates
+			double                coefIcoor;                     // coef to integer Icoor1;
+			double                MaxCornerAngle;
 
 			//Constructor/Destructors
 			~Geometry(); 
-			Geometry(const Geometry & Gh); //Copy  Operator 
-			Geometry(int nbg,const Geometry** ag); // intersection operator 
+			Geometry();
+			Geometry(const Geometry & Gh);
+			Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
 
 			//Operators
@@ -46,16 +48,14 @@
 
 			//Methods
-			int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
 			void Echo();
 			I2 toI2(const R2 & P) const {
 				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
-							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
+							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
+			}
 			double MinimalHmin() {return 2.0/coefIcoor;}
 			double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
 			void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
-			void EmptyGeometry();
-			Geometry() {EmptyGeometry();}// empty Geometry
+			void Init(void);
 			void AfterRead();
-			Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}
 			long Number(const GeometricalVertex & t) const  { return &t - vertices;}
 			long Number(const GeometricalVertex * t) const  { return t - vertices;}
@@ -64,5 +64,5 @@
 			long Number(const Curve * c) const  { return c - curves;}
 			void UnMarkEdges() {for (int i=0;i<nbe;i++) edges[i].SetUnMark();}
-			GeometricalEdge *  ProjectOnCurve(const Edge & ,double,MeshVertex &,VertexOnGeom &) const ;
+			GeometricalEdge *  ProjectOnCurve(const Edge & ,double,BamgVertex &,VertexOnGeom &) const ;
 			GeometricalEdge *  Containing(const R2 P,  GeometricalEdge * start) const;
 			void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
Index: /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 5120)
@@ -44,5 +44,5 @@
 				double ba,bb;
 				TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
-				MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
+				BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
 				NewItem(A,Metric(ba,v0,bb,v1));
 				t=edge;
@@ -50,5 +50,5 @@
 				if (det(v0.i,v1.i,b)>=0) {
 					TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
-					MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
+					BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
 					NewItem(A,Metric(ba,v0,bb,v1));
 					return;
@@ -80,5 +80,5 @@
 					long int verbose=2;
 					TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
-					MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
+					BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
 					NewItem(A,Metric(ba,v0,bb,v1));
 					return;
@@ -232,5 +232,5 @@
 			lIntTria[Size].x = x;
 			Metric m0,m1,m2;
-			register MeshVertex * v;
+			register BamgVertex * v;
 			if ((v=(*tt)(0))) m0    = v->m;
 			if ((v=(*tt)(1))) m1    = v->m;
@@ -307,5 +307,5 @@
 	/*}}}1*/
 	/*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
-	long ListofIntersectionTriangles::NewPoints(MeshVertex* vertices,long &nbv,long nbvx){
+	long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
 
@@ -354,5 +354,5 @@
 
 			si += sint;
-			if ( nbv<nbvx) {
+			if ( nbv<maxnbv) {
 				vertices[nbv].r = C;
 				vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
Index: /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h	(revision 5120)
@@ -71,5 +71,5 @@
 			void  SplitEdge(const Mesh & ,const R2 &,const R2  &,int nbegin=0); 
 			double Length(); 
-			long  NewPoints(MeshVertex *,long & nbv,long nbvx);
+			long  NewPoints(BamgVertex *,long & nbv,long maxnbv);
 			void  NewSubSeg(GeometricalEdge *e,double s0,double s1){ 
 				long int verbosity=0;
Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5120)
@@ -16,5 +16,5 @@
 
 		/*Initialize fields*/
-		PreInit(0);
+		Init(0);
 
 		/*Read Geometry if provided*/
@@ -45,5 +45,5 @@
 	Mesh::Mesh(double* index,double* x,double* y,int nods,int nels):Gh(*(new Geometry())),BTh(*this){
 
-		PreInit(0);
+		Init(0);
 		ReadMesh(index,x,y,nods,nels);
 		SetIntCoor();
@@ -100,6 +100,6 @@
 		  printf("   number of triangles %i, remove = %i\n",kt,nbInT-kt);
 		  printf("   number of New boundary edge %i\n",nbNewBedge);
-		  long inbvx =k;
-		  PreInit(inbvx);
+		  long imaxnbv =k;
+		  Init(imaxnbv);
 		  for (i=0;i<Tho.nbv;i++)
 			if (kk[i]>=0) 
@@ -110,6 +110,6 @@
 				nbv++;
 			  }
-		  if (inbvx != nbv){
-			  ISSMERROR("inbvx != nbv");
+		  if (imaxnbv != nbv){
+			  ISSMERROR("imaxnbv != nbv");
 		  }
 		  for (i=0;i<Tho.nbt;i++)
@@ -154,19 +154,18 @@
 	  }
 	/*}}}1*/
-	/*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long nbvxx) COPY{{{1*/
-	Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long nbvxx)
+	/*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) COPY{{{1*/
+	Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in)
 	  : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) {
 		  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)*/
 		  Gh.NbRef++;
-		  nbvxx = Max(nbvxx,Th.nbv); 
+		  maxnbv_in = Max(maxnbv_in,Th.nbv); 
 		  long i;
 		  // do all the allocation to be sure all the pointer existe
 
-		  PreInit(nbvxx);// to make the allocation 
+		  Init(maxnbv_in);// to make the allocation 
 		  // copy of triangles
 		  nt=Th.nt;
 		  nbv = Th.nbv;
 		  nbt = Th.nbt;
-		  nbiv = Th.nbiv;
 		  nbe = Th.nbe;
 		  NbSubDomains = Th.NbSubDomains;
@@ -245,5 +244,5 @@
 			else if (BTh.NbRef==0) delete &BTh;
 		}
-		PreInit(0); // set all to zero 
+		Init(0); // set all to zero 
 	}
 	/*}}}1*/
@@ -260,12 +259,11 @@
 
 		nbv=nods;
-		nbvx=nbv;
+		maxnbv=nbv;
 		nbt=nels;
-		nbiv=nbvx;
 
 		//Vertices
 		if (verbose) printf("Reading vertices (%i)\n",nbv);
-		vertices=(MeshVertex*)xmalloc(nbv*sizeof(MeshVertex));
-		ordre=(MeshVertex**)xmalloc(nbv*sizeof(MeshVertex*));
+		vertices=(BamgVertex*)xmalloc(nbv*sizeof(BamgVertex));
+		ordre=(BamgVertex**)xmalloc(nbv*sizeof(BamgVertex*));
 		for (i=0;i<nbv;i++){
 			vertices[i].r.x=x[i];
@@ -276,5 +274,5 @@
 			vertices[i].color=0;
 		}
-		nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals 
+		nbtx=2*maxnbv-2; // for filling The Holes and quadrilaterals 
 
 		//Triangles
@@ -311,7 +309,6 @@
 
 		nbv=bamgmesh->VerticesSize[0];
-		nbvx=nbv;
+		maxnbv=nbv;
 		nbt=bamgmesh->TrianglesSize[0];
-		nbiv=nbvx;
 
 		//Vertices
@@ -319,6 +316,6 @@
 			if(verbose>5) printf("      processing Vertices\n");
 
-			vertices=(MeshVertex*)xmalloc(nbv*sizeof(MeshVertex));
-			ordre=(MeshVertex**)xmalloc(nbv*sizeof(MeshVertex*));
+			vertices=(BamgVertex*)xmalloc(nbv*sizeof(BamgVertex));
+			ordre=(BamgVertex**)xmalloc(nbv*sizeof(BamgVertex*));
 
 			for (i=0;i<nbv;i++){
@@ -330,5 +327,5 @@
 				vertices[i].color=0;
 			}
-			nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals 
+			nbtx=2*maxnbv-2; // for filling The Holes and quadrilaterals 
 		}
 		else{
@@ -450,5 +447,5 @@
 			for (i=0;i<nbe;i++){
 				for (j=0;j<2;j++) { 
-					MeshVertex *v=edges[i].v[j];
+					BamgVertex *v=edges[i].v[j];
 					long i0=v->color,j0;
 					if(i0==-1){
@@ -726,5 +723,5 @@
 				Triangle &t =triangles[i];
 				Triangle* ta;
-				MeshVertex *v0,*v1,*v2,*v3;
+				BamgVertex *v0,*v1,*v2,*v3;
 				if (reft[i]<0) continue;
 				if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { 
@@ -775,5 +772,5 @@
 				VertexOnGeom &v=VerticesOnGeomVertex[i];
 				ISSMASSERT(v.OnGeomVertex());
-				bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((MeshVertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((BamgVertex*)v)+1; //back to Matlab indexing
 				bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number((GeometricalVertex*)v)+1; //back to Matlab indexing
 			}
@@ -791,5 +788,5 @@
 					ISSMERROR("A vertices supposed to be OnGeometricEdge is actually not");
 				}
-				bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((MeshVertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((BamgVertex*)v)+1; //back to Matlab indexing
 				bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
 				bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; //absisce
@@ -1020,5 +1017,5 @@
 			for (int j=0;j<2;j++){
 
-				MeshVertex V;
+				BamgVertex V;
 				VertexOnGeom GV;
 				Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
@@ -1066,5 +1063,5 @@
 	/*}}}1*/
 	/*FUNCTION Mesh::AddVertex{{{1*/
-	void Mesh::AddVertex( MeshVertex &s,Triangle* t, Icoor2* det3) {
+	void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/
 		// -------------------------------------------
@@ -1086,5 +1083,5 @@
 		Triangle* tt[3];
 		//three vertices of t
-		MeshVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];
+		BamgVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];
 		//three determinants
 		Icoor2 det3local[3];
@@ -1435,5 +1432,5 @@
 			for (j=0;j<2;j++){ 
 				//get current vertex
-				MeshVertex* v=edges[i].v[j];
+				BamgVertex* v=edges[i].v[j];
 				//get vertex color (i0)
 				long i0=v->color;
@@ -1575,5 +1572,5 @@
 		for (i=0;i<nbv;i++){
 			if((j=colorV[i])>=0){
-				MeshVertex & v = Gh.vertices[j];
+				BamgVertex & v = Gh.vertices[j];
 				v = vertices[i];
 				v.color =0;
@@ -2326,5 +2323,5 @@
 		delete [] Edgeflags;
 
-		//Reset MeshVertex to On
+		//Reset BamgVertex to On
 		SetVertexFieldOn();
 
@@ -2557,5 +2554,5 @@
 						ISSMERROR("&e==NULL");
 					}
-					MeshVertex * v0 =  e(0),*v1 = e(1);
+					BamgVertex * v0 =  e(0),*v1 = e(1);
 					Triangle *t  = v0->t;
 					int sens = Gh.subdomains[i].sens;
@@ -2645,5 +2642,5 @@
 
 			/*Call NearestVertex*/
-			MeshVertex *a = quadtree->NearestVertex(B.x,B.y) ;
+			BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
 
 			/*Check output (Vertex a)*/
@@ -2722,6 +2719,2047 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Mesh::GeomToTriangles0{{{1*/
-	void Mesh::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){
+	/*FUNCTION Mesh::Init{{{1*/
+	void Mesh::Init(long maxnbv_in) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/
+
+		long int verbose=0;
+
+		srand(19999999);
+		NbRef=0;
+		NbOfTriangleSearchFind =0;
+		NbOfSwapTriangle =0;
+		nbv=0;
+		maxnbv=maxnbv_in;
+		nbt=0;
+		NbOfQuad = 0;
+		nbtx=2*maxnbv_in-2;
+		NbSubDomains=0;
+		NbVertexOnBThVertex=0;
+		NbVertexOnBThEdge=0;
+		VertexOnBThVertex=NULL;
+		VertexOnBThEdge=NULL;
+
+		NbCrackedVertices=0;
+		CrackedVertices  =NULL;  
+		NbCrackedEdges   =0;
+		CrackedEdges     =NULL;  
+		nbe = 0; 
+
+		if (maxnbv_in) {
+			vertices=new BamgVertex[maxnbv];
+			ISSMASSERT(vertices);
+			ordre=new (BamgVertex* [maxnbv]);
+			ISSMASSERT(ordre);
+			triangles=new Triangle[nbtx];
+			ISSMASSERT(triangles);
+		}
+		else {
+			vertices=NULL;
+			ordre=NULL;
+			triangles=NULL;
+			nbtx=0;
+		} 
+
+		quadtree=NULL;
+		edges=NULL;
+		VerticesOnGeomVertex=NULL;
+		VerticesOnGeomEdge=NULL;
+		NbVerticesOnGeomVertex=0;
+		NbVerticesOnGeomEdge=0;
+		subdomains=NULL;
+		NbSubDomains=0;
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::Insert{{{1*/
+	void Mesh::Insert() {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
+
+		/*Insert points in the existing Geometry*/
+
+		//Intermediary
+		int i;
+
+		/*Get options*/
+		long int verbose=2;
+
+		//Display info
+		if (verbose>2) printf("   Insert initial %i vertices\n",nbv);
+
+		//Compute integer coordinates and determinants for the existing vertices (from Geometry)
+		SetIntCoor();
+
+		/*Now we want to build a list (ordre) of the vertices in a random
+		 * order. To do so, we use the following method:
+		 *
+		 * From an initial k0 in [0 nbv[ random (vertex number)
+		 * the next k (vertex number) is computed using a big
+		 * prime number (PN>>nbv) following:
+		 *
+		 * k_{i+1} = k_i + PN  [nbv]
+		 *
+		 * let's show that:
+		 *
+		 *   for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j
+		 *
+		 * Let's assume that there are 2 distinct j1 and j2 such that
+		 * k_j1 = k_j2
+		 *
+		 * This means that
+		 *  
+		 *  k0+j1*PN = k0+j2*PN [nbv]
+		 *  (j1-j2)*PN =0       [nbv]
+		 * since PN is a prime number larger than nbv, and nbv!=1
+		 *  j1-j2=0             [nbv]
+		 * BUT
+		 *  j1 and j2 are in [0 nbv[ which is impossible.
+		 *
+		 *  We hence have built a random list of nbv elements of
+		 *  [0 nbv[ all distincts*/
+		for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ;
+		const long PrimeNumber= BigPrimeNumber(nbv) ;
+		int   k0=rand()%nbv; 
+		for (int is3=0; is3<nbv; is3++){
+			ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv];
+		}
+
+		/*Modify ordre such that the first 3 vertices form a triangle*/
+
+		//get first vertex i such that [0,1,i] are not aligned
+		for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
+			//if i is higher than nbv, it means that all the determinants are 0,
+			//all vertices are aligned!
+			if  ( ++i >= nbv) {
+				ISSMERROR("all the vertices are aligned");
+			}
+		}
+		// exchange i et 2 in "ordre" so that
+		// the first 3 vertices are not aligned (real triangle)
+		Exchange(ordre[2], ordre[i]);
+
+		/*Take the first edge formed by the first two vertices and build
+		 * the initial simple mesh from this edge and 2 boundary triangles*/
+
+		BamgVertex *  v0=ordre[0], *v1=ordre[1];
+
+		nbt = 2;
+		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;
+
+		//Build adjacence
+		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 triangle: det = -1
+		triangles[1].det = -1;  //boundary triangle: det = -1
+
+		triangles[0].SetTriangleContainingTheVertex();
+		triangles[1].SetTriangleContainingTheVertex();
+
+		triangles[0].link=&triangles[1];
+		triangles[1].link=&triangles[0];
+
+		//build quadtree
+		if (!quadtree)  quadtree = new QuadTree(this,0);
+		quadtree->Add(*v0);
+		quadtree->Add(*v1);
+
+		/*Now, add the vertices One by One*/
+
+		long NbSwap=0;
+		if (verbose>3) printf("   Begining of insertion process...\n");
+
+		for (int icount=2; icount<nbv; icount++) {
+
+			//Get new vertex
+			BamgVertex *newvertex=ordre[icount];
+
+			//Find the triangle in which newvertex is located
+			Icoor2 dete[3];
+			Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
+
+			//Add newvertex to the quadtree
+			quadtree->Add(*newvertex); 
+
+			//Add newvertex to the existing mesh
+			AddVertex(*newvertex,tcvi,dete);
+
+			//Make the mesh Delaunay around newvertex by swaping the edges
+			NbSwap += newvertex->Optim(1,0);
+		}
+
+		//Display info
+		if (verbose>3) {
+			printf("      NbSwap of insertion: %i\n",NbSwap);
+			printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
+		}
+
+#ifdef NBLOOPOPTIM
+
+		k0 = rand()%nbv ; 
+		for (int is4=0; is4<nbv; is4++) 
+		 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];
+
+		for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){
+			long  NbSwap = 0;
+			for (int is1=0; is1<nbv; is1++) 
+			 NbSwap += ordre[is1]->Optim(0,0);
+			if (verbose>3) {
+				printf("      Optim Loop: %i\n",Nbloop);
+				printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
+			}
+			if(!NbSwap) break;
+		}
+		ReMakeTriangleContainingTheVertex(); 
+		// because we break the TriangleContainingTheVertex
+#endif
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::InsertNewPoints{{{1*/
+	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
+
+		long int verbose=0;
+		double seuil= 1.414/2 ;// for two close point 
+		long i;
+		long NbSwap=0;
+		Icoor2 dete[3];
+
+		//number of new points
+		const long nbvnew=nbv-nbvold;
+
+		//display info if required
+		if (verbose>5) printf("      Try to Insert %i new points\n",nbvnew);
+
+		//return if no new points
+		if (!nbvnew) return 0; 
+
+		/*construction of a random order*/
+		const long PrimeNumber= BigPrimeNumber(nbv)  ;
+		//remainder of the division of rand() by nbvnew
+		long k3 = rand()%nbvnew;
+		//loop over the new points
+		for (int is3=0; is3<nbvnew; is3++){
+			register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
+			register long i=nbvold+is3; 
+			ordre[i]= vertices + j;
+			ordre[i]->ReferenceNumber=i;
+		}
+
+		// for all the new point
+		long iv=nbvold;
+		for (i=nbvold;i<nbv;i++){
+			BamgVertex &vi=*ordre[i];
+			vi.i=toI2(vi.r);
+			vi.r=toR2(vi.i);
+			double hx,hy;
+			vi.m.Box(hx,hy);
+			Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
+			if (!quadtree->ToClose(vi,seuil,hi,hj)){
+				// a good new point 
+				BamgVertex &vj = vertices[iv];
+				long  j=vj.ReferenceNumber; 
+				if (&vj!=ordre[j]){
+					ISSMERROR("&vj!= ordre[j]");
+				}
+				if(i!=j){ 
+					Exchange(vi,vj);
+					Exchange(ordre[j],ordre[i]);
+				}
+				vj.ReferenceNumber=0; 
+				Triangle *tcvj=FindTriangleContaining(vj.i,dete);
+				if (tcvj && !tcvj->link){
+					tcvj->Echo();
+					ISSMERROR("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link);
+				}
+				quadtree->Add(vj);
+				AddVertex(vj,tcvj,dete);
+				NbSwap += vj.Optim(1);          
+				iv++;
+			}
+		} 
+		if (verbose>3) {
+			printf("         number of new points: %i\n",iv);
+			printf("         number of to close (?) points: %i\n",nbv-iv);
+			printf("         number of swap after: %i\n",NbSwap);
+		}
+		nbv = iv;
+
+		for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);  
+		if (verbose>3) printf("   NbSwap=%i\n",NbSwap);
+
+		NbTSwap +=  NbSwap ;
+		return nbv-nbvold;
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::MakeGeometricalEdgeToEdge{{{1*/
+	Edge** Mesh::MakeGeometricalEdgeToEdge() {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeometricalEdgeToEdge)*/
+
+		if (!Gh.nbe){
+			ISSMERROR("!Gh.nbe");
+		}
+		Edge **e= new (Edge* [Gh.nbe]);
+
+		long i;
+		for ( i=0;i<Gh.nbe ; i++)
+		 e[i]=NULL;
+		for ( i=0;i<nbe ; i++) 
+		  { 
+			Edge * ei = edges+i;
+			GeometricalEdge *onGeometry = ei->onGeometry; 
+			e[Gh.Number(onGeometry)] = ei;    
+		  }
+		for ( i=0;i<nbe ; i++) 
+		 for (int ii=0;ii<2;ii++) { 
+			 Edge * ei = edges+i;
+			 GeometricalEdge *onGeometry = ei->onGeometry;
+			 int j= ii;
+			 while (!(*onGeometry)[j].Required()) { 
+				 Adj(onGeometry,j); // next geom edge
+				 j=1-j;
+				 if (e[Gh.Number(onGeometry)])  break; // optimisation
+				 e[Gh.Number(onGeometry)] = ei; 
+			 }
+		 }
+
+		int kk=0;
+		for ( i=0;i<Gh.nbe ; i++){
+			if (!e[i]){
+				kk++;
+				if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);
+			}
+		}
+		if(kk) ISSMERROR("See above");
+
+		return e;
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::MakeQuadrangles{{{1*/
+	void Mesh::MakeQuadrangles(double costheta){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/
+
+		long int verbose=0;
+
+		if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta);
+
+		if (costheta >1) {
+			if (verbose>5) printf("   do nothing: costheta > 1\n");
+		}
+
+
+			long nbqq = (nbt*3)/2;
+			DoubleAndInt *qq = new DoubleAndInt[nbqq];
+
+			long i,ij;
+			int j;
+			long k=0;
+			for (i=0;i<nbt;i++)
+			 for (j=0;j<3;j++)
+			  if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)
+				qq[k++].i3j=i*3+j;
+			//  sort  qq
+			HeapSort(qq,k);
+
+			long kk=0;
+			for (ij=0;ij<k;ij++) { 
+				i=qq[ij].i3j/3;
+				j=(int) (qq[ij].i3j%3);
+				// optisamition no float computation  
+				if (triangles[i].QualityQuad(j,0) >=costheta) 
+				 triangles[i].SetHidden(j),kk++;
+			  }
+			NbOfQuad = kk;
+			if (verbose>2){
+				printf("   number of quadrilaterals    = %i\n",NbOfQuad);
+				printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
+				printf("   number of outside triangles = %i\n",NbOutT);
+			}
+			delete [] qq;
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::MakeQuadTree{{{1*/
+	void Mesh::MakeQuadTree() {  
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/
+
+		long int verbose=0;
+		if (  !quadtree )  quadtree = new QuadTree(this);
+
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::MaxSubDivision{{{1*/
+	void  Mesh::MaxSubDivision(double maxsubdiv) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
+
+		long int verbose=0;
+
+		const  double maxsubdiv2 = maxsubdiv*maxsubdiv;
+		if(verbose>1) printf("   Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);
+		// for all the edges 
+		// if the len of the edge is to long 
+		long it,nbchange=0;    
+		double lmax=0;
+		for (it=0;it<nbt;it++){
+			Triangle &t=triangles[it];
+			for (int j=0;j<3;j++){
+				Triangle &tt = *t.TriangleAdj(j);
+				if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)){
+					BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
+					BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
+					R2 AB= (R2) v1-(R2) v0;
+					Metric M = v0;
+					double l = M(AB,AB);
+					lmax = Max(lmax,l);
+					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;
+						D2xD2 Rt1(Rt.inv());
+						D2xD2 D(maxsubdiv2,0,0,lc);
+						D2xD2 MM = Rt1*D*Rt1.t();
+						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
+						double lc = M(AC,AC);
+						D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
+						D2xD2 Rt1(Rt.inv());
+						D2xD2 D(maxsubdiv2,0,0,lc);
+						D2xD2  MM = Rt1*D*Rt1.t();
+						v1.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
+						nbchange++;
+					}
+				}
+			}
+		}
+		if(verbose>3){
+			printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
+		}
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::MetricAt{{{1*/
+	Metric Mesh::MetricAt(const R2 & A) const { 
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/
+
+		I2 a = toI2(A);
+		Icoor2 deta[3];
+		Triangle * t =FindTriangleContaining(a,deta);
+		if (t->det <0) { // outside
+			double ba,bb;
+			TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;
+			return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
+		else { // inside
+			double   aa[3];
+			double s = deta[0]+deta[1]+deta[2];
+			aa[0]=deta[0]/s;
+			aa[1]=deta[1]/s;
+			aa[2]=deta[2]/s;
+			return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
+		}
+	}
+	/*}}}1*/
+/*FUNCTION Mesh::NearestVertex{{{1*/
+BamgVertex* Mesh::NearestVertex(Icoor1 i,Icoor1 j) {
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/
+	return  quadtree->NearestVertex(i,j); 
+} 
+/*}}}1*/
+	/*FUNCTION Mesh::NewPoints{{{1*/
+	void  Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
+
+		int i,j,k;
+		long NbTSwap=0;
+		long nbtold=nbt;
+		long nbvold=nbv;
+		long Headt=0;
+		long next_t;
+		long* first_np_or_next_t=new long[nbtx];
+		Triangle* t=NULL;
+
+		/*Recover options*/
+		int verbose=bamgopts->verbose;
+
+		/*First, insert old points if requested*/
+		if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){
+			if (verbose>5) printf("         Inserting initial mesh points\n");
+			for (i=0;i<Bh.nbv;i++){ 
+				BamgVertex &bv=Bh[i];
+				if (!bv.onGeometry){
+					vertices[nbv].r   = bv.r;
+					vertices[nbv++].m = bv.m;
+				}
+			}
+			Bh.ReMakeTriangleContainingTheVertex();     
+			InsertNewPoints(nbvold,NbTSwap)   ;            
+		}  
+		else Bh.ReMakeTriangleContainingTheVertex();     
+
+		// generation of the list of next Triangle 
+		for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
+		// the next traingle of i is -first_np_or_next_t[i]
+
+		// Big loop (most time consuming)
+		int iter=0;
+		if (verbose>5) printf("         Big loop\n");
+		do {
+			/*Update variables*/
+			iter++;
+			nbtold=nbt;
+			nbvold=nbv;
+
+			/*We test all triangles*/
+			i=Headt;
+			next_t=-first_np_or_next_t[i];
+			for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
+
+				//check i
+				if (i<0 || i>=nbt ){
+					ISSMERROR("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1);
+				}
+				//change first_np_or_next_t[i]
+				first_np_or_next_t[i] = iter; 
+
+				//Loop over the edges of t
+				for(j=0;j<3;j++){
+					TriangleAdjacent tj(t,j);
+					BamgVertex &vA = *tj.EdgeVertex(0);
+					BamgVertex &vB = *tj.EdgeVertex(1);
+
+					//if t is a boundary triangle, or tj locked, continue
+					if (!t->link)     continue;
+					if (t->det <0)    continue;
+					if (t->Locked(j)) continue;
+
+					TriangleAdjacent tadjj = t->Adj(j);	  
+					Triangle* ta=tadjj;
+
+					//if the adjacent triangle is a boundary triangle, continur
+					if (ta->det<0) continue;	  
+
+					R2 A=vA;
+					R2 B=vB;
+					k=Number(ta);
+
+					//if this edge has already been done, go to next edge of triangle
+					if(first_np_or_next_t[k]==iter) continue;
+
+					lIntTria.SplitEdge(Bh,A,B);
+					lIntTria.NewPoints(vertices,nbv,maxnbv);
+				  } // end loop for each edge 
+			  }// for triangle   
+
+			if (!InsertNewPoints(nbvold,NbTSwap)) break;
+			for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
+			Headt = nbt; // empty list 
+
+			// for all the triangle containing the vertex i
+			for (i=nbvold;i<nbv;i++){ 
+				BamgVertex*          s  = vertices + i;
+				TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
+				Triangle*        tbegin= (Triangle*) ta;
+				long kt;
+				do { 
+					kt = Number((Triangle*) ta);
+					if (first_np_or_next_t[kt]>0){
+						first_np_or_next_t[kt]=-Headt;
+						Headt=kt;
+					}
+					if (ta.EdgeVertex(0)!=s){
+						ISSMERROR("ta.EdgeVertex(0)!=s");
+					}
+					ta = Next(Adj(ta));
+				} while ( (tbegin != (Triangle*) ta)); 
+			}
+
+		} while (nbv!=nbvold);
+
+		delete []  first_np_or_next_t;
+
+		long NbSwapf =0,NbSwp;
+
+		NbSwp = NbSwapf;
+		for (i=0;i<nbv;i++)
+		 NbSwapf += vertices[i].Optim(0);
+		NbTSwap +=  NbSwapf ;
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::ProjectOnCurve{{{1*/
+	GeometricalEdge*   Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex &  vA, BamgVertex & vB,
+				double theta,BamgVertex & R,VertexOnEdge &  BR,VertexOnGeom & GR) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/
+
+		void *pA=0,*pB=0;
+		double tA=0,tB=0;
+		R2 A=vA,B=vB;
+		BamgVertex * pvA=&vA, * pvB=&vB;
+		if (vA.vint == IsVertexOnVertex){
+			pA=vA.onBackgroundVertex;
+		}
+		else if (vA.vint == IsVertexOnEdge){
+			pA=vA.onBackgroundEdge->be;
+			tA=vA.onBackgroundEdge->abcisse;
+		}
+		else {
+			ISSMERROR("ProjectOnCurve On BamgVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA));
+		} 
+
+		if (vB.vint == IsVertexOnVertex){
+			pB=vB.onBackgroundVertex;
+		}
+		else if(vB.vint == IsVertexOnEdge){
+			pB=vB.onBackgroundEdge->be;
+			tB=vB.onBackgroundEdge->abcisse;
+		}
+		else {
+			ISSMERROR("ProjectOnCurve On BamgVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB));
+		} 
+		Edge * e = &BhAB;
+		if (!pA || !pB || !e){
+			ISSMERROR("!pA || !pB || !e");
+		}
+		// be carefull the back ground edge e is on same geom edge 
+		// of the initiale edge def by the 2 vertex A B;
+		//check Is a background Mesh;   
+		if (e<BTh.edges || e>=BTh.edges+BTh.nbe){
+			ISSMERROR("e<BTh.edges || e>=BTh.edges+BTh.nbe");
+		}
+		// walk on BTh edge 
+		//not finish ProjectOnCurve with BackGround Mesh);
+		// 1 first find a back ground edge contening the vertex A
+		// 2 walk n back gound boundary to find the final vertex B
+
+		if( vA.vint == IsVertexOnEdge) 
+		  { // find the start edge 
+			e = vA.onBackgroundEdge->be;	 
+
+		  } 
+		else if (vB.vint == IsVertexOnEdge) 
+		  {
+			theta = 1-theta;
+			Exchange(tA,tB);
+			Exchange(pA,pB);
+			Exchange(pvA,pvB);
+			Exchange(A,B);
+			e =  vB.onBackgroundEdge->be;
+
+		  } 
+		else{ // do the search by walking 
+			ISSMERROR("case not supported yet");
+		  }
+
+		// find the direction of walking with sens of edge and pA,PB;
+		R2 AB=B-A;
+
+		double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
+		int kkk=0;
+		int sens = (cosE01AB>0) ? 1 : 0;
+
+		//   double l=0; // length of the edge AB
+		double abscisse = -1;
+
+		for (int cas=0;cas<2;cas++){
+			// 2 times algo:
+			//    1 for computing the length l
+			//    2 for find the vertex 
+			int  iii;
+			BamgVertex  *v0=pvA,*v1; 
+			Edge *neee,*eee;
+			double lg =0; // length of the curve 
+			double te0;
+			// we suppose take the curve's abcisse 
+			for ( eee=e,iii=sens,te0=tA;
+						eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
+						neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { 
+
+				kkk=kkk+1;
+				if (kkk>=100){
+					ISSMERROR("kkk>=100");
+				}
+				if (!eee){
+					ISSMERROR("!eee");
+				}
+				double lg0 = lg;
+				double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
+				lg += dp;
+				if (cas && abscisse <= lg) { // ok we find the geom edge 
+					double sss  =   (abscisse-lg0)/dp;
+					double thetab = te0*(1-sss)+ sss*iii;
+					if (thetab<0 || thetab>1){
+						ISSMERROR("thetab<0 || thetab>1");
+					}
+					BR = VertexOnEdge(&R,eee,thetab);
+					return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
+				  }
+			  }
+			// we find the end 
+			if (v1 != pvB){
+				if (( void*) v1 == pB)
+				 tB = iii;
+
+				double lg0 = lg;
+				if (!eee){
+					ISSMERROR("!eee");
+				}
+				v1 = pvB;
+				double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
+				lg += dp;	
+				abscisse = lg*theta;
+				if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
+				  { // ok we find the geom edge 
+					double sss  =   (abscisse-lg0)/dp;
+					double thetab = te0*(1-sss)+ sss*tB;
+					if (thetab<0 || thetab>1){
+						ISSMERROR("thetab<0 || thetab>1");
+					}
+					BR = VertexOnEdge(&R,eee,thetab);
+					return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
+				  }
+			  }
+			abscisse = lg*theta;
+
+		  }
+		ISSMERROR("Big bug...");
+		return 0; // just for the compiler 
+	}                  
+	/*}}}1*/
+/*FUNCTION Mesh::ReconstructExistingMesh{{{1*/
+void Mesh::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 verbose=0;
+
+	// generation of the integer coordinate
+
+	// find extrema coordinates of vertices pmin,pmax
+	long i;
+	if(verbose>2) printf("      Reconstruct mesh of %i vertices\n",nbv); 
+
+	//initialize ordre
+	ISSMASSERT(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){ 
+		ISSMERROR("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) {
+				ISSMASSERT(!triangles[i].TriangleAdj(j));
+				ISSMASSERT(!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 {
+				ISSMERROR("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(verbose>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) {
+		ISSMERROR("%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) {
+		 ISSMERROR("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*/
+	BamgVertex *  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++) {
+		BamgVertex *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
+			BamgVertex &a=edges[i][0], &b=edges[i][1];
+			if (a.t && b.t){
+				knbe++;
+				if (ForceEdge(a,b,ta)<0) nbloss++;
+			}
+		}
+	}
+	if(nbloss) {
+		ISSMERROR("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;
+					BamgVertex *v0= ta.EdgeVertex(0);
+					BamgVertex *v1= ta.EdgeVertex(1);
+					long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
+
+					ISSMASSERT(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((BamgVertex *) NULL,(BamgVertex *) NULL,(BamgVertex *) NULL);
+			triangles[i].color=-1;
+		}
+		else{
+			triangles[i].color= savenbt+ NbTfillHoll++;
+		}
+	}
+	ISSMASSERT(savenbt+NbTfillHoll<=savenbtx);
+
+	// copy of the outside triangles in saveMesh 
+	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) {
+		ISSMERROR("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 BamgVertex number %i\n",Number(edges[i][j])+1,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",Number(edges[i][j])+1,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");
+						ISSMERROR("See above (might be cryptic...)");
+					}
+				}
+			}
+		}
+	}
+}
+/*}}}1*/
+	/*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/
+	void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
+
+		long int verbose=0;
+		long *renu= new long[nbt];
+		register Triangle *t0,*t,*te=triangles+nbt;
+		register long k=0,it,i,j;
+
+		for ( it=0;it<nbt;it++) 
+		 renu[it]=-1; // outside triangle 
+		for ( i=0;i<NbSubDomains;i++)
+		  { 
+			t=t0=subdomains[i].head;
+			if (!t0){ // not empty sub domain
+				ISSMERROR("!t0");
+			}
+			do { 
+				long kt = Number(t);
+				if (kt<0 || kt >= nbt ){
+					ISSMERROR("kt<0 || kt >= nbt");
+				}
+				if (renu[kt]!=-1){
+					ISSMERROR("renu[kt]!=-1");
+				}
+				renu[kt]=k++;
+			}
+			while (t0 != (t=t->link));
+		  }
+		// take is same numbering if possible    
+		if(justcompress)
+		 for ( k=0,it=0;it<nbt;it++) 
+		  if(renu[it] >=0 ) 
+			renu[it]=k++;
+
+		// put the outside triangles at the end
+		for ( it=0;it<nbt;it++){
+			if (renu[it]==-1) renu[it]=k++;
+		}
+		if (k != nbt){
+			ISSMERROR("k != nbt");
+		}
+		// do the change on all the pointeur 
+		for ( it=0;it<nbt;it++)
+		 triangles[it].ReNumbering(triangles,te,renu);
+
+		for ( i=0;i<NbSubDomains;i++)
+		 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
+
+		// move the Triangles  without a copy of the array 
+		// be carefull not trivial code 
+		for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
+		 if (renu[it] >= 0) // a new sub cycle
+			{ 
+			 i=it;
+			 Triangle ti=triangles[i],tj;
+			 while ( (j=renu[i]) >= 0) 
+				{ // i is old, and j is new 
+				 renu[i] = -1; // mark 
+				 tj = triangles[j]; // save new
+				 triangles[j]= ti; // new <- old
+				 i=j;     // next 
+				 ti = tj;
+				}  
+			}
+		delete [] renu;
+		nt = nbt - NbOutT;
+
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::ReNumberingVertex{{{1*/
+	void Mesh::ReNumberingVertex(long * renu) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
+
+		// warning be carfull because pointer
+		// from on mesh to over mesh 
+		//  --  so do ReNumbering at the beginning
+		BamgVertex * ve = vertices+nbv;
+		long it,ie,i;
+
+		printf("renumbering triangles\n");
+		for ( it=0;it<nbt;it++) 
+		 triangles[it].ReNumbering(vertices,ve,renu);
+
+		printf("renumbering edges\n");
+		for ( ie=0;ie<nbe;ie++) 
+		 edges[ie].ReNumbering(vertices,ve,renu);
+
+		printf("renumbering vertices on geom\n");
+		for (i=0;i< NbVerticesOnGeomVertex;i++)
+		  {
+			BamgVertex *v = VerticesOnGeomVertex[i].mv;
+			if (v>=vertices && v < ve)
+			 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];
+		  }
+
+		printf("renumbering vertices on edge\n");
+		for (i=0;i< NbVerticesOnGeomEdge;i++)
+		  {
+			BamgVertex *v =VerticesOnGeomEdge[i].mv;
+			if (v>=vertices && v < ve)
+			 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];
+		  }
+
+		printf("renumbering vertices on Bth vertex\n");
+		for (i=0;i< NbVertexOnBThVertex;i++)
+		  {
+			BamgVertex *v=VertexOnBThVertex[i].v;
+			if (v>=vertices && v < ve)
+			 VertexOnBThVertex[i].v=vertices+renu[Number(v)];
+		  }
+
+		for (i=0;i< NbVertexOnBThEdge;i++)
+		  {
+			BamgVertex *v=VertexOnBThEdge[i].v;
+			if (v>=vertices && v < ve)
+			 VertexOnBThEdge[i].v=vertices+renu[Number(v)];
+		  }
+
+		// move the Vertices without a copy of the array 
+		// be carefull not trivial code 
+		long j;
+		for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
+		 if (renu[it] >= 0) // a new sub cycle
+			{ 
+			 i=it;
+			 BamgVertex ti=vertices[i],tj;
+			 while ( (j=renu[i]) >= 0){
+				 // i is old, and j is new 
+				 renu[i] = -1-renu[i]; // mark 
+				 tj = vertices[j];     // save new
+				 vertices[j]= ti;      // new <- old
+				 i=j;     // next 
+				 ti = tj;
+				}  
+			}
+		if (quadtree){
+			delete quadtree;
+			quadtree = new QuadTree(this);
+		}
+
+		for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
+	}
+	/*}}}1*/
+/*FUNCTION Mesh::SetIntCoor{{{1*/
+void Mesh::SetIntCoor(const char * strfrom) {
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/
+
+	/*Set integer coordinate for existing vertices*/
+
+	//Get extrema coordinates of the existing vertices
+	pmin =  vertices[0].r;
+	pmax =  vertices[0].r;
+	long i;
+	for (i=0;i<nbv;i++) {
+		pmin.x = Min(pmin.x,vertices[i].r.x);
+		pmin.y = Min(pmin.y,vertices[i].r.y);
+		pmax.x = Max(pmax.x,vertices[i].r.x);
+		pmax.y = Max(pmax.y,vertices[i].r.y);
+	}
+	R2 DD = (pmax-pmin)*0.05;
+	pmin = pmin-DD;
+	pmax = pmax+DD; 
+
+	//Compute coefIcoor
+	coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
+	if (coefIcoor<=0){
+		ISSMERROR("coefIcoor should be positive, a problem in the geometry is likely");
+	}
+
+	// generation of integer coord  
+	for (i=0;i<nbv;i++) {
+		vertices[i].i = toI2(vertices[i].r);    
+	}
+
+	// computation of the det 
+	int number_of_errors=0;
+	for (i=0;i<nbt;i++) {
+		BamgVertex & v0 = triangles[i][0];
+		BamgVertex & v1 = triangles[i][1];
+		BamgVertex & v2 = triangles[i][2];
+
+		//If this is not a boundary triangle
+		if ( &v0 && &v1 &&  &v2 ){
+
+			/*Compute determinant*/
+			triangles[i].det= det(v0,v1,v2);
+
+			/*Check that determinant is positive*/
+			if (triangles[i].det <=0){
+
+				/*increase number_of_errors and print error only for the first 20 triangles*/
+				number_of_errors++;
+				if (number_of_errors<20){
+					printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det);
+				}
+			}
+		}
+
+		//else, set as -1
+		else triangles[i].det=-1;
+	}
+
+	if (number_of_errors) ISSMERROR("Fatal error: some triangles have negative areas, see above");
+}
+/*}}}1*/
+/*FUNCTION Mesh::ShowRegulaty{{{1*/
+void  Mesh::ShowRegulaty() const {
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
+
+	const  double  sqrt32=sqrt(3.)*0.5; 
+	const double  aireKh=sqrt32*0.5;
+	D2  Beq(1,0),Heq(0.5,sqrt32);
+	D2xD2 Br(D2xD2(Beq,Heq).t());
+	D2xD2 B1r(Br.inv());
+	double gammamn=1e100,hmin=1e100;
+	double gammamx=0,hmax=0;
+	double beta=1e100;
+	double beta0=0;
+	double  alpha2=0;
+	double area=0,Marea=0;
+	// double cf= double(coefIcoor);
+	// double cf2= 6.*cf*cf;
+	int nt=0;
+	for (int it=0;it<nbt;it++)
+	 if ( triangles[it].link) 
+		{
+		 nt++;
+		 Triangle &K=triangles[it];
+		 double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
+		 area+= area3;
+		 D2xD2 B_Kt(K[0],K[1],K[2]);
+		 D2xD2 B_K(B_Kt.t());
+		 D2xD2 B1K = Br*B_K.inv();
+		 D2xD2 BK =  B_K*B1r;
+		 D2xD2 B1B1 = B1K.t()*B1K;
+		 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
+		 MatVVP2x2 VMK(MK);
+		 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
+		 double betaK=0;
+
+		 for (int j=0;j<3;j++)
+			{
+			 double he= Norme2(R2(K[j],K[(j+1)%3]));
+			 hmin=Min(hmin,he);
+			 hmax=Max(hmax,he);
+			 BamgVertex & v=K[j];
+			 D2xD2 M((Metric)v);
+			 betaK += sqrt(M.det());
+
+			 D2xD2 BMB = BK.t()*M*BK;
+			 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
+			 MatVVP2x2 VM1(M1);
+			 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
+			 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);		
+			}
+		 betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
+		 Marea+= betaK;
+		 beta=min(beta,betaK);
+		 beta0=max(beta0,betaK);
+		}   
+	area*=3; 
+	gammamn=sqrt(gammamn);
+	gammamx=sqrt(gammamx);    
+	printf("   Adaptmesh info:\n");
+	printf("      number of triangles = %i\n",nt);
+	printf("      hmin = %g, hmax=%g\n",hmin,hmax);
+	printf("      area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt));
+	printf("      infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx);
+	printf("      anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5));
+}
+/*}}}1*/
+/*FUNCTION Mesh::ShowHistogram{{{1*/
+void  Mesh::ShowHistogram() const {
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
+
+	const long kmax=10;
+	const double llmin = 0.5,llmax=2;
+	const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
+	long histo[kmax+1];
+	long i,it,k, nbedges =0;
+	for (i=0;i<=kmax;i++) histo[i]=0;
+	for (it=0;it<nbt;it++)
+	 if ( triangles[it].link) 
+		{
+
+		 for (int j=0;j<3;j++)
+			{
+			 Triangle *ta = triangles[it].TriangleAdj(j);
+			 if ( !ta || !ta->link || Number(ta) >= it) 
+				{ 
+				 BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
+				 BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
+				 if ( !& vP || !&vQ) continue;
+				 R2 PQ = vQ.r - vP.r;
+				 double l = log(LengthInterpole(vP,vQ,PQ));
+				 nbedges++;
+				 k = (int) ((l - lmin)*delta);
+				 k = Min(Max(k,0L),kmax);
+				 histo[k]++;
+				}
+			}
+		}  
+	printf(" --- Histogram of the unit mesh,  nb of edges = %i\n",nbedges);
+	printf("      length of edge in   | %% of edge  | Nb of edges \n"); 
+	printf("      --------------------+-------------+-------------\n"); 
+	for   (i=0;i<=kmax;i++){ 
+		if (i==0) printf("      %10i",0);
+		else      printf("      %10g",exp(lmin+i/delta));
+		if (i==kmax) printf("          +inf   ");
+		else printf("      %10g",exp(lmin+(i+1)/delta));
+		printf("|  %10g |\n",((long)  ((10000.0 * histo[i])/ nbedges))/100.0);
+		printf("  %i\n",histo[i]);
+	}
+	printf("      --------------------+-------------+-------------\n"); 
+}
+/*}}}1*/
+/*FUNCTION Mesh::SmoothingVertex{{{1*/
+void Mesh::SmoothingVertex(int nbiter,double omega ) { 
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
+
+	long int verbose=0;
+	//  if quatree exist remove it end reconstruct
+	if (quadtree) delete quadtree;
+	quadtree=0;
+	ReMakeTriangleContainingTheVertex();
+	Triangle vide; // a triangle to mark the boundary vertex
+	Triangle   ** tstart= new Triangle* [nbv];
+	long i,j,k;
+	//   attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide 
+	if ( this == & BTh)
+	 for ( i=0;i<nbv;i++)
+	  tstart[i]=vertices[i].t;     
+	else 
+	 for ( i=0;i<nbv;i++)
+	  tstart[i]=0;
+	for ( j=0;j<NbVerticesOnGeomVertex;j++ ) 
+	 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;
+	for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 
+	 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
+	if(verbose>2) printf("   SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);
+	for (k=0;k<nbiter;k++)
+	  {
+		long i,NbSwap =0;
+		double delta =0;
+		for ( i=0;i<nbv;i++)
+		 if (tstart[i] != &vide) // not a boundary vertex 
+		  delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
+		if (!NbOfQuad)
+		 for ( i=0;i<nbv;i++)
+		  if (tstart[i] != &vide) // not a boundary vertex 
+			NbSwap += vertices[i].Optim(1);
+		if (verbose>3) printf("      move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);
+	  }
+
+	delete [] tstart;
+	if (quadtree) quadtree= new QuadTree(this);
+}
+/*}}}1*/
+/*FUNCTION Mesh::SmoothMetric{{{1*/
+void Mesh::SmoothMetric(double raisonmax) { 
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
+
+	long int verbose=0;
+
+	if(raisonmax<1.1) return;
+	if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
+	ReMakeTriangleContainingTheVertex();
+	long i,j,kch,kk,ip;
+	long *first_np_or_next_t0 = new long[nbv];
+	long *first_np_or_next_t1 = new long[nbv];
+	long Head0 =0,Head1=-1;
+	double logseuil= log(raisonmax);
+
+	for(i=0;i<nbv-1;i++)
+	 first_np_or_next_t0[i]=i+1; 
+	first_np_or_next_t0[nbv-1]=-1;// end;
+	for(i=0;i<nbv;i++)
+	 first_np_or_next_t1[i]=-1;
+	kk=0;
+	while (Head0>=0&& kk++<100) {
+		kch=0;
+		for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
+			//  pour tous les triangles autour du sommet s
+			register Triangle* t= vertices[i].t;
+			if (!t){
+				ISSMERROR("!t");
+			}
+			BamgVertex & vi = vertices[i];
+			TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
+			BamgVertex *pvj0 = ta.EdgeVertex(0);
+			while (1) {
+				ta=Previous(Adj(ta));
+				if (vertices+i != ta.EdgeVertex(1)){
+					ISSMERROR("vertices+i != ta.EdgeVertex(1)");
+				}
+				BamgVertex & vj = *(ta.EdgeVertex(0));
+				if ( &vj ) {
+					j= &vj-vertices;
+					if (j<0 || j >= nbv){
+						ISSMERROR("j<0 || j >= nbv");
+					}
+					R2 Aij = (R2) vj - (R2) vi;
+					double ll =  Norme2(Aij);
+					if (0) {  
+						double hi = ll/vi.m(Aij);
+						double hj = ll/vj.m(Aij);
+						if(hi < hj)
+						  {
+							double dh=(hj-hi)/ll;
+							if (dh>logseuil) {
+								vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
+								if(first_np_or_next_t1[j]<0)
+								 kch++,first_np_or_next_t1[j]=Head1,Head1=j;
+							}
+						  }
+					} 
+					else
+					  {
+						double li = vi.m(Aij);
+						if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
+						 if(first_np_or_next_t1[j]<0) // if the metrix change 
+						  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
+					  }
+				}
+				if  ( &vj ==  pvj0 ) break;
+			}
+		}
+		Head0 = Head1;
+		Head1 = -1;
+		Exchange(first_np_or_next_t0,first_np_or_next_t1);
+	}
+	if(verbose>2) printf("      number of iterations = %i\n",kch); 
+	delete [] first_np_or_next_t0;
+	delete [] first_np_or_next_t1;
+}
+/*}}}1*/
+	/*FUNCTION Mesh::SplitElement{{{1*/
+	int  Mesh::SplitElement(int choice){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
+
+		long int verbose=0;
+
+		Direction NoDirOfSearch;
+		const  int withBackground = &BTh != this && &BTh;
+
+		ReNumberingTheTriangleBySubDomain();
+		//int nswap =0;
+		const long nfortria( choice ? 4 : 6);
+		if(withBackground) 
+		  {
+			BTh.SetVertexFieldOn();
+			SetVertexFieldOnBTh();
+		  }
+		else
+		 BTh.SetVertexFieldOn();
+
+		long newnbt=0,newnbv=0;
+		long * kedge = 0;
+		long newNbOfQuad=NbOfQuad;
+		long * ksplit= 0, * ksplitarray=0;
+		long kkk=0;
+		int ret =0;
+		if (maxnbv<nbv+nbe) return 1;//   
+		// 1) create  the new points by spliting the internal edges 
+		// set the 
+		long nbvold = nbv;
+		long nbtold = nbt;
+		long NbOutTold  = NbOutT;
+		long  NbEdgeOnGeom=0;
+		long i;
+
+		nbt = nbt - NbOutT; // remove all the  the ouside triangles 
+		long nbtsave = nbt;
+		Triangle * lastT = triangles + nbt;
+		for (i=0;i<nbe;i++)
+		 if(edges[i].onGeometry) NbEdgeOnGeom++;
+		long newnbe=nbe+nbe;
+		//  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
+		long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
+		// long newNbVertexOnBThVertex=NbVertexOnBThVertex;
+		long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
+
+		// do allocation for pointeur to the geometry and background
+		VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
+		VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
+		if (NbVerticesOnGeomEdge)
+		 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
+		if (NbVertexOnBThEdge)
+		 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
+		Edge *newedges = new Edge [newnbe];
+		//  memcpy(newedges,edges,sizeof(Edge)*nbe);
+		SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
+		long k=nbv;
+		long kk=0;
+		long kvb = NbVertexOnBThEdge;
+		long kvg = NbVerticesOnGeomEdge;
+		long ie =0;
+		Edge ** edgesGtoB=0;
+		if (withBackground)
+		 edgesGtoB= BTh.MakeGeometricalEdgeToEdge();
+		long ferr=0;
+		for (i=0;i<nbe;i++)
+		 newedges[ie].onGeometry=0;
+
+		for (i=0;i<nbe;i++)
+		  {
+			GeometricalEdge *ong =  edges[i].onGeometry;
+
+			newedges[ie]=edges[i];
+			newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
+			newedges[ie].adj[1]=newedges + ie +1;
+			R2 A = edges[i][0],B = edges[i][1];
+
+
+			kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
+			if (ong) // a geometrical edges 
+			  { 
+				if (withBackground){
+					// walk on back ground mesh 
+					//  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 
+					// a faire -- difficile 
+					// the first PB is to now a background edge between the 2 vertices
+					if (!edgesGtoB){
+						ISSMERROR("!edgesGtoB");
+					}
+					ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)],
+								edges[i][0],edges[i][1],0.5,vertices[k],
+								newVertexOnBThEdge[kvb],
+								newVerticesOnGeomEdge[kvg++]);
+					vertices[k].ReferenceNumber= edges[i].ref;
+					vertices[k].DirOfSearch =   NoDirOfSearch;        
+					;
+					// get the Info on background mesh 
+					double s =        newVertexOnBThEdge[kvb];
+					BamgVertex &  bv0  = newVertexOnBThEdge[kvb][0];
+					BamgVertex &  bv1  = newVertexOnBThEdge[kvb][1];
+					// compute the metrix of the new points 
+					vertices[k].m =  Metric(1-s,bv0,s,bv1); 
+					kvb++;
+				  }
+				else 
+				  {
+					ong=Gh.ProjectOnCurve(edges[i],
+								0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
+					// vertices[k].i = toI2( vertices[k].r);
+					vertices[k].ReferenceNumber = edges[i].ref;
+					vertices[k].DirOfSearch = NoDirOfSearch;
+					vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);	      
+				  }  
+			  }
+			else // straigth line edge ---
+			  { 
+				vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
+				vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
+				vertices[k].onGeometry = 0;
+			  }
+			//vertices[k].i = toI2( vertices[k].r);
+			R2 AB =  vertices[k].r;
+			R2 AA = (A+AB)*0.5;
+			R2 BB = (AB+B)*0.5;
+			vertices[k].ReferenceNumber = edges[i].ref;
+			vertices[k].DirOfSearch = NoDirOfSearch;
+
+			newedges[ie].onGeometry = Gh.Containing(AA,ong);
+			newedges[ie++].v[1]=vertices+k;
+
+			newedges[ie]=edges[i];
+			newedges[ie].adj[0]=newedges + ie -1;
+			newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
+			newedges[ie].onGeometry =  Gh.Containing(BB,ong);
+			newedges[ie++].v[0]=vertices+k;
+			k++;
+		  }
+		if (edgesGtoB) delete [] edgesGtoB;
+		edgesGtoB=0;
+
+		newnbv=k;
+		newNbVerticesOnGeomEdge=kvg;
+		if (newnbv> maxnbv) goto Error;// bug 
+
+		nbv = k;
+
+
+		kedge = new long[3*nbt+1];
+		ksplitarray = new long[nbt+1];
+		ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
+
+		for (i=0;i<3*nbt;i++)
+		 kedge[i]=-1;
+
+		//  
+
+		for (i=0;i<nbt;i++) {
+			Triangle & t = triangles[i];
+			if (!t.link){
+				ISSMERROR("!t.link");
+			}
+			for(int j=0;j<3;j++)
+			  {
+				const TriangleAdjacent ta = t.Adj(j);
+				const Triangle & tt = ta;
+				if (&tt >= lastT)
+				 t.SetAdj2(j,0,0);// unset adj
+				const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
+				const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
+				long  ke =edge4->SortAndFind(Number(v0),Number(v1));
+				if (ke>0) 
+				  {
+					long ii = Number(tt);
+					int  jj = ta;
+					long ks = ke + nbvold;
+					kedge[3*i+j] = ks;
+					if (ii<nbt) // good triangle
+					 kedge[3*ii+jj] = ks;
+					BamgVertex &A=vertices[ks];
+					double aa,bb,cc,dd;
+					if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
+						// warning PB roundoff error 
+						if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0 
+										||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0  
+										||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
+							printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd);
+							ISSMERROR("Number of triangles with P2 interpolation Problem");
+						}
+					}
+					else {
+						if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0 
+										||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0 
+										||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
+							printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd);
+							ISSMERROR("Number of triangles with P2 interpolation Problem");
+						}
+					} 
+				  }
+			  }
+		}
+
+		for (i=0;i<nbt;i++){
+			ksplit[i]=1; // no split by default
+			const Triangle & t = triangles[ i];
+			int nbsplitedge =0;
+			int nbinvisible =0;
+			int invisibleedge=0;
+			int kkk[3];      
+			for (int j=0;j<3;j++)
+			  {
+				if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
+
+				const TriangleAdjacent ta = t.Adj(j);
+				const Triangle & tt = ta;
+
+
+				const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
+				const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
+				if ( kedge[3*i+j] < 0) 
+				  {
+					long  ke =edge4->SortAndFind(Number(v0),Number(v1));
+					if (ke<0) // new 
+					  {
+						if (&tt) // internal triangles all the boundary 
+						  { // new internal edges 
+							long ii = Number(tt);
+							int  jj = ta;
+
+							kedge[3*i+j]=k;// save the vertex number 
+							kedge[3*ii+jj]=k;
+							if (k<maxnbv) 
+							  {
+								vertices[k].r = ((R2) v0+(R2) v1 )/2;
+								//vertices[k].i = toI2( vertices[k].r);
+								vertices[k].ReferenceNumber=0;
+								vertices[k].DirOfSearch =NoDirOfSearch;
+								vertices[k].m =  Metric(0.5,v0,0.5,v1);
+							  }
+							k++;
+							kkk[nbsplitedge++]=j;		      
+						  } // tt 
+						else
+						 ISSMERROR("Bug...");
+					  } // ke<0	       
+					else
+					  { // ke >=0
+						kedge[3*i+j]=nbvold+ke;
+						kkk[nbsplitedge++]=j;// previously splited
+					  }
+				  }
+				else 
+				 kkk[nbsplitedge++]=j;// previously splited
+
+			  } 
+			if (nbinvisible>=2){
+				ISSMERROR("nbinvisible>=2");
+			}
+			switch (nbsplitedge) {
+				case 0: ksplit[i]=10; newnbt++; break;   // nosplit
+				case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 
+				case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 
+				case 3:
+						  if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
+						  else   ksplit[i]=10*nfortria,newnbt+=nfortria;
+						  break;
+			} 
+			if (ksplit[i]<40){
+				ISSMERROR("ksplit[i]<40");
+			}
+		  }
+		//  now do the element split
+		newNbOfQuad = 4*NbOfQuad;
+		nbv = k;
+		kkk = nbt;
+		ksplit[-1] = nbt;
+		// look on  old true  triangles 
+
+		for (i=0;i<nbtsave;i++){
+			int  nbmkadj=0;
+			long mkadj [100];
+			mkadj[0]=i;
+			long kk=ksplit[i]/10;
+			int  ke=(int) (ksplit[i]%10);
+			if (kk>=7 || kk<=0){
+				ISSMERROR("kk>=7 || kk<=0");
+			}
+
+			// def the numbering   k (edge) i vertex 
+			int k0 = ke;
+			int k1 = NextEdge[k0];
+			int k2 = PreviousEdge[k0];
+			int i0 = OppositeVertex[k0];
+			int i1 = OppositeVertex[k1];
+			int i2 = OppositeVertex[k2];
+
+			Triangle &t0=triangles[i];
+			BamgVertex * v0=t0(i0);           
+			BamgVertex * v1=t0(i1);           
+			BamgVertex * v2=t0(i2);
+
+			if (nbmkadj>=10){
+				ISSMERROR("nbmkadj>=10");
+			}
+			// --------------------------
+			TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
+			// save the flag Hidden
+			int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
+			// un set all adj -- save Hidden flag --
+			t0.SetAdj2(0,0,hid[0]);
+			t0.SetAdj2(1,0,hid[1]);
+			t0.SetAdj2(2,0,hid[2]);
+			// --  remake 
+			switch  (kk) {
+				case 1: break;// nothing 
+				case 2: // 
+						  {
+							Triangle &t1=triangles[kkk++];
+							t1=t0;
+							if (kedge[3*i+i0]<0){
+								ISSMERROR("kedge[3*i+i0]<0");
+							}
+							BamgVertex * v3 = vertices + kedge[3*i+k0];
+
+							t0(i2) = v3;
+							t1(i1) = v3;
+							t0.SetAllFlag(k2,0);
+							t1.SetAllFlag(k1,0);
+						  } 
+						break; 
+				case 3: //
+						  {
+							Triangle &t1=triangles[kkk++];
+							Triangle &t2=triangles[kkk++];
+							t2=t1=t0;
+							if (kedge[3*i+k1]<0){
+								ISSMERROR("kedge[3*i+k1]<0");
+							}
+							if (kedge[3*i+k2]<0){
+								ISSMERROR("kedge[3*i+k2]<0");
+							}
+
+							BamgVertex * v01 = vertices + kedge[3*i+k2];
+							BamgVertex * v02 = vertices + kedge[3*i+k1]; 
+							t0(i1) = v01; 
+							t0(i2) = v02; 
+							t1(i2) = v02;
+							t1(i0) = v01; 
+							t2(i0) = v02; 
+							t0.SetAllFlag(k0,0);
+							t1.SetAllFlag(k1,0);
+							t1.SetAllFlag(k0,0);
+							t2.SetAllFlag(k2,0);
+						  } 
+						break;
+				case 4: // 
+				case 6: // split in 4 
+						  {
+							Triangle &t1=triangles[kkk++];
+							Triangle &t2=triangles[kkk++];
+							Triangle &t3=triangles[kkk++];
+							t3=t2=t1=t0;
+							if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
+								ISSMERROR("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
+							}
+							BamgVertex * v12 = vertices + kedge[3*i+k0];
+							BamgVertex * v02 = vertices + kedge[3*i+k1]; 
+							BamgVertex * v01 = vertices + kedge[3*i+k2];
+							t0(i1) = v01;
+							t0(i2) = v02;
+							t0.SetAllFlag(k0,hid[k0]);
+
+							t1(i0) = v01;
+							t1(i2) = v12;
+							t0.SetAllFlag(k1,hid[k1]);
+
+							t2(i0) = v02;
+							t2(i1) = v12;
+							t2.SetAllFlag(k2,hid[k2]);
+
+							t3(i0) = v12;
+							t3(i1) = v02;
+							t3(i2) = v01;
+
+							t3.SetAllFlag(0,hid[0]);	   
+							t3.SetAllFlag(1,hid[1]);	   
+							t3.SetAllFlag(2,hid[2]);
+
+							if ( kk == 6)
+							  {
+
+								Triangle &t4=triangles[kkk++];
+								Triangle &t5=triangles[kkk++];
+
+								t4 = t3;
+								t5 = t3;
+
+								t0.SetHidden(k0);
+								t1.SetHidden(k1);
+								t2.SetHidden(k2);
+								t3.SetHidden(0);
+								t4.SetHidden(1);
+								t5.SetHidden(2);
+
+								if (nbv < maxnbv ) 
+								  {
+									vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
+									vertices[nbv].ReferenceNumber =0;
+									vertices[nbv].DirOfSearch =NoDirOfSearch;
+									//vertices[nbv].i = toI2(vertices[nbv].r);
+									double a3[]={1./3.,1./3.,1./3.};
+									vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
+									BamgVertex * vc =  vertices +nbv++;
+									t3(i0) = vc;
+									t4(i1) = vc;
+									t5(i2) = vc;
+
+								  }
+								else
+								 goto Error; 
+							  }
+
+						  } 
+						break;         
+			}
+
+			// save all the new triangles
+			mkadj[nbmkadj++]=i;
+			long jj;
+			if (t0.link) 
+			 for (jj=nbt;jj<kkk;jj++)
+				{
+				 triangles[jj].link=t0.link;
+				 t0.link= triangles+jj;
+				 mkadj[nbmkadj++]=jj;
+				}
+			if (nbmkadj>13){// 13 = 6 + 4 +
+				ISSMERROR("nbmkadj>13");
+			}
+
+			if (kk==6)  newNbOfQuad+=3;
+			for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
+			ksplit[i]= nbt; // save last adresse of the new triangles
+			kkk = nbt;
+		  }
+
+		for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
+
+		if(withBackground)
+		 for (i=0;i<BTh.nbv;i++)
+		  BTh.vertices[i].m =  BTh.vertices[i].m*2.;
+
+
+		ret = 2;
+		if (nbt>= nbtx) goto Error; // bug 
+		if (nbv>= maxnbv) goto Error; // bug 
+		// generation of the new triangles 
+
+		SetIntCoor("In SplitElement"); 
+
+		ReMakeTriangleContainingTheVertex();
+		if(withBackground)
+		 BTh.ReMakeTriangleContainingTheVertex();
+
+		delete [] edges;
+		edges = newedges;
+		nbe = newnbe;
+		NbOfQuad = newNbOfQuad;
+
+		for (i=0;i<NbSubDomains;i++)
+		  { 
+			long k = subdomains[i].edge- edges;
+			subdomains[i].edge =  edges+2*k; // spilt all edge in 2 
+		  }
+
+		if (ksplitarray) delete [] ksplitarray;
+		if (kedge) delete [] kedge;
+		if (edge4) delete edge4;
+		if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
+		VerticesOnGeomEdge= newVerticesOnGeomEdge;
+		if(VertexOnBThEdge) delete []  VertexOnBThEdge;
+		VertexOnBThEdge = newVertexOnBThEdge;
+		NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
+		NbVertexOnBThEdge=newNbVertexOnBThEdge;
+		//  ReMakeTriangleContainingTheVertex();
+
+		ReconstructExistingMesh();
+
+		if (verbose>2){
+			printf("   number of quadrilaterals    = %i\n",NbOfQuad);
+			printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
+			printf("   number of outside triangles = %i\n",NbOutT);
+		}
+
+		return 0; //ok
+
+Error:
+		nbv = nbvold;
+		nbt = nbtold;
+		NbOutT = NbOutTold;
+		// cleaning memory ---
+		delete newedges;
+		if (ksplitarray) delete [] ksplitarray;
+		if (kedge) delete [] kedge;
+		if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
+		if (edge4) delete edge4;
+		if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
+
+		return ret; // ok 
+	}
+	/*}}}1*/
+/*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{1*/
+long  Mesh::SplitInternalEdgeWithBorderVertices(){
+	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/
+
+	long NbSplitEdge=0;
+	SetVertexFieldOn();  
+	long it;
+	long nbvold=nbv;
+	long int verbose=2;
+	for (it=0;it<nbt;it++){
+		Triangle &t=triangles[it];
+		if (t.link)
+		 for (int j=0;j<3;j++)
+		  if(!t.Locked(j) && !t.Hidden(j)){
+			  Triangle &tt = *t.TriangleAdj(j);
+			  if ( &tt && tt.link && it < Number(tt)) 
+				 { // an internal edge 
+				  BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
+				  BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
+				  if (v0.onGeometry && v1.onGeometry){
+					  R2 P= ((R2) v0 + (R2) v1)*0.5;
+					  if ( nbv<maxnbv) {
+						  vertices[nbv].r = P;
+						  vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
+						  vertices[nbv].ReferenceNumber=0;
+						  vertices[nbv].DirOfSearch = NoDirOfSearch ;
+					  }
+					  NbSplitEdge++;
+				  }
+				 }
+		  }
+	}
+	ReMakeTriangleContainingTheVertex();    
+	if (nbvold!=nbv){
+		long  iv = nbvold;
+		long NbSwap = 0;
+		Icoor2 dete[3];  
+		for (int i=nbvold;i<nbv;i++) {// for all the new point
+			BamgVertex & vi = vertices[i];
+			vi.i = toI2(vi.r);
+			vi.r = toR2(vi.i);
+
+			// a good new point 
+			vi.ReferenceNumber=0; 
+			vi.DirOfSearch =NoDirOfSearch;
+			Triangle *tcvi = FindTriangleContaining(vi.i,dete);
+			if (tcvi && !tcvi->link) {
+				printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
+			}
+
+			quadtree->Add(vi);
+			if (!tcvi || tcvi->det<0){// internal
+				ISSMERROR("!tcvi || tcvi->det < 0");
+			}
+			AddVertex(vi,tcvi,dete);
+			NbSwap += vi.Optim(1);          
+			iv++;
+			//      }
+	}
+	if (verbose>3) {
+		printf("   number of points: %i\n",iv);
+		printf("   number of swap to  split internal edges with border vertices: %i\n",NbSwap);
+		nbv = iv;
+	}
+}
+if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices  to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold));
+if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge);
+
+return  NbSplitEdge;
+}
+/*}}}1*/
+	/*FUNCTION Mesh::TriangulateFromGeom0{{{1*/
+	void Mesh::TriangulateFromGeom0(long imaxnbv,BamgOpts* bamgopts){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)*/
 
@@ -2737,5 +4775,5 @@
 		R2 AB;
 		GeometricalVertex *a,*b;
-		MeshVertex *va,*vb;
+		BamgVertex *va,*vb;
 		GeometricalEdge *e;
 
@@ -2745,5 +4783,5 @@
 		//initialize this
 		if (verbose>3) printf("      Generating Boundary vertices\n");
-		PreInit(inbvx);
+		Init(imaxnbv);
 		nbv=0;
 		NbVerticesOnGeomVertex=0;
@@ -2760,6 +4798,6 @@
 		//initialize VerticesOnGeomVertex
 		VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];  
-		if( NbVerticesOnGeomVertex >= nbvx) {
-			ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);
+		if( NbVerticesOnGeomVertex >= maxnbv) {
+			ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);
 		}
 
@@ -2773,10 +4811,10 @@
 			if (Gh[i].Required() && Gh[i].IsThe()) {//Gh  vertices Required
 
-				//Add the vertex (provided that nbv<nbvx)
-				if (nbv<nbvx){
+				//Add the vertex (provided that nbv<maxnbv)
+				if (nbv<maxnbv){
 					vertices[nbv]=Gh[i];
 				}
 				else{
-					ISSMERROR("Maximum number of vertices (nbvx = %i) too small",nbvx);
+					ISSMERROR("Maximum number of vertices (maxnbv = %i) too small",maxnbv);
 				}
 				
@@ -2860,5 +4898,5 @@
 								NbNewPoints=0;
 								NbEdgeCurve=0;
-								if (nbvend>=nbvx) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase nbvx");
+								if (nbvend>=maxnbv) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase maxnbv");
 								lcurve =0;
 								s = lstep;
@@ -3033,6 +5071,6 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Mesh::GeomToTriangles1{{{1*/
-	void Mesh::GeomToTriangles1(long inbvx,BamgOpts* bamgopts,int KeepVertices){ 
+	/*FUNCTION Mesh::TriangulateFromGeom1{{{1*/
+	void Mesh::TriangulateFromGeom1(long imaxnbv,BamgOpts* bamgopts,int KeepVertices){ 
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles1)*/
 
@@ -3064,5 +5102,5 @@
 
 		//Initialize new mesh
-		this->PreInit(inbvx);
+		this->Init(imaxnbv);
 		BTh.SetVertexFieldOn();
 		int* bcurve = new int[Gh.NbOfCurves]; // 
@@ -3082,10 +5120,10 @@
 		int i; 
 		for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
-		if(NbVerticesOnGeomVertex >= nbvx) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);}
+		if(NbVerticesOnGeomVertex >= maxnbv) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);}
 
 		VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
 		VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
 
-		//At this point there is NO vertex but vertices should have been allocated by PreInit
+		//At this point there is NO vertex but vertices should have been allocated by Init
 		ISSMASSERT(vertices);
 		for (i=0;i<Gh.nbv;i++){
@@ -3103,5 +5141,5 @@
 			if (vog.IsRequiredVertex()){
 				GeometricalVertex* gv=vog;
-				MeshVertex *bv = vog;
+				BamgVertex *bv = vog;
 				ISSMASSERT(gv->to); // use of Geom -> Th
 				VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->to,bv);
@@ -3208,7 +5246,7 @@
 						long i=0;// index of new points on the curve
 						register GeometricalVertex * GA0 = *(*peequi)[k0equi].onGeometry;
-						MeshVertex *A0;
+						BamgVertex *A0;
 						A0 = GA0->to;  // the vertex in new mesh
-						MeshVertex *A1;
+						BamgVertex *A1;
 						VertexOnGeom *GA1;
 						Edge* PreviousNewEdge = 0;
@@ -3228,5 +5266,5 @@
 								ISSMASSERT(pe && ee.onGeometry);
 								ee.onGeometry->SetMark();
-								MeshVertex & v0=ee[0], & v1=ee[1];
+								BamgVertex & v0=ee[0], & v1=ee[1];
 								R2 AB=(R2)v1-(R2)v0;
 								double L0=L,LAB;
@@ -3241,5 +5279,5 @@
 										ISSMASSERT(sNew>=L0);
 										ISSMASSERT(LAB);
-										ISSMASSERT(vertices && nbv<nbvx);
+										ISSMASSERT(vertices && nbv<maxnbv);
 										ISSMASSERT(edges && nbe<nbex);
 										ISSMASSERT(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
@@ -3328,6 +5366,6 @@
 			//Allocate memory
 			if(step==0){
-				if(nbv+NbOfNewPoints > nbvx) {
-					ISSMERROR("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx);
+				if(nbv+NbOfNewPoints > maxnbv) {
+					ISSMERROR("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,maxnbv);
 				}
 				edges = new Edge[NbOfNewEdge];
@@ -3365,2053 +5403,4 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Mesh::Insert{{{1*/
-	void Mesh::Insert() {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
-
-		/*Insert points in the existing Geometry*/
-
-		//Intermediary
-		int i;
-
-		/*Get options*/
-		long int verbose=2;
-
-		//Display info
-		if (verbose>2) printf("   Insert initial %i vertices\n",nbv);
-
-		//Compute integer coordinates and determinants for the existing vertices (from Geometry)
-		SetIntCoor();
-
-		/*Now we want to build a list (ordre) of the vertices in a random
-		 * order. To do so, we use the following method:
-		 *
-		 * From an initial k0 in [0 nbv[ random (vertex number)
-		 * the next k (vertex number) is computed using a big
-		 * prime number (PN>>nbv) following:
-		 *
-		 * k_{i+1} = k_i + PN  [nbv]
-		 *
-		 * let's show that:
-		 *
-		 *   for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j
-		 *
-		 * Let's assume that there are 2 distinct j1 and j2 such that
-		 * k_j1 = k_j2
-		 *
-		 * This means that
-		 *  
-		 *  k0+j1*PN = k0+j2*PN [nbv]
-		 *  (j1-j2)*PN =0       [nbv]
-		 * since PN is a prime number larger than nbv, and nbv!=1
-		 *  j1-j2=0             [nbv]
-		 * BUT
-		 *  j1 and j2 are in [0 nbv[ which is impossible.
-		 *
-		 *  We hence have built a random list of nbv elements of
-		 *  [0 nbv[ all distincts*/
-		for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ;
-		const long PrimeNumber= BigPrimeNumber(nbv) ;
-		int   k0=rand()%nbv; 
-		for (int is3=0; is3<nbv; is3++){
-			ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv];
-		}
-
-		/*Modify ordre such that the first 3 vertices form a triangle*/
-
-		//get first vertex i such that [0,1,i] are not aligned
-		for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
-			//if i is higher than nbv, it means that all the determinants are 0,
-			//all vertices are aligned!
-			if  ( ++i >= nbv) {
-				ISSMERROR("all the vertices are aligned");
-			}
-		}
-		// exchange i et 2 in "ordre" so that
-		// the first 3 vertices are not aligned (real triangle)
-		Exchange(ordre[2], ordre[i]);
-
-		/*Take the first edge formed by the first two vertices and build
-		 * the initial simple mesh from this edge and 2 boundary triangles*/
-
-		MeshVertex *  v0=ordre[0], *v1=ordre[1];
-
-		nbt = 2;
-		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;
-
-		//Build adjacence
-		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 triangle: det = -1
-		triangles[1].det = -1;  //boundary triangle: det = -1
-
-		triangles[0].SetTriangleContainingTheVertex();
-		triangles[1].SetTriangleContainingTheVertex();
-
-		triangles[0].link=&triangles[1];
-		triangles[1].link=&triangles[0];
-
-		//build quadtree
-		if (!quadtree)  quadtree = new QuadTree(this,0);
-		quadtree->Add(*v0);
-		quadtree->Add(*v1);
-
-		/*Now, add the vertices One by One*/
-
-		long NbSwap=0;
-		if (verbose>3) printf("   Begining of insertion process...\n");
-
-		for (int icount=2; icount<nbv; icount++) {
-
-			//Get new vertex
-			MeshVertex *newvertex=ordre[icount];
-
-			//Find the triangle in which newvertex is located
-			Icoor2 dete[3];
-			Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
-
-			//Add newvertex to the quadtree
-			quadtree->Add(*newvertex); 
-
-			//Add newvertex to the existing mesh
-			AddVertex(*newvertex,tcvi,dete);
-
-			//Make the mesh Delaunay around newvertex by swaping the edges
-			NbSwap += newvertex->Optim(1,0);
-		}
-
-		//Display info
-		if (verbose>3) {
-			printf("      NbSwap of insertion: %i\n",NbSwap);
-			printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
-		}
-
-#ifdef NBLOOPOPTIM
-
-		k0 = rand()%nbv ; 
-		for (int is4=0; is4<nbv; is4++) 
-		 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];
-
-		for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){
-			long  NbSwap = 0;
-			for (int is1=0; is1<nbv; is1++) 
-			 NbSwap += ordre[is1]->Optim(0,0);
-			if (verbose>3) {
-				printf("      Optim Loop: %i\n",Nbloop);
-				printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
-			}
-			if(!NbSwap) break;
-		}
-		ReMakeTriangleContainingTheVertex(); 
-		// because we break the TriangleContainingTheVertex
-#endif
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::InsertNewPoints{{{1*/
-	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
-
-		long int verbose=0;
-		double seuil= 1.414/2 ;// for two close point 
-		long i;
-		long NbSwap=0;
-		Icoor2 dete[3];
-
-		//number of new points
-		const long nbvnew=nbv-nbvold;
-
-		//display info if required
-		if (verbose>5) printf("      Try to Insert %i new points\n",nbvnew);
-
-		//return if no new points
-		if (!nbvnew) return 0; 
-
-		/*construction of a random order*/
-		const long PrimeNumber= BigPrimeNumber(nbv)  ;
-		//remainder of the division of rand() by nbvnew
-		long k3 = rand()%nbvnew;
-		//loop over the new points
-		for (int is3=0; is3<nbvnew; is3++){
-			register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
-			register long i=nbvold+is3; 
-			ordre[i]= vertices + j;
-			ordre[i]->ReferenceNumber=i;
-		}
-
-		// for all the new point
-		long iv=nbvold;
-		for (i=nbvold;i<nbv;i++){
-			MeshVertex &vi=*ordre[i];
-			vi.i=toI2(vi.r);
-			vi.r=toR2(vi.i);
-			double hx,hy;
-			vi.m.Box(hx,hy);
-			Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
-			if (!quadtree->ToClose(vi,seuil,hi,hj)){
-				// a good new point 
-				MeshVertex &vj = vertices[iv];
-				long  j=vj.ReferenceNumber; 
-				if (&vj!=ordre[j]){
-					ISSMERROR("&vj!= ordre[j]");
-				}
-				if(i!=j){ 
-					Exchange(vi,vj);
-					Exchange(ordre[j],ordre[i]);
-				}
-				vj.ReferenceNumber=0; 
-				Triangle *tcvj=FindTriangleContaining(vj.i,dete);
-				if (tcvj && !tcvj->link){
-					tcvj->Echo();
-					ISSMERROR("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link);
-				}
-				quadtree->Add(vj);
-				AddVertex(vj,tcvj,dete);
-				NbSwap += vj.Optim(1);          
-				iv++;
-			}
-		} 
-		if (verbose>3) {
-			printf("         number of new points: %i\n",iv);
-			printf("         number of to close (?) points: %i\n",nbv-iv);
-			printf("         number of swap after: %i\n",NbSwap);
-		}
-		nbv = iv;
-
-		for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);  
-		if (verbose>3) printf("   NbSwap=%i\n",NbSwap);
-
-		NbTSwap +=  NbSwap ;
-		return nbv-nbvold;
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::MakeGeometricalEdgeToEdge{{{1*/
-	Edge** Mesh::MakeGeometricalEdgeToEdge() {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeometricalEdgeToEdge)*/
-
-		if (!Gh.nbe){
-			ISSMERROR("!Gh.nbe");
-		}
-		Edge **e= new (Edge* [Gh.nbe]);
-
-		long i;
-		for ( i=0;i<Gh.nbe ; i++)
-		 e[i]=NULL;
-		for ( i=0;i<nbe ; i++) 
-		  { 
-			Edge * ei = edges+i;
-			GeometricalEdge *onGeometry = ei->onGeometry; 
-			e[Gh.Number(onGeometry)] = ei;    
-		  }
-		for ( i=0;i<nbe ; i++) 
-		 for (int ii=0;ii<2;ii++) { 
-			 Edge * ei = edges+i;
-			 GeometricalEdge *onGeometry = ei->onGeometry;
-			 int j= ii;
-			 while (!(*onGeometry)[j].Required()) { 
-				 Adj(onGeometry,j); // next geom edge
-				 j=1-j;
-				 if (e[Gh.Number(onGeometry)])  break; // optimisation
-				 e[Gh.Number(onGeometry)] = ei; 
-			 }
-		 }
-
-		int kk=0;
-		for ( i=0;i<Gh.nbe ; i++){
-			if (!e[i]){
-				kk++;
-				if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);
-			}
-		}
-		if(kk) ISSMERROR("See above");
-
-		return e;
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::MakeQuadrangles{{{1*/
-	void Mesh::MakeQuadrangles(double costheta){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/
-
-		long int verbose=0;
-
-		if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta);
-
-		if (costheta >1) {
-			if (verbose>5) printf("   do nothing: costheta > 1\n");
-		}
-
-
-			long nbqq = (nbt*3)/2;
-			DoubleAndInt *qq = new DoubleAndInt[nbqq];
-
-			long i,ij;
-			int j;
-			long k=0;
-			for (i=0;i<nbt;i++)
-			 for (j=0;j<3;j++)
-			  if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)
-				qq[k++].i3j=i*3+j;
-			//  sort  qq
-			HeapSort(qq,k);
-
-			long kk=0;
-			for (ij=0;ij<k;ij++) { 
-				i=qq[ij].i3j/3;
-				j=(int) (qq[ij].i3j%3);
-				// optisamition no float computation  
-				if (triangles[i].QualityQuad(j,0) >=costheta) 
-				 triangles[i].SetHidden(j),kk++;
-			  }
-			NbOfQuad = kk;
-			if (verbose>2){
-				printf("   number of quadrilaterals    = %i\n",NbOfQuad);
-				printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
-				printf("   number of outside triangles = %i\n",NbOutT);
-			}
-			delete [] qq;
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::MakeQuadTree{{{1*/
-	void Mesh::MakeQuadTree() {  
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/
-
-		long int verbose=0;
-		if (  !quadtree )  quadtree = new QuadTree(this);
-
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::MaxSubDivision{{{1*/
-	void  Mesh::MaxSubDivision(double maxsubdiv) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
-
-		long int verbose=0;
-
-		const  double maxsubdiv2 = maxsubdiv*maxsubdiv;
-		if(verbose>1) printf("   Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);
-		// for all the edges 
-		// if the len of the edge is to long 
-		long it,nbchange=0;    
-		double lmax=0;
-		for (it=0;it<nbt;it++){
-			Triangle &t=triangles[it];
-			for (int j=0;j<3;j++){
-				Triangle &tt = *t.TriangleAdj(j);
-				if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)){
-					MeshVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
-					MeshVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
-					R2 AB= (R2) v1-(R2) v0;
-					Metric M = v0;
-					double l = M(AB,AB);
-					lmax = Max(lmax,l);
-					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;
-						D2xD2 Rt1(Rt.inv());
-						D2xD2 D(maxsubdiv2,0,0,lc);
-						D2xD2 MM = Rt1*D*Rt1.t();
-						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
-						double lc = M(AC,AC);
-						D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
-						D2xD2 Rt1(Rt.inv());
-						D2xD2 D(maxsubdiv2,0,0,lc);
-						D2xD2  MM = Rt1*D*Rt1.t();
-						v1.m =  M = Metric(MM.x.x,MM.y.x,MM.y.y);
-						nbchange++;
-					}
-				}
-			}
-		}
-		if(verbose>3){
-			printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
-		}
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::MetricAt{{{1*/
-	Metric Mesh::MetricAt(const R2 & A) const { 
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/
-
-		I2 a = toI2(A);
-		Icoor2 deta[3];
-		Triangle * t =FindTriangleContaining(a,deta);
-		if (t->det <0) { // outside
-			double ba,bb;
-			TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;
-			return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
-		else { // inside
-			double   aa[3];
-			double s = deta[0]+deta[1]+deta[2];
-			aa[0]=deta[0]/s;
-			aa[1]=deta[1]/s;
-			aa[2]=deta[2]/s;
-			return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
-		}
-	}
-	/*}}}1*/
-/*FUNCTION Mesh::NearestVertex{{{1*/
-MeshVertex* Mesh::NearestVertex(Icoor1 i,Icoor1 j) {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/
-	return  quadtree->NearestVertex(i,j); 
-} 
-/*}}}1*/
-	/*FUNCTION Mesh::NewPoints{{{1*/
-	void  Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
-
-		int i,j,k;
-		long NbTSwap=0;
-		long nbtold=nbt;
-		long nbvold=nbv;
-		long Headt=0;
-		long next_t;
-		long* first_np_or_next_t=new long[nbtx];
-		Triangle* t=NULL;
-
-		/*Recover options*/
-		int verbose=bamgopts->verbose;
-
-		/*First, insert old points if requested*/
-		if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){
-			if (verbose>5) printf("         Inserting initial mesh points\n");
-			for (i=0;i<Bh.nbv;i++){ 
-				MeshVertex &bv=Bh[i];
-				if (!bv.onGeometry){
-					vertices[nbv].r   = bv.r;
-					vertices[nbv++].m = bv.m;
-				}
-			}
-			Bh.ReMakeTriangleContainingTheVertex();     
-			InsertNewPoints(nbvold,NbTSwap)   ;            
-		}  
-		else Bh.ReMakeTriangleContainingTheVertex();     
-
-		// generation of the list of next Triangle 
-		for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
-		// the next traingle of i is -first_np_or_next_t[i]
-
-		// Big loop (most time consuming)
-		int iter=0;
-		if (verbose>5) printf("         Big loop\n");
-		do {
-			/*Update variables*/
-			iter++;
-			nbtold=nbt;
-			nbvold=nbv;
-
-			/*We test all triangles*/
-			i=Headt;
-			next_t=-first_np_or_next_t[i];
-			for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
-
-				//check i
-				if (i<0 || i>=nbt ){
-					ISSMERROR("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1);
-				}
-				//change first_np_or_next_t[i]
-				first_np_or_next_t[i] = iter; 
-
-				//Loop over the edges of t
-				for(j=0;j<3;j++){
-					TriangleAdjacent tj(t,j);
-					MeshVertex &vA = *tj.EdgeVertex(0);
-					MeshVertex &vB = *tj.EdgeVertex(1);
-
-					//if t is a boundary triangle, or tj locked, continue
-					if (!t->link)     continue;
-					if (t->det <0)    continue;
-					if (t->Locked(j)) continue;
-
-					TriangleAdjacent tadjj = t->Adj(j);	  
-					Triangle* ta=tadjj;
-
-					//if the adjacent triangle is a boundary triangle, continur
-					if (ta->det<0) continue;	  
-
-					R2 A=vA;
-					R2 B=vB;
-					k=Number(ta);
-
-					//if this edge has already been done, go to next edge of triangle
-					if(first_np_or_next_t[k]==iter) continue;
-
-					lIntTria.SplitEdge(Bh,A,B);
-					lIntTria.NewPoints(vertices,nbv,nbvx);
-				  } // end loop for each edge 
-			  }// for triangle   
-
-			if (!InsertNewPoints(nbvold,NbTSwap)) break;
-			for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
-			Headt = nbt; // empty list 
-
-			// for all the triangle containing the vertex i
-			for (i=nbvold;i<nbv;i++){ 
-				MeshVertex*          s  = vertices + i;
-				TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
-				Triangle*        tbegin= (Triangle*) ta;
-				long kt;
-				do { 
-					kt = Number((Triangle*) ta);
-					if (first_np_or_next_t[kt]>0){
-						first_np_or_next_t[kt]=-Headt;
-						Headt=kt;
-					}
-					if (ta.EdgeVertex(0)!=s){
-						ISSMERROR("ta.EdgeVertex(0)!=s");
-					}
-					ta = Next(Adj(ta));
-				} while ( (tbegin != (Triangle*) ta)); 
-			}
-
-		} while (nbv!=nbvold);
-
-		delete []  first_np_or_next_t;
-
-		long NbSwapf =0,NbSwp;
-
-		NbSwp = NbSwapf;
-		for (i=0;i<nbv;i++)
-		 NbSwapf += vertices[i].Optim(0);
-		NbTSwap +=  NbSwapf ;
-	}
-	/*}}}1*/
-/*FUNCTION Mesh::PreInit{{{1*/
-void Mesh::PreInit(long inbvx) {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/
-
-	long int verbose=0;
-
-	srand(19999999);
-	NbRef=0;
-	//  allocGeometry=0;
-	NbOfTriangleSearchFind =0;
-	NbOfSwapTriangle =0;
-	nbiv=0;
-	nbv=0;
-	nbvx=inbvx;
-	nbt=0;
-	NbOfQuad = 0;
-	nbtx=2*inbvx-2;
-	NbSubDomains=0;
-	NbVertexOnBThVertex=0;
-	NbVertexOnBThEdge=0;
-	VertexOnBThVertex=NULL;
-	VertexOnBThEdge=NULL;
-
-	NbCrackedVertices=0;
-	CrackedVertices  =NULL;  
-	NbCrackedEdges   =0;
-	CrackedEdges     =NULL;  
-	nbe = 0; 
-
-	if (inbvx) {
-		vertices=new MeshVertex[nbvx];
-		if (!vertices){
-			ISSMERROR("!vertices");
-		}
-		ordre=new (MeshVertex* [nbvx]);
-		if (!ordre){
-			ISSMERROR("!ordre");
-		}
-		triangles=new Triangle[nbtx];
-		if (!triangles){
-			ISSMERROR("!triangles");
-		}
-	}
-	else {
-		vertices=NULL;
-		ordre=NULL;
-		triangles=NULL;
-		nbtx=0;
-	} 
-
-	quadtree=NULL;
-	edges=NULL;
-	VerticesOnGeomVertex=NULL;
-	VerticesOnGeomEdge=NULL;
-	NbVerticesOnGeomVertex=0;
-	NbVerticesOnGeomEdge=0;
-	subdomains=NULL;
-	NbSubDomains=0;
-}
-/*}}}1*/
-	/*FUNCTION Mesh::ProjectOnCurve{{{1*/
-	GeometricalEdge*   Mesh::ProjectOnCurve( Edge & BhAB, MeshVertex &  vA, MeshVertex & vB,
-				double theta,MeshVertex & R,VertexOnEdge &  BR,VertexOnGeom & GR) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/
-
-		void *pA=0,*pB=0;
-		double tA=0,tB=0;
-		R2 A=vA,B=vB;
-		MeshVertex * pvA=&vA, * pvB=&vB;
-		if (vA.vint == IsVertexOnVertex){
-			pA=vA.onBackgroundVertex;
-		}
-		else if (vA.vint == IsVertexOnEdge){
-			pA=vA.onBackgroundEdge->be;
-			tA=vA.onBackgroundEdge->abcisse;
-		}
-		else {
-			ISSMERROR("ProjectOnCurve On MeshVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA));
-		} 
-
-		if (vB.vint == IsVertexOnVertex){
-			pB=vB.onBackgroundVertex;
-		}
-		else if(vB.vint == IsVertexOnEdge){
-			pB=vB.onBackgroundEdge->be;
-			tB=vB.onBackgroundEdge->abcisse;
-		}
-		else {
-			ISSMERROR("ProjectOnCurve On MeshVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB));
-		} 
-		Edge * e = &BhAB;
-		if (!pA || !pB || !e){
-			ISSMERROR("!pA || !pB || !e");
-		}
-		// be carefull the back ground edge e is on same geom edge 
-		// of the initiale edge def by the 2 vertex A B;
-		//check Is a background Mesh;   
-		if (e<BTh.edges || e>=BTh.edges+BTh.nbe){
-			ISSMERROR("e<BTh.edges || e>=BTh.edges+BTh.nbe");
-		}
-		// walk on BTh edge 
-		//not finish ProjectOnCurve with BackGround Mesh);
-		// 1 first find a back ground edge contening the vertex A
-		// 2 walk n back gound boundary to find the final vertex B
-
-		if( vA.vint == IsVertexOnEdge) 
-		  { // find the start edge 
-			e = vA.onBackgroundEdge->be;	 
-
-		  } 
-		else if (vB.vint == IsVertexOnEdge) 
-		  {
-			theta = 1-theta;
-			Exchange(tA,tB);
-			Exchange(pA,pB);
-			Exchange(pvA,pvB);
-			Exchange(A,B);
-			e =  vB.onBackgroundEdge->be;
-
-		  } 
-		else{ // do the search by walking 
-			ISSMERROR("case not supported yet");
-		  }
-
-		// find the direction of walking with sens of edge and pA,PB;
-		R2 AB=B-A;
-
-		double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
-		int kkk=0;
-		int sens = (cosE01AB>0) ? 1 : 0;
-
-		//   double l=0; // length of the edge AB
-		double abscisse = -1;
-
-		for (int cas=0;cas<2;cas++){
-			// 2 times algo:
-			//    1 for computing the length l
-			//    2 for find the vertex 
-			int  iii;
-			MeshVertex  *v0=pvA,*v1; 
-			Edge *neee,*eee;
-			double lg =0; // length of the curve 
-			double te0;
-			// we suppose take the curve's abcisse 
-			for ( eee=e,iii=sens,te0=tA;
-						eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
-						neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { 
-
-				kkk=kkk+1;
-				if (kkk>=100){
-					ISSMERROR("kkk>=100");
-				}
-				if (!eee){
-					ISSMERROR("!eee");
-				}
-				double lg0 = lg;
-				double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
-				lg += dp;
-				if (cas && abscisse <= lg) { // ok we find the geom edge 
-					double sss  =   (abscisse-lg0)/dp;
-					double thetab = te0*(1-sss)+ sss*iii;
-					if (thetab<0 || thetab>1){
-						ISSMERROR("thetab<0 || thetab>1");
-					}
-					BR = VertexOnEdge(&R,eee,thetab);
-					return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
-				  }
-			  }
-			// we find the end 
-			if (v1 != pvB){
-				if (( void*) v1 == pB)
-				 tB = iii;
-
-				double lg0 = lg;
-				if (!eee){
-					ISSMERROR("!eee");
-				}
-				v1 = pvB;
-				double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
-				lg += dp;	
-				abscisse = lg*theta;
-				if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
-				  { // ok we find the geom edge 
-					double sss  =   (abscisse-lg0)/dp;
-					double thetab = te0*(1-sss)+ sss*tB;
-					if (thetab<0 || thetab>1){
-						ISSMERROR("thetab<0 || thetab>1");
-					}
-					BR = VertexOnEdge(&R,eee,thetab);
-					return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
-				  }
-			  }
-			abscisse = lg*theta;
-
-		  }
-		ISSMERROR("Big bug...");
-		return 0; // just for the compiler 
-	}                  
-	/*}}}1*/
-/*FUNCTION Mesh::ReconstructExistingMesh{{{1*/
-void Mesh::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 verbose=0;
-
-	// generation of the integer coordinate
-
-	// find extrema coordinates of vertices pmin,pmax
-	long i;
-	if(verbose>2) printf("      Reconstruct mesh of %i vertices\n",nbv); 
-
-	//initialize ordre
-	ISSMASSERT(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){ 
-		ISSMERROR("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) {
-				ISSMASSERT(!triangles[i].TriangleAdj(j));
-				ISSMASSERT(!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 {
-				ISSMERROR("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(verbose>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) {
-		ISSMERROR("%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) {
-		 ISSMERROR("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*/
-	MeshVertex *  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++) {
-		MeshVertex *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
-			MeshVertex &a=edges[i][0], &b=edges[i][1];
-			if (a.t && b.t){
-				knbe++;
-				if (ForceEdge(a,b,ta)<0) nbloss++;
-			}
-		}
-	}
-	if(nbloss) {
-		ISSMERROR("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;
-					MeshVertex *v0= ta.EdgeVertex(0);
-					MeshVertex *v1= ta.EdgeVertex(1);
-					long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
-
-					ISSMASSERT(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((MeshVertex *) NULL,(MeshVertex *) NULL,(MeshVertex *) NULL);
-			triangles[i].color=-1;
-		}
-		else{
-			triangles[i].color= savenbt+ NbTfillHoll++;
-		}
-	}
-	ISSMASSERT(savenbt+NbTfillHoll<=savenbtx);
-
-	// copy of the outside triangles in saveMesh 
-	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) {
-		ISSMERROR("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 MeshVertex number %i\n",Number(edges[i][j])+1,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",Number(edges[i][j])+1,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");
-						ISSMERROR("See above (might be cryptic...)");
-					}
-				}
-			}
-		}
-	}
-}
-/*}}}1*/
-	/*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/
-	void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
-
-		long int verbose=0;
-		long *renu= new long[nbt];
-		register Triangle *t0,*t,*te=triangles+nbt;
-		register long k=0,it,i,j;
-
-		for ( it=0;it<nbt;it++) 
-		 renu[it]=-1; // outside triangle 
-		for ( i=0;i<NbSubDomains;i++)
-		  { 
-			t=t0=subdomains[i].head;
-			if (!t0){ // not empty sub domain
-				ISSMERROR("!t0");
-			}
-			do { 
-				long kt = Number(t);
-				if (kt<0 || kt >= nbt ){
-					ISSMERROR("kt<0 || kt >= nbt");
-				}
-				if (renu[kt]!=-1){
-					ISSMERROR("renu[kt]!=-1");
-				}
-				renu[kt]=k++;
-			}
-			while (t0 != (t=t->link));
-		  }
-		// take is same numbering if possible    
-		if(justcompress)
-		 for ( k=0,it=0;it<nbt;it++) 
-		  if(renu[it] >=0 ) 
-			renu[it]=k++;
-
-		// put the outside triangles at the end
-		for ( it=0;it<nbt;it++){
-			if (renu[it]==-1) renu[it]=k++;
-		}
-		if (k != nbt){
-			ISSMERROR("k != nbt");
-		}
-		// do the change on all the pointeur 
-		for ( it=0;it<nbt;it++)
-		 triangles[it].ReNumbering(triangles,te,renu);
-
-		for ( i=0;i<NbSubDomains;i++)
-		 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
-
-		// move the Triangles  without a copy of the array 
-		// be carefull not trivial code 
-		for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
-		 if (renu[it] >= 0) // a new sub cycle
-			{ 
-			 i=it;
-			 Triangle ti=triangles[i],tj;
-			 while ( (j=renu[i]) >= 0) 
-				{ // i is old, and j is new 
-				 renu[i] = -1; // mark 
-				 tj = triangles[j]; // save new
-				 triangles[j]= ti; // new <- old
-				 i=j;     // next 
-				 ti = tj;
-				}  
-			}
-		delete [] renu;
-		nt = nbt - NbOutT;
-
-	}
-	/*}}}1*/
-	/*FUNCTION Mesh::ReNumberingVertex{{{1*/
-	void Mesh::ReNumberingVertex(long * renu) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
-
-		// warning be carfull because pointer
-		// from on mesh to over mesh 
-		//  --  so do ReNumbering at the beginning
-		MeshVertex * ve = vertices+nbv;
-		long it,ie,i;
-
-		printf("renumbering triangles\n");
-		for ( it=0;it<nbt;it++) 
-		 triangles[it].ReNumbering(vertices,ve,renu);
-
-		printf("renumbering edges\n");
-		for ( ie=0;ie<nbe;ie++) 
-		 edges[ie].ReNumbering(vertices,ve,renu);
-
-		printf("renumbering vertices on geom\n");
-		for (i=0;i< NbVerticesOnGeomVertex;i++)
-		  {
-			MeshVertex *v = VerticesOnGeomVertex[i].mv;
-			if (v>=vertices && v < ve)
-			 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];
-		  }
-
-		printf("renumbering vertices on edge\n");
-		for (i=0;i< NbVerticesOnGeomEdge;i++)
-		  {
-			MeshVertex *v =VerticesOnGeomEdge[i].mv;
-			if (v>=vertices && v < ve)
-			 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];
-		  }
-
-		printf("renumbering vertices on Bth vertex\n");
-		for (i=0;i< NbVertexOnBThVertex;i++)
-		  {
-			MeshVertex *v=VertexOnBThVertex[i].v;
-			if (v>=vertices && v < ve)
-			 VertexOnBThVertex[i].v=vertices+renu[Number(v)];
-		  }
-
-		for (i=0;i< NbVertexOnBThEdge;i++)
-		  {
-			MeshVertex *v=VertexOnBThEdge[i].v;
-			if (v>=vertices && v < ve)
-			 VertexOnBThEdge[i].v=vertices+renu[Number(v)];
-		  }
-
-		// move the Vertices without a copy of the array 
-		// be carefull not trivial code 
-		long j;
-		for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
-		 if (renu[it] >= 0) // a new sub cycle
-			{ 
-			 i=it;
-			 MeshVertex ti=vertices[i],tj;
-			 while ( (j=renu[i]) >= 0){
-				 // i is old, and j is new 
-				 renu[i] = -1-renu[i]; // mark 
-				 tj = vertices[j];     // save new
-				 vertices[j]= ti;      // new <- old
-				 i=j;     // next 
-				 ti = tj;
-				}  
-			}
-		if (quadtree){
-			delete quadtree;
-			quadtree = new QuadTree(this);
-		}
-
-		for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
-	}
-	/*}}}1*/
-/*FUNCTION Mesh::SetIntCoor{{{1*/
-void Mesh::SetIntCoor(const char * strfrom) {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/
-
-	/*Set integer coordinate for existing vertices*/
-
-	//Get extrema coordinates of the existing vertices
-	pmin =  vertices[0].r;
-	pmax =  vertices[0].r;
-	long i;
-	for (i=0;i<nbv;i++) {
-		pmin.x = Min(pmin.x,vertices[i].r.x);
-		pmin.y = Min(pmin.y,vertices[i].r.y);
-		pmax.x = Max(pmax.x,vertices[i].r.x);
-		pmax.y = Max(pmax.y,vertices[i].r.y);
-	}
-	R2 DD = (pmax-pmin)*0.05;
-	pmin = pmin-DD;
-	pmax = pmax+DD; 
-
-	//Compute coefIcoor
-	coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
-	if (coefIcoor<=0){
-		ISSMERROR("coefIcoor should be positive, a problem in the geometry is likely");
-	}
-
-	// generation of integer coord  
-	for (i=0;i<nbv;i++) {
-		vertices[i].i = toI2(vertices[i].r);    
-	}
-
-	// computation of the det 
-	int number_of_errors=0;
-	for (i=0;i<nbt;i++) {
-		MeshVertex & v0 = triangles[i][0];
-		MeshVertex & v1 = triangles[i][1];
-		MeshVertex & v2 = triangles[i][2];
-
-		//If this is not a boundary triangle
-		if ( &v0 && &v1 &&  &v2 ){
-
-			/*Compute determinant*/
-			triangles[i].det= det(v0,v1,v2);
-
-			/*Check that determinant is positive*/
-			if (triangles[i].det <=0){
-
-				/*increase number_of_errors and print error only for the first 20 triangles*/
-				number_of_errors++;
-				if (number_of_errors<20){
-					printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det);
-				}
-			}
-		}
-
-		//else, set as -1
-		else triangles[i].det=-1;
-	}
-
-	if (number_of_errors) ISSMERROR("Fatal error: some triangles have negative areas, see above");
-}
-/*}}}1*/
-/*FUNCTION Mesh::ShowRegulaty{{{1*/
-void  Mesh::ShowRegulaty() const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
-
-	const  double  sqrt32=sqrt(3.)*0.5; 
-	const double  aireKh=sqrt32*0.5;
-	D2  Beq(1,0),Heq(0.5,sqrt32);
-	D2xD2 Br(D2xD2(Beq,Heq).t());
-	D2xD2 B1r(Br.inv());
-	double gammamn=1e100,hmin=1e100;
-	double gammamx=0,hmax=0;
-	double beta=1e100;
-	double beta0=0;
-	double  alpha2=0;
-	double area=0,Marea=0;
-	// double cf= double(coefIcoor);
-	// double cf2= 6.*cf*cf;
-	int nt=0;
-	for (int it=0;it<nbt;it++)
-	 if ( triangles[it].link) 
-		{
-		 nt++;
-		 Triangle &K=triangles[it];
-		 double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
-		 area+= area3;
-		 D2xD2 B_Kt(K[0],K[1],K[2]);
-		 D2xD2 B_K(B_Kt.t());
-		 D2xD2 B1K = Br*B_K.inv();
-		 D2xD2 BK =  B_K*B1r;
-		 D2xD2 B1B1 = B1K.t()*B1K;
-		 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
-		 MatVVP2x2 VMK(MK);
-		 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
-		 double betaK=0;
-
-		 for (int j=0;j<3;j++)
-			{
-			 double he= Norme2(R2(K[j],K[(j+1)%3]));
-			 hmin=Min(hmin,he);
-			 hmax=Max(hmax,he);
-			 MeshVertex & v=K[j];
-			 D2xD2 M((Metric)v);
-			 betaK += sqrt(M.det());
-
-			 D2xD2 BMB = BK.t()*M*BK;
-			 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
-			 MatVVP2x2 VM1(M1);
-			 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
-			 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);		
-			}
-		 betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
-		 Marea+= betaK;
-		 beta=min(beta,betaK);
-		 beta0=max(beta0,betaK);
-		}   
-	area*=3; 
-	gammamn=sqrt(gammamn);
-	gammamx=sqrt(gammamx);    
-	printf("   Adaptmesh info:\n");
-	printf("      number of triangles = %i\n",nt);
-	printf("      hmin = %g, hmax=%g\n",hmin,hmax);
-	printf("      area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt));
-	printf("      infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx);
-	printf("      anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5));
-}
-/*}}}1*/
-/*FUNCTION Mesh::ShowHistogram{{{1*/
-void  Mesh::ShowHistogram() const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
-
-	const long kmax=10;
-	const double llmin = 0.5,llmax=2;
-	const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
-	long histo[kmax+1];
-	long i,it,k, nbedges =0;
-	for (i=0;i<=kmax;i++) histo[i]=0;
-	for (it=0;it<nbt;it++)
-	 if ( triangles[it].link) 
-		{
-
-		 for (int j=0;j<3;j++)
-			{
-			 Triangle *ta = triangles[it].TriangleAdj(j);
-			 if ( !ta || !ta->link || Number(ta) >= it) 
-				{ 
-				 MeshVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
-				 MeshVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
-				 if ( !& vP || !&vQ) continue;
-				 R2 PQ = vQ.r - vP.r;
-				 double l = log(LengthInterpole(vP,vQ,PQ));
-				 nbedges++;
-				 k = (int) ((l - lmin)*delta);
-				 k = Min(Max(k,0L),kmax);
-				 histo[k]++;
-				}
-			}
-		}  
-	printf(" --- Histogram of the unit mesh,  nb of edges = %i\n",nbedges);
-	printf("      length of edge in   | %% of edge  | Nb of edges \n"); 
-	printf("      --------------------+-------------+-------------\n"); 
-	for   (i=0;i<=kmax;i++){ 
-		if (i==0) printf("      %10i",0);
-		else      printf("      %10g",exp(lmin+i/delta));
-		if (i==kmax) printf("          +inf   ");
-		else printf("      %10g",exp(lmin+(i+1)/delta));
-		printf("|  %10g |\n",((long)  ((10000.0 * histo[i])/ nbedges))/100.0);
-		printf("  %i\n",histo[i]);
-	}
-	printf("      --------------------+-------------+-------------\n"); 
-}
-/*}}}1*/
-/*FUNCTION Mesh::SmoothingVertex{{{1*/
-void Mesh::SmoothingVertex(int nbiter,double omega ) { 
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
-
-	long int verbose=0;
-	//  if quatree exist remove it end reconstruct
-	if (quadtree) delete quadtree;
-	quadtree=0;
-	ReMakeTriangleContainingTheVertex();
-	Triangle vide; // a triangle to mark the boundary vertex
-	Triangle   ** tstart= new Triangle* [nbv];
-	long i,j,k;
-	//   attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide 
-	if ( this == & BTh)
-	 for ( i=0;i<nbv;i++)
-	  tstart[i]=vertices[i].t;     
-	else 
-	 for ( i=0;i<nbv;i++)
-	  tstart[i]=0;
-	for ( j=0;j<NbVerticesOnGeomVertex;j++ ) 
-	 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;
-	for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 
-	 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
-	if(verbose>2) printf("   SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);
-	for (k=0;k<nbiter;k++)
-	  {
-		long i,NbSwap =0;
-		double delta =0;
-		for ( i=0;i<nbv;i++)
-		 if (tstart[i] != &vide) // not a boundary vertex 
-		  delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
-		if (!NbOfQuad)
-		 for ( i=0;i<nbv;i++)
-		  if (tstart[i] != &vide) // not a boundary vertex 
-			NbSwap += vertices[i].Optim(1);
-		if (verbose>3) printf("      move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);
-	  }
-
-	delete [] tstart;
-	if (quadtree) quadtree= new QuadTree(this);
-}
-/*}}}1*/
-/*FUNCTION Mesh::SmoothMetric{{{1*/
-void Mesh::SmoothMetric(double raisonmax) { 
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
-
-	long int verbose=0;
-
-	if(raisonmax<1.1) return;
-	if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
-	ReMakeTriangleContainingTheVertex();
-	long i,j,kch,kk,ip;
-	long *first_np_or_next_t0 = new long[nbv];
-	long *first_np_or_next_t1 = new long[nbv];
-	long Head0 =0,Head1=-1;
-	double logseuil= log(raisonmax);
-
-	for(i=0;i<nbv-1;i++)
-	 first_np_or_next_t0[i]=i+1; 
-	first_np_or_next_t0[nbv-1]=-1;// end;
-	for(i=0;i<nbv;i++)
-	 first_np_or_next_t1[i]=-1;
-	kk=0;
-	while (Head0>=0&& kk++<100) {
-		kch=0;
-		for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
-			//  pour tous les triangles autour du sommet s
-			register Triangle* t= vertices[i].t;
-			if (!t){
-				ISSMERROR("!t");
-			}
-			MeshVertex & vi = vertices[i];
-			TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
-			MeshVertex *pvj0 = ta.EdgeVertex(0);
-			while (1) {
-				ta=Previous(Adj(ta));
-				if (vertices+i != ta.EdgeVertex(1)){
-					ISSMERROR("vertices+i != ta.EdgeVertex(1)");
-				}
-				MeshVertex & vj = *(ta.EdgeVertex(0));
-				if ( &vj ) {
-					j= &vj-vertices;
-					if (j<0 || j >= nbv){
-						ISSMERROR("j<0 || j >= nbv");
-					}
-					R2 Aij = (R2) vj - (R2) vi;
-					double ll =  Norme2(Aij);
-					if (0) {  
-						double hi = ll/vi.m(Aij);
-						double hj = ll/vj.m(Aij);
-						if(hi < hj)
-						  {
-							double dh=(hj-hi)/ll;
-							if (dh>logseuil) {
-								vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
-								if(first_np_or_next_t1[j]<0)
-								 kch++,first_np_or_next_t1[j]=Head1,Head1=j;
-							}
-						  }
-					} 
-					else
-					  {
-						double li = vi.m(Aij);
-						if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
-						 if(first_np_or_next_t1[j]<0) // if the metrix change 
-						  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
-					  }
-				}
-				if  ( &vj ==  pvj0 ) break;
-			}
-		}
-		Head0 = Head1;
-		Head1 = -1;
-		Exchange(first_np_or_next_t0,first_np_or_next_t1);
-	}
-	if(verbose>2) printf("      number of iterations = %i\n",kch); 
-	delete [] first_np_or_next_t0;
-	delete [] first_np_or_next_t1;
-}
-/*}}}1*/
-	/*FUNCTION Mesh::SplitElement{{{1*/
-	int  Mesh::SplitElement(int choice){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
-
-		long int verbose=0;
-
-		Direction NoDirOfSearch;
-		const  int withBackground = &BTh != this && &BTh;
-
-		ReNumberingTheTriangleBySubDomain();
-		//int nswap =0;
-		const long nfortria( choice ? 4 : 6);
-		if(withBackground) 
-		  {
-			BTh.SetVertexFieldOn();
-			SetVertexFieldOnBTh();
-		  }
-		else
-		 BTh.SetVertexFieldOn();
-
-		long newnbt=0,newnbv=0;
-		long * kedge = 0;
-		long newNbOfQuad=NbOfQuad;
-		long * ksplit= 0, * ksplitarray=0;
-		long kkk=0;
-		int ret =0;
-		if (nbvx<nbv+nbe) return 1;//   
-		// 1) create  the new points by spliting the internal edges 
-		// set the 
-		long nbvold = nbv;
-		long nbtold = nbt;
-		long NbOutTold  = NbOutT;
-		long  NbEdgeOnGeom=0;
-		long i;
-
-		nbt = nbt - NbOutT; // remove all the  the ouside triangles 
-		long nbtsave = nbt;
-		Triangle * lastT = triangles + nbt;
-		for (i=0;i<nbe;i++)
-		 if(edges[i].onGeometry) NbEdgeOnGeom++;
-		long newnbe=nbe+nbe;
-		//  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
-		long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
-		// long newNbVertexOnBThVertex=NbVertexOnBThVertex;
-		long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
-
-		// do allocation for pointeur to the geometry and background
-		VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
-		VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
-		if (NbVerticesOnGeomEdge)
-		 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
-		if (NbVertexOnBThEdge)
-		 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
-		Edge *newedges = new Edge [newnbe];
-		//  memcpy(newedges,edges,sizeof(Edge)*nbe);
-		SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
-		long k=nbv;
-		long kk=0;
-		long kvb = NbVertexOnBThEdge;
-		long kvg = NbVerticesOnGeomEdge;
-		long ie =0;
-		Edge ** edgesGtoB=0;
-		if (withBackground)
-		 edgesGtoB= BTh.MakeGeometricalEdgeToEdge();
-		long ferr=0;
-		for (i=0;i<nbe;i++)
-		 newedges[ie].onGeometry=0;
-
-		for (i=0;i<nbe;i++)
-		  {
-			GeometricalEdge *ong =  edges[i].onGeometry;
-
-			newedges[ie]=edges[i];
-			newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
-			newedges[ie].adj[1]=newedges + ie +1;
-			R2 A = edges[i][0],B = edges[i][1];
-
-
-			kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
-			if (ong) // a geometrical edges 
-			  { 
-				if (withBackground){
-					// walk on back ground mesh 
-					//  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 
-					// a faire -- difficile 
-					// the first PB is to now a background edge between the 2 vertices
-					if (!edgesGtoB){
-						ISSMERROR("!edgesGtoB");
-					}
-					ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)],
-								edges[i][0],edges[i][1],0.5,vertices[k],
-								newVertexOnBThEdge[kvb],
-								newVerticesOnGeomEdge[kvg++]);
-					vertices[k].ReferenceNumber= edges[i].ref;
-					vertices[k].DirOfSearch =   NoDirOfSearch;        
-					;
-					// get the Info on background mesh 
-					double s =        newVertexOnBThEdge[kvb];
-					MeshVertex &  bv0  = newVertexOnBThEdge[kvb][0];
-					MeshVertex &  bv1  = newVertexOnBThEdge[kvb][1];
-					// compute the metrix of the new points 
-					vertices[k].m =  Metric(1-s,bv0,s,bv1); 
-					kvb++;
-				  }
-				else 
-				  {
-					ong=Gh.ProjectOnCurve(edges[i],
-								0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
-					// vertices[k].i = toI2( vertices[k].r);
-					vertices[k].ReferenceNumber = edges[i].ref;
-					vertices[k].DirOfSearch = NoDirOfSearch;
-					vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);	      
-				  }  
-			  }
-			else // straigth line edge ---
-			  { 
-				vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
-				vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
-				vertices[k].onGeometry = 0;
-			  }
-			//vertices[k].i = toI2( vertices[k].r);
-			R2 AB =  vertices[k].r;
-			R2 AA = (A+AB)*0.5;
-			R2 BB = (AB+B)*0.5;
-			vertices[k].ReferenceNumber = edges[i].ref;
-			vertices[k].DirOfSearch = NoDirOfSearch;
-
-			newedges[ie].onGeometry = Gh.Containing(AA,ong);
-			newedges[ie++].v[1]=vertices+k;
-
-			newedges[ie]=edges[i];
-			newedges[ie].adj[0]=newedges + ie -1;
-			newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
-			newedges[ie].onGeometry =  Gh.Containing(BB,ong);
-			newedges[ie++].v[0]=vertices+k;
-			k++;
-		  }
-		if (edgesGtoB) delete [] edgesGtoB;
-		edgesGtoB=0;
-
-		newnbv=k;
-		newNbVerticesOnGeomEdge=kvg;
-		if (newnbv> nbvx) goto Error;// bug 
-
-		nbv = k;
-
-
-		kedge = new long[3*nbt+1];
-		ksplitarray = new long[nbt+1];
-		ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
-
-		for (i=0;i<3*nbt;i++)
-		 kedge[i]=-1;
-
-		//  
-
-		for (i=0;i<nbt;i++) {
-			Triangle & t = triangles[i];
-			if (!t.link){
-				ISSMERROR("!t.link");
-			}
-			for(int j=0;j<3;j++)
-			  {
-				const TriangleAdjacent ta = t.Adj(j);
-				const Triangle & tt = ta;
-				if (&tt >= lastT)
-				 t.SetAdj2(j,0,0);// unset adj
-				const MeshVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
-				const MeshVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				long  ke =edge4->SortAndFind(Number(v0),Number(v1));
-				if (ke>0) 
-				  {
-					long ii = Number(tt);
-					int  jj = ta;
-					long ks = ke + nbvold;
-					kedge[3*i+j] = ks;
-					if (ii<nbt) // good triangle
-					 kedge[3*ii+jj] = ks;
-					MeshVertex &A=vertices[ks];
-					double aa,bb,cc,dd;
-					if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
-						// warning PB roundoff error 
-						if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0 
-										||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0  
-										||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
-							printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd);
-							ISSMERROR("Number of triangles with P2 interpolation Problem");
-						}
-					}
-					else {
-						if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0 
-										||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0 
-										||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
-							printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd);
-							ISSMERROR("Number of triangles with P2 interpolation Problem");
-						}
-					} 
-				  }
-			  }
-		}
-
-		for (i=0;i<nbt;i++){
-			ksplit[i]=1; // no split by default
-			const Triangle & t = triangles[ i];
-			int nbsplitedge =0;
-			int nbinvisible =0;
-			int invisibleedge=0;
-			int kkk[3];      
-			for (int j=0;j<3;j++)
-			  {
-				if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
-
-				const TriangleAdjacent ta = t.Adj(j);
-				const Triangle & tt = ta;
-
-
-				const MeshVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
-				const MeshVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				if ( kedge[3*i+j] < 0) 
-				  {
-					long  ke =edge4->SortAndFind(Number(v0),Number(v1));
-					if (ke<0) // new 
-					  {
-						if (&tt) // internal triangles all the boundary 
-						  { // new internal edges 
-							long ii = Number(tt);
-							int  jj = ta;
-
-							kedge[3*i+j]=k;// save the vertex number 
-							kedge[3*ii+jj]=k;
-							if (k<nbvx) 
-							  {
-								vertices[k].r = ((R2) v0+(R2) v1 )/2;
-								//vertices[k].i = toI2( vertices[k].r);
-								vertices[k].ReferenceNumber=0;
-								vertices[k].DirOfSearch =NoDirOfSearch;
-								vertices[k].m =  Metric(0.5,v0,0.5,v1);
-							  }
-							k++;
-							kkk[nbsplitedge++]=j;		      
-						  } // tt 
-						else
-						 ISSMERROR("Bug...");
-					  } // ke<0	       
-					else
-					  { // ke >=0
-						kedge[3*i+j]=nbvold+ke;
-						kkk[nbsplitedge++]=j;// previously splited
-					  }
-				  }
-				else 
-				 kkk[nbsplitedge++]=j;// previously splited
-
-			  } 
-			if (nbinvisible>=2){
-				ISSMERROR("nbinvisible>=2");
-			}
-			switch (nbsplitedge) {
-				case 0: ksplit[i]=10; newnbt++; break;   // nosplit
-				case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 
-				case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 
-				case 3:
-						  if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
-						  else   ksplit[i]=10*nfortria,newnbt+=nfortria;
-						  break;
-			} 
-			if (ksplit[i]<40){
-				ISSMERROR("ksplit[i]<40");
-			}
-		  }
-		//  now do the element split
-		newNbOfQuad = 4*NbOfQuad;
-		nbv = k;
-		kkk = nbt;
-		ksplit[-1] = nbt;
-		// look on  old true  triangles 
-
-		for (i=0;i<nbtsave;i++){
-			int  nbmkadj=0;
-			long mkadj [100];
-			mkadj[0]=i;
-			long kk=ksplit[i]/10;
-			int  ke=(int) (ksplit[i]%10);
-			if (kk>=7 || kk<=0){
-				ISSMERROR("kk>=7 || kk<=0");
-			}
-
-			// def the numbering   k (edge) i vertex 
-			int k0 = ke;
-			int k1 = NextEdge[k0];
-			int k2 = PreviousEdge[k0];
-			int i0 = OppositeVertex[k0];
-			int i1 = OppositeVertex[k1];
-			int i2 = OppositeVertex[k2];
-
-			Triangle &t0=triangles[i];
-			MeshVertex * v0=t0(i0);           
-			MeshVertex * v1=t0(i1);           
-			MeshVertex * v2=t0(i2);
-
-			if (nbmkadj>=10){
-				ISSMERROR("nbmkadj>=10");
-			}
-			// --------------------------
-			TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
-			// save the flag Hidden
-			int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
-			// un set all adj -- save Hidden flag --
-			t0.SetAdj2(0,0,hid[0]);
-			t0.SetAdj2(1,0,hid[1]);
-			t0.SetAdj2(2,0,hid[2]);
-			// --  remake 
-			switch  (kk) {
-				case 1: break;// nothing 
-				case 2: // 
-						  {
-							Triangle &t1=triangles[kkk++];
-							t1=t0;
-							if (kedge[3*i+i0]<0){
-								ISSMERROR("kedge[3*i+i0]<0");
-							}
-							MeshVertex * v3 = vertices + kedge[3*i+k0];
-
-							t0(i2) = v3;
-							t1(i1) = v3;
-							t0.SetAllFlag(k2,0);
-							t1.SetAllFlag(k1,0);
-						  } 
-						break; 
-				case 3: //
-						  {
-							Triangle &t1=triangles[kkk++];
-							Triangle &t2=triangles[kkk++];
-							t2=t1=t0;
-							if (kedge[3*i+k1]<0){
-								ISSMERROR("kedge[3*i+k1]<0");
-							}
-							if (kedge[3*i+k2]<0){
-								ISSMERROR("kedge[3*i+k2]<0");
-							}
-
-							MeshVertex * v01 = vertices + kedge[3*i+k2];
-							MeshVertex * v02 = vertices + kedge[3*i+k1]; 
-							t0(i1) = v01; 
-							t0(i2) = v02; 
-							t1(i2) = v02;
-							t1(i0) = v01; 
-							t2(i0) = v02; 
-							t0.SetAllFlag(k0,0);
-							t1.SetAllFlag(k1,0);
-							t1.SetAllFlag(k0,0);
-							t2.SetAllFlag(k2,0);
-						  } 
-						break;
-				case 4: // 
-				case 6: // split in 4 
-						  {
-							Triangle &t1=triangles[kkk++];
-							Triangle &t2=triangles[kkk++];
-							Triangle &t3=triangles[kkk++];
-							t3=t2=t1=t0;
-							if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
-								ISSMERROR("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
-							}
-							MeshVertex * v12 = vertices + kedge[3*i+k0];
-							MeshVertex * v02 = vertices + kedge[3*i+k1]; 
-							MeshVertex * v01 = vertices + kedge[3*i+k2];
-							t0(i1) = v01;
-							t0(i2) = v02;
-							t0.SetAllFlag(k0,hid[k0]);
-
-							t1(i0) = v01;
-							t1(i2) = v12;
-							t0.SetAllFlag(k1,hid[k1]);
-
-							t2(i0) = v02;
-							t2(i1) = v12;
-							t2.SetAllFlag(k2,hid[k2]);
-
-							t3(i0) = v12;
-							t3(i1) = v02;
-							t3(i2) = v01;
-
-							t3.SetAllFlag(0,hid[0]);	   
-							t3.SetAllFlag(1,hid[1]);	   
-							t3.SetAllFlag(2,hid[2]);
-
-							if ( kk == 6)
-							  {
-
-								Triangle &t4=triangles[kkk++];
-								Triangle &t5=triangles[kkk++];
-
-								t4 = t3;
-								t5 = t3;
-
-								t0.SetHidden(k0);
-								t1.SetHidden(k1);
-								t2.SetHidden(k2);
-								t3.SetHidden(0);
-								t4.SetHidden(1);
-								t5.SetHidden(2);
-
-								if (nbv < nbvx ) 
-								  {
-									vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
-									vertices[nbv].ReferenceNumber =0;
-									vertices[nbv].DirOfSearch =NoDirOfSearch;
-									//vertices[nbv].i = toI2(vertices[nbv].r);
-									double a3[]={1./3.,1./3.,1./3.};
-									vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
-									MeshVertex * vc =  vertices +nbv++;
-									t3(i0) = vc;
-									t4(i1) = vc;
-									t5(i2) = vc;
-
-								  }
-								else
-								 goto Error; 
-							  }
-
-						  } 
-						break;         
-			}
-
-			// save all the new triangles
-			mkadj[nbmkadj++]=i;
-			long jj;
-			if (t0.link) 
-			 for (jj=nbt;jj<kkk;jj++)
-				{
-				 triangles[jj].link=t0.link;
-				 t0.link= triangles+jj;
-				 mkadj[nbmkadj++]=jj;
-				}
-			if (nbmkadj>13){// 13 = 6 + 4 +
-				ISSMERROR("nbmkadj>13");
-			}
-
-			if (kk==6)  newNbOfQuad+=3;
-			for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
-			ksplit[i]= nbt; // save last adresse of the new triangles
-			kkk = nbt;
-		  }
-
-		for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
-
-		if(withBackground)
-		 for (i=0;i<BTh.nbv;i++)
-		  BTh.vertices[i].m =  BTh.vertices[i].m*2.;
-
-
-		ret = 2;
-		if (nbt>= nbtx) goto Error; // bug 
-		if (nbv>= nbvx) goto Error; // bug 
-		// generation of the new triangles 
-
-		SetIntCoor("In SplitElement"); 
-
-		ReMakeTriangleContainingTheVertex();
-		if(withBackground)
-		 BTh.ReMakeTriangleContainingTheVertex();
-
-		delete [] edges;
-		edges = newedges;
-		nbe = newnbe;
-		NbOfQuad = newNbOfQuad;
-
-		for (i=0;i<NbSubDomains;i++)
-		  { 
-			long k = subdomains[i].edge- edges;
-			subdomains[i].edge =  edges+2*k; // spilt all edge in 2 
-		  }
-
-		if (ksplitarray) delete [] ksplitarray;
-		if (kedge) delete [] kedge;
-		if (edge4) delete edge4;
-		if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
-		VerticesOnGeomEdge= newVerticesOnGeomEdge;
-		if(VertexOnBThEdge) delete []  VertexOnBThEdge;
-		VertexOnBThEdge = newVertexOnBThEdge;
-		NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
-		NbVertexOnBThEdge=newNbVertexOnBThEdge;
-		//  ReMakeTriangleContainingTheVertex();
-
-		ReconstructExistingMesh();
-
-		if (verbose>2){
-			printf("   number of quadrilaterals    = %i\n",NbOfQuad);
-			printf("   number of triangles         = %i\n",nbt-NbOutT- NbOfQuad*2);
-			printf("   number of outside triangles = %i\n",NbOutT);
-		}
-
-		return 0; //ok
-
-Error:
-		nbv = nbvold;
-		nbt = nbtold;
-		NbOutT = NbOutTold;
-		// cleaning memory ---
-		delete newedges;
-		if (ksplitarray) delete [] ksplitarray;
-		if (kedge) delete [] kedge;
-		if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
-		if (edge4) delete edge4;
-		if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
-
-		return ret; // ok 
-	}
-	/*}}}1*/
-/*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{1*/
-long  Mesh::SplitInternalEdgeWithBorderVertices(){
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/
-
-	long NbSplitEdge=0;
-	SetVertexFieldOn();  
-	long it;
-	long nbvold=nbv;
-	long int verbose=2;
-	for (it=0;it<nbt;it++){
-		Triangle &t=triangles[it];
-		if (t.link)
-		 for (int j=0;j<3;j++)
-		  if(!t.Locked(j) && !t.Hidden(j)){
-			  Triangle &tt = *t.TriangleAdj(j);
-			  if ( &tt && tt.link && it < Number(tt)) 
-				 { // an internal edge 
-				  MeshVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
-				  MeshVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
-				  if (v0.onGeometry && v1.onGeometry){
-					  R2 P= ((R2) v0 + (R2) v1)*0.5;
-					  if ( nbv<nbvx) {
-						  vertices[nbv].r = P;
-						  vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
-						  vertices[nbv].ReferenceNumber=0;
-						  vertices[nbv].DirOfSearch = NoDirOfSearch ;
-					  }
-					  NbSplitEdge++;
-				  }
-				 }
-		  }
-	}
-	ReMakeTriangleContainingTheVertex();    
-	if (nbvold!=nbv){
-		long  iv = nbvold;
-		long NbSwap = 0;
-		Icoor2 dete[3];  
-		for (int i=nbvold;i<nbv;i++) {// for all the new point
-			MeshVertex & vi = vertices[i];
-			vi.i = toI2(vi.r);
-			vi.r = toR2(vi.i);
-
-			// a good new point 
-			vi.ReferenceNumber=0; 
-			vi.DirOfSearch =NoDirOfSearch;
-			Triangle *tcvi = FindTriangleContaining(vi.i,dete);
-			if (tcvi && !tcvi->link) {
-				printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
-			}
-
-			quadtree->Add(vi);
-			if (!tcvi || tcvi->det<0){// internal
-				ISSMERROR("!tcvi || tcvi->det < 0");
-			}
-			AddVertex(vi,tcvi,dete);
-			NbSwap += vi.Optim(1);          
-			iv++;
-			//      }
-	}
-	if (verbose>3) {
-		printf("   number of points: %i\n",iv);
-		printf("   number of swap to  split internal edges with border vertices: %i\n",NbSwap);
-		nbv = iv;
-	}
-}
-if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices  to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold));
-if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge);
-
-return  NbSplitEdge;
-}
-/*}}}1*/
 /*FUNCTION Mesh::TriangleReferenceList{{{1*/
 long  Mesh::TriangleReferenceList(long* reft) const {
@@ -5484,6 +5473,6 @@
 				ISSMERROR("kkk>=1000");
 			}
-			MeshVertex  &vI =  *edge.EdgeVertex(0);
-			MeshVertex  &vJ =  *edge.EdgeVertex(1);
+			BamgVertex  &vI =  *edge.EdgeVertex(0);
+			BamgVertex  &vJ =  *edge.EdgeVertex(1);
 			I2 I=vI, J=vJ, IJ= J-I;
 			IJ_IA = (IJ ,(A-I));
@@ -5525,5 +5514,5 @@
 		// the probleme is in case of  the fine and long internal hole
 		// for exemple neart the training edge of a wing
-		MeshVertex * s=0,*s1=0, *s0=0;
+		BamgVertex * s=0,*s1=0, *s0=0;
 		Icoor2 imax = MaxICoor22;
 		Icoor2 l0 = imax,l1 = imax;
@@ -5612,5 +5601,5 @@
 				if ((link + linkp) == 1) 
 				  { // a boundary edge 
-					MeshVertex * st = edge.EdgeVertex(1);
+					BamgVertex * st = edge.EdgeVertex(1);
 					I2 I=*st;
 					Icoor2  ll = Norme2_2 (C-I);
@@ -5657,5 +5646,5 @@
 	/*}}}1*/
 /*FUNCTION ForceEdge{{{1*/
-int ForceEdge(MeshVertex &a, MeshVertex & b,TriangleAdjacent & taret)  { 
+int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret)  { 
 	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/
 
@@ -5668,5 +5657,5 @@
 
 	TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]);
-	MeshVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
+	BamgVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
 	// we turn around a in the  direct sens  
 
@@ -5693,5 +5682,5 @@
 		if((det1 < 0) && (det2 >0)) { 
 			// try to force the edge 
-			MeshVertex * va = &a, *vb = &b;
+			BamgVertex * va = &a, *vb = &b;
 			tc = Previous(tc);
 			if (!v1 || !v2){
@@ -5703,5 +5692,5 @@
 				 ISSMERROR("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i) too big",NbSwap,l);
 			 }
-			MeshVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
+			BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
 			if (( aa == &a ) && (bb == &b) ||  (bb ==  &a ) && (aa == &b)) {
 				tc.SetLock();
@@ -5733,5 +5722,5 @@
 /*}}}1*/
 /*FUNCTION swap{{{1*/
-void  swap(Triangle *t1,short a1, Triangle *t2,short a2, MeshVertex *s1,MeshVertex *s2,Icoor2 det1,Icoor2 det2){ 
+void  swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2){ 
 	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/
 	// --------------------------------------------------------------
@@ -5778,5 +5767,5 @@
 /*}}}1*/
 	/*FUNCTION SwapForForcingEdge{{{1*/
-	int SwapForForcingEdge(MeshVertex   *  & pva ,MeshVertex  * &   pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {
+	int SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/
 		// l'arete ta coupe l'arete pva pvb
@@ -5795,7 +5784,7 @@
 		}
 
-		MeshVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
-		MeshVertex & s1= (*t1)[OppositeVertex[a1]];
-		MeshVertex & s2= (*t2)[OppositeVertex[a2]];
+		BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
+		BamgVertex & s1= (*t1)[OppositeVertex[a1]];
+		BamgVertex & s2= (*t2)[OppositeVertex[a2]];
 
 
Index: /issm/trunk/src/c/objects/Bamg/Mesh.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.h	(revision 5120)
@@ -22,49 +22,50 @@
 
 	class Mesh {
+
 		public:
 
-			Geometry & Gh;   // Geometry
-			Mesh & BTh; // Background Mesh Bth==*this =>no  background 
-			long NbRef;      // counter of ref on the this class if 0 we can delete
-			long nbvx,nbtx;  // nombre max  de sommets , de  triangles
-			long nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,
-			long NbOfQuad;  // nb of quadrangle 
-			long NbSubDomains;
-			long NbOutT;    // Nb of oudeside triangle
-			long NbOfTriangleSearchFind;
-			long NbOfSwapTriangle;
-			MeshVertex* vertices;
-			long NbVerticesOnGeomVertex;
-			VertexOnGeom * VerticesOnGeomVertex;
-			long NbVerticesOnGeomEdge;
-			VertexOnGeom * VerticesOnGeomEdge;
-			long NbVertexOnBThVertex;
-			VertexOnVertex *VertexOnBThVertex;
-			long NbVertexOnBThEdge;
-			VertexOnEdge *VertexOnBThEdge;
-			long NbCrackedVertices;
-			long* CrackedVertices;
-			long NbCrackedEdges;
-			CrackedEdge* CrackedEdges;
-			R2 pmin,pmax;    // extrema
-			double coefIcoor; // coef to integer Icoor1;
-			Triangle* triangles;
-			Edge* edges; 
-			QuadTree* quadtree;
-			MeshVertex** ordre;
-			SubDomain* subdomains;
-			ListofIntersectionTriangles lIntTria;
+			Geometry                    & Gh;                    // Geometry
+			Mesh                        & BTh;                   // Background Mesh Bth== *this =>no background
+			long                          NbRef;                 // counter of ref on the this class if 0 we can delete
+			long                          maxnbv,nbtx;           // nombre max de sommets , de triangles
+			long                          nt,nbv,nbt,nbe;        // nb of legal triangles, nb of vertex, of triangles, of edges with reference,
+			long                          NbOfQuad;              // nb of quadrangle
+			long                          NbSubDomains;
+			long                          NbOutT;                // Nb of oudeside triangle
+			long                          NbOfTriangleSearchFind;
+			long                          NbOfSwapTriangle;
+			BamgVertex                   *vertices;
+			long                          NbVerticesOnGeomVertex;
+			VertexOnGeom                 *VerticesOnGeomVertex;
+			long                          NbVerticesOnGeomEdge;
+			VertexOnGeom                 *VerticesOnGeomEdge;
+			long                          NbVertexOnBThVertex;
+			VertexOnVertex               *VertexOnBThVertex;
+			long                          NbVertexOnBThEdge;
+			VertexOnEdge                 *VertexOnBThEdge;
+			long                          NbCrackedVertices;
+			long                         *CrackedVertices;
+			long                          NbCrackedEdges;
+			CrackedEdge                  *CrackedEdges;
+			R2                            pmin,pmax;             // extrema
+			double                        coefIcoor;             // coef to integer Icoor1;
+			Triangle                     *triangles;
+			Edge                         *edges;
+			QuadTree                     *quadtree;
+			BamgVertex                  **ordre;
+			SubDomain                    *subdomains;
+			ListofIntersectionTriangles   lIntTria;
 
 			//Constructors/Destructors
 			Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts);
 			Mesh(double* index,double* x,double* y,int nods,int nels);
-			Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long nbvxx=0 ); //copy operator
+			Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator
 			Mesh(const Mesh &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
-			Mesh(long nbvx,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
-				try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
+			Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
+				try {TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);}
 				catch(...) { this->~Mesh(); throw; }
 			}
-			Mesh(long nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
-				try { GeomToTriangles0(nbvx,bamgopts);}
+			Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
+				try { TriangulateFromGeom0(maxnbv,bamgopts);}
 				catch(...) { this->~Mesh(); throw; }
 			}
@@ -72,6 +73,6 @@
 
 			//Operators
-			const MeshVertex & operator[]  (long i) const { return vertices[i];};
-			MeshVertex & operator[](long i) { return vertices[i];};
+			const BamgVertex & operator[]  (long i) const { return vertices[i];};
+			BamgVertex & operator[](long i) { return vertices[i];};
 			const Triangle & operator()  (long i) const { return triangles[i];};
 			Triangle & operator()(long i) { return triangles[i];};
@@ -87,5 +88,5 @@
 				return  R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
 			}
-			void AddVertex(MeshVertex & s,Triangle * t,Icoor2 *  =0) ;
+			void AddVertex(BamgVertex & s,Triangle * t,Icoor2 *  =0) ;
 			void Insert();
 			void ForceBoundary();
@@ -110,13 +111,13 @@
 			void SmoothingVertex(int =3,double=0.3);
 			Metric MetricAt (const R2 &) const;
-			GeometricalEdge* ProjectOnCurve( Edge & AB, MeshVertex &  A, MeshVertex & B,double theta, MeshVertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
+			GeometricalEdge* ProjectOnCurve( Edge & AB, BamgVertex &  A, BamgVertex & B,double theta, BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
 			long Number(const Triangle & t) const  { return &t - triangles;}
 			long Number(const Triangle * t) const  { return t - triangles;}
-			long Number(const MeshVertex & t) const  { return &t - vertices;}
-			long Number(const MeshVertex * t) const  { return t - vertices;}
+			long Number(const BamgVertex & t) const  { return &t - vertices;}
+			long Number(const BamgVertex * t) const  { return t - vertices;}
 			long Number(const Edge & t) const  { return &t - edges;}
 			long Number(const Edge * t) const  { return t - edges;}
 			long Number2(const Triangle * t) const  { return t - triangles; }
-			MeshVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
+			BamgVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
 			Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
 			void ReadMesh(double* index,double* x,double* y,int nods,int nels);
@@ -155,7 +156,7 @@
 
 		private:
-			void GeomToTriangles1(long nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
-			void GeomToTriangles0(long nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
-			void PreInit(long);
+			void TriangulateFromGeom1(long maxnbv,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
+			void TriangulateFromGeom0(long maxnbv,BamgOpts* bamgopts);// the real constructor mesh generator
+			void Init(long);
 
 	};
@@ -166,9 +167,9 @@
 	void  swap(Triangle *t1,short a1,
 				Triangle *t2,short a2,
-				MeshVertex *s1,MeshVertex *s2,Icoor2 det1,Icoor2 det2);
-	int SwapForForcingEdge(MeshVertex   *  & pva ,MeshVertex  * &   pvb ,
+				BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2);
+	int SwapForForcingEdge(BamgVertex   *  & pva ,BamgVertex  * &   pvb ,
 				TriangleAdjacent & tt1,Icoor2 & dets1,
 				Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
-	int ForceEdge(MeshVertex &a, MeshVertex & b,TriangleAdjacent & taret) ;
+	int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret) ;
 	inline TriangleAdjacent Previous(const TriangleAdjacent & ta){
 		return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);
@@ -183,5 +184,5 @@
 		int j=i;i=on->DirAdj[i];on=on->Adj[j];
 	}
-	inline double qualite(const MeshVertex &va,const MeshVertex &vb,const MeshVertex &vc){
+	inline double qualite(const BamgVertex &va,const BamgVertex &vb,const BamgVertex &vc){
 		double ret; 
 		I2 ia=va,ib=vb,ic=vc;
Index: sm/trunk/src/c/objects/Bamg/MeshVertex.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/MeshVertex.cpp	(revision 5119)
+++ 	(revision )
@@ -1,251 +1,0 @@
-#include <cstdio>
-#include <cstring>
-#include <cmath>
-#include <ctime>
-
-#include "../objects.h"
-
-namespace bamg {
-
-	/*Methods*/
-	/*FUNCTION MeshVertex::Smoothing{{{1*/
-	double  MeshVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/
-
-		register MeshVertex* s=this;
-		MeshVertex &vP = *s,vPsave=vP;
-
-		register Triangle* tbegin= t , *tria = t , *ttc;
-
-		register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc;
-		R2 P(s->r),PNew(0,0);
-		do {
-			k++; 
-
-			if (!tria->Hidden(j)){
-				MeshVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]]; 
-
-				R2 Q = vQ,QP(P-Q);
-				double lQP = LengthInterpole(vP,vQ,QP);
-				PNew += Q+QP/Max(lQP,1e-20);
-				kk ++;
-			}
-			ttc =  tria->TriangleAdj(j);
-			jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
-			tria = ttc;
-			j = NextEdge[jc];
-			if (k>=2000){
-				ISSMERROR("k>=2000 (Maximum number of iterations reached)");
-			}
-		} while ( tbegin != tria); 
-		if (kk<4) return 0;
-		PNew = PNew/(double)kk;
-		R2 Xmove((PNew-P)*omega);
-		PNew = P+Xmove;
-		double delta=Norme2_2(Xmove); 
-
-		Icoor2 deta[3];
-		I2 IBTh  = BTh.toI2(PNew);
-
-		tstart=BTh.FindTriangleContaining(IBTh,deta,tstart);  
-
-		if (tstart->det <0){ // outside
-			double ba,bb;
-			TriangleAdjacent edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ;
-			tstart = edge;
-			vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));
-		}
-		else { // inside
-			double   aa[3];
-			double s = deta[0]+deta[1]+deta[2];
-			aa[0]=deta[0]/s;
-			aa[1]=deta[1]/s;
-			aa[2]=deta[2]/s;
-			vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]);
-		}
-
-		// recompute the det of the triangle
-		vP.r = PNew;
-
-		vP.i = Th.toI2(PNew);
-
-		MeshVertex vPnew = vP;
-
-		int ok=1;
-		int loop=1;
-		k=0;
-		while (ok){
-			ok =0;
-			do {
-				k++; 
-				double detold = tria->det;
-				tria->det =  bamg::det( (*tria)[0],(*tria)[1]  ,(*tria)[2]);
-				if (loop) {
-					MeshVertex *v0,*v1,*v2,*v3;
-					if (tria->det<0) ok =1;			       
-					else if (tria->Quadrangle(v0,v1,v2,v3))
-					  {
-						vP = vPsave;
-						double qold =QuadQuality(*v0,*v1,*v2,*v3);
-						vP = vPnew;
-						double qnew =QuadQuality(*v0,*v1,*v2,*v3);
-						if (qnew<qold) ok = 1;
-					  }
-					else if ( (double)tria->det < detold/2 ) ok=1;
-
-				}
-				tria->SetUnMarkUnSwap(0);
-				tria->SetUnMarkUnSwap(1);
-				tria->SetUnMarkUnSwap(2);
-				ttc =  tria->TriangleAdj(j);
-				jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
-				tria = ttc;
-				j = NextEdge[jc];
-				if (k>=2000){
-					ISSMERROR("k>=2000");
-				}
-			}while ( tbegin != tria); 
-
-			if (ok && loop) vP=vPsave; // no move 
-			loop=0;
-		}
-		return delta;
-	}
-	/*}}}1*/
-	/*FUNCTION MeshVertex::MetricFromHessian{{{1*/
-	void MeshVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){
-		/*Compute Metric from Hessian*/
-
-		/*get options*/
-		double power=(bamgopts->power);
-		double anisomax=(bamgopts->anisomax);
-		double CutOff=bamgopts->cutoff;
-		double hmin=(bamgopts->hmin);
-		double hmax=(bamgopts->hmax);
-		double coef=bamgopts->coeff;
-		int    Metrictype=(bamgopts->Metrictype);
-
-		/*Intermediary*/
-		double ci;
-
-		/*compute multiplicative coefficient depending on Metric Type (2/9 because it is 2d)*/
-
-		//Absolute Error
-		/*
-		 *            2         1       
-		 *Metric M = ---  ------------   Abs(Hessian)
-		 *            9   err * coeff^2  
-		 */
-		if (Metrictype==0){
-			ci= 2.0/9.0 * 1/(err*coef*coef);
-		}
-
-		//Relative Error
-		/*
-		 *            2         1            Abs(Hessian)
-		 *Metric M = ---  ------------  ----------------------
-		 *            9   err * coeff^2  max( |s| , cutoff*max(|s|) )
-		 *
-		 */
-		else if (Metrictype==1){
-			ci= 2.0/9.0 * 1/(err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
-		}
-
-		//Rescaled absolute error
-		/*
-		 *            2         1            Abs(Hessian)
-		 *Metric M = ---  ------------  ---------------------- 
-		 *            9   err * coeff^2       (smax-smin)
-		 */
-		else if (Metrictype==2){
-			ci= 2.0/9.0 * 1/(err*coef*coef) * 1/(smax-smin);
-		}
-		else{
-			ISSMERROR("Metrictype %i not supported yet (use 0,1 or 2(default))",Metrictype);
-		}
-
-		//initialize metric Miv with ci*H
-		Metric Miv(Hxx*ci,Hyx*ci,Hyy*ci);
-
-		//Get eigen values and vectors of Miv
-		MatVVP2x2 Vp(Miv);
-
-		//move eigen valuse to their absolute values
-		Vp.Abs();
-
-		//Apply a power if requested by user
-		if(power!=1.0) Vp.pow(power);
-
-		//modify eigen values according to hmin and hmax
-		Vp.Maxh(hmax);
-		Vp.Minh(hmin);
-
-		//Bound anisotropy by 1/(anisomax)^2
-		Vp.BoundAniso2(1/(anisomax*anisomax));
-
-		//rebuild Metric from Vp
-		Metric MVp(Vp);
-
-		//Apply Metric to vertex
-		m.IntersectWith(MVp);
-
-	}
-	/*}}}1*/
-	/*FUNCTION MeshVertex::Echo {{{1*/
-
-	void MeshVertex::Echo(void){
-
-		printf("Vertex:\n");
-		printf("  integer   coordinates i.x: %i, i.y: %i\n",i.x,i.y);
-		printf("  Euclidean coordinates r.x: %g, r.y: %g\n",r.x,r.y);
-		printf("  ReferenceNumber = %i\n",ReferenceNumber);
-		m.Echo();
-
-		return;
-	}
-	/*}}}*/
-	/*FUNCTION MeshVertex::Optim {{{1*/
-	long MeshVertex::Optim(int i,int koption){ 
-		long ret=0;
-		if ( t && (vint >= 0 ) && (vint <3) ){
-			ret = t->Optim(vint,koption);
-			if(!i){
-				t =0; // for no future optime 
-				vint= 0;
-			}
-		}
-		return ret;
-	}
-	/*}}}*/
-
-	/*Intermediary*/
-	/*FUNCTION QuadQuality{{{1*/
-	double QuadQuality(const MeshVertex & a,const MeshVertex &b,const MeshVertex &c,const MeshVertex &d) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/
-
-		// calcul de 4 angles --
-		R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d);
-		R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D);
-		//  Move(A),Line(B),Line(C),Line(D),Line(A);
-		const Metric & Ma  = a;
-		const Metric & Mb  = b;
-		const Metric & Mc  = c;
-		const Metric & Md  = d;
-
-		double lAB=Norme2(AB);
-		double lBC=Norme2(BC);
-		double lCD=Norme2(CD);
-		double lDA=Norme2(DA);
-		AB /= lAB;  BC /= lBC;  CD /= lCD;  DA /= lDA;
-		// version aniso 
-		double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB);
-		double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC);
-		double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD);
-		double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA);
-		double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA));
-		if (sinmin<=0) return sinmin;
-		return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA)));
-	}
-	/*}}}1*/
-
-} 
Index: sm/trunk/src/c/objects/Bamg/MeshVertex.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/MeshVertex.h	(revision 5119)
+++ 	(revision )
@@ -1,55 +1,0 @@
-#ifndef _MESHVERTEX_H_
-#define _MESHVERTEX_H_
-
-#include "./include.h"
-#include "./Metric.h"
-#include "./Direction.h"
-#include "./BamgOpts.h"
-
-namespace bamg {
-
-	//classes
-	class Triangle;
-	class Mesh;
-	class VertexOnGeom;
-	class VertexOnEdge;
-
-	class MeshVertex {
-
-		public:
-			I2 i;  // integer coordinates
-			R2 r;  // real coordinates
-			Metric m;
-			long ReferenceNumber;
-			Direction DirOfSearch;
-			short vint;  // the vertex number in triangle; varies between 0 and 2 in t
-			union {
-				Triangle* t;   // one triangle which is containing the vertex
-				long      color;
-				MeshVertex*   to;  // used in geometry MeshVertex to know the Mesh MeshVertex associated 
-				VertexOnGeom* onGeometry;        // if vint == 8; // set with Mesh::SetVertexFieldOn()
-				MeshVertex*       onBackgroundVertex;// if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()
-				VertexOnEdge* onBackgroundEdge;  // if vint == 32 on Background edge
-			};
-
-			//Operators
-			operator I2() const {return i;}             // Cast operator
-			operator const R2 & () const {return r;}    // Cast operator
-			operator Metric () const {return m;}        // Cast operator
-			double operator()(R2 x) const { return m(x);} // Get x in the metric m
-
-			//methods (No constructor and no destructors...)
-			double Smoothing(Mesh & ,const Mesh & ,Triangle  * & ,double =1);
-			void   MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
-			void   Echo();
-			int    ref() const { return ReferenceNumber;}
-			long   Optim(int =1,int =0); 
-
-			//inline functions
-			inline void Set(const MeshVertex &rec,const Mesh & ,Mesh & ){*this=rec;}
-	};
-
-	//Intermediary
-	double QuadQuality(const MeshVertex &,const MeshVertex &,const MeshVertex &,const MeshVertex &);
-}
-#endif
Index: /issm/trunk/src/c/objects/Bamg/QuadTree.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/QuadTree.cpp	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/QuadTree.cpp	(revision 5120)
@@ -24,7 +24,7 @@
 	/*}}}*/
 	/*DOCUMENTATION What is a QuadTree? {{{1
-	 * A Quadtree is a very simple way to group the vertices according
-	 * to their location. A square that holds all the points of the mesh
-	 * (or the geometry) is divided into 4 boxes. As soom as one box
+	 * A Quadtree is a very simple way to group vertices according
+	 * to their locations. A square that holds all the points of the mesh
+	 * (or the geometry) is divided into 4 boxes. As soon as one box
 	 * hold more than 4 vertices, it is divided into 4 new boxes, etc...
 	 * There cannot be more than MAXDEEP (=30) subdivision.
@@ -67,5 +67,5 @@
 	 * 0 1 1 1 .... 1    0 1 1 1 .... 1 in binary
 	 *  \--   29  --/     \--   29  --/
-	 * Using the binaries is therefor very easy to locate a vertex in a box:
+	 * Using binaries is therefore very easy to locate a vertex in a box:
 	 * we just need to look at the bits from the left to the right (See ::Add)
 	 }}}1*/
@@ -76,5 +76,5 @@
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/QuadTree)*/
 
-		lenStorageQuadTreeBox(t->nbvx/8+10),
+		lenStorageQuadTreeBox(t->maxnbv/8+10),
 		th(t),
 		NbQuadTreeBox(0),
@@ -120,5 +120,5 @@
 	/*Methods*/
 	/*FUNCTION QuadTree::Add{{{1*/
-	void  QuadTree::Add(MeshVertex &w){
+	void  QuadTree::Add(BamgVertex &w){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/Add)*/
 
@@ -163,5 +163,5 @@
 
 			//Copy the 4 vertices in the current QuadTreebox
-			MeshVertex* v4[4];
+			BamgVertex* v4[4];
 			v4[0]= b->v[0];
 			v4[1]= b->v[1];
@@ -205,5 +205,5 @@
 	/*}}}1*/
 	/*FUNCTION QuadTree::NearestVertex{{{1*/
-	MeshVertex*  QuadTree::NearestVertex(Icoor1 i,Icoor1 j) {
+	BamgVertex*  QuadTree::NearestVertex(Icoor1 i,Icoor1 j) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertex)*/
 
@@ -224,5 +224,5 @@
 		Icoor1   jplus( j<MaxISize?(j<0?0:j):MaxISize-1);
 		//initial nearest vertex pointer
-		MeshVertex*  vn=NULL;
+		BamgVertex*  vn=NULL;
 
 		//Get initial Quadtree box (largest)
@@ -293,5 +293,5 @@
 
 				//if the current subbox is holding vertices,
-				if (b->n>0){ // MeshVertex QuadTreeBox not empty
+				if (b->n>0){ // BamgVertex QuadTreeBox not empty
 					NbVerticesSearch++;
 					I2 i2 =  b->v[k]->i;
@@ -344,5 +344,5 @@
 	/*}}}1*/
 	/*FUNCTION QuadTree::NearestVertexWithNormal{{{1*/
-	MeshVertex*  QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) {
+	BamgVertex*  QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertexWithNormal)*/
 
@@ -358,5 +358,5 @@
 		Icoor1  jplus( j<MaxISize?(j<0?0:j):MaxISize-1);
 
-		MeshVertex *vn=0;
+		BamgVertex *vn=0;
 
 		// init for optimisation ---
@@ -411,5 +411,5 @@
 				int k = pi[l];
 
-				if (b->n>0) // MeshVertex QuadTreeBox none empty
+				if (b->n>0) // BamgVertex QuadTreeBox none empty
 				  { 
 					NbVerticesSearch++;
@@ -472,5 +472,5 @@
 	/*}}}1*/
 	/*FUNCTION QuadTree::ToClose {{{1*/
-	MeshVertex *   QuadTree::ToClose(MeshVertex & v,double seuil,Icoor1 hx,Icoor1 hy){
+	BamgVertex *   QuadTree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/ToClose)*/
 
@@ -489,5 +489,5 @@
 		Icoor1 i0=0,j0=0;
 
-		//  MeshVertex *vn=0;
+		//  BamgVertex *vn=0;
 
 		if (!root->n)
@@ -506,5 +506,5 @@
 				register int k = pi[l];
 
-				if (b->n>0) // MeshVertex QuadTreeBox none empty
+				if (b->n>0) // BamgVertex QuadTreeBox none empty
 				  { 
 					NbVerticesSearch++;
Index: /issm/trunk/src/c/objects/Bamg/QuadTree.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/QuadTree.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/QuadTree.h	(revision 5120)
@@ -11,5 +11,5 @@
 
 	class Mesh;
-	class MeshVertex;
+	class BamgVertex;
 
 	class QuadTree{
@@ -18,9 +18,9 @@
 				public:
 					long n; 
-					//contains only one object form the list (either MeshVertex or QuadTreeBox)
-					// if n < 4 => MeshVertex else =>  QuadTreeBox;
+					//contains only one object form the list (either BamgVertex or QuadTreeBox)
+					// if n < 4 => BamgVertex else =>  QuadTreeBox;
 					union{
 						QuadTreeBox* b[4];
-						MeshVertex* v[4];
+						BamgVertex* v[4];
 					};
 			};
@@ -50,9 +50,9 @@
 			long    NbQuadTreeBox,NbVertices;
 			long    NbQuadTreeBoxSearch,NbVerticesSearch;
-			MeshVertex* NearestVertex(Icoor1 i,Icoor1 j);
-			MeshVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);
-			MeshVertex* ToClose(MeshVertex & ,double ,Icoor1,Icoor1);
+			BamgVertex* NearestVertex(Icoor1 i,Icoor1 j);
+			BamgVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);
+			BamgVertex* ToClose(BamgVertex & ,double ,Icoor1,Icoor1);
 			long    SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();}
-			void    Add( MeshVertex & w);
+			void    Add( BamgVertex & w);
 			QuadTreeBox* NewQuadTreeBox(){
 				if(! (sb->bc<sb->be)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
Index: /issm/trunk/src/c/objects/Bamg/R2.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/R2.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/R2.h	(revision 5120)
@@ -7,10 +7,11 @@
 namespace bamg {
 
-	template <class R,class RR> class P2xP2;
-	template <class R,class RR>
-	  class P2 {
+	template <class R,class RR> class P2 {
+
 		  public:  
+
 			  //fields
 			  R x,y;
+
 			  //functions
 			  P2 () :x(0),y(0) {};
@@ -19,6 +20,6 @@
 			  void Echo(){
 				  printf("Member of P2:\n");
-				  printf("   x: %g\n",x);
-				  printf("   y: %g\n",y);
+				  printf("   x: %g or %i\n",x,x);
+				  printf("   y: %g or %i\n",y,y);
 			  }
 			  //operators
@@ -33,15 +34,19 @@
 			  P2<R,RR> operator*=(const  R r) {x *= r;y *= r;return *this;}
 			  P2<R,RR> operator-=(const  P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;}
+
 	  };
 
-	template <class R,class RR>
-	  class P2xP2 {
+	template <class R,class RR> class P2xP2 {
+
 		  private:
+
 			  friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc){
 				  return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y);
 			  } 
 		  public:
+
 			  //fields
 			  P2<R,RR> x,y; 
+
 			  //functions
 			  P2xP2 (): x(),y()  {}
Index: /issm/trunk/src/c/objects/Bamg/Triangle.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Triangle.cpp	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Triangle.cpp	(revision 5120)
@@ -11,5 +11,5 @@
 	/*FUNCTION Triangle(Mesh *Th,long i,long j,long k) {{{1*/
 	Triangle::Triangle(Mesh *Th,long i,long j,long k) {
-		MeshVertex *v=Th->vertices;
+		BamgVertex *v=Th->vertices;
 		long nbv = Th->nbv;
 		if (i<0 || j<0 || k<0){
@@ -27,6 +27,6 @@
 	}
 	/*}}}*/
-	/*FUNCTION Triangle(MeshVertex *v0,MeshVertex *v1,MeshVertex *v2) {{{1*/
-	Triangle::Triangle(MeshVertex *v0,MeshVertex *v1,MeshVertex *v2){
+	/*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{1*/
+	Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){
 		TriaVertices[0]=v0;
 		TriaVertices[1]=v1;
@@ -154,5 +154,5 @@
 	/*}}}1*/
 	/*FUNCTION Triangle::Quadrangle {{{1*/
-	Triangle* Triangle::Quadrangle(MeshVertex * & v0,MeshVertex * & v1,MeshVertex * & v2,MeshVertex * & v3) const{
+	Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{
 		// return the other triangle of the quad if a quad or 0 if not a quat
 		Triangle * t =0;
@@ -186,8 +186,8 @@
 			 q= -1;
 			else if(option){ 
-				const MeshVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
-				const MeshVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
-				const MeshVertex & v1 = *TriaVertices[OppositeEdge[a]];
-				const MeshVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
+				const BamgVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
+				const BamgVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
+				const BamgVertex & v1 = *TriaVertices[OppositeEdge[a]];
+				const BamgVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
 				q =  QuadQuality(v0,v1,v2,v3); // do the float part
 			}
@@ -220,8 +220,8 @@
 		if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
 
-		register MeshVertex  *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]];
-		register MeshVertex  *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]];
-		register MeshVertex  *s1=t1->TriaVertices[OppositeVertex[a1]];
-		register MeshVertex  *s2=t2->TriaVertices[OppositeVertex[a2]];
+		register BamgVertex  *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]];
+		register BamgVertex  *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]];
+		register BamgVertex  *s1=t1->TriaVertices[OppositeVertex[a1]];
+		register BamgVertex  *s2=t2->TriaVertices[OppositeVertex[a2]];
 
 		Icoor2 det1=t1->det , det2=t2->det ;
Index: /issm/trunk/src/c/objects/Bamg/Triangle.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Triangle.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/Triangle.h	(revision 5120)
@@ -9,5 +9,5 @@
 	//classes
 	class Mesh;
-	class MeshVertex;
+	class BamgVertex;
 	class Triangle;
 
@@ -17,5 +17,5 @@
 
 		private:
-			MeshVertex*   TriaVertices[3];      // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
+			BamgVertex*   TriaVertices[3];      // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
 			Triangle* TriaAdjTriangles[3];  // 3 pointers toward the adjacent triangles
 			short     TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1
@@ -31,11 +31,11 @@
 			Triangle() {}
 			Triangle(Mesh *Th,long i,long j,long k);
-			Triangle(MeshVertex *v0,MeshVertex *v1,MeshVertex *v2);
+			Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2);
 
 			//Operators
-			const MeshVertex & operator[](int i) const {return *TriaVertices[i];};
-			MeshVertex & operator[](int i)  {return *TriaVertices[i];};
-			const MeshVertex * operator()(int i) const {return TriaVertices[i];};
-			MeshVertex * & operator()(int i)  {return TriaVertices[i];};
+			const BamgVertex & operator[](int i) const {return *TriaVertices[i];};
+			BamgVertex & operator[](int i)  {return *TriaVertices[i];};
+			const BamgVertex * operator()(int i) const {return TriaVertices[i];};
+			BamgVertex * & operator()(int i)  {return TriaVertices[i];};
 
 			//Methods
@@ -52,5 +52,5 @@
 			TriangleAdjacent Adj(int i)  const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
 			Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
-			Triangle* Quadrangle(MeshVertex * & v0,MeshVertex * & v1,MeshVertex * & v2,MeshVertex * & v3) const ;
+			Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ;
 			void  ReNumbering(Triangle *tb,Triangle *te, long *renu){
 				if (link  >=tb && link  <te) link  = tb + renu[link -tb];
@@ -59,5 +59,5 @@
 				if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];    
 			}
-			void ReNumbering(MeshVertex *vb,MeshVertex *ve, long *renu){
+			void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
 				if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
 				if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
@@ -118,5 +118,5 @@
 			double qualite() ;
 			void  Set(const Triangle &,const Mesh &,Mesh &);
-			int   In(MeshVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
+			int   In(BamgVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
 
 	};
Index: /issm/trunk/src/c/objects/Bamg/TriangleAdjacent.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/TriangleAdjacent.cpp	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/TriangleAdjacent.cpp	(revision 5120)
@@ -41,10 +41,10 @@
 	/*}}}*/
 	/*FUNCTION TriangleAdjacent::EdgeVertex {{{1*/
-	MeshVertex* TriangleAdjacent::EdgeVertex(const int & i) const {
+	BamgVertex* TriangleAdjacent::EdgeVertex(const int & i) const {
 		return t->TriaVertices[VerticesOfTriangularEdge[a][i]];
 	}
 	/*}}}*/
 	/*FUNCTION TriangleAdjacent::OppositeVertex {{{1*/
-	MeshVertex* TriangleAdjacent::OppositeVertex() const {
+	BamgVertex* TriangleAdjacent::OppositeVertex() const {
 		return t->TriaVertices[bamg::OppositeVertex[a]]; 
 	}
Index: /issm/trunk/src/c/objects/Bamg/TriangleAdjacent.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/TriangleAdjacent.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/TriangleAdjacent.h	(revision 5120)
@@ -3,5 +3,5 @@
 
 #include "./include.h"
-#include "./MeshVertex.h"
+#include "./BamgVertex.h"
 
 namespace bamg {
@@ -38,6 +38,6 @@
 			int  swap();
 			TriangleAdjacent Adj() const;
-			MeshVertex* EdgeVertex(const int & i) const;
-			MeshVertex* OppositeVertex() const;
+			BamgVertex* EdgeVertex(const int & i) const;
+			BamgVertex* OppositeVertex() const;
 			Icoor2& det() const;
 	};
Index: /issm/trunk/src/c/objects/Bamg/VertexOnEdge.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/VertexOnEdge.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/VertexOnEdge.h	(revision 5120)
@@ -9,22 +9,22 @@
 	//classes
 	class Mesh;
-	class MeshVertex;
+	class BamgVertex;
 
 	class VertexOnEdge {
 
 		public:
-			MeshVertex* v;
+			BamgVertex* v;
 			Edge*   be;
 			double abcisse;
 
 			//Constructors
-			VertexOnEdge( MeshVertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {}
+			VertexOnEdge( BamgVertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {}
 			VertexOnEdge(){}
 
 			//Operators
 			operator double () const { return abcisse;}
-			operator MeshVertex* () const { return v;}  
+			operator BamgVertex* () const { return v;}  
 			operator Edge* () const { return be;}  
-			MeshVertex & operator[](int i) const { return (*be)[i];}
+			BamgVertex & operator[](int i) const { return (*be)[i];}
 
 			//Methods
Index: /issm/trunk/src/c/objects/Bamg/VertexOnGeom.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/VertexOnGeom.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/VertexOnGeom.h	(revision 5120)
@@ -9,5 +9,5 @@
 	//classes
 	class Mesh;
-	class MeshVertex;
+	class BamgVertex;
 	class GeometricalEdge;
 
@@ -16,5 +16,5 @@
 		public:
 
-			MeshVertex* mv;
+			BamgVertex* mv;
 			double abscisse;  
 			union{ 
@@ -25,9 +25,9 @@
 			//Constructors/Destructors
 			VertexOnGeom(): mv(0),abscisse(0){gv=0;} 
-			VertexOnGeom(MeshVertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
-			VertexOnGeom(MeshVertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;}
+			VertexOnGeom(BamgVertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
+			VertexOnGeom(BamgVertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;}
 
 			//Operators
-			operator MeshVertex*() const  {return mv;}
+			operator BamgVertex*() const  {return mv;}
 			operator GeometricalVertex * () const  {return gv;}
 			operator GeometricalEdge * () const  {return ge;}
Index: /issm/trunk/src/c/objects/Bamg/VertexOnVertex.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/VertexOnVertex.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/VertexOnVertex.h	(revision 5120)
@@ -3,5 +3,5 @@
 
 #include "./include.h"
-#include "./MeshVertex.h"
+#include "./BamgVertex.h"
 
 namespace bamg {
@@ -13,9 +13,9 @@
 
 		public:
-			MeshVertex* v;
-			MeshVertex* bv;
+			BamgVertex* v;
+			BamgVertex* bv;
 
 			//Constructors
-			VertexOnVertex(MeshVertex * w,MeshVertex *bw) :v(w),bv(bw){}
+			VertexOnVertex(BamgVertex * w,BamgVertex *bw) :v(w),bv(bw){}
 			VertexOnVertex() {};
 
Index: /issm/trunk/src/c/objects/Bamg/typedefs.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/typedefs.h	(revision 5119)
+++ /issm/trunk/src/c/objects/Bamg/typedefs.h	(revision 5120)
@@ -6,10 +6,13 @@
 namespace bamg {
 
-	typedef int  Icoor1;  
-#if LONG_BIT > 63 //64 bits or more
+	/*Integer coordinates types*/
+	typedef int  Icoor1; 
+	#if LONG_BIT > 63 //64 bits or more
 	typedef long Icoor2;
-#else //32 bits
+	#else //32 bits
 	typedef double Icoor2;
-#endif
+	#endif
+
+	/*I2 and R2*/
 	typedef P2<Icoor1,Icoor2>  I2;
 	typedef P2xP2<short,long>  I2xI2;
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 5119)
+++ /issm/trunk/src/c/objects/objects.h	(revision 5120)
@@ -99,5 +99,5 @@
 #include "./Bamg/DoubleAndInt.h"
 #include "./Bamg/Direction.h"
-#include "./Bamg/MeshVertex.h"
+#include "./Bamg/BamgVertex.h"
 #include "./Bamg/TriangleAdjacent.h"
 #include "./Bamg/Edge.h"