Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2956)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2957)
@@ -660,5 +660,4 @@
 			Int4 NbOfTriangleSearchFind;
 			Int4 NbOfSwapTriangle;
-			char * name, *identity;
 			Vertex * vertices;   // data of vertices des sommets
 
@@ -710,5 +709,5 @@
 
 			Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR
-			Triangles(const Triangles &,const int *flag,const int *bb); // truncature
+			Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
 
 			void SetIntCoor(const char * from =0);
@@ -786,5 +785,5 @@
 			int UnCrack();
 
-			void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo 
+			void BuildGeometryFromMesh(BamgOpts* bamgopts);
 			void FillHoleInMesh() ;
 			int  CrackMesh();
@@ -792,5 +791,5 @@
 			void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption
 			void GeomToTriangles0(Int4 nbvx);                   // the real constructor mesh generator
-			void PreInit(Int4,char* =NULL );
+			void PreInit(Int4);
 
 	};
@@ -801,5 +800,4 @@
 			Int4 NbRef; // counter of ref on the this class if 0 we can delete
 
-			char *name;
 			Int4 nbvx,nbtx; // nombre max  de sommets , de  Geometry
 			Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux,
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2956)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2957)
@@ -26,7 +26,4 @@
 		NbRef =0;
 		quadtree=0;
-		name = new char[strlen(Gh.name)+4];
-		strcpy(name,"cp:");
-		strcat(name,Gh.name);
 		vertices = nbv ? new GeometricalVertex[nbv] : NULL;
 		triangles = nbt ? new  Triangle[nbt]:NULL;
@@ -60,5 +57,4 @@
 		if(quadtree)  delete  quadtree;quadtree=0;
 		if(curves)  delete  []curves;curves=0;NbOfCurves=0;
-		if(name) delete [] name;name=0;
 		if(subdomains) delete [] subdomains;subdomains=0;
 		//  if(ordre)     delete [] ordre;
@@ -877,5 +873,4 @@
 
 		printf("Geometry:\n");
-		printf("   name: %s\n",name);
 		printf("   nbvx (maximum number of vertices) : %i\n",nbvx);
 		printf("   nbtx (maximum number of triangles): %i\n",nbtx);
@@ -905,5 +900,4 @@
 	void Geometry::EmptyGeometry() {
 		NbRef=0;
-		name =0;
 		quadtree=0;
 		curves=0;
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2956)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2957)
@@ -26,5 +26,5 @@
 	Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
 
-		PreInit(0,"none");
+		PreInit(0);
 		ReadMesh(bamgmesh,bamgopts);
 		SetIntCoor();
@@ -33,7 +33,5 @@
 	/*}}}1*/
 	/*FUNCTION Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb){{{1*/
-	Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb) : Gh(*(new Geometry())), BTh(*this) {
-
-		  char cname[] = "trunc";
+	Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb,BamgOpts* bamgopts) : Gh(*(new Geometry())), BTh(*this) {
 
 		  int i,k,itadj;
@@ -83,5 +81,5 @@
 		  printf("   number of New boundary edge %i\n",nbNewBedge);
 		  Int4 inbvx =k;
-		  PreInit(inbvx,cname);
+		  PreInit(inbvx);
 		  for (i=0;i<Tho.nbv;i++)
 			if (kk[i]>=0) 
@@ -121,6 +119,6 @@
 		  delete [] reft;
 		  delete [] refv;
-		  double cutoffradian = 10.0/180.0*Pi;
-		  ConsGeometry(cutoffradian);
+		  //double cutoffradian = 10.0/180.0*Pi;
+		  BuildGeometryFromMesh(bamgopts);
 		  Gh.AfterRead(); 
 		  SetIntCoor();
@@ -148,6 +146,4 @@
 		if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
 		if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex;
-		//if (name) delete [] name; //TEST crash if not commented
-		if (identity) delete [] identity;
 		if (VertexOnBThVertex) delete [] VertexOnBThVertex;
 		if (VertexOnBThEdge) delete [] VertexOnBThEdge;
@@ -172,11 +168,5 @@
 		  // do all the allocation to be sure all the pointer existe
 
-		  char * cname = 0;
-		  if (Th.name) 
-			 {
-			  cname = new char[strlen(Th.name)+1];
-			  strcpy(cname,Th.name);
-			 }
-		  PreInit(nbvxx,cname);// to make the allocation 
+		  PreInit(nbvxx);// to make the allocation 
 		  // copy of triangles
 		  nt=Th.nt;
@@ -302,4 +292,5 @@
 			triangles =new Triangle[nbt];
 			for (i=0;i<bamgmesh->NumQuadrilaterals;i++){
+				//divide the quad into two triangles
 				Triangle & t1 = triangles[2*i];
 				Triangle & t2 = triangles[2*i+1];
@@ -466,5 +457,5 @@
 			printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
 			/*Recreate geometry: */
-			ConsGeometry(bamgopts->MaximalAngleOfCorner*Pi/180);
+			BuildGeometryFromMesh(bamgopts);
 			Gh.AfterRead();
 		}
@@ -885,523 +876,30 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::BuildMetric0 (double P2 projection){{{1*/
-	void Triangles::BuildMetric0(BamgOpts* bamgopts){
-
-		/*Options*/
-		const int dim = 2;
-		double* s=NULL;
-		Int4    nbsol;
-		int     verbosity;
-
-		int   i,j,k,iA,iB,iC;
-		int   iv;
+	/*FUNCTION Triangles::BuildGeometryFromMesh{{{1*/
+	void Triangles::BuildGeometryFromMesh(BamgOpts* bamgopts){
+		/*Reconstruct Geometry from Mesh*/
+
+		/*Intermediary*/
+		int i,j,k;
 
 		/*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->err[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*/
-		const int dim = 2;
-		double* s=NULL;
-		Int4 nbsol;
-		int NbJacobi;
-		int verbosity;
-
-		/*Recover options*/
-		verbosity=bamgopts->verbose;
-		NbJacobi=bamgopts->nbjacobi;
-
-		/*Get and process fields*/
-		s=bamgopts->field;
-		nbsol=bamgopts->numfields;
-
-		//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",nbsol,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
-				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[nusol],smax=ss[nusol];
-			Real8 h1=1.e30,h2=1e-30,rx=0;
-			Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-
-			//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;
-			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,nbsol);
-
-			//skip constant field
-			if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
-				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; 
-
-				//initialize the hessian matrix
-				for ( iv=0,k=0; iv<nbv; iv++) 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*nbsol+nusol];
-						sB = ss[iB*nbsol+nusol];
-						sC = ss[iC*nbsol+nusol];
-
-						/*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];
-
-						//Use Green to compute Hessian Matrix
-
-						// if edge on boundary no contribution  => normal = 0
-						if ( !tBC || !tBC->link ) nBC=O;
-						if ( !tCA || !tCA->link ) nCA=O;
-						if ( !tAB || !tAB->link ) nAB=O;
-
-						// remark we forgot a 1/2 because
-						//       int_{edge} w_i = 1/2 if i is in edge 
-						//                         0  if not
-						// if we don't take the  boundary 
-						dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
-						dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
-						dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
-
-						//warning optimization (1) the division by 2 is done on the metric construction
-						dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
-						dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
-						dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 
-
-						dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
-						dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
-						dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
-
-					} // for real all triangles 
-				}
-
-				Int4 kk=0;
-				for ( iv=0,k=0 ; iv<nbv; iv++){
-					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 leat 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[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
-				}
-
-		}// 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::ConsGeometry{{{1*/
-	void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges){
-		//  if equiedges existe taille nbe 
-		//   equiedges[i]/2 == i  original
-		//   equiedges[i]/2 = j  =>   equivalence entre i et j => meme maillage
-		//   equiedges[i]%2   : 0 meme sens , 1 pas meme sens 
-		//       
-		// --------------------------
-		long int verbosity=0;
+		int    verbosity=bamgopts->verbose;
+		double cutoffradian=bamgopts->MaximalAngleOfCorner*Pi/180;
+		CurrentTh=this;
+
+		//display info
 		if (verbosity>1) printf("   construction of the geometry from the 2d mesh\n");
+
+		//check that the mesh is not empty
 		if (nbt<=0 || nbv <=0 ) {
-			throw ErrorException(__FUNCT__,exprintf("nbt or nbv is negative"));
-		}
-
-		// construction of the edges 
-		//  Triangles * OldCurrentTh =CurrentTh;
-		CurrentTh=this;
-		//  Int4 NbTold = nbt;
-		// generation of the integer coor
-		// generation of the adjacence of the triangles
-		if (cutoffradian>=0)
-		 Gh.MaximalAngleOfCorner = cutoffradian;
-		SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
-		Int4 * st = new Int4[nbt*3];
-		Int4 i,k;
-		int j; 
-		if (Gh.name) delete Gh.name;
-		Gh.name = new char [ name ? strlen(name) + 15 : 50 ];
-		Gh.name[0]=0;
-		strcat(Gh.name,"cons from: ");
-		if (name) strcat(Gh.name,name);
-		else strcat(Gh.name," a mesh with no name");
+			throw ErrorException(__FUNCT__,exprintf("nbt or nbv is negative (Mesh empty?)"));
+		}
+
+		//Gh is the geometry of the mesh (this), initialize MaximalAngleOfCorner
+		if (cutoffradian>=0) Gh.MaximalAngleOfCorner = cutoffradian;
+
+		//Construction of the edges 
+		SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
+		Int4*        st   = new Int4[nbt*3];
 		for (i=0;i<nbt*3;i++)
 		 st[i]=-1;
@@ -1445,5 +943,5 @@
 
 		if(verbosity>5) {
-			printf("         info on Mesh %s:\n",name);
+			printf("         info on Meshs:\n");
 			printf("            - number of vertices    = %i \n",nbv); 
 			printf("            - number of triangles   = %i \n",nbt); 
@@ -1697,16 +1195,4 @@
 			if(requis) kreq++;
 			edges[i].onGeometry =  Gh.edges + i;
-			if(equiedges && i < nbeold ) {
-				int j=equiedges[i]/2;
-				int sens=equiedges[i]%2;
-				if(i!=j && equiedges[i]>=0) {
-					if( sens==0)
-					 Gh.edges[i].SetEqui();
-					else 
-					 Gh.edges[i].SetReverseEqui();
-					Gh.edges[i].link= & Gh.edges[j];
-				}
-
-			}
 			if(requis)  {  // correction fevr 2009 JYU ...
 				Gh.edges[i].v[0]->SetRequired();
@@ -1766,5 +1252,492 @@
 		  triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j));
 
-	  }
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::BuildMetric0 (double P2 projection){{{1*/
+	void Triangles::BuildMetric0(BamgOpts* bamgopts){
+
+		/*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->err[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*/
+		const int dim = 2;
+		double* s=NULL;
+		Int4 nbsol;
+		int NbJacobi;
+		int verbosity;
+
+		/*Recover options*/
+		verbosity=bamgopts->verbose;
+		NbJacobi=bamgopts->nbjacobi;
+
+		/*Get and process fields*/
+		s=bamgopts->field;
+		nbsol=bamgopts->numfields;
+
+		//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",nbsol,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
+				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[nusol],smax=ss[nusol];
+			Real8 h1=1.e30,h2=1e-30,rx=0;
+			Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
+
+			//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;
+			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,nbsol);
+
+			//skip constant field
+			if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
+				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; 
+
+				//initialize the hessian matrix
+				for ( iv=0,k=0; iv<nbv; iv++) 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*nbsol+nusol];
+						sB = ss[iB*nbsol+nusol];
+						sC = ss[iC*nbsol+nusol];
+
+						/*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];
+
+						//Use Green to compute Hessian Matrix
+
+						// if edge on boundary no contribution  => normal = 0
+						if ( !tBC || !tBC->link ) nBC=O;
+						if ( !tCA || !tCA->link ) nCA=O;
+						if ( !tAB || !tAB->link ) nAB=O;
+
+						// remark we forgot a 1/2 because
+						//       int_{edge} w_i = 1/2 if i is in edge 
+						//                         0  if not
+						// if we don't take the  boundary 
+						dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
+						dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
+						dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
+
+						//warning optimization (1) the division by 2 is done on the metric construction
+						dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
+						dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
+						dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 
+
+						dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
+						dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
+						dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
+
+					} // for real all triangles 
+				}
+
+				Int4 kk=0;
+				for ( iv=0,k=0 ; iv<nbv; iv++){
+					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 leat 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[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
+				}
+
+		}// 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::ConsRefTriangle{{{1*/
@@ -2067,5 +2040,5 @@
 			}
 			if(verbosity>5) {
-				printf("         info on Mesh %s:\n",name);
+				printf("         info on Mesh\n");
 				printf("            - number of vertices    = %i \n",nbv); 
 				printf("            - number of triangles   = %i \n",nbt); 
@@ -3851,5 +3824,5 @@
 /*}}}1*/
 /*FUNCTION Triangles::PreInit{{{1*/
-void Triangles::PreInit(Int4 inbvx,char *fname) {
+void Triangles::PreInit(Int4 inbvx) {
 	long int verbosity=0;
 
@@ -3857,5 +3830,4 @@
 	NbRef=0;
 	//  allocGeometry=0;
-	identity=0;
 	NbOfTriangleSearchFind =0;
 	NbOfSwapTriangle =0;
@@ -3876,5 +3848,4 @@
 	CrackedEdges  =0;  
 	nbe = 0; 
-	name = fname ;
 
 	if (inbvx) {