Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2939)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2940)
@@ -138,4 +138,5 @@
 		Real8  Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,Real8 =1);
 		void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,BamgOpts* bamgopts);
+		void Echo();
 		int ref() const { return ReferenceNumber;}
 
@@ -304,7 +305,7 @@
 
 		private: // les arete sont opposes a un sommet
-		Vertex * ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
-		Triangle * at [3]; // nu triangle adjacent  
-		Int1  aa[3];  // les nu des arete dans le triangles (mod 4)
+		Vertex* ns [3];    // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
+		Triangle* at [3];  // 3 adjacent triangles (nu)
+		Int1  aa[3];       // les nu des arete dans le triangles (mod 4)
 		public: 
 		Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
@@ -313,4 +314,5 @@
 			Int4 color;
 		};
+		void Echo();
 		void SetDet() {
 			if(ns[0] && ns[1] && ns[2])    det = bamg::det(*ns[0],*ns[1],*ns[2]);
@@ -323,18 +325,15 @@
 		TriangleAdjacent FindBoundaryEdge(int ) const;
 
-		void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu) 
-		  {
+		void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
 			if (link  >=tb && link  <te) link  = tb + renu[link -tb];
 			if (at[0] >=tb && at[0] <te) at[0] = tb + renu[at[0]-tb];
 			if (at[1] >=tb && at[1] <te) at[1] = tb + renu[at[1]-tb];
 			if (at[2] >=tb && at[2] <te) at[2] = tb + renu[at[2]-tb];    
-		  }
-		void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu) 
-		  {
+		}
+		void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
 			if (ns[0] >=vb && ns[0] <ve) ns[0] = vb + renu[ns[0]-vb];
 			if (ns[1] >=vb && ns[1] <ve) ns[1] = vb + renu[ns[1]-vb];
 			if (ns[2] >=vb && ns[2] <ve) ns[2] = vb + renu[ns[2]-vb];    
-		  }
-
+		}
 
 		const Vertex & operator[](int i) const {return *ns[i];};
@@ -353,5 +352,4 @@
 
 		inline Real4 qualite() ;
-
 
 		void SetAdjAdj(Int1 a) 
@@ -390,28 +388,33 @@
 			register Triangle * t = at[a];
 			if(t) t->aa[aa[a] & 3] |=16;
-			aa[a] |= 16;}
-			void SetCracked(int a){
-				register Triangle * t = at[a];
-				if(t) t->aa[aa[a] & 3] |=32;
-				aa[a] |= 32;}
-
-				double   QualityQuad(int a,int option=1) const;
-				Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
-
-				void SetLocked(int a){
-					register Triangle * t = at[a];
-					t->aa[aa[a] & 3] |=4;
-					aa[a] |= 4;}
-
-					void SetMarkUnSwap(int a){
-						register Triangle * t = at[a];
-						t->aa[aa[a] & 3] |=8;
-						aa[a] |=8 ;}
-
-
-						void SetUnMarkUnSwap(int a){ 
-							register Triangle * t = at[a];
-							t->aa[aa[a] & 3] &=55; // 23 + 32 
-							aa[a] &=55 ;}
+			aa[a] |= 16;
+		}
+		void SetCracked(int a){
+			register Triangle * t = at[a];
+			if(t) t->aa[aa[a] & 3] |=32;
+			aa[a] |= 32;
+		}
+
+		double   QualityQuad(int a,int option=1) const;
+		Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
+
+		void SetLocked(int a){
+			register Triangle * t = at[a];
+			t->aa[aa[a] & 3] |=4;
+			aa[a] |= 4;
+		}
+
+		void SetMarkUnSwap(int a){
+			register Triangle * t = at[a];
+			t->aa[aa[a] & 3] |=8;
+			aa[a] |=8 ;
+		}
+
+
+		void SetUnMarkUnSwap(int a){ 
+			register Triangle * t = at[a];
+			t->aa[aa[a] & 3] &=55; // 23 + 32 
+			aa[a] &=55 ;
+		}
 
 	};  // end of Triangle class 
@@ -799,5 +802,4 @@
 						 void BuildMetric0(BamgOpts* bamgopts);
 						 void BuildMetric1(BamgOpts* bamgopts);
-						 void BuildMetric2(BamgOpts* bamgopts);
 						 void IntersectGeomMetric(const Real8 err,const int iso);
 
Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 2939)
+++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 2940)
@@ -44,4 +44,36 @@
 	  }
 	/*}}}1*/
+	/*FUNCTION Triangle::Echo {{{1*/
+
+	void Triangle::Echo(void){
+
+		int i;
+
+		printf("Triangle:\n");
+		printf("   ns pointer towards three vertices\n");
+		printf("      ns[0] ns[1] ns[2] = %p %p %p\n",ns[0],ns[1],ns[2]);
+		printf("   at pointer towards three adjacent triangles\n");
+		printf("      at[0] at[1] at[2] = %p %p %p\n",at[0],at[1],at[2]);
+		printf("   det (integer triangle determinant) = %i\n",det);
+		if (link){
+			printf("   link (pointer toward duplicate triangle)= %p\n",link);
+		}
+		else{
+			printf("   color = %i\n",color);
+		}
+
+		printf("\nThree vertices:\n");
+		for(i=0;i<3;i++){
+			if (ns[i]){
+				ns[i]->Echo();
+			}
+			else{
+				printf("   vertex %i does not exist\n",i+1);
+			}
+		}
+
+		return;
+	}
+	/*}}}*/
 
 }
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2939)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2940)
@@ -891,10 +891,205 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::BuildMetric0 (Green formula){{{1*/
+	/*FUNCTION Triangles::BuildMetric0 (double P2 projection){{{1*/
 	void Triangles::BuildMetric0(BamgOpts* bamgopts){
-		//  the array of solution s is store    
-		// sol0,sol1,...,soln    on vertex 0
-		//  sol0,sol1,...,soln   on vertex 1
-		//  etc.
+
+		/*Options*/
+		const int dim = 2;
+		double* s=NULL;
+		Int4    nbsol;
+		int     verbosity;
+
+		int   i,j,k,iA,iB,iC;
+		int   iv;
+
+		/*Recover options*/
+		verbosity=bamgopts->verbose;
+
+		/*Get and process fields*/
+		s=bamgopts->field;
+		nbsol=bamgopts->numfields;
+
+		//initialization of some variables
+		double* ss=(double*)s;
+		double  sA,sB,sC;
+		Real8*  detT = new Real8[nbt];
+		Real8*  sumareas = new Real8[nbv];
+		Real8*  alpha= new Real8[nbt*3];
+		Real8*  beta = new Real8[nbt*3];
+		Real8*  dx_elem    = new Real8[nbt];
+		Real8*  dy_elem    = new Real8[nbt];
+		Real8*  dx_vertex  = new Real8[nbv];
+		Real8*  dy_vertex  = new Real8[nbv];
+		Real8*  dxdx_elem  = new Real8[nbt];
+		Real8*  dxdy_elem  = new Real8[nbt];
+		Real8*  dydy_elem  = new Real8[nbt];
+		Real8*  dxdx_vertex= new Real8[nbv];
+		Real8*  dxdy_vertex= new Real8[nbv];
+		Real8*  dydy_vertex= new Real8[nbv];
+
+		//display infos
+		if(verbosity>1) {
+			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv);
+		}
+
+		//first, build the chains that will be used for the Hessian computation, as weel as the area of each element
+		int head_s[nbv];
+		int next_p[3*nbt];
+		int  p=0;
+		//initialization
+		for(i=0;i<nbv;i++){
+			sumareas[i]=0;
+			head_s[i]=-1;
+		}
+		for(i=0;i<nbt;i++){
+
+			//lopp over the real triangles (no boundary elements)
+			if(triangles[i].link){ 
+
+				//get current triangle t
+				const Triangle &t=triangles[i];
+
+				// coor of 3 vertices 
+				R2 A=t[0];
+				R2 B=t[1];
+				R2 C=t[2];
+
+				//compute triangle determinant (2*Area)
+				Real8 dett = bamg::Area2(A,B,C);
+				detT[i]=dett;
+
+				/*The nodal functions are such that for a vertex A:
+				 *    N_A(x,y)=alphaA x + beta_A y +gamma_A
+				 *    N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
+				 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
+				 * leads to:
+				 *    N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
+				 * and this gives:
+				 *    alpha_A = (yB-yC)/(2*Area(T))*/
+				alpha[i*3+0]=(B.y-C.y)/dett;
+				alpha[i*3+1]=(C.y-A.y)/dett;
+				alpha[i*3+2]=(A.y-B.y)/dett;
+				beta[ i*3+0]=(C.x-B.x)/dett;
+				beta[ i*3+1]=(A.x-C.x)/dett;
+				beta[ i*3+2]=(B.x-A.x)/dett;
+
+				//compute chains
+				for(j=0;j<3;j++){
+					k=Number(triangles[i][j]);
+					next_p[p]=head_s[k];
+					head_s[k]=p++;
+
+					//add area to sumareas
+					sumareas[k]+=dett;
+				}
+
+			}
+		}
+
+		//for all Solutions
+		for (Int4 nusol=0;nusol<nbsol;nusol++) {
+			Real8 smin=ss[nusol],smax=ss[nusol];
+
+			//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
+			for ( iv=0,k=0; iv<nbv; iv++){
+				smin=Min(smin,ss[iv*nbsol+nusol]);
+				smax=Max(smax,ss[iv*nbsol+nusol]);
+			}
+			Real8 sdelta=smax-smin;
+			Real8 absmax=Max(Abs(smin),Abs(smax));
+
+			//display info
+			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g\n",nusol,smin,smax,sdelta);
+
+			//skip constant field
+			if (sdelta < 1.0e-10*Max(absmax,1e-20)){
+				printf("      Solution %i is constant, skipping...\n",nusol);
+				continue;
+			}
+
+			//initialize the hessian and gradient matrices
+			for ( iv=0,k=0; iv<nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
+
+			//1: Compute gradient for each element (exact)
+			for (i=0;i<nbt;i++){
+				if(triangles[i].link){
+					// number of the 3 vertices
+					iA = Number(triangles[i][0]);
+					iB = Number(triangles[i][1]);
+					iC = Number(triangles[i][2]);
+
+					// value of the P1 fonction on 3 vertices 
+					sA = ss[iA*nbsol+nusol];
+					sB = ss[iB*nbsol+nusol];
+					sC = ss[iC*nbsol+nusol];
+
+					//gradient = (sum alpha_i s_i, sum_i beta_i s_i)
+					dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
+					dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
+				}
+			}
+
+			//2: then compute a gradient for each vertex using a P2 projection
+			for(i=0;i<nbv;i++){
+				for(p=head_s[i];p!=-1;p=next_p[p]){
+					//Get triangle number
+					k=(Int4)(p/3);
+					dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
+					dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
+				}
+			}
+
+			//3: compute Hessian matrix on each element
+			for (i=0;i<nbt;i++){
+				if(triangles[i].link){
+					// number of the 3 vertices
+					iA = Number(triangles[i][0]);
+					iB = Number(triangles[i][1]);
+					iC = Number(triangles[i][2]);
+
+					//Hessian
+					dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
+					dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
+					dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
+				}
+			}
+
+			//4: finaly compute Hessian on each vertex using the second P2 projection
+			for(i=0;i<nbv;i++){
+				for(p=head_s[i];p!=-1;p=next_p[p]){
+					//Get triangle number
+					k=(Int4)(p/3);
+					dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
+					dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
+					dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
+				}
+			}
+
+			/*Compute Metric from Hessian*/
+			for ( iv=0;iv<nbv;iv++){
+				vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts);
+			}
+
+		}//for all solutions
+
+		//clean up
+		delete [] detT;
+		delete [] alpha;
+		delete [] beta;
+		delete [] sumareas;
+		delete [] dx_elem;
+		delete [] dy_elem;
+		delete [] dx_vertex;
+		delete [] dy_vertex;
+		delete [] dxdx_elem;
+		delete [] dxdy_elem;
+		delete [] dydy_elem;
+		delete [] dxdx_vertex;
+		delete [] dxdy_vertex;
+		delete [] dydy_vertex;
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::BuildMetric1 (Green formula){{{1*/
+	void Triangles::BuildMetric1(BamgOpts* bamgopts){
 
 		/*Options*/
@@ -902,5 +1097,4 @@
 		double* s=NULL;
 		Int4 nbsol;
-		int* typsols=NULL;
 		int NbJacobi;
 		int verbosity;
@@ -913,22 +1107,4 @@
 		s=bamgopts->field;
 		nbsol=bamgopts->numfields;
-		typsols=(int*)xmalloc(1*sizeof(int));
-		typsols[0]=0; // only one dof per node
-
-		//sizeoftype = {1,2,3,4}
-		int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 
-
-		// computation of the number of fields
-		Int4 ntmp = 0;
-		if (typsols){
-			//if there is more than one dof per vertex for one field
-			//increase ntmp to take into account all the fields
-			//if nbsol=1
-			for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
-		}
-		else ntmp = nbsol;
-
-		//n is the total number of fields
-		const Int4 n = ntmp;
 
 		//initialization of some variables
@@ -949,5 +1125,5 @@
 		//display infos
 		if(verbosity>1) {
-			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
+			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv);
 		}
 
@@ -1018,31 +1194,13 @@
 		for (Int4 nusol=0;nusol<nbsol;nusol++) {
 
-			Real8 smin=ss[0],smax=ss[0];
+			Real8 smin=ss[nusol],smax=ss[nusol];
 			Real8 h1=1.e30,h2=1e-30,rx=0;
 			Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-			int   nbfield=typsols?sizeoftype[typsols[nusol]]:1; 
-
-			//only one field
-			if (nbfield == 1) {
-				//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
-				for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
-					dxdx[iv]=dxdy[iv]=dydy[iv]=0;
-					smin=Min(smin,ss[k]);
-					smax=Max(smax,ss[k]);
-				}
-			}
-
-			//vectorial case
-			else{
-				for (iv=0,k=0;iv<nbv;iv++,k+=n ){
-					//compute v = √sum(s[i]^2)
-					double v=0;		     
-					for (int i=0;i<nbfield;i++){
-						v += ss[k+i]*ss[k+i];
-					}
-					v = sqrt(v);
-					smin=Min(smin,v);
-					smax=Max(smax,v);
-				}
+
+			//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
+			for ( iv=0,k=0; iv<nbv; iv++ ){
+				dxdx[iv]=dxdy[iv]=dydy[iv]=0;
+				smin=Min(smin,ss[iv*nbsol+nusol]);
+				smax=Max(smax,ss[iv*nbsol+nusol]);
 			}
 			Real8 sdelta=smax-smin;
@@ -1050,8 +1208,8 @@
 
 			//display info
-			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
+			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbsol);
 
 			//skip constant field
-			if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){
+			if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
 				if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
 				continue;
@@ -1061,10 +1219,6 @@
 			double* sf=ss; 
 
-			//loop over all the fields of the solution
-			for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){
-				//ss++ so that for each iteration ss points toward the right field
-
 				//initialize the hessian matrix
-				for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
+				for ( iv=0,k=0; iv<nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
 
 				//loop over the triangles
@@ -1097,7 +1251,7 @@
 
 						// value of the P1 fonction on 3 vertices 
-						sA = ss[iA*n];
-						sB = ss[iB*n];
-						sC = ss[iC*n];
+						sA = ss[iA*nbsol+nusol];
+						sB = ss[iB*nbsol+nusol];
+						sC = ss[iC*nbsol+nusol];
 
 						/*The nodal functions are such that for a vertex A:
@@ -1145,5 +1299,5 @@
 
 				Int4 kk=0;
-				for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
+				for ( iv=0,k=0 ; iv<nbv; iv++){
 					if(Mmassxx[iv]>0){
 						dxdx[iv] /= 2*Mmassxx[iv];
@@ -1207,8 +1361,7 @@
 				/*Compute Metric from Hessian*/
 				for ( iv=0;iv<nbv;iv++){
-					vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts);
-				}
-
-			} //  end of for all field
+					vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts);
+				}
+
 		}// end for all solution 
 
@@ -1223,653 +1376,4 @@
 		delete []  OnBoundary;
 
-	}
-	/*}}}1*/
-	/*FUNCTION Triangles::BuildMetric1 (from P2 on 4T){{{1*/
-	void Triangles::BuildMetric1(BamgOpts* bamgopts){
-		//  the array of solution s is store    
-		// sol0,sol1,...,soln    on vertex 0
-		//  sol0,sol1,...,soln   on vertex 1
-		//  etc.
-
-		/*Options*/
-		const int dim = 2;
-		double* s=NULL;
-		Int4 nbsol;
-		int* typsols=NULL;
-		int NbJacobi;
-		int verbosity;
-
-		/*Recover options*/
-		verbosity=bamgopts->verbose;
-		NbJacobi=bamgopts->nbjacobi;
-
-		/*Get and process fields*/
-		s=bamgopts->field;
-		nbsol=bamgopts->numfields;
-		typsols=(int*)xmalloc(1*sizeof(int));
-		typsols[0]=0; // only one dof per node
-
-		//sizeoftype = {1,2,3,4}
-		int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 
-
-		// computation of the number of fields
-		Int4 ntmp = 0;
-		if (typsols){
-			//if there is more than one dof per vertex for one field
-			//increase ntmp to take into account all the fields
-			//if nbsol=1
-			for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
-		}
-		else ntmp = nbsol;
-
-		//n is the total number of fields
-		const Int4 n = ntmp;
-
-		//initialization of some variables
-		Int4    i,k,iA,iB,iC,iv;
-		R2      O(0,0);
-		double* ss=(double*)s;
-		double  sA,sB,sC;
-		Real8*  detT = new Real8[nbt];
-		Real8*  Mmass= new Real8[nbv];
-		Real8*  Mmassxx= new Real8[nbv];
-		Real8*  dxdx= new Real8[nbv];
-		Real8*  dxdy= new Real8[nbv];
-		Real8*  dydy= new Real8[nbv];
-		Real8*  workT= new Real8[nbt];
-		Real8*  workV= new Real8[nbv];
-		int*    OnBoundary = new int[nbv];
-
-		//display infos
-		if(verbosity>1) {
-			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
-		}
-
-		//initialize Mmass, OnBoundary and Massxx by zero
-		for (iv=0;iv<nbv;iv++){
-			Mmass[iv]=0;
-			OnBoundary[iv]=0;
-			Mmassxx[iv]=0;
-		}
-
-		//Build detT Mmas Mmassxx workT and OnBoundary
-		for (i=0;i<nbt;i++){ 
-
-			//lopp over the real triangles (no boundary elements)
-			if(triangles[i].link){ 
-
-				//get current triangle t
-				const Triangle &t=triangles[i];
-
-				// coor of 3 vertices 
-				R2 A=t[0];
-				R2 B=t[1];
-				R2 C=t[2];
-
-				// number of the 3 vertices
-				iA = Number(t[0]);
-				iB = Number(t[1]);
-				iC = Number(t[2]);
-
-				//compute triangle determinant (2*Area)
-				Real8 dett = bamg::Area2(A,B,C);
-				detT[i]=dett;
-				dett /= 6;
-
-				// construction of OnBoundary (flag=1 if on boundary, else 0)
-				int nbb=0;
-				for(int j=0;j<3;j++){
-					//get adjacent triangle
-					Triangle *ta=t.Adj(j);
-					//if there is no adjacent triangle, the edge of the triangle t is on boundary
-					if ( !ta || !ta->link){
-						//mark the two vertices of the edge as OnBoundary
-						OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1;
-						OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1;
-						nbb++;
-					}
-				}
-
-				//number of vertices on boundary for current triangle t
-				workT[i] = nbb;
-
-				//Build Mmass Mmass[i] = Mmass[i] + Area/3
-				Mmass[iA] += dett;
-				Mmass[iB] += dett;
-				Mmass[iC] += dett;
-
-				//Build Massxx = Mmass
-				if(nbb==0){
-					Mmassxx[iA] += dett;
-					Mmassxx[iB] += dett;
-					Mmassxx[iC] += dett;
-				}
-			}
-
-			//else: the triangle is a boundary triangle -> workT=-1
-			else workT[i]=-1;
-		}
-
-		//for all Solution  
-		for (Int4 nusol=0;nusol<nbsol;nusol++) {
-
-			Real8 smin=ss[0],smax=ss[0];
-			Real8 h1=1.e30,h2=1e-30,rx=0;
-			Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-			int   nbfield=typsols?sizeoftype[typsols[nusol]]:1; 
-
-			//only one field
-			if (nbfield == 1) {
-				//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
-				for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
-					dxdx[iv]=dxdy[iv]=dydy[iv]=0;
-					smin=Min(smin,ss[k]);
-					smax=Max(smax,ss[k]);
-				}
-			}
-
-			//vectorial case
-			else{
-				for (iv=0,k=0;iv<nbv;iv++,k+=n ){
-					//compute v = √sum(s[i]^2)
-					double v=0;		     
-					for (int i=0;i<nbfield;i++){
-						v += ss[k+i]*ss[k+i];
-					}
-					v = sqrt(v);
-					smin=Min(smin,v);
-					smax=Max(smax,v);
-				}
-			}
-			Real8 sdelta=smax-smin;
-			Real8 absmax=Max(Abs(smin),Abs(smax));
-
-			//display info
-			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
-
-			//skip constant field
-			if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){
-				if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
-				continue;
-			}
-
-			//pointer toward ss that is also a pointer toward s (solutions)
-			double* sf=ss; 
-
-			//loop over all the fields of the solution
-			for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){
-				//ss++ so that for each iteration ss points toward the right field
-
-				//initialize the hessian matrix
-				for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
-
-				//loop over the triangles
-				for (i=0;i<nbt;i++){
-
-					//for real all triangles 
-					if(triangles[i].link){
-
-						// coor of 3 vertices 
-						R2 A=triangles[i][0];
-						R2 B=triangles[i][1];
-						R2 C=triangles[i][2];
-
-						//warning: the normal is internal and the size is the length of the edge
-						R2 nAB = Orthogonal(B-A);
-						R2 nBC = Orthogonal(C-B);
-						R2 nCA = Orthogonal(A-C);
-						//note that :  nAB + nBC + nCA == 0 
-
-						// number of the 3 vertices
-						iA = Number(triangles[i][0]);
-						iB = Number(triangles[i][1]);
-						iC = Number(triangles[i][2]);
-
-						// for the test of  boundary edge
-						// the 3 adj triangles 
-						Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
-						Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
-						Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);
-
-						// value of the P1 fonction on 3 vertices 
-						sA = ss[iA*n];
-						sB = ss[iB*n];
-						sC = ss[iC*n];
-
-						/*The nodal functions are such that for a vertex A:
-						  N_A(x,y)=alphaA x + beta_A y +gamma_A
-						  N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
-						  solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
-						  leads to:
-						  N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
-						  and this gives:
-						  alpha_A = (yB-yC)/(2*Area(T))
-						  beta_A = (xC-xB)/(2*Area(T))
-						  and therefore:
-						  grad N_A = nA / detT
-						  for an interpolation of a solution s:
-						  grad(s) = s * sum_{i=A,B,C} grad(N_i) */
-
-						R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i];
-
-						//from P2 on 4T to compute Hessian
-
-						int   nbb=0;
-						Real8 dd = detT[i];
-						Real8 taa[3][3],bb[3];
-
-						// construction of the transpose of lin system
-						for (int j=0;j<3;j++){
-							int              ie = OppositeEdge[j];
-							TriangleAdjacent ta = triangles[i].Adj(ie);
-							Triangle*        tt = ta;
-
-
-							/* V is the vertex opposed to the edge shared by the
-							 * current triangle i and its neighbor ta
-							 *         C       
-							 *       / | \     
-							 *      /  |  \    
-							 *     /   |   \   
-							 * V  / ta | i  \ B
-							 *    \    |    /  
-							 *     \   |   /   
-							 *      \  |  /
-							 *       \ | /     
-							 *         A
-							 * lA=area(VBC)/Area(ABC)
-							 * lB=area(AVC)/Area(ABC)
-							 * lC=area(ABV)/Area(ABC)
-							 */
-
-							//first, get V
-							Vertex &v = *ta.OppositeVertex();
-							Int4   iV = Number(v);
-
-							//if the adjacent triangle is not a boundary triangle:
-							if (tt && tt->link){
-								Vertex &v = *ta.OppositeVertex();
-								R2      V = v;
-								Int4   iV = Number(v);
-
-								//get lA lB and lC
-								Real8  lA = bamg::Area2(V,B,C)/dd;
-								Real8  lB = bamg::Area2(A,V,C)/dd;
-								Real8  lC = bamg::Area2(A,B,V)/dd;
-
-								//fill taa
-								taa[0][j] = lB*lC;
-								taa[1][j] = lC*lA;
-								taa[2][j] = lA*lB;
-
-								//fill bb
-								bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ; //value of the solution on V minus...
-							}
-							else{
-								//there might be a problem here: if 2 nodes are inline lA, lB or lC
-								//is 0. This means that one column has only one term.
-								//if instead of 0.1 we use 0, and taa[j][j]=1, we might have det(taa)=0
-								nbb++;
-								taa[0][j]=0.1;
-								taa[1][j]=0.1;
-								taa[2][j]=0.1;
-								taa[j][j]=1;
-								bb[j]=0;
-							}
-						}
-
-						// resolution of 3x3 linear system transpose
-						Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);	
-						Real8 cBC   =  det3x3(bb,taa[1],taa[2]);
-						Real8 cCA   =  det3x3(taa[0],bb,taa[2]);
-						Real8 cAB   =  det3x3(taa[0],taa[1],bb);
-						if (!det33){
-							throw ErrorException(__FUNCT__,exprintf("for triangle %i, det33==0! cannot compute theHessian matrix using Hessiantype=1",i+1));
-						}
-
-						// computation of the Hessian in the element 
-						// H( li*lj) = grad li grad lj + grad lj grad lj
-						// grad li = njk  / detT ; with i j k =(A,B,C)
-						Real8 Hxx = cAB * ( nBC.x*nCA.x) +  cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);
-						Real8 Hyy = cAB * ( nBC.y*nCA.y) +  cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);
-						Real8 Hxy = cAB * ( nBC.y*nCA.x) +  cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 
-						  + cAB * ( nBC.x*nCA.y) +  cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);
-						if(nbb==0){
-							dxdx[iA] += Hxx;
-							dydy[iA] += Hyy;
-							dxdy[iA] += Hxy;
-
-							dxdx[iB] += Hxx;
-							dydy[iB] += Hyy;
-							dxdy[iB] += Hxy;
-
-							dxdx[iC] += Hxx;
-							dydy[iC] += Hyy;
-							dxdy[iC] += Hxy;
-						}
-					} // for real all triangles 
-				}
-
-				Int4 kk=0;
-				for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){
-					if(Mmassxx[iv]>0){
-						dxdx[iv] /= 2*Mmassxx[iv];
-						// warning optimization (1) on term dxdy[iv]*ci/2 
-						dxdy[iv] /= 4*Mmassxx[iv];
-						dydy[iv] /= 2*Mmassxx[iv];
-						// Compute the matrix with abs(eigen value)
-						Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
-						MatVVP2x2 Vp(M);
-						Vp.Abs();
-						M = Vp;
-						dxdx[iv] = M.a11;
-						dxdy[iv] = M.a21;
-						dydy[iv] = M.a22;
-					}
-					else kk++;
-				}
-
-				// correction of second derivative
-				// by a laplacien
-				Real8* d2[3] = {dxdx, dxdy, dydy};
-				Real8* dd;
-				for (int xy = 0;xy<3;xy++) {
-					dd = d2[xy];
-					// do least 2 iteration for boundary problem
-					for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
-						for (i=0;i<nbt;i++) 
-						 if(triangles[i].link){// the real triangles 
-							 // number of the 3 vertices
-							 iA = Number(triangles[i][0]);
-							 iB = Number(triangles[i][1]);
-							 iC = Number(triangles[i][2]);
-							 Real8 cc=3;
-							 if(ijacobi==0)
-							  cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
-							 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
-						 }
-						for (iv=0;iv<nbv;iv++) workV[iv]=0;
-
-						for (i=0;i<nbt;i++){ 
-							if(triangles[i].link){ // the real triangles 
-								// number of the 3 vertices
-								iA = Number(triangles[i][0]);
-								iB = Number(triangles[i][1]);
-								iC = Number(triangles[i][2]);
-								Real8 cc =  workT[i]*detT[i];
-								workV[iA] += cc;
-								workV[iB] += cc;
-								workV[iC] += cc;
-							}
-						}
-
-						for (iv=0;iv<nbv;iv++){
-							if( ijacobi<NbJacobi || OnBoundary[iv]){
-								dd[iv] = workV[iv]/(Mmass[iv]*6);
-							}
-						}
-					}
-				}
-
-				/*Compute Metric from Hessian*/
-				for ( iv=0;iv<nbv;iv++){
-					vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[n*iv],bamgopts);
-				}
-
-			} //  end of for all field
-		}// end for all solution 
-
-		delete [] detT;
-		delete [] Mmass;
-		delete [] dxdx;
-		delete [] dxdy;
-		delete [] dydy;
-		delete []  workT;
-		delete [] workV;
-		delete [] Mmassxx;
-		delete []  OnBoundary;
-
-	}
-	/*}}}1*/
-	/*FUNCTION Triangles::BuildMetric2 (double P2 projection){{{1*/
-	void Triangles::BuildMetric2(BamgOpts* bamgopts){
-
-		/*Options*/
-		const int dim = 2;
-		double* s=NULL;
-		Int4   nbsol;
-		int*   typsols=NULL;
-		int    verbosity;
-
-		int   i,j,k,iA,iB,iC;
-		int   iv,nbfield;
-
-		/*Recover options*/
-		verbosity=bamgopts->verbose;
-
-		/*Get and process fields*/
-		s=bamgopts->field;
-		nbsol=bamgopts->numfields;
-		typsols=(int*)xmalloc(1*sizeof(int));
-		typsols[0]=0; // only one dof per node
-
-		//sizeoftype = {1,2,3,4}
-		int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 
-
-		// computation of the number of fields
-		Int4 ntmp = 0;
-		if (typsols){
-			//if there is more than one dof per vertex for one field
-			//increase ntmp to take into account all the fields
-			//if nbsol=1
-			for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
-		}
-		else ntmp = nbsol;
-
-		//n is the total number of fields
-		const Int4 n = ntmp;
-
-		//initialization of some variables
-		double* ss=(double*)s;
-		double  sA,sB,sC;
-		Real8*  detT = new Real8[nbt];
-		Real8*  sumareas = new Real8[nbv];
-		Real8*  alpha= new Real8[nbt*3];
-		Real8*  beta = new Real8[nbt*3];
-		Real8*  dx_elem    = new Real8[nbt];
-		Real8*  dy_elem    = new Real8[nbt];
-		Real8*  dx_vertex  = new Real8[nbv];
-		Real8*  dy_vertex  = new Real8[nbv];
-		Real8*  dxdx_elem  = new Real8[nbt];
-		Real8*  dxdy_elem  = new Real8[nbt];
-		Real8*  dydy_elem  = new Real8[nbt];
-		Real8*  dxdx_vertex= new Real8[nbv];
-		Real8*  dxdy_vertex= new Real8[nbv];
-		Real8*  dydy_vertex= new Real8[nbv];
-
-		//display infos
-		if(verbosity>1) {
-			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
-		}
-
-		//first, build the chains that will be used for the Hessian computation, as weel as the area of each element
-		int head_s[nbv];
-		int next_p[3*nbt];
-		int  p=0;
-		//initialization
-		for(i=0;i<nbv;i++){
-			sumareas[i]=0;
-			head_s[i]=-1;
-		}
-		for(i=0;i<nbt;i++){
-
-			//lopp over the real triangles (no boundary elements)
-			if(triangles[i].link){ 
-
-				//get current triangle t
-				const Triangle &t=triangles[i];
-
-				// coor of 3 vertices 
-				R2 A=t[0];
-				R2 B=t[1];
-				R2 C=t[2];
-
-				//compute triangle determinant (2*Area)
-				Real8 dett = bamg::Area2(A,B,C);
-				detT[i]=dett;
-
-				/*The nodal functions are such that for a vertex A:
-				 *    N_A(x,y)=alphaA x + beta_A y +gamma_A
-				 *    N_A(A) = 1,   N_A(B) = 0,   N_A(C) = 0
-				 * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
-				 * leads to:
-				 *    N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
-				 * and this gives:
-				 *    alpha_A = (yB-yC)/(2*Area(T))*/
-				alpha[i*3+0]=(B.y-C.y)/dett;
-				alpha[i*3+1]=(C.y-A.y)/dett;
-				alpha[i*3+2]=(A.y-B.y)/dett;
-				beta[ i*3+0]=(C.x-B.x)/dett;
-				beta[ i*3+1]=(A.x-C.x)/dett;
-				beta[ i*3+2]=(B.x-A.x)/dett;
-
-				//compute chains
-				for(j=0;j<3;j++){
-					k=Number(triangles[i][j]);
-					next_p[p]=head_s[k];
-					head_s[k]=p++;
-
-					//add area to sumareas
-					sumareas[k]+=dett;
-				}
-
-			}
-		}
-
-		//for all Solutions
-		for (Int4 nusol=0;nusol<nbsol;nusol++) {
-			int   nbfield=typsols?sizeoftype[typsols[nusol]]:1; 
-			Real8 smin=ss[0],smax=ss[0];
-
-			//only one field
-			if (nbfield == 1) {
-				//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
-				for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
-					smin=Min(smin,ss[k]);
-					smax=Max(smax,ss[k]);
-				}
-			}
-
-			//vectorial case
-			else{
-				for (iv=0,k=0;iv<nbv;iv++,k+=n ){
-					//compute v = √sum(s[i]^2)
-					double v=0;		     
-					for (int i=0;i<nbfield;i++){
-						v += ss[k+i]*ss[k+i];
-					}
-					v = sqrt(v);
-					smin=Min(smin,v);
-					smax=Max(smax,v);
-				}
-			}
-			Real8 sdelta=smax-smin;
-			Real8 absmax=Max(Abs(smin),Abs(smax));
-
-			//display info
-			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbfield);
-
-			//skip constant field
-			if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){
-				if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
-				continue;
-			}
-
-
-			//loop over all the fields of the solution
-			for (Int4 nufield=0;nufield<nbfield;nufield++,s++){
-				//ss++ so that for each iteration ss points toward the right field
-
-				//initialize the hessian and gradient matrices
-				for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
-
-				//1: Compute gradient for each element (exact)
-				for (i=0;i<nbt;i++){
-					if(triangles[i].link){
-						// number of the 3 vertices
-						iA = Number(triangles[i][0]);
-						iB = Number(triangles[i][1]);
-						iC = Number(triangles[i][2]);
-
-						// value of the P1 fonction on 3 vertices 
-						sA = ss[iA*n];
-						sB = ss[iB*n];
-						sC = ss[iC*n];
-
-						//gradient = (sum alpha_i s_i, sum_i beta_i s_i)
-						dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
-						dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
-					}
-				}
-
-				//2: then compute a gradient for each vertex using a P2 projection
-				for(i=0;i<nbv;i++){
-					for(p=head_s[i];p!=-1;p=next_p[p]){
-						//Get triangle number
-						k=(Int4)(p/3);
-						dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
-						dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
-					}
-				}
-
-				//3: compute Hessian matrix on each element
-				for (i=0;i<nbt;i++){
-					if(triangles[i].link){
-						// number of the 3 vertices
-						iA = Number(triangles[i][0]);
-						iB = Number(triangles[i][1]);
-						iC = Number(triangles[i][2]);
-
-						//Hessian
-						dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
-						dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
-						dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
-					}
-				}
-
-				//4: finaly compute Hessian on each vertex using the second P2 projection
-				for(i=0;i<nbv;i++){
-					for(p=head_s[i];p!=-1;p=next_p[p]){
-						//Get triangle number
-						k=(Int4)(p/3);
-						dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
-						dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
-						dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
-					}
-				}
-
-				/*Compute Metric from Hessian*/
-				for ( iv=0;iv<nbv;iv++){
-					vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[n*iv],bamgopts);
-				}
-
-			}//for all fields
-		}//for all solutions
-
-		//clean up
-		delete [] detT;
-		delete [] alpha;
-		delete [] beta;
-		delete [] sumareas;
-		delete [] dx_elem;
-		delete [] dy_elem;
-		delete [] dx_vertex;
-		delete [] dy_vertex;
-		delete [] dxdx_elem;
-		delete [] dxdy_elem;
-		delete [] dydy_elem;
-		delete [] dxdx_vertex;
-		delete [] dxdy_vertex;
-		delete [] dydy_vertex;
 	}
 	/*}}}1*/
@@ -3032,9 +2536,9 @@
 		else {
 			if (!quadtree){
-				throw ErrorException(__FUNCT__,exprintf("!quadtree"));
+				throw ErrorException(__FUNCT__,exprintf("no starting triangle provided and no quadtree available"));
 			}
 			Vertex *a = quadtree->NearestVertex(B.x,B.y) ;
 
-			if (! a || !a->t ) {
+			if (!a || !a->t ) {
 				if (a) {
 					printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a));
@@ -3051,9 +2555,12 @@
 		}
 		Icoor2  detop ;
-		int kkkk =0; // number of test triangle 
-
-		while ( t->det < 0){ // the initial triangles is outside  
+
+		// number of test triangle 
+		int kkkk =0; 
+
+		// if the initial triangles is outside  
+		while ( t->det < 0){ 
 			int k0=(*t)(0) ?  ((  (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1  )) : 0;
-			if (k0<0){ // k0 the NULL  vertex
+			if (k0<0){ // k0 the NULL vertex
 				throw ErrorException(__FUNCT__,exprintf("k0<0"));
 			}
@@ -3151,12 +2658,26 @@
 		nbv=0;
 		for (i=0;i<Gh.nbv;i++){
-			/*Add vertex only if required and not duplicate 
+			/* Add vertex only if required and not duplicate 
 			 * (IsThe is a method of GeometricalVertex that checks
 			 * that we are taking into acount only one vertex is
 			 * 2 vertices are superimposed) */
 			if (Gh[i].Required() && Gh[i].IsThe()) {//Gh  vertices Required
-				if (nbv<nbvx) vertices[nbv]=Gh[i];
-				Gh[i].to = vertices + nbv;// save Geom -> Th
+
+				//Add the vertex (provided that nbv<nbvx)
+				if (nbv<nbvx){
+					vertices[nbv]=Gh[i];
+				}
+				else{
+					throw ErrorException(__FUNCT__,exprintf("Maximum number of vertices (nbvx = %i) too small",nbvx));
+				}
+				
+				//Add pointer from geometry (Gh) to vertex from mesh (Th)
+				Gh[i].to=vertices+nbv;
+
+
+				//Build VerticesOnGeomVertex for current point
 				VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);
+
+				//nbv increment
 				nbv++;
 			}
@@ -3165,5 +2686,5 @@
 		//check that edges is empty
 		if (edges){
-			throw ErrorException(__FUNCT__,exprintf("edges"));
+			throw ErrorException(__FUNCT__,exprintf("edges is empty"));
 		}
 
@@ -3184,6 +2705,9 @@
 			//go through the edges of the geometry
 			for (i=0;i<Gh.nbe;i++) {
+
 				GeometricalEdge &ei=Gh.edges[i];   
-				if (!ei.Dup()){ // a good curve (not dup)
+
+				// a good curve (not dup)
+				if (!ei.Dup()){ 
 					for(int j=0;j<2;j++) {
 						if (!ei.Mark() && ei[j].Required()) { 
@@ -3734,9 +3258,6 @@
 		BuildMetric1(bamgopts);
 	}
-	else if (Hessiantype==2){
-		BuildMetric2(bamgopts);
-	}
 	else{
-		throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (0->use Green formula, 1-> from P2 on 4T, 2-> double P2 projection)",Hessiantype));
+		throw ErrorException(__FUNCT__,exprintf("Hessiantype %i not supported yet (1->use Green formula, 0-> double P2 projection)",Hessiantype));
 	}
 }
@@ -3898,67 +3419,70 @@
 	/*FUNCTION Triangles::InsertNewPoints{{{1*/
 	Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) {
-		long int verbosity=2;
+
+		long int verbosity=0;
 		Real8 seuil= 1.414/2 ;// for two close point 
 		Int4 i;
-		// insertion part --- 
-
-		const Int4 nbvnew = nbv-nbvold;
-		if (verbosity>5) printf("      Try to Insert the %i new points\n",nbvnew);
 		Int4 NbSwap=0;
 		Icoor2 dete[3];
 
-		// construction d'un ordre aleatoire 
-		if (! nbvnew) 
-		 return 0; 
-		if (nbvnew) {
-			const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv)  ;
-			Int4 k3 = rand()%nbvnew ; 
-			for (Int4 is3=0; is3<nbvnew; is3++) {
-				register Int4 j = nbvold +(k3 = (k3 + PrimeNumber)% nbvnew);
-				register Int4 i = nbvold+is3; 
-				ordre[i]= vertices + j;
-				ordre[i]->ReferenceNumber=i;
-			}
-			// be carefull 
-			Int4  iv = nbvold;
-			for (i=nbvold;i<nbv;i++) 
-			  {// for all the new point
-				Vertex & vi = *ordre[i];
-				vi.i = toI2(vi.r);
-				vi.r = toR2(vi.i);
-				Real4 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 
-					Vertex & vj = vertices[iv];
-					Int4 j = vj.ReferenceNumber; 
-					if ( &vj!= ordre[j]){
-						throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]"));
-					}
-					if(i!=j)
-					  { //  for valgring 
-						Exchange(vi,vj);
-						Exchange(ordre[j],ordre[i]);
-					  }
-					vj.ReferenceNumber=0; 
-					Triangle *tcvj = FindTriangleContening(vj.i,dete);
-					if (tcvj && !tcvj->link){
-						throw ErrorException(__FUNCT__,exprintf("problem inserting point"));
-					}
-					quadtree->Add(vj);
-					Add(vj,tcvj,dete);
-					NbSwap += vj.Optim(1);          
-					iv++;
-				  }
-			  } 
-			if (verbosity>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;
-		}
+		//number of new points
+		const Int4 nbvnew=nbv-nbvold;
+
+		//display info if required
+		if (verbosity>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 Int4 PrimeNumber= AGoodNumberPrimeWith(nbv)  ;
+		//remainder of the division of rand() by nbvnew
+		Int4 k3 = rand()%nbvnew;
+		//loop over the new points
+		for (Int4 is3=0; is3<nbvnew; is3++){
+			register Int4 j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
+			register Int4 i=nbvold+is3; 
+			ordre[i]= vertices + j;
+			ordre[i]->ReferenceNumber=i;
+		}
+
+		// for all the new point
+		Int4 iv=nbvold;
+		for (i=nbvold;i<nbv;i++){
+			Vertex &vi=*ordre[i];
+			vi.i=toI2(vi.r);
+			vi.r=toR2(vi.i);
+			Real4 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 
+				Vertex &vj = vertices[iv];
+				Int4  j=vj.ReferenceNumber; 
+				if (&vj!=ordre[j]){
+					throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]"));
+				}
+				if(i!=j){ 
+					Exchange(vi,vj);
+					Exchange(ordre[j],ordre[i]);
+				}
+				vj.ReferenceNumber=0; 
+				Triangle *tcvj=FindTriangleContening(vj.i,dete);
+				if (tcvj && !tcvj->link){
+					tcvj->Echo();
+					throw ErrorException(__FUNCT__,exprintf("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link));
+				}
+				quadtree->Add(vj);
+				Add(vj,tcvj,dete);
+				NbSwap += vj.Optim(1);          
+				iv++;
+			}
+		} 
+		if (verbosity>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);  
@@ -5488,5 +5012,5 @@
 			Triangle *tcvi = FindTriangleContening(vi.i,dete);
 			if (tcvi && !tcvi->link) {
-				printf("problem inserting point\n");
+				printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
 			}
 
@@ -5909,4 +5433,6 @@
 	/*FUNCTION AGoodNumberPrimeWith{{{1*/
 	Int4 AGoodNumberPrimeWith(Int4 n){
+
+		//list of big prime numbers
 		const Int4 BigPrimeNumber[] ={ 567890359L,
 			567890431L,  567890437L,  567890461L,  567890471L,
@@ -5914,12 +5440,22 @@
 			567890591L,  567890599L,  567890621L,  567890629L , 0};
 
-		Int4 o = 0;
-		Int4 pi = BigPrimeNumber[1];
-		for (int i=0; BigPrimeNumber[i]; i++) {
+		//initialize o and pi
+		Int4 o =0;
+		Int4 pi=BigPrimeNumber[1];
+
+		//loop until BigPrimeNumber[i]==0 (end of BigPrimeNumber)
+		for (int i=0; BigPrimeNumber[i]; i++){
+
+			//compute r, rest of the remainder of the division of BigPrimeNumber[i] by n
 			Int4 r = BigPrimeNumber[i] % n;
+
+			/*compute oo = min ( r , n-r , |n - 2r|, |n-3r|)*/
 			Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));
-			if ( o < oo) 
-			 o=oo,pi=BigPrimeNumber[i];}
-			return pi; 
+			if ( o < oo){
+				o=oo;
+				pi=BigPrimeNumber[i];
+			}
+		}
+		return pi; 
 	}
 	/*}}}1*/
Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2939)
+++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2940)
@@ -206,7 +206,19 @@
 	}
 	/*}}}1*/
+	/*FUNCTION Vertex::Echo {{{1*/
+
+	void Vertex::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);
+
+		return;
+	}
+	/*}}}*/
 
 	/*Intermediary*/
-	/*FUNCTION QuadQuality{{{{1*/
+	/*FUNCTION QuadQuality{{{1*/
 	double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) {
 		// calcul de 4 angles --
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 2939)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 2940)
@@ -83,5 +83,5 @@
 bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1);
 bamg_options.AbsError=getfieldvalue(options,'AbsError',0);
-bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',2);
+bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
 bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
 bamg_options.NbSmooth=getfieldvalue(options,'NbSmooth',3);
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2939)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2940)
@@ -152,5 +152,4 @@
 	FetchData(&field,NULL,&numfields,mxGetField(BAMGOPTIONS,0,"field"));
 	bamgopts.numfields=numfields;
-	printf("numfields = %i\n",numfields);
 	bamgopts.field=field;