Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.cpp	(revision 2812)
+++ /issm/trunk/src/c/Bamgx/Mesh2.cpp	(revision 2813)
@@ -1074,92 +1074,4 @@
 		}
 
-		Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j) 
-		{ return  quadtree->NearestVertex(i,j); } 
-
-
-		void Triangles::PreInit(Int4 inbvx,char *fname)
-		{
-			long int verbosity=0;
-
-			srand(19999999);
-			OnDisk =0;
-			NbRef=0;
-			//  allocGeometry=0;
-			identity=0;
-			NbOfTriangleSearchFind =0;
-			NbOfSwapTriangle =0;
-			nbiv=0;
-			nbv=0;
-			nbvx=inbvx;
-			nbt=0;
-			NbOfQuad = 0;
-			nbtx=2*inbvx-2;
-			NbSubDomains=0;
-			NbVertexOnBThVertex=0;
-			NbVertexOnBThEdge=0;
-			VertexOnBThVertex=0;
-			VertexOnBThEdge=0;
-
-			NbCrackedVertices=0;
-			NbCrackedEdges =0;
-			CrackedEdges  =0;  
-			nbe = 0; 
-			name = fname ;
-
-			if (inbvx) {
-				vertices=new Vertex[nbvx];
-				assert(vertices);
-				ordre=new (Vertex* [nbvx]);
-				assert(ordre);
-				triangles=new Triangle[nbtx];
-				assert(triangles);}
-			else {
-				vertices=0;
-				ordre=0;
-				triangles=0;
-				nbtx=0;
-			}
-			if ( name || inbvx)
-			{
-				time_t timer =time(0);
-				char buf[70];     
-				strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer));
-				counter++; 
-				char countbuf[30];   
-				sprintf(countbuf,"%d",counter);
-				int lg =0 ;
-				if (&BTh != this && BTh.name)
-					lg = strlen(BTh.name)+4;
-				identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2  + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)];
-				identity[0]=0;
-				if (lg)
-					strcat(strcat(strcat(identity,"B="),BTh.name),", ");
-
-				if (Gh.name)
-					strcat(strcat(identity,"G="),Gh.name);
-				strcat(strcat(identity,";"),countbuf);
-				strcat(identity,buf);
-				// cout << "New MAILLAGE "<< identity << endl;
-			} 
-
-			quadtree=0;
-			//  edgescomponante=0;
-			edges=0;
-			VerticesOnGeomVertex=0;
-			VerticesOnGeomEdge=0;
-			NbVerticesOnGeomVertex=0;
-			NbVerticesOnGeomEdge=0;
-			//  nbMaxIntersectionTriangles=0;
-			//  lIntTria;
-			subdomains=0;
-			NbSubDomains=0;
-			//  Meshbegin = vertices;
-			//  Meshend  = vertices + nbvx;
-			if (verbosity>98) 
-				cout << "Triangles::PreInit() " << nbvx << " " << nbtx 
-					<< " " << vertices 
-					<< " " << ordre << " " <<  triangles <<endl;
-		}
-
 
 double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) {
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2812)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2813)
@@ -961,4 +961,11 @@
 }; // End Class Geometry
 
+//From Metric.cpp
+inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3])
+{ return    A[0] * ( B[1]*C[2]-B[2]*C[1])
+          - A[1] * ( B[0]*C[2]-B[2]*C[0])
+          + A[2] * ( B[0]*C[1]-B[1]*C[0]);
+}
+
 /////////////////////////////////////////////////////////////////////////////////////
 /////////////////////////////////////////////////////////////////////////////////////
Index: /issm/trunk/src/c/Bamgx/Metric.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.cpp	(revision 2812)
+++ /issm/trunk/src/c/Bamgx/Metric.cpp	(revision 2813)
@@ -31,11 +31,4 @@
 
 namespace bamg {
-
-
-inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3])
-{ return    A[0] * ( B[1]*C[2]-B[2]*C[1])
-          - A[1] * ( B[0]*C[2]-B[2]*C[0])
-          + A[2] * ( B[0]*C[1]-B[1]*C[0]);
-}
 
 SaveMetricInterpole  LastMetricInterpole;
@@ -206,793 +199,4 @@
       return r;
 }
-void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0)
-
-{
-	long int verbosity=0;
-
-  if(verbosity>1)
-    cout << "  -- IntersectGeomMetric geometric err=" << err << (iso ? " iso " : " aniso "  ) << endl;
-  Real8 ss[2]={0.00001,0.99999};
-  Real8 errC = 2*sqrt(2*err);
-  Real8 hmax = Gh.MaximalHmax();
-  Real8 hmin = Gh.MinimalHmin();
-  Real8 maxaniso = 1e6;
-  assert(hmax>0);
-  SetVertexFieldOn();
-  if (errC > 1) errC = 1;
-  for (Int4  i=0;i<nbe;i++)
-   for (int j=0;j<2;j++)
-    {
-      
-      Vertex V;
-      VertexOnGeom GV;
-      // cerr << Number(edges[i]) << " " << ss[j] << endl;
-      Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
-	{
-	  GeometricalEdge * eg = GV;
-	  Real8 s = GV;
-	  R2 tg;
-	  //	   cerr << i << " " << j << " " << Number(V) << " on = " 
-	  //	<< Gh.Number(eg) << " at s = " << s << " " << endl;
-	  Real8  R1= eg->R1tg(s,tg);
-	  // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x " 
-	  //    << V.r << errC/ Max(R1,1e-20) <<  " hold=" <<V.m(tg) << " "  << endl;
-	  Real8 ht = hmax;
-          if (R1>1.0e-20) 
-	    {  // err relative to the length of the edge
-	      ht = Min(Max(errC/R1,hmin),hmax);
-	    }
-	  Real8 hn = iso? ht : Min(hmax,ht*maxaniso);
-	  //cerr << ht << " " << hn << "m=" << edges[i][j].m <<  endl;
-	  assert(ht>0 && hn>0);
-	  MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
-	  //cerr << " : " ;
-	  Metric MVp(Vp);
-	  // cerr << " : "  << MVp  << endl;
-	  edges[i][j].m.IntersectWith(MVp);
-	  //cerr << " . " << endl;
-	}
-
-    }
-  // the problem is for the vertex on vertex 
-     
-}
-/*
-void  Triangles::BoundAnisotropy(Real8 anisomax)
-{
-  if (verbosity > 1) 
-    cout << "  -- BoundAnisotropy by  " << anisomax << endl; 
-  Real8 h1=1.e30,h2=1e-30,rx=0;
-  Real8 coef = 1./(anisomax*anisomax);
-  Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-  for (Int4 i=0;i<nbv;i++)
-    {
-
-      MatVVP2x2 Vp(vertices[i]);
-      
-      h1=Min(h1,Vp.lmin());
-      h2=Max(h2,Vp.lmax());
-      rx = Max(rx,Vp.Aniso2());
-      
-      Vp.BoundAniso2(coef);
-      
-      hn1=Min(hn1,Vp.lmin());
-      hn2=Max(hn2,Vp.lmax());
-      rnx = Max(rnx,Vp.Aniso2());
-
-      
-      vertices[i].m = Vp;
-
-    }
-
-  if (verbosity>2)
-    {
-      cout << "     input :  Hmin = " << sqrt(1/h2)  << " Hmax = " << sqrt(1/h1) 
-	   << " factor of anisotropy max  = " << sqrt(rx) << endl;
-      cout << "     output:  Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 
-	   << " factor of anisotropy max  = " << sqrt(rnx) << endl;
-    }
-}
-*/
-void  Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso)
-{
-	long int verbosity=0;
-
-  double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
-  if (verbosity > 1) 
-    cout << "  -- BoundAnisotropy by  " << anisomax << endl; 
-  Real8 h1=1.e30,h2=1e-30,rx=0;
-  Real8 coef = 1./(anisomax*anisomax);
-  Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-  for (Int4 i=0;i<nbv;i++)
-    {
-
-      MatVVP2x2 Vp(vertices[i]);
-      double lmax=Vp.lmax();
-      h1=Min(h1,Vp.lmin());
-      h2=Max(h2,Vp.lmax());
-      rx = Max(rx,Vp.Aniso2());
-
-      Vp *= Min(lminaniso,lmax)/lmax;
-      
-      Vp.BoundAniso2(coef);
-      
-      hn1=Min(hn1,Vp.lmin());
-      hn2=Max(hn2,Vp.lmax());
-      rnx = Max(rnx,Vp.Aniso2());
-
-      
-      vertices[i].m = Vp;
-
-    }
-
-  if (verbosity>2)
-    {
-      cout << "     input :  Hmin = " << sqrt(1/h2)  << " Hmax = " << sqrt(1/h1) 
-	   << " factor of anisotropy max  = " << sqrt(rx) << endl;
-      cout << "     output:  Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 
-	   << " factor of anisotropy max  = " << sqrt(rnx) << endl;
-    }
-}
-void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols,
-				    const  Real8 hmin1,const Real8 hmax1,const Real8 coef,
-				    const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,
-				    const int DoNormalisation,const double power,const int choice)
-{ //  the array of solution s is store    
-  // sol0,sol1,...,soln    on vertex 0
-  //  sol0,sol1,...,soln   on vertex 1
-  //  etc.
-  //  choise = 0 =>  H is computed with green formule
-  //   otherwise  => H is computed from P2 on 4T 
-  const int dim = 2;
-  
-  long int verbosity=0;
-
-  int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 
-
-  // computation of the nb of field 
-  Int4 ntmp = 0;
-  if (typsols)
-    {
-      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;
-
-  Int4 i,k,iA,iB,iC,iv;
-  R2 O(0,0);
-  int RelativeMetric = CutOff>1e-30;
-  Real8 hmin = Max(hmin1,MinimalHmin());
-  Real8 hmax = Min(hmax1,MaximalHmax());
-  Real8 coef2 = 1/(coef*coef);
-
-  if(verbosity>1) 
-    {
-      cout << "  -- Construction of Metric: Nb of field. " << n << " nbt = " 
-	   << nbt << " nbv= " << nbv 
-	   << " coef = " << coef << endl
-	   << "     hmin = " << hmin << " hmax=" << hmax 
-	   << " anisomax = " << anisomax <<  " Nb Jacobi " << NbJacobi << " Power = " << power ;
-      if (RelativeMetric)
-	cout << " RelativeErr with CutOff= "  <<  CutOff << endl;
-      else
-	cout << " Absolute Err" <<endl;
-    }
-  double *ss=(double*)s;//, *ssiii = ss;
-
-  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];
-  for (iv=0;iv<nbv;iv++)
-    {
-      Mmass[iv]=0;
-      OnBoundary[iv]=0;
-      Mmassxx[iv]=0;
-    }
-
-  for (i=0;i<nbt;i++) 
-    if(triangles[i].link) // the real triangles 
-      {
-	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]);
-	
-	Real8 dett = bamg::Area2(A,B,C);
-	detT[i]=dett;
-	dett /= 6;
-
-	// construction of on boundary 
-	int nbb =0;
-	for(int j=0;j<3;j++)
-          {
-	    Triangle *ta=t.Adj(j);
-	    if ( ! ta || !ta->link) // no adj triangle => edge on boundary
-	      OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1,
-		OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1,
-		nbb++;
-	  }
-	
-	workT[i] = nbb;
-	Mmass[iA] += dett;
-	Mmass[iB] += dett;
-	Mmass[iC] += dett;
-	
-	if((nbb==0)|| !choice)
-	  {
-	    Mmassxx[iA] += dett;
-	    Mmassxx[iB] += dett;
-	    Mmassxx[iC] += dett;
-	  }
-      }
-  else
-    workT[i]=-1;
-
-//  for (Int4 kcount=0;kcount<n;kcount++,ss++)
-    for (Int4 nusol=0;nusol<nbsol;nusol++)
-    { //for all Solution  
-
-      Real8 smin=ss[0],smax=ss[0];
-      
-      Real8 h1=1.e30,h2=1e-30,rx=0;
-      Real8 coef = 1./(anisomax*anisomax);
-      Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
-      int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 
-      if (nbfield == 1) 
-       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]);
-				 }
-			  else
-			   {
-         //  cas vectoriel 
-          for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-          {	
-           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));
-      Real8 cnorm = DoNormalisation ? coef2/sdelta : coef2;
-      
-      if(verbosity>2) 
-	     cout << "    Solution " << nusol <<  " Min = " << smin << " Max = " 
-	       << smax << " Delta =" << sdelta << " cnorm = " << cnorm <<  " Nb of fields =" << nbfield << endl;
-
-      
-      if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) 
-				{
-				  if (verbosity>2)
-				    cout << "      Solution " << nusol << " is constant. We skip. " 
-					 << " Min = " << smin << " Max = " << smax << endl;
-				continue;
-				}
-				
-	 double *sf  = ss; 
-	 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++) 
-	   {
-	     for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-		       dxdx[iv]=dxdy[iv]=dydy[iv]=0;
-       for (i=0;i<nbt;i++) 
-	      if(triangles[i].link)
-	  {// for real all triangles 
-	    // 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);
-	    // remark :  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];
-
-	    R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ;
-	    if(choice) 
-	      {
-		int nbb = 0;
-		Real8 dd = detT[i];
-		Real8 lla,llb,llc,llf;
-		Real8  taa[3][3],bb[3];
-		// construction of the trans of lin system
-		for (int j=0;j<3;j++)
-		  {
-		    int ie = OppositeEdge[j];
-		    TriangleAdjacent ta = triangles[i].Adj(ie);
-		    Triangle *tt = ta;
-		    if (tt && tt->link)
-		      {
-			Vertex &v = *ta.OppositeVertex();
-			R2 V = v;
-			Int4 iV = Number(v);
-			Real8 lA  = bamg::Area2(V,B,C)/dd;
-			Real8 lB  = bamg::Area2(A,V,C)/dd;
-			Real8 lC  = bamg::Area2(A,B,V)/dd;
-			taa[0][j] =  lB*lC;
-			taa[1][j] =  lC*lA;
-			taa[2][j] =  lA*lB;
-			//Real8 xx = V.x-V.y;
-			//Real8 yy = V.x + V.y;
-			//cout << " iv " << ss[iV*n] << " == " << (8*xx*xx+yy*yy)
-			//     << " l = " << lA << " " << lB << " " << lC 
-			//     << " = " << lA+lB+lC << " " <<  V << " == " << A*lA+B*lB+C*lC << endl;
-			
-			lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;
-
-			bb[j]     =  ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ;
-		      }
-		    else
-		      {
-			nbb++;
-			taa[0][j]=0;
-			taa[1][j]=0;
-			taa[2][j]=0;
-			taa[j][j]=1;
-			bb[j]=0;
-		      }
-		  }
-
-		// resolution of 3x3 lineaire 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);
-		
-		assert(det33);
-		//	det33=1;
-		// verif
-		//	cout << " " << (taa[0][0]*cBC +  taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ;
-		//	cout << " " << (taa[0][1]*cBC +  taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1];
-		//	cout << " " << (taa[0][2]*cBC +  taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2] 
-		//	     << "  -- " ;
-		//cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB +  llb*llc* cBC + llc*lla*cCA)/det33 
-		//   << " == " << llf <<  endl;
-		// computation of the gradient 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);
-		Real8 coef = 1.0/(3*dd*det33);
-		Real8 coef2 = 2*coef;
-		//	cout << " H = " << Hxx << " " << Hyy << " " <<  Hxy/2 << " coef2 = " << coef2 << endl;
-		Hxx *= coef2;
-		Hyy *= coef2;
-		Hxy *= coef2;
-		//cout << i  << " H = " << 3*Hxx/dd << " " << 3*Hyy/dd << " " <<  3*Hxy/(dd*2) << " nbb = " << nbb << endl;
-		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;
-		  }
-		
-	      }
-	    else
-	      {
-		
-		// 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[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 divide by 2 is done on the metrix 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++,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);
-	    //cout <<iv <<  "  M  = " <<  M <<  " aniso= " << Vp.Aniso() ;
-	    Vp.Abs();
-	    M = Vp;
-	      dxdx[iv] = M.a11;
-	      dxdy[iv] = M.a21;
-	      dydy[iv] = M.a22;
-	      //  cout << " (abs)  iv M  = " <<  M <<  " aniso= " << Vp.Aniso() <<endl;
-	  }
-	else kk++;
-      
-      
-      // correction of second derivate
-      // 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);
-	      
-
-	    }
-
-	  
-	}
-
-      // constuction  of the metrix from the Hessian dxdx. dxdy,dydy
-
-      Real8 rCutOff=CutOff*absmax;// relative cut off 
-
-      for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )
-	{ // for all vertices 
-	  //{
-	  //Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
-	  // MatVVP2x2 Vp(M);	  
-	  //cout << " iv M="<<  M << "  Vp = " << Vp << " aniso  " << Vp.Aniso() << endl;
-	  //}
-	  MetricIso Miso;
-// new code to compute ci ---	  
-	  Real8 ci ;
-	  if (RelativeMetric)
-	    { //   compute the norm of the solution
-	       double xx =0,*sfk=sf+k; 
-	       for (int ifield=0;ifield<nbfield;ifield++,sfk++)
-	          xx += *sfk* *sfk;	       
-	       xx=sqrt(xx);
-	       ci = coef2/Max(xx,rCutOff);
-	    }
-	  else ci = cnorm;
-	  
- // old 
-//	  Real8 ci = RelativeMetric ? coef2/(Max(Abs(ss[k]),rCutOff)) : cnorm ;
- //   modif F Hecht 101099
-	  Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci,  dydy[iv]*ci);
-	  MatVVP2x2 Vp(Miv);
-
-	  Vp.Abs();
-	 if(power!=1.0) 
-	      Vp.pow(power);
-	  
-
-
-	  h1=Min(h1,Vp.lmin());
-	  h2=Max(h2,Vp.lmax());
-
-	  Vp.Maxh(hmin);
-	  Vp.Minh(hmax);
-
-	  rx = Max(rx,Vp.Aniso2());
-
-	  Vp.BoundAniso2(coef);
-
-	  hn1=Min(hn1,Vp.lmin());
-	  hn2=Max(hn2,Vp.lmax());
-	  rnx = Max(rnx,Vp.Aniso2());
-
-	  Metric MVp(Vp);
-	  vertices[iv].m.IntersectWith(MVp);
-	}// for all vertices 
-      if (verbosity>2)
-	{ 
-	  cout << "              Field " << nufield << " of solution " << nusol  << endl;
-	  cout << "              before bounding :  Hmin = " << sqrt(1/h2) << " Hmax = " 
-	       << sqrt(1/h1)  << " factor of anisotropy max  = " << sqrt(rx) << endl;
-	  cout << "              after  bounding :  Hmin = " << sqrt(1/hn2) << " Hmax = " 
-	       << sqrt(1/hn1)  << " factor of anisotropy max  = " << sqrt(rnx) << endl;
-	}
-	 } //  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;
- 
-}
-
-
-void Triangles::ReadMetric(BamgOpts* bamgopts,const Real8 hmin1=1.0e-30,const Real8 hmax1=1.0e30,const Real8 coef=1)
-{
-	int  i,j;
-
-	if(bamgopts->verbose>3) printf("      processing metric\n");
-
-	Real8 hmin = Max(hmin1,MinimalHmin());
-	Real8 hmax = Min(hmax1,MaximalHmax());
-
-	//for now we only use j==3
-	j=3;
-
-	for (i=0;i<nbv;i++){
-		Real8 h;
-		if (j == 1){
-			h=bamgopts->metric[i];
-			vertices[i].m=Metric(Max(hmin,Min(hmax, h*coef)));
-		}
-		else if (j==3){
-			Real8 a,b,c;	     
-			a=bamgopts->metric[i*3+0];
-			b=bamgopts->metric[i*3+1];
-			c=bamgopts->metric[i*3+2];
-			MetricAnIso M(a,b,c);
-			MatVVP2x2 Vp(M/coef);
-
-			Vp.Maxh(hmin);
-			Vp.Minh(hmax);
-			vertices[i].m = Vp;
-		}
-	}
-}
-
-void Triangles::WriteMetric(ostream & f,int iso)
-{
-  if (iso)
-    {
-      f <<  nbv <<" " << 1 << endl ;
-      for (Int4 iv=0;iv<nbv;iv++)
-	{
-	  MatVVP2x2 V=vertices[iv].m;
-	  f <<  V.hmin()  << endl;
-	}
-    }
-else
-  {
-    f <<  nbv <<" " << 3 << endl ;
-    for (Int4 iv=0;iv<nbv;iv++)
-      f <<  vertices[iv].m.a11 << " " 
-	<<  vertices[iv].m.a21 << " " 
-	<<  vertices[iv].m.a22 << endl;
-  }
-}
-void  Triangles::MaxSubDivision(Real8 maxsubdiv)
-{
-	long int verbosity=0;
-
-const  Real8 maxsubdiv2 = maxsubdiv*maxsubdiv;
-  if(verbosity>1)
-    cout << "  -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv <<   endl  ;
-  // for all the edges 
-  // if the len of the edge is to long 
-  Int4 it,nbchange=0;    
-  Real8 lmax=0;
-  for (it=0;it<nbt;it++)
-    {
-      Triangle &t=triangles[it];
-      for (int j=0;j<3;j++)
-	{
-	  Triangle &tt = *t.TriangleAdj(j);
-	  if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)) 
-	    {
-		Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];
-		Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];
-		R2 AB= (R2) v1-(R2) v0;
-		Metric M = v0;
-		Real8 l = M(AB,AB);
-		lmax = Max(lmax,l);
-		if(l> maxsubdiv2)
-		  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
-		    Real8 lc = M(AC,AC);
-		    D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
-		    D2xD2 Rt1(Rt.inv());
-		    D2xD2 D(maxsubdiv2,0,0,lc);
-		    D2xD2 MM = Rt1*D*Rt1.t();
-		    v0.m =  M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);
-		    nbchange++;
-		  }
-		M = v1;
-		l = M(AB,AB);
-		lmax = Max(lmax,l);
-		if(l> maxsubdiv2)
-		  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
-		    Real8 lc = M(AC,AC);
-		    D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
-		    D2xD2 Rt1(Rt.inv());
-		    D2xD2 D(maxsubdiv2,0,0,lc);
-		    D2xD2  MM = Rt1*D*Rt1.t();
-		    v1.m =  M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);
-		    nbchange++;
-		  }
-		
-		
-	    }
-	}
-    }
-  if(verbosity>3)
-  cout << "    Nb of metric change = " << nbchange 
-       << " Max  of the subdivision of a edges before change  = " << sqrt(lmax) << endl;
-
-}
-
-void Triangles::SmoothMetric(Real8 raisonmax) 
-{ 
-	long int verbosity=0;
-
-  if(raisonmax<1.1) return;
-  if(verbosity > 1)
-     cout << "  -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl;
-  ReMakeTriangleContainingTheVertex();
-  Int4 i,j,kch,kk,ip;
-  Int4 *first_np_or_next_t0 = new Int4[nbv];
-  Int4 *first_np_or_next_t1 = new Int4[nbv];
-  Int4 Head0 =0,Head1=-1;
-  Real8 logseuil= log(raisonmax);
-
-  for(i=0;i<nbv-1;i++)
-    first_np_or_next_t0[i]=i+1; 
-  first_np_or_next_t0[nbv-1]=-1;// end;
-  for(i=0;i<nbv;i++)
-    first_np_or_next_t1[i]=-1;
-  kk=0;
-  while (Head0>=0&& kk++<100) {
-    kch=0;
-    for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1)
-      {  //  pour tous les triangles autour du sommet s
-	// 	cout << kk << " i = " << i << " " << ip << endl;
-	register Triangle * t= vertices[i].t;
-	assert(t);
-	Vertex & vi = vertices[i];
-	TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
-	Vertex *pvj0 = ta.EdgeVertex(0);
-	while (1) {
-	  //	  cout << i << " " <<  Number(ta.EdgeVertex(0)) << " "
-	  //      << Number(ta.EdgeVertex(1)) << "  ---> " ;
-	  ta=Previous(Adj(ta));
-	  // cout <<  Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl;
-	  assert(vertices+i == ta.EdgeVertex(1));
-	  Vertex & vj = *(ta.EdgeVertex(0));
-	  if ( &vj ) {
-	    j= &vj-vertices;
-	    assert(j>=0 && j < nbv);
-	    R2 Aij = (R2) vj - (R2) vi;
-	    Real8 ll =  Norme2(Aij);
-	    if (0) {  
-	      Real8 hi = ll/vi.m(Aij);
-	      Real8 hj = ll/vj.m(Aij);
-	      if(hi < hj)
-		{
-		  Real8 dh=(hj-hi)/ll;
-		  //cout << " dh = " << dh << endl;
-		  if (dh>logseuil) {
-		    vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
-		    if(first_np_or_next_t1[j]<0)
-		      kch++,first_np_or_next_t1[j]=Head1,Head1=j;
-		  }
-		}
-	    } 
-	    else
-	      {
-		Real8 li = vi.m(Aij);
-		//Real8 lj = vj.m(Aij);
-		//		if ( i == 2 || j == 2)
-		//  cout << " inter " << i << " " << j << " " << ((1 +logseuil*li)) <<  endl;
-	      	if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
-		  //if( vj.m.IntersectWith(vi.m*(lj/li/(1 +logseuil*lj))) )
-		  if(first_np_or_next_t1[j]<0) // if the metrix change 
-		    kch++,first_np_or_next_t1[j]=Head1,Head1=j;
-	      }
-	  }
-	  if  ( &vj ==  pvj0 ) break;
-	}
-      }
-    Head0 = Head1;
-    Head1 = -1;
-    Exchange(first_np_or_next_t0,first_np_or_next_t1);
-    if(verbosity>5)
-    cout << "     Iteration " << kk << " Nb de  vertices with change  " << kch << endl;
-  }
-  if(verbosity>2 && verbosity < 5) 
-    cout << "    Nb of Loop " << kch << endl;
-  delete [] first_np_or_next_t0;
-  delete [] first_np_or_next_t1;
-}
 
 Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB)
Index: /issm/trunk/src/c/Bamgx/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Triangles.cpp	(revision 2812)
+++ /issm/trunk/src/c/Bamgx/Triangles.cpp	(revision 2813)
@@ -233,205 +233,4 @@
 
 	/*IO*/
-	/*FUNCTION Triangles::WriteMesh {{{1*/
-	void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
-
-		int i,j;
-		int verbose;
-
-		verbose=bamgopts->verbose;
-		Int4 *reft = new Int4[nbt];
-		Int4 nbInT = ConsRefTriangle(reft);
-
-		//Vertices
-		if(verbose>3) printf("      writing Vertices\n");
-		bamgmesh->NumVertices=nbv;
-		xfree((void**)&bamgmesh->Vertices);
-		if (nbv){
-			bamgmesh->Vertices=(double*)xmalloc(3*nbv*sizeof(double));
-			for (i=0;i<nbv;i++){
-				bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
-				bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
-				bamgmesh->Vertices[i*3+2]=vertices[i].ref();
-			}
-		}
-		else{
-			bamgmesh->Vertices=NULL;
-		}
-
-		//Edges
-		if(verbose>3) printf("      writing Edges\n");
-		bamgmesh->NumEdges=nbe;
-		xfree((void**)&bamgmesh->Edges);
-		if (nbe){
-			bamgmesh->Edges=(double*)xmalloc(3*nbe*sizeof(double));
-			for (i=0;i<nbe;i++){
-				bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing
-				bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing
-				bamgmesh->Edges[i*3+2]=edges[i].ref;
-			}
-		}
-		else{
-			bamgmesh->Edges=NULL;
-		}
-
-		//CrackedEdges
-		if(verbose>3) printf("      writing CrackedEdges\n");
-		bamgmesh->NumCrackedEdges=NbCrackedEdges;
-		xfree((void**)&bamgmesh->CrackedEdges);
-		if (NbCrackedEdges){
-			bamgmesh->CrackedEdges=(double*)xmalloc(2*NbCrackedEdges*sizeof(double));
-			for (i=0;i<NbCrackedEdges;i++){
-				bamgmesh->CrackedEdges[i*2+0]=Number(CrackedEdges[i].a.edge)+1; //back to M indexing
-				bamgmesh->CrackedEdges[i*2+1]=Number(CrackedEdges[i].b.edge)+1; //back to M indexing
-			}
-		}
-		else{
-			bamgmesh->CrackedEdges=NULL;
-		}
-
-		//Triangles
-		if(verbose>3) printf("      writing Triangles\n");
-		Int4 k=nbInT-NbOfQuad*2;
-		Int4 num=0;
-		bamgmesh->NumTriangles=k;
-		xfree((void**)&bamgmesh->Triangles);
-		if (k){
-			bamgmesh->Triangles=(double*)xmalloc(4*k*sizeof(double));
-			for (i=0;i<nbt;i++){
-				Triangle &t=triangles[i];
-				if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
-					bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing
-					bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing
-					bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing
-					bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
-					num=num+1;
-				}
-			}
-		}
-		else{
-			bamgmesh->Triangles=NULL;
-		}
-
-		//Quadrilaterals
-		if(verbose>3) printf("      writing Quadrilaterals\n");
-		bamgmesh->NumQuadrilaterals=NbOfQuad;
-		xfree((void**)&bamgmesh->Quadrilaterals);
-		if (NbOfQuad){
-			bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double));
-			for (i=0;i<nbt;i++){
-				Triangle &t =triangles[i];
-				Triangle* ta;
-				Vertex *v0,*v1,*v2,*v3;
-				if (reft[i]<0) continue;
-				if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { 
-					bamgmesh->Quadrilaterals[i*5+0]=Number(v0)+1; //back to M indexing
-					bamgmesh->Quadrilaterals[i*5+1]=Number(v1)+1; //back to M indexing
-					bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing
-					bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing
-					bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref;
-				}
-			}
-		}
-		else{
-			bamgmesh->Quadrilaterals=NULL;
-		}
-
-		//SubDomains
-		if(verbose>3) printf("      writing SubDomains\n");
-		bamgmesh->NumSubDomains=NbSubDomains;
-		xfree((void**)&bamgmesh->SubDomains);
-		if (NbSubDomains){
-			bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
-			for (i=0;i<NbSubDomains;i++){
-				bamgmesh->SubDomains[i*4+0]=3;
-				bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)];
-				bamgmesh->SubDomains[i*4+2]=1;
-				bamgmesh->SubDomains[i*4+3]=subdomains[i].ref;
-			}
-		}
-		else{
-			bamgmesh->SubDomains=NULL;
-		}
-
-		//SubDomainsFromGeom
-		if(verbose>3) printf("      writing SubDomainsFromGeom\n");
-		bamgmesh->NumSubDomainsFromGeom=Gh.NbSubDomains;
-		xfree((void**)&bamgmesh->SubDomainsFromGeom);
-		if (Gh.NbSubDomains){
-			bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double));
-			for (i=0;i<Gh.NbSubDomains;i++){
-				bamgmesh->SubDomainsFromGeom[i*4+0]=2;
-				bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
-				bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens;
-				bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref;
-			}
-		}
-		else{
-			bamgmesh->SubDomainsFromGeom=NULL;
-		}
-
-		//VerticesOnGeomVertex
-		if(verbose>3) printf("      writing VerticesOnGeometricVertex\n");
-		bamgmesh->NumVerticesOnGeometricVertex=NbVerticesOnGeomVertex;
-		xfree((void**)&bamgmesh->VerticesOnGeometricVertex);
-		if (NbVerticesOnGeomVertex){
-			bamgmesh->VerticesOnGeometricVertex=(double*)xmalloc(2*NbVerticesOnGeomVertex*sizeof(double));
-			for (i=0;i<NbVerticesOnGeomVertex;i++){
-				VertexOnGeom &v=VerticesOnGeomVertex[i];
-				if (!v.OnGeomVertex()){
-					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not"));
-				}
-				bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing
-				bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing
-			}
-		}
-		else{
-			bamgmesh->VerticesOnGeometricVertex=NULL;
-		}
-
-		//VertexOnGeometricEdge
-		if(verbose>3) printf("      writing VerticesOnGeometricEdge\n");
-		bamgmesh->NumVerticesOnGeometricEdge=NbVerticesOnGeomEdge;
-		xfree((void**)&bamgmesh->VerticesOnGeometricEdge);
-		if (NbVerticesOnGeomEdge){
-			bamgmesh->VerticesOnGeometricEdge=(double*)xmalloc(3*NbVerticesOnGeomEdge*sizeof(double));
-			for (i=0;i<NbVerticesOnGeomEdge;i++){
-				const VertexOnGeom &v=VerticesOnGeomEdge[i];
-				if (!v.OnGeomEdge()){
-					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricEdge is actually not"));
-				}
-				bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing
-				bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
-				bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v;
-			}
-		}
-		else{
-			bamgmesh->VerticesOnGeometricEdge=NULL;
-		}
-
-		//EdgesOnGeometricEdge
-		if(verbose>3) printf("      writing EdgesOnGeometricEdge\n");
-		k=0;
-		for (i=0;i<nbe;i++){
-			if (edges[i].on) k=k+1;
-		}
-		bamgmesh->NumEdgesOnGeometricEdge=k;
-		xfree((void**)&bamgmesh->EdgesOnGeometricEdge);
-		if (k){
-			bamgmesh->EdgesOnGeometricEdge=(double*)xmalloc(2*(int)k*sizeof(double));
-			int count=0;
-			for (i=0;i<nbe;i++){
-				if (edges[i].on){
-					bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing
-					bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].on)+1; //back to Matlab indexing
-					count=count+1;
-				}
-			}
-		}
-		else{
-			bamgmesh->EdgesOnGeometricEdge=NULL;
-		}
-	}
-	/*}}}1*/
 	/*FUNCTION Triangles::ReadMesh{{{1*/
 	void Triangles::ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){
@@ -662,4 +461,259 @@
 			Gh.AfterRead();
 		}
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::WriteMesh {{{1*/
+	void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
+
+		int i,j;
+		int verbose;
+
+		verbose=bamgopts->verbose;
+		Int4 *reft = new Int4[nbt];
+		Int4 nbInT = ConsRefTriangle(reft);
+
+		//Vertices
+		if(verbose>3) printf("      writing Vertices\n");
+		bamgmesh->NumVertices=nbv;
+		xfree((void**)&bamgmesh->Vertices);
+		if (nbv){
+			bamgmesh->Vertices=(double*)xmalloc(3*nbv*sizeof(double));
+			for (i=0;i<nbv;i++){
+				bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
+				bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
+				bamgmesh->Vertices[i*3+2]=vertices[i].ref();
+			}
+		}
+		else{
+			bamgmesh->Vertices=NULL;
+		}
+
+		//Edges
+		if(verbose>3) printf("      writing Edges\n");
+		bamgmesh->NumEdges=nbe;
+		xfree((void**)&bamgmesh->Edges);
+		if (nbe){
+			bamgmesh->Edges=(double*)xmalloc(3*nbe*sizeof(double));
+			for (i=0;i<nbe;i++){
+				bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing
+				bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing
+				bamgmesh->Edges[i*3+2]=edges[i].ref;
+			}
+		}
+		else{
+			bamgmesh->Edges=NULL;
+		}
+
+		//CrackedEdges
+		if(verbose>3) printf("      writing CrackedEdges\n");
+		bamgmesh->NumCrackedEdges=NbCrackedEdges;
+		xfree((void**)&bamgmesh->CrackedEdges);
+		if (NbCrackedEdges){
+			bamgmesh->CrackedEdges=(double*)xmalloc(2*NbCrackedEdges*sizeof(double));
+			for (i=0;i<NbCrackedEdges;i++){
+				bamgmesh->CrackedEdges[i*2+0]=Number(CrackedEdges[i].a.edge)+1; //back to M indexing
+				bamgmesh->CrackedEdges[i*2+1]=Number(CrackedEdges[i].b.edge)+1; //back to M indexing
+			}
+		}
+		else{
+			bamgmesh->CrackedEdges=NULL;
+		}
+
+		//Triangles
+		if(verbose>3) printf("      writing Triangles\n");
+		Int4 k=nbInT-NbOfQuad*2;
+		Int4 num=0;
+		bamgmesh->NumTriangles=k;
+		xfree((void**)&bamgmesh->Triangles);
+		if (k){
+			bamgmesh->Triangles=(double*)xmalloc(4*k*sizeof(double));
+			for (i=0;i<nbt;i++){
+				Triangle &t=triangles[i];
+				if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
+					bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing
+					bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing
+					bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing
+					bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
+					num=num+1;
+				}
+			}
+		}
+		else{
+			bamgmesh->Triangles=NULL;
+		}
+
+		//Quadrilaterals
+		if(verbose>3) printf("      writing Quadrilaterals\n");
+		bamgmesh->NumQuadrilaterals=NbOfQuad;
+		xfree((void**)&bamgmesh->Quadrilaterals);
+		if (NbOfQuad){
+			bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double));
+			for (i=0;i<nbt;i++){
+				Triangle &t =triangles[i];
+				Triangle* ta;
+				Vertex *v0,*v1,*v2,*v3;
+				if (reft[i]<0) continue;
+				if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { 
+					bamgmesh->Quadrilaterals[i*5+0]=Number(v0)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+1]=Number(v1)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref;
+				}
+			}
+		}
+		else{
+			bamgmesh->Quadrilaterals=NULL;
+		}
+
+		//SubDomains
+		if(verbose>3) printf("      writing SubDomains\n");
+		bamgmesh->NumSubDomains=NbSubDomains;
+		xfree((void**)&bamgmesh->SubDomains);
+		if (NbSubDomains){
+			bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
+			for (i=0;i<NbSubDomains;i++){
+				bamgmesh->SubDomains[i*4+0]=3;
+				bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)];
+				bamgmesh->SubDomains[i*4+2]=1;
+				bamgmesh->SubDomains[i*4+3]=subdomains[i].ref;
+			}
+		}
+		else{
+			bamgmesh->SubDomains=NULL;
+		}
+
+		//SubDomainsFromGeom
+		if(verbose>3) printf("      writing SubDomainsFromGeom\n");
+		bamgmesh->NumSubDomainsFromGeom=Gh.NbSubDomains;
+		xfree((void**)&bamgmesh->SubDomainsFromGeom);
+		if (Gh.NbSubDomains){
+			bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double));
+			for (i=0;i<Gh.NbSubDomains;i++){
+				bamgmesh->SubDomainsFromGeom[i*4+0]=2;
+				bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
+				bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens;
+				bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref;
+			}
+		}
+		else{
+			bamgmesh->SubDomainsFromGeom=NULL;
+		}
+
+		//VerticesOnGeomVertex
+		if(verbose>3) printf("      writing VerticesOnGeometricVertex\n");
+		bamgmesh->NumVerticesOnGeometricVertex=NbVerticesOnGeomVertex;
+		xfree((void**)&bamgmesh->VerticesOnGeometricVertex);
+		if (NbVerticesOnGeomVertex){
+			bamgmesh->VerticesOnGeometricVertex=(double*)xmalloc(2*NbVerticesOnGeomVertex*sizeof(double));
+			for (i=0;i<NbVerticesOnGeomVertex;i++){
+				VertexOnGeom &v=VerticesOnGeomVertex[i];
+				if (!v.OnGeomVertex()){
+					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not"));
+				}
+				bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing
+			}
+		}
+		else{
+			bamgmesh->VerticesOnGeometricVertex=NULL;
+		}
+
+		//VertexOnGeometricEdge
+		if(verbose>3) printf("      writing VerticesOnGeometricEdge\n");
+		bamgmesh->NumVerticesOnGeometricEdge=NbVerticesOnGeomEdge;
+		xfree((void**)&bamgmesh->VerticesOnGeometricEdge);
+		if (NbVerticesOnGeomEdge){
+			bamgmesh->VerticesOnGeometricEdge=(double*)xmalloc(3*NbVerticesOnGeomEdge*sizeof(double));
+			for (i=0;i<NbVerticesOnGeomEdge;i++){
+				const VertexOnGeom &v=VerticesOnGeomEdge[i];
+				if (!v.OnGeomEdge()){
+					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricEdge is actually not"));
+				}
+				bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v;
+			}
+		}
+		else{
+			bamgmesh->VerticesOnGeometricEdge=NULL;
+		}
+
+		//EdgesOnGeometricEdge
+		if(verbose>3) printf("      writing EdgesOnGeometricEdge\n");
+		k=0;
+		for (i=0;i<nbe;i++){
+			if (edges[i].on) k=k+1;
+		}
+		bamgmesh->NumEdgesOnGeometricEdge=k;
+		xfree((void**)&bamgmesh->EdgesOnGeometricEdge);
+		if (k){
+			bamgmesh->EdgesOnGeometricEdge=(double*)xmalloc(2*(int)k*sizeof(double));
+			int count=0;
+			for (i=0;i<nbe;i++){
+				if (edges[i].on){
+					bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing
+					bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].on)+1; //back to Matlab indexing
+					count=count+1;
+				}
+			}
+		}
+		else{
+			bamgmesh->EdgesOnGeometricEdge=NULL;
+		}
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::ReadMetric{{{1*/
+	void Triangles::ReadMetric(BamgOpts* bamgopts,const Real8 hmin1=1.0e-30,const Real8 hmax1=1.0e30,const Real8 coef=1) {
+		int  i,j;
+
+		if(bamgopts->verbose>3) printf("      processing metric\n");
+
+		Real8 hmin = Max(hmin1,MinimalHmin());
+		Real8 hmax = Min(hmax1,MaximalHmax());
+
+		//for now we only use j==3
+		j=3;
+
+		for (i=0;i<nbv;i++){
+			Real8 h;
+			if (j == 1){
+				h=bamgopts->metric[i];
+				vertices[i].m=Metric(Max(hmin,Min(hmax, h*coef)));
+			}
+			else if (j==3){
+				Real8 a,b,c;	     
+				a=bamgopts->metric[i*3+0];
+				b=bamgopts->metric[i*3+1];
+				c=bamgopts->metric[i*3+2];
+				MetricAnIso M(a,b,c);
+				MatVVP2x2 Vp(M/coef);
+
+				Vp.Maxh(hmin);
+				Vp.Minh(hmax);
+				vertices[i].m = Vp;
+			}
+		}
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::WriteMetric TO UPDATE{{{1*/
+	void Triangles::WriteMetric(ostream & f,int iso) {
+		if (iso)
+		  {
+			f <<  nbv <<" " << 1 << endl ;
+			for (Int4 iv=0;iv<nbv;iv++)
+			  {
+				MatVVP2x2 V=vertices[iv].m;
+				f <<  V.hmin()  << endl;
+			  }
+		  }
+		else
+		  {
+			f <<  nbv <<" " << 3 << endl ;
+			for (Int4 iv=0;iv<nbv;iv++)
+			 f <<  vertices[iv].m.a11 << " " 
+				<<  vertices[iv].m.a21 << " " 
+				<<  vertices[iv].m.a22 << endl;
+		  }
 	}
 	/*}}}1*/
@@ -4918,4 +4972,792 @@
 	}
 	/*}}}1*/
+	/*FUNCTION Triangles::IntersectGeomMetric{{{1*/
+	void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0){
+		long int verbosity=0;
+
+		if(verbosity>1)
+		 cout << "  -- IntersectGeomMetric geometric err=" << err << (iso ? " iso " : " aniso "  ) << endl;
+		Real8 ss[2]={0.00001,0.99999};
+		Real8 errC = 2*sqrt(2*err);
+		Real8 hmax = Gh.MaximalHmax();
+		Real8 hmin = Gh.MinimalHmin();
+		Real8 maxaniso = 1e6;
+		assert(hmax>0);
+		SetVertexFieldOn();
+		if (errC > 1) errC = 1;
+		for (Int4  i=0;i<nbe;i++)
+		 for (int j=0;j<2;j++)
+			{
+
+			 Vertex V;
+			 VertexOnGeom GV;
+			 // cerr << Number(edges[i]) << " " << ss[j] << endl;
+			 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
+				{
+				 GeometricalEdge * eg = GV;
+				 Real8 s = GV;
+				 R2 tg;
+				 //	   cerr << i << " " << j << " " << Number(V) << " on = " 
+				 //	<< Gh.Number(eg) << " at s = " << s << " " << endl;
+				 Real8  R1= eg->R1tg(s,tg);
+				 // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x " 
+				 //    << V.r << errC/ Max(R1,1e-20) <<  " hold=" <<V.m(tg) << " "  << endl;
+				 Real8 ht = hmax;
+				 if (R1>1.0e-20) 
+					{  // err relative to the length of the edge
+					 ht = Min(Max(errC/R1,hmin),hmax);
+					}
+				 Real8 hn = iso? ht : Min(hmax,ht*maxaniso);
+				 //cerr << ht << " " << hn << "m=" << edges[i][j].m <<  endl;
+				 assert(ht>0 && hn>0);
+				 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
+				 //cerr << " : " ;
+				 Metric MVp(Vp);
+				 // cerr << " : "  << MVp  << endl;
+				 edges[i][j].m.IntersectWith(MVp);
+				 //cerr << " . " << endl;
+				}
+
+			}
+		// the problem is for the vertex on vertex 
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::BoundAnisotropy{{{1*/
+	void  Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso) {
+		long int verbosity=0;
+
+		double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
+		if (verbosity > 1) 
+		 cout << "  -- BoundAnisotropy by  " << anisomax << endl; 
+		Real8 h1=1.e30,h2=1e-30,rx=0;
+		Real8 coef = 1./(anisomax*anisomax);
+		Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
+		for (Int4 i=0;i<nbv;i++)
+		  {
+
+			MatVVP2x2 Vp(vertices[i]);
+			double lmax=Vp.lmax();
+			h1=Min(h1,Vp.lmin());
+			h2=Max(h2,Vp.lmax());
+			rx = Max(rx,Vp.Aniso2());
+
+			Vp *= Min(lminaniso,lmax)/lmax;
+
+			Vp.BoundAniso2(coef);
+
+			hn1=Min(hn1,Vp.lmin());
+			hn2=Max(hn2,Vp.lmax());
+			rnx = Max(rnx,Vp.Aniso2());
+
+
+			vertices[i].m = Vp;
+
+		  }
+
+		if (verbosity>2)
+		  {
+			cout << "     input :  Hmin = " << sqrt(1/h2)  << " Hmax = " << sqrt(1/h1) 
+			  << " factor of anisotropy max  = " << sqrt(rx) << endl;
+			cout << "     output:  Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 
+			  << " factor of anisotropy max  = " << sqrt(rnx) << endl;
+		  }
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::IntersectConsMetric{{{1*/
+	void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols,
+				const  Real8 hmin1,const Real8 hmax1,const Real8 coef,
+				const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,
+				const int DoNormalisation,const double power,const int choice)
+	  { //  the array of solution s is store    
+		// sol0,sol1,...,soln    on vertex 0
+		//  sol0,sol1,...,soln   on vertex 1
+		//  etc.
+		//  choise = 0 =>  H is computed with green formule
+		//   otherwise  => H is computed from P2 on 4T 
+		const int dim = 2;
+
+		long int verbosity=0;
+
+		int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 
+
+		// computation of the nb of field 
+		Int4 ntmp = 0;
+		if (typsols)
+		  {
+			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;
+
+		Int4 i,k,iA,iB,iC,iv;
+		R2 O(0,0);
+		int RelativeMetric = CutOff>1e-30;
+		Real8 hmin = Max(hmin1,MinimalHmin());
+		Real8 hmax = Min(hmax1,MaximalHmax());
+		Real8 coef2 = 1/(coef*coef);
+
+		if(verbosity>1) 
+		  {
+			cout << "  -- Construction of Metric: Nb of field. " << n << " nbt = " 
+			  << nbt << " nbv= " << nbv 
+			  << " coef = " << coef << endl
+			  << "     hmin = " << hmin << " hmax=" << hmax 
+			  << " anisomax = " << anisomax <<  " Nb Jacobi " << NbJacobi << " Power = " << power ;
+			if (RelativeMetric)
+			 cout << " RelativeErr with CutOff= "  <<  CutOff << endl;
+			else
+			 cout << " Absolute Err" <<endl;
+		  }
+		double *ss=(double*)s;//, *ssiii = ss;
+
+		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];
+		for (iv=0;iv<nbv;iv++)
+		  {
+			Mmass[iv]=0;
+			OnBoundary[iv]=0;
+			Mmassxx[iv]=0;
+		  }
+
+		for (i=0;i<nbt;i++) 
+		 if(triangles[i].link) // the real triangles 
+			{
+			 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]);
+
+			 Real8 dett = bamg::Area2(A,B,C);
+			 detT[i]=dett;
+			 dett /= 6;
+
+			 // construction of on boundary 
+			 int nbb =0;
+			 for(int j=0;j<3;j++)
+				{
+				 Triangle *ta=t.Adj(j);
+				 if ( ! ta || !ta->link) // no adj triangle => edge on boundary
+				  OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1,
+					 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1,
+					 nbb++;
+				}
+
+			 workT[i] = nbb;
+			 Mmass[iA] += dett;
+			 Mmass[iB] += dett;
+			 Mmass[iC] += dett;
+
+			 if((nbb==0)|| !choice)
+				{
+				 Mmassxx[iA] += dett;
+				 Mmassxx[iB] += dett;
+				 Mmassxx[iC] += dett;
+				}
+			}
+		 else
+		  workT[i]=-1;
+
+		for (Int4 nusol=0;nusol<nbsol;nusol++)
+		  { //for all Solution  
+
+			Real8 smin=ss[0],smax=ss[0];
+
+			Real8 h1=1.e30,h2=1e-30,rx=0;
+			Real8 coef = 1./(anisomax*anisomax);
+			Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
+			int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 
+			if (nbfield == 1) 
+			 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]);
+				}
+			else
+			  {
+				//  cas vectoriel 
+				for ( iv=0,k=0; iv<nbv; iv++,k+=n )
+				  {	
+					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));
+			Real8 cnorm = DoNormalisation ? coef2/sdelta : coef2;
+
+			if(verbosity>2) 
+			 cout << "    Solution " << nusol <<  " Min = " << smin << " Max = " 
+				<< smax << " Delta =" << sdelta << " cnorm = " << cnorm <<  " Nb of fields =" << nbfield << endl;
+
+
+			if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) 
+			  {
+				if (verbosity>2)
+				 cout << "      Solution " << nusol << " is constant. We skip. " 
+					<< " Min = " << smin << " Max = " << smax << endl;
+				continue;
+			  }
+
+			double *sf  = ss; 
+			for (Int4 nufield=0;nufield<nbfield;nufield++,ss++) 
+			  {
+				for ( iv=0,k=0; iv<nbv; iv++,k+=n )
+				 dxdx[iv]=dxdy[iv]=dydy[iv]=0;
+				for (i=0;i<nbt;i++) 
+				 if(triangles[i].link)
+					{// for real all triangles 
+					 // 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);
+					 // remark :  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];
+
+					 R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ;
+					 if(choice) 
+						{
+						 int nbb = 0;
+						 Real8 dd = detT[i];
+						 Real8 lla,llb,llc,llf;
+						 Real8  taa[3][3],bb[3];
+						 // construction of the trans of lin system
+						 for (int j=0;j<3;j++)
+							{
+							 int ie = OppositeEdge[j];
+							 TriangleAdjacent ta = triangles[i].Adj(ie);
+							 Triangle *tt = ta;
+							 if (tt && tt->link)
+								{
+								 Vertex &v = *ta.OppositeVertex();
+								 R2 V = v;
+								 Int4 iV = Number(v);
+								 Real8 lA  = bamg::Area2(V,B,C)/dd;
+								 Real8 lB  = bamg::Area2(A,V,C)/dd;
+								 Real8 lC  = bamg::Area2(A,B,V)/dd;
+								 taa[0][j] =  lB*lC;
+								 taa[1][j] =  lC*lA;
+								 taa[2][j] =  lA*lB;
+								 //Real8 xx = V.x-V.y;
+								 //Real8 yy = V.x + V.y;
+								 //cout << " iv " << ss[iV*n] << " == " << (8*xx*xx+yy*yy)
+								 //     << " l = " << lA << " " << lB << " " << lC 
+								 //     << " = " << lA+lB+lC << " " <<  V << " == " << A*lA+B*lB+C*lC << endl;
+
+								 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;
+
+								 bb[j]     =  ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ;
+								}
+							 else
+								{
+								 nbb++;
+								 taa[0][j]=0;
+								 taa[1][j]=0;
+								 taa[2][j]=0;
+								 taa[j][j]=1;
+								 bb[j]=0;
+								}
+							}
+
+						 // resolution of 3x3 lineaire 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);
+
+						 assert(det33);
+						 //	det33=1;
+						 // verif
+						 //	cout << " " << (taa[0][0]*cBC +  taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ;
+						 //	cout << " " << (taa[0][1]*cBC +  taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1];
+						 //	cout << " " << (taa[0][2]*cBC +  taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2] 
+						 //	     << "  -- " ;
+						 //cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB +  llb*llc* cBC + llc*lla*cCA)/det33 
+						 //   << " == " << llf <<  endl;
+						 // computation of the gradient 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);
+						 Real8 coef = 1.0/(3*dd*det33);
+						 Real8 coef2 = 2*coef;
+						 //	cout << " H = " << Hxx << " " << Hyy << " " <<  Hxy/2 << " coef2 = " << coef2 << endl;
+						 Hxx *= coef2;
+						 Hyy *= coef2;
+						 Hxy *= coef2;
+						 //cout << i  << " H = " << 3*Hxx/dd << " " << 3*Hyy/dd << " " <<  3*Hxy/(dd*2) << " nbb = " << nbb << endl;
+						 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;
+							}
+
+						}
+					 else
+						{
+
+						 // 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[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 divide by 2 is done on the metrix 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++,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);
+					 //cout <<iv <<  "  M  = " <<  M <<  " aniso= " << Vp.Aniso() ;
+					 Vp.Abs();
+					 M = Vp;
+					 dxdx[iv] = M.a11;
+					 dxdy[iv] = M.a21;
+					 dydy[iv] = M.a22;
+					 //  cout << " (abs)  iv M  = " <<  M <<  " aniso= " << Vp.Aniso() <<endl;
+					}
+				 else kk++;
+
+
+				// correction of second derivate
+				// 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);
+
+
+					  }
+
+
+				  }
+
+				// constuction  of the metrix from the Hessian dxdx. dxdy,dydy
+
+				Real8 rCutOff=CutOff*absmax;// relative cut off 
+
+				for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )
+				  { // for all vertices 
+					//{
+					//Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
+					// MatVVP2x2 Vp(M);	  
+					//cout << " iv M="<<  M << "  Vp = " << Vp << " aniso  " << Vp.Aniso() << endl;
+					//}
+					MetricIso Miso;
+					// new code to compute ci ---	  
+					Real8 ci ;
+					if (RelativeMetric)
+					  { //   compute the norm of the solution
+						double xx =0,*sfk=sf+k; 
+						for (int ifield=0;ifield<nbfield;ifield++,sfk++)
+						 xx += *sfk* *sfk;	       
+						xx=sqrt(xx);
+						ci = coef2/Max(xx,rCutOff);
+					  }
+					else ci = cnorm;
+
+					// old 
+					//	  Real8 ci = RelativeMetric ? coef2/(Max(Abs(ss[k]),rCutOff)) : cnorm ;
+					//   modif F Hecht 101099
+					Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci,  dydy[iv]*ci);
+					MatVVP2x2 Vp(Miv);
+
+					Vp.Abs();
+					if(power!=1.0) 
+					 Vp.pow(power);
+
+
+
+					h1=Min(h1,Vp.lmin());
+					h2=Max(h2,Vp.lmax());
+
+					Vp.Maxh(hmin);
+					Vp.Minh(hmax);
+
+					rx = Max(rx,Vp.Aniso2());
+
+					Vp.BoundAniso2(coef);
+
+					hn1=Min(hn1,Vp.lmin());
+					hn2=Max(hn2,Vp.lmax());
+					rnx = Max(rnx,Vp.Aniso2());
+
+					Metric MVp(Vp);
+					vertices[iv].m.IntersectWith(MVp);
+				  }// for all vertices 
+				if (verbosity>2)
+				  { 
+					cout << "              Field " << nufield << " of solution " << nusol  << endl;
+					cout << "              before bounding :  Hmin = " << sqrt(1/h2) << " Hmax = " 
+					  << sqrt(1/h1)  << " factor of anisotropy max  = " << sqrt(rx) << endl;
+					cout << "              after  bounding :  Hmin = " << sqrt(1/hn2) << " Hmax = " 
+					  << sqrt(1/hn1)  << " factor of anisotropy max  = " << sqrt(rnx) << endl;
+				  }
+			  } //  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::MaxSubDivision{{{1*/
+	void  Triangles::MaxSubDivision(Real8 maxsubdiv) {
+		long int verbosity=0;
+
+		const  Real8 maxsubdiv2 = maxsubdiv*maxsubdiv;
+		if(verbosity>1)
+		 cout << "  -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv <<   endl  ;
+		// for all the edges 
+		// if the len of the edge is to long 
+		Int4 it,nbchange=0;    
+		Real8 lmax=0;
+		for (it=0;it<nbt;it++)
+		  {
+			Triangle &t=triangles[it];
+			for (int j=0;j<3;j++)
+			  {
+				Triangle &tt = *t.TriangleAdj(j);
+				if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)) 
+				  {
+					Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];
+					Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];
+					R2 AB= (R2) v1-(R2) v0;
+					Metric M = v0;
+					Real8 l = M(AB,AB);
+					lmax = Max(lmax,l);
+					if(l> maxsubdiv2)
+					  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
+						Real8 lc = M(AC,AC);
+						D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
+						D2xD2 Rt1(Rt.inv());
+						D2xD2 D(maxsubdiv2,0,0,lc);
+						D2xD2 MM = Rt1*D*Rt1.t();
+						v0.m =  M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);
+						nbchange++;
+					  }
+					M = v1;
+					l = M(AB,AB);
+					lmax = Max(lmax,l);
+					if(l> maxsubdiv2)
+					  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
+						Real8 lc = M(AC,AC);
+						D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
+						D2xD2 Rt1(Rt.inv());
+						D2xD2 D(maxsubdiv2,0,0,lc);
+						D2xD2  MM = Rt1*D*Rt1.t();
+						v1.m =  M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);
+						nbchange++;
+					  }
+
+
+				  }
+			  }
+		  }
+		if(verbosity>3)
+		 cout << "    Nb of metric change = " << nbchange 
+			<< " Max  of the subdivision of a edges before change  = " << sqrt(lmax) << endl;
+
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::SmoothMetric{{{1*/
+	void Triangles::SmoothMetric(Real8 raisonmax) { 
+		long int verbosity=0;
+
+		if(raisonmax<1.1) return;
+		if(verbosity > 1)
+		 cout << "  -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl;
+		ReMakeTriangleContainingTheVertex();
+		Int4 i,j,kch,kk,ip;
+		Int4 *first_np_or_next_t0 = new Int4[nbv];
+		Int4 *first_np_or_next_t1 = new Int4[nbv];
+		Int4 Head0 =0,Head1=-1;
+		Real8 logseuil= log(raisonmax);
+
+		for(i=0;i<nbv-1;i++)
+		 first_np_or_next_t0[i]=i+1; 
+		first_np_or_next_t0[nbv-1]=-1;// end;
+		for(i=0;i<nbv;i++)
+		 first_np_or_next_t1[i]=-1;
+		kk=0;
+		while (Head0>=0&& kk++<100) {
+			kch=0;
+			for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1)
+			  {  //  pour tous les triangles autour du sommet s
+				// 	cout << kk << " i = " << i << " " << ip << endl;
+				register Triangle * t= vertices[i].t;
+				assert(t);
+				Vertex & vi = vertices[i];
+				TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
+				Vertex *pvj0 = ta.EdgeVertex(0);
+				while (1) {
+					//	  cout << i << " " <<  Number(ta.EdgeVertex(0)) << " "
+					//      << Number(ta.EdgeVertex(1)) << "  ---> " ;
+					ta=Previous(Adj(ta));
+					// cout <<  Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl;
+					assert(vertices+i == ta.EdgeVertex(1));
+					Vertex & vj = *(ta.EdgeVertex(0));
+					if ( &vj ) {
+						j= &vj-vertices;
+						assert(j>=0 && j < nbv);
+						R2 Aij = (R2) vj - (R2) vi;
+						Real8 ll =  Norme2(Aij);
+						if (0) {  
+							Real8 hi = ll/vi.m(Aij);
+							Real8 hj = ll/vj.m(Aij);
+							if(hi < hj)
+							  {
+								Real8 dh=(hj-hi)/ll;
+								//cout << " dh = " << dh << endl;
+								if (dh>logseuil) {
+									vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
+									if(first_np_or_next_t1[j]<0)
+									 kch++,first_np_or_next_t1[j]=Head1,Head1=j;
+								}
+							  }
+						} 
+						else
+						  {
+							Real8 li = vi.m(Aij);
+							//Real8 lj = vj.m(Aij);
+							//		if ( i == 2 || j == 2)
+							//  cout << " inter " << i << " " << j << " " << ((1 +logseuil*li)) <<  endl;
+							if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
+							 //if( vj.m.IntersectWith(vi.m*(lj/li/(1 +logseuil*lj))) )
+							 if(first_np_or_next_t1[j]<0) // if the metrix change 
+							  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
+						  }
+					}
+					if  ( &vj ==  pvj0 ) break;
+				}
+			  }
+			Head0 = Head1;
+			Head1 = -1;
+			Exchange(first_np_or_next_t0,first_np_or_next_t1);
+			if(verbosity>5)
+			 cout << "     Iteration " << kk << " Nb de  vertices with change  " << kch << endl;
+		}
+		if(verbosity>2 && verbosity < 5) 
+		 cout << "    Nb of Loop " << kch << endl;
+		delete [] first_np_or_next_t0;
+		delete [] first_np_or_next_t1;
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::NearestVertex{{{1*/
+	Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j) {
+		return  quadtree->NearestVertex(i,j); 
+	} 
+	/*}}}1*/
+	/*FUNCTION Triangles::PreInit{{{1*/
+	void Triangles::PreInit(Int4 inbvx,char *fname) {
+		long int verbosity=0;
+
+		srand(19999999);
+		OnDisk =0;
+		NbRef=0;
+		//  allocGeometry=0;
+		identity=0;
+		NbOfTriangleSearchFind =0;
+		NbOfSwapTriangle =0;
+		nbiv=0;
+		nbv=0;
+		nbvx=inbvx;
+		nbt=0;
+		NbOfQuad = 0;
+		nbtx=2*inbvx-2;
+		NbSubDomains=0;
+		NbVertexOnBThVertex=0;
+		NbVertexOnBThEdge=0;
+		VertexOnBThVertex=0;
+		VertexOnBThEdge=0;
+
+		NbCrackedVertices=0;
+		NbCrackedEdges =0;
+		CrackedEdges  =0;  
+		nbe = 0; 
+		name = fname ;
+
+		if (inbvx) {
+			vertices=new Vertex[nbvx];
+			assert(vertices);
+			ordre=new (Vertex* [nbvx]);
+			assert(ordre);
+			triangles=new Triangle[nbtx];
+			assert(triangles);}
+		else {
+			vertices=0;
+			ordre=0;
+			triangles=0;
+			nbtx=0;
+		}
+		if ( name || inbvx)
+		  {
+			time_t timer =time(0);
+			char buf[70];     
+			strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer));
+			counter++; 
+			char countbuf[30];   
+			sprintf(countbuf,"%d",counter);
+			int lg =0 ;
+			if (&BTh != this && BTh.name)
+			 lg = strlen(BTh.name)+4;
+			identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2  + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)];
+			identity[0]=0;
+			if (lg)
+			 strcat(strcat(strcat(identity,"B="),BTh.name),", ");
+
+			if (Gh.name)
+			 strcat(strcat(identity,"G="),Gh.name);
+			strcat(strcat(identity,";"),countbuf);
+			strcat(identity,buf);
+			// cout << "New MAILLAGE "<< identity << endl;
+		  } 
+
+		quadtree=0;
+		//  edgescomponante=0;
+		edges=0;
+		VerticesOnGeomVertex=0;
+		VerticesOnGeomEdge=0;
+		NbVerticesOnGeomVertex=0;
+		NbVerticesOnGeomEdge=0;
+		//  nbMaxIntersectionTriangles=0;
+		//  lIntTria;
+		subdomains=0;
+		NbSubDomains=0;
+		//  Meshbegin = vertices;
+		//  Meshend  = vertices + nbvx;
+		if (verbosity>98) 
+		 cout << "Triangles::PreInit() " << nbvx << " " << nbtx 
+			<< " " << vertices 
+			<< " " << ordre << " " <<  triangles <<endl;
+	}
+	/*}}}1*/
 
 } // end of namespace bamg