Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3230)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3231)
@@ -143,9 +143,9 @@
 			if (bamgopts->field){
 				if (verbosity>1) printf("   Generating Metric from solution field...\n");
-				BTh.IntersectConsMetric(bamgopts);
+				BTh.AddMetric(bamgopts);
 			}
 		}
 
-		BTh.IntersectGeomMetric(bamgopts);
+		BTh.AddGeometryMetric(bamgopts);
 		if (verbosity>1) printf("   Smoothing Metric...");
 		if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation);
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3230)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3231)
@@ -705,10 +705,10 @@
 				return  R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
 			}
-			void Add( Vertex & s,Triangle * t,Icoor2 *  =0) ;
+			void AddVertex(Vertex & s,Triangle * t,Icoor2 *  =0) ;
 			void Insert();
 			void ForceBoundary();
 			void Renumber(BamgOpts* bamgopts);
 			void FindSubDomain(int OutSide=0);
-			Int4 ConsRefTriangle(Int4 *) const;
+			Int4 TriangleReferenceList(Int4 *) const;
 			void ShowHistogram() const;
 			void ShowRegulaty() const;
@@ -748,13 +748,13 @@
 			void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
 			void WriteMetric(BamgOpts* bamgopts);
-			void IntersectConsMetric(BamgOpts* bamgopts);
+			void AddMetric(BamgOpts* bamgopts);
 			void BuildMetric0(BamgOpts* bamgopts);
 			void BuildMetric1(BamgOpts* bamgopts);
-			void IntersectGeomMetric(BamgOpts* bamgopts);
+			void AddGeometryMetric(BamgOpts* bamgopts);
 			int  isCracked() const {return NbCrackedVertices != 0;}
 			int  Crack();
 			int  UnCrack();
 			void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
-			void FillHoleInMesh() ;
+			void GenerateMeshProperties() ;
 			int  CrackMesh();
 
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3230)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3231)
@@ -27,5 +27,5 @@
 		ReadMesh(bamgmesh,bamgopts);
 		SetIntCoor();
-		FillHoleInMesh();
+		GenerateMeshProperties();
 	}
 	/*}}}1*/
@@ -36,5 +36,5 @@
 		ReadMesh(index,x,y,nods,nels);
 		SetIntCoor();
-		FillHoleInMesh();
+		GenerateMeshProperties();
 	}
 	/*}}}1*/
@@ -46,5 +46,5 @@
 		  int * kk    = new int [Tho.nbv];
 		  Int4 * reft = new Int4[Tho.nbt];
-		  Int4 nbInT =    Tho.ConsRefTriangle(reft);
+		  Int4 nbInT =    Tho.TriangleReferenceList(reft);
 		  Int4 * refv = new Int4[Tho.nbv];
 
@@ -130,5 +130,5 @@
 		  Gh.AfterRead(); 
 		  SetIntCoor();
-		  FillHoleInMesh();
+		  GenerateMeshProperties();
 
 		  if (!NbSubDomains){
@@ -552,5 +552,5 @@
 		//Build reft that holds the number the subdomain number of each triangle
 		Int4* reft = new Int4[nbt];
-		Int4 nbInT = ConsRefTriangle(reft);
+		Int4 nbInT = TriangleReferenceList(reft);
 
 		//Vertices
@@ -862,6 +862,78 @@
 
 	/*Methods*/
-	/*FUNCTION Triangles::Add{{{1*/
-	void Triangles::Add( Vertex &s,Triangle* t, Icoor2* det3) {
+	/*FUNCTION Triangles::AddGeometryMetric{{{1*/
+	void Triangles::AddGeometryMetric(BamgOpts* bamgopts){
+
+		/*Get options*/
+		int    verbosity=bamgopts->verbose;
+		double anisomax =bamgopts->anisomax;
+		double errg     =bamgopts->errg;
+
+		Real8 ss[2]={0.00001,0.99999};
+		Real8 errC = 2*sqrt(2*errg);
+		Real8 hmax = Gh.MaximalHmax();
+		Real8 hmin = Gh.MinimalHmin();
+
+		//check that hmax is positive
+		if (hmax<=0){
+			throw ErrorException(__FUNCT__,exprintf("hmax<=0"));
+		}
+
+		//errC cannot be higher than 1
+		if (errC > 1) errC = 1;
+
+		//Set all vertices to "on"
+		SetVertexFieldOn();
+
+		//loop over all the vertices on edges
+		for (Int4  i=0;i<nbe;i++){
+			for (int j=0;j<2;j++){
+
+				Vertex V;
+				VertexOnGeom GV;
+				Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
+
+				GeometricalEdge * eg = GV;
+				Real8 s = GV;
+				R2 tg;
+				Real8  R1= eg->R1tg(s,tg);
+				Real8  ht=hmax;
+				// err relative to the length of the edge
+				if (R1>1.0e-20) {  
+					ht = Min(Max(errC/R1,hmin),hmax);
+				}
+				Real8 hn=Min(hmax,ht*anisomax);
+				if (ht<=0 || hn<=0){
+					throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
+				}
+				MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
+				Metric MVp(Vp);
+				edges[i][j].m.IntersectWith(MVp);
+			}
+		}
+		// the problem is for the vertex on vertex 
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::AddMetric{{{1*/
+	void Triangles::AddMetric(BamgOpts* bamgopts){
+		//  Hessiantype = 0 =>  H is computed using double P2 projection
+		//  Hessiantype = 1 =>  H is computed with green formula
+
+		/*Options*/
+		int Hessiantype=bamgopts->Hessiantype;
+
+		if (Hessiantype==0){
+			BuildMetric0(bamgopts);
+		}
+		else if (Hessiantype==1){
+			BuildMetric1(bamgopts);
+		}
+		else{
+			throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (1->use Green formula, 0-> double P2 projection)",Hessiantype));
+		}
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::AddVertex{{{1*/
+	void Triangles::AddVertex( Vertex &s,Triangle* t, Icoor2* det3) {
 		// -------------------------------------------
 		//             s2
@@ -935,5 +1007,5 @@
 					if ((( Triangle *) ta)->det < 0 ) {
 						// add in outside triangle 
-						Add(s,( Triangle *) ta);
+						AddVertex(s,( Triangle *) ta);
 						return;
 					}
@@ -1978,45 +2050,4 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::ConsRefTriangle{{{1*/
-	Int4  Triangles::ConsRefTriangle(Int4* reft) const {
-		long int verbosity=0;
-		register Triangle *t0,*t;
-		register Int4 k=0, num;   
-
-		//initialize all triangles as -1 (outside)
-		for (Int4 it=0;it<nbt;it++) reft[it]=-1;
-
-		//loop over all subdomains
-		for (Int4 i=0;i<NbSubDomains;i++){ 
-
-			//first triangle of the subdomain i
-			t=t0=subdomains[i].head;
-
-			//check that the subdomain is not empty
-			if (!t0){ // no empty sub domai{
-				throw ErrorException(__FUNCT__,exprintf("At least one subdomain is empty"));
-			}
-
-			//loop
-			do{
-				k++;
-
-				//get current triangle number
-				num = Number(t);
-
-				//check that num is in [0 nbt[
-				if (num<0 || num>=nbt){
-					throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt"));
-				}
-
-				//reft of this triangle is the subdomain number
-				reft[num]=i;
-			} while (t0 != (t=t->link));
-			//stop when all triangles of subdomains have been tagged
-
-			}
-			return k;   
-		}
-		/*}}}1*/
 		/*FUNCTION Triangles::Crack{{{1*/
 		int  Triangles::Crack() { 
@@ -2215,291 +2246,4 @@
 			if (verbosity > 3) printf("      number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 
 		}
-	/*}}}1*/
-	/*FUNCTION Triangles::FillHoleInMesh{{{1*/
-	void Triangles::FillHoleInMesh() {
-
-		int verbosity=0;
-
-		// generation of the integer coordinate
-		  {
-
-			// find extrema coordinates of vertices pmin,pmax
-			Int4 i;
-			if(verbosity>2) printf("      Filling holes in mesh of %i vertices\n",nbv); 
-
-			//initialize ordre
-			if (!ordre){
-				throw ErrorException(__FUNCT__,exprintf("!ordre"));
-			}
-			for (i=0;i<nbv;i++) ordre[i]=0;
-
-			NbSubDomains =0;
-
-			/* generation of the adjacence of the triangles*/
-
-			SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
-
-			//initialize st
-			Int4* st = new Int4[nbt*3];
-			for (i=0;i<nbt*3;i++) st[i]=-1;
-
-			//check number of edges
-			Int4 kk =0;
-			for (i=0;i<nbe;i++){
-				kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
-			}
-			if (kk != nbe) { 
-				throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe));
-			}
-
-			//
-			for (i=0;i<nbt;i++){
-				for (int j=0;j<3;j++) {
-					Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
-								Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
-					Int4 invisible = triangles[i].Hidden(j);
-					if(st[k]==-1){
-						st[k]=3*i+j;
-					}
-					else if(st[k]>=0) {
-						if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){
-							throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || 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);
-						}
-						st[k]=-2-st[k]; 
-					}
-					else {
-						throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",
-										Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));
-					}
-				}
-			}
-			if(verbosity>5) {
-				printf("         info on 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 consistant of edge[].adj and the geometrical required  vertex
-			Int4 k=0;
-			for (i=0;i<edge4->nb();i++){
-				if (st[i] >=0){ // edge alone 
-					if (i < nbe) {
-						Int4 i0=edge4->i(i);
-						ordre[i0] = vertices+i0;
-						Int4 i1=edge4->j(i);
-						ordre[i1] = vertices+i1;
-					}
-					else {
-						k=k+1;
-						if (k <20) {
-							throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));
-						}
-					}
-				}
-			}
-			if(k != 0) {
-				throw ErrorException(__FUNCT__,exprintf("%i boundary edges  are not defined as edges",k));
-			}
-
-			/* mesh generation with boundary points*/
-			Int4 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;
-			Int4 savenbt=nbt;
-			Int4 savenbtx=nbtx;
-			SubDomain * savesubdomains = subdomains;
-			subdomains = 0;
-
-			Int4  Nbtriafillhole = 2*nbvb;
-			Triangle* triafillhole =new Triangle[Nbtriafillhole];
-			triangles =  triafillhole;
-
-			nbt=2;
-			nbtx= Nbtriafillhole;
-
-			for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 
-			 if  ( ++i >= nbvb) {
-				 throw ErrorException(__FUNCT__,exprintf("FillHoleInMesh: All the vertices are aligned"));
-			 }
-			Exchange( ordre[2], ordre[i]);
-
-			Vertex *  v0=ordre[0], *v1=ordre[1];
-
-
-			triangles[0](0) = 0; // sommet pour infini 
-			triangles[0](1) = v0;
-			triangles[0](2) = v1;
-
-			triangles[1](0) = 0;// sommet pour infini 
-			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;  // faux triangles
-			triangles[1].det = -1;  // faux 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);
-
-			// We add the vertices one by one
-			Int4 NbSwap=0;
-			for (Int4 icount=2; icount<nbvb; icount++) {
-				Vertex *vi  = ordre[icount];
-				Icoor2 dete[3];
-				Triangle *tcvi = FindTriangleContaining(vi->i,dete);
-				quadtree->Add(*vi); 
-				Add(*vi,tcvi,dete);
-				NbSwap += vi->Optim(1,1);
-			}
-
-			// inforce the boundary 
-			TriangleAdjacent ta(0,0);
-			Int4 nbloss = 0,knbe=0;
-			for ( i = 0; i < nbe; i++){
-				if (st[i] >=0){  // edge alone => on border ...  FH oct 2009
-					Vertex & a=edges[i][0], & b =    edges[i][1];
-					if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
-					  {
-						knbe++;
-						if (ForceEdge(a,b,ta)<0)
-						 nbloss++;
-					  }
-				}
-			}
-			if(nbloss) {
-				throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe));
-			}
-
-			FindSubDomain(1);
-			// remove all the hole 
-			// remove all the good sub domain
-			Int4 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(! tta.link) // edge between remove and not remove 
-						  { // change the link of ta;
-							int ja = ta;
-							Vertex *v0= ta.EdgeVertex(0);
-							Vertex *v1= ta.EdgeVertex(1);
-							Int4 k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
-							if (st[k]<0){
-								throw ErrorException(__FUNCT__,exprintf("st[k]<0"));
-							}
-							tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
-							ta.SetLock();
-							st[k]=-2-st[k]; 
-						  }
-					  }
-				}
-			}
-			Int4 NbTfillHoll =0;
-			for (i=0;i<nbt;i++){
-				if (triangles[i].link) {
-					triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
-					triangles[i].color=-1;
-				}
-				else
-				  {
-					triangles[i].color= savenbt+ NbTfillHoll++;
-				  }
-			}
-			if (savenbt+NbTfillHoll>savenbtx ){
-				throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx"));
-			}
-			// copy of the outside triangles in saveTriangles 
-			for (i=0;i<nbt;i++){
-				if(triangles[i].color>=0) {
-					savetriangles[savenbt]=triangles[i];
-					savetriangles[savenbt].link=0;
-					savenbt++;
-				}
-			}
-			// gestion of the adj
-			k =0;
-			Triangle * tmax = triangles + nbt;
-			for (i=0;i<savenbt;i++)  
-			  { 
-				Triangle & ti = savetriangles[i];
-				for (int j=0;j<3;j++)
-				  {
-					Triangle * ta = ti.TriangleAdj(j);
-					int aa = ti.NuEdgeTriangleAdj(j);
-					int lck = ti.Locked(j);
-					if (!ta) k++; // bug 
-					else if ( ta >= triangles && ta < tmax) 
-					  {
-						ta= savetriangles + ta->color;
-						ti.SetAdj2(j,ta,aa);
-						if(lck) ti.SetLocked(j);
-					  }
-				  }
-			  }
-			//	 OutSidesTriangles = triangles;
-			//	Int4 NbOutSidesTriangles = nbt;
-
-			// restore triangles;
-			nbt=savenbt;
-			nbtx=savenbtx;
-			delete [] triangles;
-			delete [] subdomains;
-			triangles = savetriangles;
-			subdomains = savesubdomains;
-			if (k) {
-				throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));
-			}
-			FindSubDomain();
-
-			delete edge4;
-			delete [] st;
-			for (i=0;i<nbv;i++)
-			 quadtree->Add(vertices[i]);
-
-			SetVertexFieldOn();
-
-			for (i=0;i<nbe;i++)
-			 if(edges[i].onGeometry) 
-			  for(int j=0;j<2;j++)
-				if (!edges[i].adj[j])
-				 if(!edges[i][j].onGeometry->IsRequiredVertex()) {
-					 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)"));
-				 }
-		  }
-	}
 	/*}}}1*/
 	/*FUNCTION Triangles::FindSubDomain{{{1*/
@@ -2846,4 +2590,289 @@
 	}
 	/*}}}1*/
+	/*FUNCTION Triangles::GenerateMeshProperties{{{1*/
+	void Triangles::GenerateMeshProperties() {
+
+		int verbosity=0;
+
+		// generation of the integer coordinate
+
+		// find extrema coordinates of vertices pmin,pmax
+		Int4 i;
+		if(verbosity>2) printf("      Filling holes in mesh of %i vertices\n",nbv); 
+
+		//initialize ordre
+		if (!ordre){
+			throw ErrorException(__FUNCT__,exprintf("!ordre"));
+		}
+		for (i=0;i<nbv;i++) ordre[i]=0;
+
+		NbSubDomains =0;
+
+		/* generation of the adjacence of the triangles*/
+
+		SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
+
+		//initialize st
+		Int4* st = new Int4[nbt*3];
+		for (i=0;i<nbt*3;i++) st[i]=-1;
+
+		//check number of edges
+		Int4 kk =0;
+		for (i=0;i<nbe;i++){
+			kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
+		}
+		if (kk != nbe) { 
+			throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe));
+		}
+
+		//
+		for (i=0;i<nbt;i++){
+			for (int j=0;j<3;j++) {
+				Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
+							Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
+				Int4 invisible = triangles[i].Hidden(j);
+				if(st[k]==-1){
+					st[k]=3*i+j;
+				}
+				else if(st[k]>=0) {
+					if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){
+						throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || 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);
+					}
+					st[k]=-2-st[k]; 
+				}
+				else {
+					throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",
+									Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));
+				}
+			}
+		}
+		if(verbosity>5) {
+			printf("         info on 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 consistant of edge[].adj and the geometrical required  vertex
+		Int4 k=0;
+		for (i=0;i<edge4->nb();i++){
+			if (st[i] >=0){ // edge alone 
+				if (i < nbe) {
+					Int4 i0=edge4->i(i);
+					ordre[i0] = vertices+i0;
+					Int4 i1=edge4->j(i);
+					ordre[i1] = vertices+i1;
+				}
+				else {
+					k=k+1;
+					if (k <20) {
+						throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));
+					}
+				}
+			}
+		}
+		if(k != 0) {
+			throw ErrorException(__FUNCT__,exprintf("%i boundary edges  are not defined as edges",k));
+		}
+
+		/* mesh generation with boundary points*/
+		Int4 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;
+		Int4 savenbt=nbt;
+		Int4 savenbtx=nbtx;
+		SubDomain * savesubdomains = subdomains;
+		subdomains = 0;
+
+		Int4  Nbtriafillhole = 2*nbvb;
+		Triangle* triafillhole =new Triangle[Nbtriafillhole];
+		triangles =  triafillhole;
+
+		nbt=2;
+		nbtx= Nbtriafillhole;
+
+		for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 
+		 if  ( ++i >= nbvb) {
+			 throw ErrorException(__FUNCT__,exprintf("GenerateMeshProperties: All the vertices are aligned"));
+		 }
+		Exchange( ordre[2], ordre[i]);
+
+		Vertex *  v0=ordre[0], *v1=ordre[1];
+
+
+		triangles[0](0) = 0; // sommet pour infini 
+		triangles[0](1) = v0;
+		triangles[0](2) = v1;
+
+		triangles[1](0) = 0;// sommet pour infini 
+		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;  // faux triangles
+		triangles[1].det = -1;  // faux 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);
+
+		// We add the vertices one by one
+		Int4 NbSwap=0;
+		for (Int4 icount=2; icount<nbvb; icount++) {
+			Vertex *vi  = ordre[icount];
+			Icoor2 dete[3];
+			Triangle *tcvi = FindTriangleContaining(vi->i,dete);
+			quadtree->Add(*vi); 
+			AddVertex(*vi,tcvi,dete);
+			NbSwap += vi->Optim(1,1);
+		}
+
+		// inforce the boundary 
+		TriangleAdjacent ta(0,0);
+		Int4 nbloss = 0,knbe=0;
+		for ( i = 0; i < nbe; i++){
+			if (st[i] >=0){  // edge alone => on border ...  FH oct 2009
+				Vertex & a=edges[i][0], & b =    edges[i][1];
+				if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
+				  {
+					knbe++;
+					if (ForceEdge(a,b,ta)<0)
+					 nbloss++;
+				  }
+			}
+		}
+		if(nbloss) {
+			throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe));
+		}
+
+		FindSubDomain(1);
+		// remove all the hole 
+		// remove all the good sub domain
+		Int4 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(! tta.link) // edge between remove and not remove 
+					  { // change the link of ta;
+						int ja = ta;
+						Vertex *v0= ta.EdgeVertex(0);
+						Vertex *v1= ta.EdgeVertex(1);
+						Int4 k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
+						if (st[k]<0){
+							throw ErrorException(__FUNCT__,exprintf("st[k]<0"));
+						}
+						tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
+						ta.SetLock();
+						st[k]=-2-st[k]; 
+					  }
+				  }
+			}
+		}
+		Int4 NbTfillHoll =0;
+		for (i=0;i<nbt;i++){
+			if (triangles[i].link) {
+				triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
+				triangles[i].color=-1;
+			}
+			else
+			  {
+				triangles[i].color= savenbt+ NbTfillHoll++;
+			  }
+		}
+		if (savenbt+NbTfillHoll>savenbtx ){
+			throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx"));
+		}
+		// copy of the outside triangles in saveTriangles 
+		for (i=0;i<nbt;i++){
+			if(triangles[i].color>=0) {
+				savetriangles[savenbt]=triangles[i];
+				savetriangles[savenbt].link=0;
+				savenbt++;
+			}
+		}
+		// gestion of the adj
+		k =0;
+		Triangle * tmax = triangles + nbt;
+		for (i=0;i<savenbt;i++)  
+		  { 
+			Triangle & ti = savetriangles[i];
+			for (int j=0;j<3;j++)
+			  {
+				Triangle * ta = ti.TriangleAdj(j);
+				int aa = ti.NuEdgeTriangleAdj(j);
+				int lck = ti.Locked(j);
+				if (!ta) k++; // bug 
+				else if ( ta >= triangles && ta < tmax) 
+				  {
+					ta= savetriangles + ta->color;
+					ti.SetAdj2(j,ta,aa);
+					if(lck) ti.SetLocked(j);
+				  }
+			  }
+		  }
+		//	 OutSidesTriangles = triangles;
+		//	Int4 NbOutSidesTriangles = nbt;
+
+		// restore triangles;
+		nbt=savenbt;
+		nbtx=savenbtx;
+		delete [] triangles;
+		delete [] subdomains;
+		triangles = savetriangles;
+		subdomains = savesubdomains;
+		if (k) {
+			throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));
+		}
+		FindSubDomain();
+
+		delete edge4;
+		delete [] st;
+		for (i=0;i<nbv;i++)
+		 quadtree->Add(vertices[i]);
+
+		SetVertexFieldOn();
+
+		for (i=0;i<nbe;i++)
+		 if(edges[i].onGeometry) 
+		  for(int j=0;j<2;j++)
+			if (!edges[i].adj[j])
+			 if(!edges[i][j].onGeometry->IsRequiredVertex()) {
+				 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)"));
+			 }
+	}
+	/*}}}1*/
 	/*FUNCTION Triangles::GeomToTriangles0{{{1*/
 	void Triangles::GeomToTriangles0(Int4 inbvx,BamgOpts* bamgopts) {
@@ -3503,76 +3532,4 @@
 	}
 	/*}}}1*/
-/*FUNCTION Triangles::IntersectConsMetric{{{1*/
-void Triangles::IntersectConsMetric(BamgOpts* bamgopts){
-	//  Hessiantype = 0 =>  H is computed using double P2 projection
-	//  Hessiantype = 1 =>  H is computed with green formula
-
-	/*Options*/
-	int Hessiantype=bamgopts->Hessiantype;
-
-	if (Hessiantype==0){
-		BuildMetric0(bamgopts);
-	}
-	else if (Hessiantype==1){
-		BuildMetric1(bamgopts);
-	}
-	else{
-		throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (1->use Green formula, 0-> double P2 projection)",Hessiantype));
-	}
-}
-/*}}}1*/
-/*FUNCTION Triangles::IntersectGeomMetric{{{1*/
-void Triangles::IntersectGeomMetric(BamgOpts* bamgopts){
-
-	/*Get options*/
-	int    verbosity=bamgopts->verbose;
-	double anisomax =bamgopts->anisomax;
-	double errg     =bamgopts->errg;
-
-	Real8 ss[2]={0.00001,0.99999};
-	Real8 errC = 2*sqrt(2*errg);
-	Real8 hmax = Gh.MaximalHmax();
-	Real8 hmin = Gh.MinimalHmin();
-
-	//check that hmax is positive
-	if (hmax<=0){
-		throw ErrorException(__FUNCT__,exprintf("hmax<=0"));
-	}
-
-	//errC cannot be higher than 1
-	if (errC > 1) errC = 1;
-
-	//Set all vertices to "on"
-	SetVertexFieldOn();
-
-	//loop over all the vertices on edges
-	for (Int4  i=0;i<nbe;i++){
-		for (int j=0;j<2;j++){
-
-			Vertex V;
-			VertexOnGeom GV;
-			Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
-
-			GeometricalEdge * eg = GV;
-			Real8 s = GV;
-			R2 tg;
-			Real8  R1= eg->R1tg(s,tg);
-			Real8  ht=hmax;
-			// err relative to the length of the edge
-			if (R1>1.0e-20) {  
-				ht = Min(Max(errC/R1,hmin),hmax);
-			}
-			Real8 hn=Min(hmax,ht*anisomax);
-			if (ht<=0 || hn<=0){
-				throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
-			}
-			MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
-			Metric MVp(Vp);
-			edges[i][j].m.IntersectWith(MVp);
-		}
-	}
-	// the problem is for the vertex on vertex 
-}
-/*}}}1*/
 	/*FUNCTION Triangles::Insert{{{1*/
 	void Triangles::Insert() {
@@ -3692,5 +3649,5 @@
 
 			//Add newvertex to the existing mesh
-			Add(*newvertex,tcvi,dete);
+			AddVertex(*newvertex,tcvi,dete);
 
 			//Make the mesh Delaunay around newvertex by swaping the edges
@@ -3782,5 +3739,5 @@
 				}
 				quadtree->Add(vj);
-				Add(vj,tcvj,dete);
+				AddVertex(vj,tcvj,dete);
 				NbSwap += vj.Optim(1);          
 				iv++;
@@ -5262,5 +5219,5 @@
 		//  ReMakeTriangleContainingTheVertex();
 
-		FillHoleInMesh();
+		GenerateMeshProperties();
 
 		if (verbosity>2){
@@ -5339,5 +5296,5 @@
 				throw ErrorException(__FUNCT__,exprintf("!tcvi || tcvi->det < 0"));
 			}
-			Add(vi,tcvi,dete);
+			AddVertex(vi,tcvi,dete);
 			NbSwap += vi.Optim(1);          
 			iv++;
@@ -5354,4 +5311,42 @@
 
 return  NbSplitEdge;
+}
+/*}}}1*/
+/*FUNCTION Triangles::TriangleReferenceList{{{1*/
+Int4  Triangles::TriangleReferenceList(Int4* reft) const {
+	long int verbosity=0;
+	register Triangle *t0,*t;
+	register Int4 k=0, num;   
+
+	//initialize all triangles as -1 (outside)
+	for (Int4 it=0;it<nbt;it++) reft[it]=-1;
+
+	//loop over all subdomains
+	for (Int4 i=0;i<NbSubDomains;i++){ 
+
+		//first triangle of the subdomain i
+		t=t0=subdomains[i].head;
+
+		//check that the subdomain is not empty
+		if (!t0){ throw ErrorException(__FUNCT__,exprintf("At least one subdomain is empty"));}
+
+		//loop
+		do{
+			k++;
+
+			//get current triangle number
+			num = Number(t);
+
+			//check that num is in [0 nbt[
+			if (num<0 || num>=nbt){ throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt"));}
+
+			//reft of this triangle is the subdomain number
+			reft[num]=i;
+
+		} while (t0 != (t=t->link));
+		//stop when all triangles of subdomains have been tagged
+
+	}
+	return k;   
 }
 /*}}}1*/