Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2857)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2858)
@@ -69,6 +69,5 @@
 	// some verification
 	if ( maxsubdiv > boundmaxsubdiv || maxsubdiv <= 1.0){
-		cerr << " -maxsubdiv " << maxsubdiv << " is not in ]1,"<< boundmaxsubdiv << "]" << endl;
-		exit(3);
+		throw ErrorException(__FUNCT__,exprintf("maxsubdiv %g is should be in ]1,%g]",maxsubdiv,boundmaxsubdiv));
 	}
 	if (iso) anisomax=1;
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2857)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2858)
@@ -815,5 +815,4 @@
   void NewPoints( Triangles &,int KeepBackVertex =1 );
   Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; 
-  void NewPointsOld( Triangles & );
   void NewPoints(int KeepBackVertex=1){ NewPoints(*this,KeepBackVertex);}
   void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2857)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2858)
@@ -179,10 +179,8 @@
 
 
-inline   void  MetricAnIso::Box(Real4 &hx,Real4 &hy) const 
-{
+inline   void  MetricAnIso::Box(Real4 &hx,Real4 &hy) const {
   Real8 d=  a11*a22-a21*a21;
   hx = sqrt(a22/d);
   hy = sqrt(a11/d);
-  //  cerr << " hx = " << hx << " hy =  " << hy << endl;
 }
 
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2857)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2858)
@@ -708,5 +708,4 @@
 				GeometricalEdge & e1=edges[i];
 				GeometricalEdge & e2=*e1.link;
-				cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " <<  e1[1].The() << " " << e2[1].The() << endl;
 				if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() )
 				  {
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2857)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2858)
@@ -1546,6 +1546,5 @@
 						  } // tt 
 						else
-						 cerr <<endl <<  " Bug " <<i<< " " << j << " t=" << t << endl;
-
+						 throw ErrorException(__FUNCT__,exprintf("Bug..."));
 					  } // ke<0	       
 					else
@@ -2349,336 +2348,11 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::NewPointsOld{{{1*/
-	void  Triangles::NewPointsOld(Triangles & Bh) {
-		long int verbosity=2;
-		Real8 seuil= 0.7 ;// for two neart point 
-		if (verbosity>1) 
-		 cout << " begin :  Triangles::NewPointsOld " << endl;
-		Int4 i,k;
-		int j;
-		Int4 BeginNewPoint[3];
-		Int4 EndNewPoint[3];
-		int step[3];
-		Int4 *first_np_or_next_t = new Int4[nbtx];
-		Int4 ColorEdge[3];
-		Int4 color=-1;
-		Triangle *t;
-		// generation of the list of next Triangle 
-		// at 1 time we test all the triangles
-		Int4 Headt =0,next_t;
-		for(i=0;i<nbt;i++)
-		 first_np_or_next_t[i]=-(i+1);
-		// end list i >= nbt 
-		// the list of test triangle is 
-		// the next Triangle on i is  -first_np_or_next_t[i]
-		Int4 nbtold,nbvold;
-
-		// Big loop 
-		do {
-			nbtold = nbt;
-			nbvold = nbv;
-			// default size of  IntersectionTriangle
-
-			i=Headt;
-			next_t=-first_np_or_next_t[i];
-			for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) { // for each triangle  t
-				// we can change first_np_or_next_t[i]
-				if (i<0 || i >= nbt ){
-					throw ErrorException(__FUNCT__,exprintf("i<0 || i >= nbt"));
-				}
-				first_np_or_next_t[i] = nbv; // to save the fist new point of triangle
-				for(j=0;j<3;j++)
-				  { // for each edge 
-					TriangleAdjacent tj(t,j);
-					// color++;// the color is 3i+j
-					color = 3*i + j ;;
-					ColorEdge[j]=color;
-					BeginNewPoint[j]=nbv;
-					EndNewPoint[j]=nbv-1;
-					step[j]=1;// right sens 
-					Vertex & vA = * tj.EdgeVertex(0);
-					Vertex & vB = * tj.EdgeVertex(1);
-
-					if (!t->link) continue;// boundary
-					if (t->det <0) continue;
-					if (t->Locked(j)) continue;
-
-					TriangleAdjacent tadjj = t->Adj(j);	  
-					Triangle * ta= tadjj;
-
-					if (ta->det <0) continue;	  
-
-					R2 A = vA;
-					R2 B = vB;
-
-					k=Number(ta);
-					// the 2 opposite vertices 
-					const Vertex & vC1  =  *tj.OppositeVertex();
-					const Vertex & vC2 = *tadjj.OppositeVertex();
-
-					if(first_np_or_next_t[k]>0) { // this edge is done before 
-						// find the color of the edge and begin , end of newpoint
-						register int kk = t->NuEdgeTriangleAdj(j);
-						if ((*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])){
-							throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])"));
-						}
-						if ((*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])){
-							throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])"));
-						}
-						register Int4 kolor =3*k + kk;
-						ColorEdge[j]=kolor;
-						register Int4 kkk= 1;
-						step[j]=-1;// other sens 
-						BeginNewPoint[j]=0;
-						EndNewPoint[j]=-1; // empty list          
-						for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++) 
-						 if (vertices[iv].color > kolor) 
-						  break; // the color is passed            
-						 else if (vertices[iv].color == kolor) {
-							 EndNewPoint[j]=iv; 
-							 if (kkk) // one time test
-							  kkk=0,BeginNewPoint[j]=iv;}
-							 continue; // next edge of the triangle 
-					} // end  if( k < i) 
-
-
-					const Int4 NbvOld = nbv;
-					lIntTria.SplitEdge(Bh,A,B);
-					// Int4 NbvNp =
-					lIntTria.NewPoints(vertices,nbv,nbvx);
-					Int4 nbvNew=nbv;
-					nbv = NbvOld;
-					for (Int4 iv=NbvOld;iv<nbvNew;iv++) {
-						vertices[nbv].color = color;
-						vertices[nbv].ReferenceNumber=nbv;// circular link
-						R2 C =  vertices[iv].r;
-						vertices[nbv].r =  C;
-						vertices[nbv].m =  vertices[iv].m ;
-						// test if the new point is not to close to the 2 opposite vertex
-						R2 CC1 = C-vC1 , CC2 = C-vC2;
-						if (   (  (vC1.m(CC1) + vertices[nbv].m(CC1)) >  seuil)
-									&& (  (vC2.m(CC2) + vertices[nbv].m(CC2)) >  seuil) )
-						 nbv++;
-					}
-
-					EndNewPoint[j] = nbv-1;
-				  } // end loop for each edge 
-
-				if (!t->link) continue;// boundary
-				if (t->det<=0) continue;// outside 
-				//      continue;
-				for(int j0=0;j0<3;j0++)
-				 for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++)
-					{
-					 // find the neart  point one the opposite edge 
-					 // to compute i1
-					 Vertex & vi0 = vertices[i0];
-					 int kstack = 0;
-					 Int4 stack[10];
-					 //   Int4 savRef[10];
-					 int j1 = j0; 
-					 while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge
-						 // computation of the intersection of edge j1 and DOrto
-						 // take the good sens 
-
-						 if (BeginNewPoint[j1]> EndNewPoint[j1])
-						  continue; // 
-						 else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) {
-							 for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++)
-							  stack[kstack++] = ii1;
-							 continue;}
-
-
-							 int k0,k1;
-							 if (step[j1]<0) k0=1,k1=0; // reverse
-							 else k0=0,k1=1; 
-							 R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]];
-							 R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]];
-							 R2 D = V11-V10;
-							 Real8 c0 =  vi0.m(D,(R2) vi0);
-
-							 Real8 c10 =  vi0.m(D,V10);
-							 Real8 c11 =  vi0.m(D,V11);
-
-							 Real8 s;
-							 //cout << " --i0 = " << i0  << D  << V10 << V11 << endl ;
-							 //cout << "   c10 " <<  c10 << " c0 " << c0 << " c11 " << c11 << endl;
-							 if (( c10 < c0 ) && (c0 < c11)) 
-							  s = (c11-c0)/(c11-c10);
-							 else if  (( c11 < c0 ) && (c0 < c10)) 
-							  s = (c11-c0) /(c11-c10);
-							 else break;
-							 R2 VP = V10*s + V11*(1-s);
-							 int sss = (c11-c10) >0 ? 1 : -1;
-							 // find the 2 point by dichotomie
-							 //cout << "   t =" << Number(t) << " c0 " << c0  ;
-							 Int4 ii0 =  BeginNewPoint[j1];
-							 Int4 ii1 =  EndNewPoint[j1];	     
-							 Real8 ciii=-1,cii0=-1,cii1=-1  ;
-							 if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 )  
-							  stack[kstack++] = ii0;//cout << " add+0 " << ii0;
-							 else if ( sss * ((cii1=  vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 )  
-							  stack[kstack++] = ii1;//cout << " add+1 " << ii1;
-							 else {
-								 while ((ii1-ii0)> 1) {
-									 Int4 iii = (ii0+ii1)/2;
-									 ciii = vi0.m( D ,(R2) vertices[iii]);
-									 //cout << " (iii = " << iii << " " <<  ciii << ") ";
-									 if ( sss * (ciii - c0) <0 )  ii0 = iii;
-									 else ii1 = iii;}	        	      
-									 stack[kstack++] = ii0;// cout << " add0 " << ii0;
-									 if (ii1 != ii0)  stack[kstack++] = ii1;//cout << " add1 " << ii1;
-							 }
-							 if (kstack >5) // bug ?
-							  cout << "NewPoints: bug????? " << kstack << " stack  " << stack[kstack]<< endl;
-					 }
-
-					 stack[kstack++] = -1; // to stop
-					 Int4 i1;
-					 kstack =0; 
-					 while( (i1=stack[kstack++]) >= 0){ // the two parameter is i0 and i1 
-						 if (i1 >= nbv || i1 < 0){
-							 throw ErrorException(__FUNCT__,exprintf("i1 >= nbv || i1 < 0"));
-						 }
-						 if (i0 >= nbv || i0 < 0){
-							 throw ErrorException(__FUNCT__,exprintf("i0 >= nbv || i0 < 0"));
-						 }
-						 if (i1 == i0){
-							 throw ErrorException(__FUNCT__,exprintf("i1 == i0"));
-						 }
-						 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0];
-						 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01));
-						 if (i0< nbvold || i1<nbvold){
-							 throw ErrorException(__FUNCT__,exprintf("i0< nbvold || i1<nbvold"));
-						 }
-						 if (d01 == 0) 
-						  break; 
-						 if ( d01 < seuil) 
-						  if (i1<nbvold) {
-							  // remove all the points i0;
-							  register Int4 ip,ipp;
-							  for  (ip=i0;i0 != (ipp = vertices[ip].ReferenceNumber);ip=ipp)
-								vertices[ip].ReferenceNumber = -1;// mark remove
-							  vertices[ip].ReferenceNumber = -1;// mark remove
-						  }
-						  else {
-							  // remove on of two points
-							  register Int4 ip0, ip1, ipp0,ipp1;
-							  register int kk0=1,kk1=1;
-							  // count the number of common points to compute weight w0,w1
-							  for  (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) kk0++;
-							  for  (ip1=i1;i1 != (ipp1 = vertices[ip1].ReferenceNumber);ip1=ipp1) kk1++;
-
-							  register Real8 w0 = ((Real8) kk0)/(kk0+kk1);
-							  register Real8 w1 = ((Real8) kk1)/(kk0+kk1);
-
-							  // make a circular link
-							  Exchange(vertices[i0].ReferenceNumber,vertices[i1].ReferenceNumber);
-							  // the new coordinate 
-							  R2 C //= vertices[i0] ;
-							  =  vertices[i0].r *w0 + vertices[i1].r* w1;
-
-							  // update the new point points of the list 
-							  for  (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0)
-								vertices[ip0].r = C;	      
-							  vertices[ip0].r = C;
-						  }
-						}
-					} // for (i0= ....
-			  }// for triangle   
-
-			// remove of all the double points
-
-			Int4 ip,ipp,kkk=nbvold;
-			for (i=nbvold;i<nbv;i++) 
-			 if (vertices[i].ReferenceNumber>=0) {// good points
-				 //  cout <<" i = " << i ;
-				 for  (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp)
-				  vertices[ip].ReferenceNumber = -1;// mark remove
-				 vertices[ip].ReferenceNumber = -1;// mark remove
-				 // cout << i << " ---> " << kkk << endl;        
-				 vertices[kkk] = vertices[i];
-				 vertices[kkk].i = toI2(vertices[kkk].r);
-				 vertices[kkk++].ReferenceNumber = 0;
-
-			 } 
-
-			// insertion part --- 
-
-			const Int4 nbvnew = kkk-nbvold;
-
-			cout << "    Remove " << nbv - kkk  << " to close  vertex " ;
-			cout << " and   Insert the " <<nbvnew<< " new points " << endl;  
-			nbv = kkk;
-			Int4 NbSwap=0;
-			Icoor2 dete[3];
-
-			// construction d'un ordre aleatoire 
-			if (! nbvnew) 
-			 break; 
-			if (nbvnew) {
-				const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv)  ;
-				Int4 k3 = rand()%nbvnew ; 
-				for (Int4 is3=0; is3<nbvnew; is3++) 
-				 ordre[nbvold+is3]= &vertices[nbvold +(k3 = (k3 + PrimeNumber)% nbvnew)];
-
-				for (i=nbvold;i<nbv;i++) 
-				  { Vertex * vi = ordre[i];
-					Triangle *tcvi = FindTriangleContening(vi->i,dete);
-					//     Vertex * nv =  quadtree->NearestVertex(vi->i.x,vi->i.y);
-					//      cout << " Neart Vertex of "  << Number(vi)<< vi->i << " is " 
-					//   << Number(nv) << nv->i  << endl;
-					// Int4  kt = Number(tcvi);
-					// 
-
-					quadtree->Add(*vi); //
-					if (tcvi->det<0){// internal
-						throw ErrorException(__FUNCT__,exprintf("tcvi->det<0"));
-					}
-					Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);          
-				  }  
-			}
-			cout << " Nb swap = " << NbSwap << " after " ;
-
-			for (i=nbvold;i<nbv;i++) 
-			 NbSwap += vertices[i].Optim(1);  
-			cout << NbSwap << endl;
-
-			for (i=nbtold;i<nbt;i++)
-			 first_np_or_next_t[i]=1;
-
-			Headt = nbt; // empty list 
-			for (i=nbvold;i<nbv;i++) 
-			  { // for all the triangle contening the vertex i
-				Vertex * s  = vertices + i;
-				TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
-				Triangle * tbegin= (Triangle*) ta;
-				Int4 kt;
-				do { 
-					kt = Number((Triangle*) ta);
-					if (first_np_or_next_t[kt]>0) 
-					 first_np_or_next_t[kt]=-Headt,Headt=kt;
-					if (ta.EdgeVertex(0) != s){
-						throw ErrorException(__FUNCT__,exprintf("ta.EdgeVertex(0) != s"));
-					}
-					ta = Next(Adj(ta));
-				} while ( (tbegin != (Triangle*) ta)); 
-			  }
-
-
-		} while (nbv!=nbvold);
-		delete []  first_np_or_next_t;
-		cout << " end :  Triangles::NewPoints old  nbv=" << nbv << endl;
-
-	}
-	/*}}}1*/
 	/*FUNCTION Triangles::Insert{{{1*/
 	void Triangles::Insert() {
 		long int verbosity=2;
-		if (verbosity>2) cout << "  -- Insert initial " << nbv << " vertices " << endl ;
+		if (verbosity>2) printf("   Insert initial %i vertices\n",nbv);
 		Triangles * OldCurrentTh =CurrentTh;
 
 		CurrentTh=this;
-		double time0=CPUtime(),time1,time2,time3;
 		SetIntCoor();
 		Int4 i;
@@ -2737,14 +2411,12 @@
 		quadtree->Add(*v1);
 
-		// on ajoute les sommets un Ò un 
+		//the vertices are added one by one
 		Int4 NbSwap=0;
 
-		time1=CPUtime();
-
-		if (verbosity>3) cout << "  -- Begin of insertion process " << endl;
+
+		if (verbosity>3) printf("   Begining of insertion process...\n");
 
 		for (Int4 icount=2; icount<nbv; icount++) {
 			Vertex *vi  = ordre[icount];
-			//    cout << " Insert " << Number(vi) << endl;
 			Icoor2 dete[3];
 			Triangle *tcvi = FindTriangleContening(vi->i,dete);
@@ -2752,13 +2424,11 @@
 			Add(*vi,tcvi,dete);
 			NbSwap += vi->Optim(1,0);
-
-		}// fin de boucle en icount
-		time2=CPUtime();
-		if (verbosity>3) 
-		 cout << "    NbSwap of insertion " <<    NbSwap 
-			<< " NbSwap/Nbv " <<  (float) NbSwap / (float) nbv 
-			<< " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv " 
-			<< (float)NbUnSwap /(float) nbv 
-			<<endl;
+		}
+		if (verbosity>3) {
+			printf("      NbSwap of insertion: %i\n",NbSwap);
+			printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
+			printf("      NbUnSwap:            %i\n",NbUnSwap);
+			printf("      NbUnSwap/nbv         %i\n",NbUnSwap/nbv);
+		}
 		NbUnSwap = 0;
 		// construction d'un ordre aleatoire 
@@ -2770,18 +2440,15 @@
 		 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
 
-		double timeloop = time2 ;
 		for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++) 
 		  {
-			double time000 = timeloop;
 			Int4  NbSwap = 0;
 			for (int is1=0; is1<nbv; is1++) 
 			 NbSwap += ordre[is1]->Optim(0,0);
-			timeloop = CPUtime();
-			if (verbosity>3) 
-			 cout << "    Optim Loop "<<Nbloop<<" NbSwap: " <<  NbSwap 
-				<< " NbSwap/Nbv " 	   <<  (float) NbSwap / (float) nbv 
-				<< " CPU=" << timeloop - time000 << "  s, " 
-				<< " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv  
-				<<  endl;
+			if (verbosity>3) {
+				printf("      Optim Loop: %i\n",Nbloop);
+				printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
+				printf("      NbUnSwap:            %i\n",NbUnSwap);
+				printf("      NbUnSwap/nbv         %i\n",NbUnSwap/nbv);
+			}
 			NbUnSwap = 0;
 			if(!NbSwap) break;
@@ -2790,11 +2457,4 @@
 		// because we break the TriangleContainingTheVertex
 #endif
-		time3=CPUtime();
-		if (verbosity>4) 
-		 cout << "    init " << time1 - time0 << " initialisation,  " 
-			<< time2 - time1 << "s, insert point  " 
-			<< time3 -time2 << "s, optim " << endl
-			<< "     Init Total Cpu Time = " << time3 - time0 << "s " << endl;
-
 		CurrentTh=OldCurrentTh;
 	}
@@ -2803,11 +2463,10 @@
 	void Triangles::ForceBoundary() {
 		long int verbosity=2;
-		if (verbosity > 2)
-		 cout << "  -- ForceBoundary  nb of edge " << nbe << endl;
+		if (verbosity > 2) printf("   ForceBoundary  nb of edge: %i\n",nbe);
 		int k=0;
 		Int4  nbfe=0,nbswp=0,Nbswap=0;
-		for (Int4 t = 0; t < nbt; t++)  
-		 if (!triangles[t].det)
-		  k++,cerr << " det T" << t << " = " << 0 << endl;
+		for (Int4 t = 0; t < nbt; t++){
+			if (!triangles[t].det) k++;
+		}
 		if (k!=0) {
 			throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k));
@@ -2822,15 +2481,6 @@
 				else Nbswap += nbswp;
 				if (nbswp) nbfe++;
-				if ( nbswp < 0 && k < 5)
-				  {
-					cerr << " Missing  Edge " << i << " v0 =  " << Number(edges[i][0]) << edges[i][0].r
-					  <<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << "  " << (Triangle *) ta ;
-					if(ta.t)
-					  {
-						Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1);  
-						cerr << " crossing with  [" << Number(aa) << ", " << Number(bb) << "]\n";
-					  }
-					else cerr << endl;
-
+				if ( nbswp < 0 && k < 5){
+					throw ErrorException(__FUNCT__,exprintf("Missing  Edge %i, v0=%i,v1=%i",i ,Number(edges[i][0]),Number(edges[i][1])));
 				  }
 				if ( nbswp >=0 && edges[i].on->Cracked())
@@ -2844,7 +2494,5 @@
 			for (Int4 j=0;j<nbv;j++)
 			 Nbswap +=  vertices[j].Optim(1,0);
-			if (verbosity > 3)
-			 cout << "     Nb of inforced edge = " << nbfe << " Nb of Swap " << Nbswap << endl;
-
+			if (verbosity > 3) printf("      number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 
 	}
 	/*}}}1*/
@@ -2853,19 +2501,8 @@
 		long int verbosity=0;
 
-		if (verbosity >2)
-		  {
-			if (OutSide)
-			 cout << "  -- Find all external sub-domain ";	
-			else
-			 cout << "  -- Find all internal sub-domain ";
-			if(verbosity>99)
-			  {
-
-				for(int i=0;i<nbt;++i)
-				 cout << i<< " " << triangles[i] << endl;
-			  }
-
+		if (verbosity >2){
+			if (OutSide) printf("   Find all external sub-domain\n"); 
+			else printf("   Find all internal sub-domain\n");
 		  }
-		// if (verbosity > 4) cout << " OutSide=" << OutSide << endl;
 		short * HeapArete = new short[nbt];
 		Triangle  **  HeapTriangle = new Triangle*  [nbt];
@@ -2919,12 +2556,12 @@
 					// infini triangle 
 					NbSubDomTot --;
-					//  cout << " triangle infini " << it << triangles[it] << endl;
 					t=&triangles[it];
 					NbOutT--;  // on fait un coup de trop. 
-					while  (t){ // cout << Number(t) << " " << endl;
+					while  (t){
 						NbOutT++;
 						t1=t;
 						t=t->link;
-						t1->link=0;}//while (t)
+						t1->link=0;
+					}
 				  }
 			  }   
@@ -2960,6 +2597,5 @@
 						mark[it]=k;
 						subdomains[k].head = t1;
-						// cout << " New -- " << Number(t1) << " " << it << endl;
-						do {// cout << " k " << k << " " << Number(t) << endl;
+						do {
 							mark[Number(t)]=k;
 							t=t->link;
@@ -2985,6 +2621,4 @@
 								  Int4 kr = mark[it];
 								  if(kr !=kl) {
-									  //cout << kl << " " << kr << " rl "  << subdomains[kl].ref
-									  // << " rr " << subdomains[kr].ref ;
 									  if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0)
 										nbk--,subdomains[kr].ref=subdomains[kl].ref-1;
@@ -2995,30 +2629,22 @@
 									  if(kl<0 && kr >=0 && subdomains[kr].ref>=0)
 										nbk--,subdomains[kr].ref=-1;
-									  //   cout << " after \t "   
-									  //	 << kl << subdomains[kl].ref << " rr " << kr 
-									  // << subdomains[kr].ref << endl;
 								  }
 								 }
-							//  cout << subdomains[0].ref << subdomains[1].ref << endl;
 							Int4  j=0;
 							for ( i=0;i<NbSubDomains;i++)
 							 if((-subdomains[i].ref) %2) { // good 
-								 //cout << " sudom ok  = " << i << " " << subdomains[i].ref
-								 // << " " << (-subdomains[i].ref) %2 << endl;
 								 if(i != j) 
 								  Exchange(subdomains[i],subdomains[j]);
 								 j++;}
-							 else
-								{ //cout << " remove sub domain " << i << endl;
+							 else{ 
 								 t= subdomains[i].head;
-								 while  (t){// cout << Number(t) << " " << endl;
+								 while (t){
 									 NbOutT++;
 									 t1=t;
 									 t=t->link;
-									 t1->link=0;}//while (t)
+									 t1->link=0;
+								 }//while (t)
 								}
-
-							if(verbosity>4)
-							 cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl;
+							if(verbosity>4) printf("      Number of removes subdomains (OutSideMesh) = %i\n",NbSubDomains-j);
 							NbSubDomains=j;
 						  }
@@ -3034,6 +2660,4 @@
 				 subdomains = new SubDomain[ Gh.NbSubDomains];
 				NbSubDomains =Gh.NbSubDomains;
-				if(verbosity>4)
-				 cout << "     find the " << NbSubDomains << " sub domain " << endl;
 				Int4 err=0;
 				ReMakeTriangleContainingTheVertex();
@@ -3058,6 +2682,4 @@
 					int sens = Gh.subdomains[i].sens;
 					// test if ge and e is in the same sens 
-					//	cout << " geom edge = " <<  Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref 
-					//     << " ref edge =" << eg.ref << " sens " << sens ;
 					if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) sens = -sens ;
 					subdomains[i].sens = sens;
@@ -3073,5 +2695,4 @@
 							throw ErrorException(__FUNCT__,exprintf("v0 != ta.EdgeVertex(1)"));
 						}
-						//	 cout << " recherche " << Number( ta.EdgeVertex(0)) << endl;
 						if (ta.EdgeVertex(0) == v1) { // ok we find the edge
 							if (sens>0)  
@@ -3079,5 +2700,4 @@
 							else 
 							 subdomains[i].head=t=ta;
-							//cout << "      triangle  =" << Number(t) << " = " << (*t)[0].r <<  (*t)[1].r <<  (*t)[2].r << endl;
 							if(t<triangles || t >= triangles+nbt || t->det < 0 || t->link == 0) {
 								throw ErrorException(__FUNCT__,exprintf("bad definition of SubSomain %i",i));
@@ -3085,9 +2705,4 @@
 							Int4 it = Number(t);
 							if (mark[it] >=0) {
-								if(verbosity>10){
-									cerr << "     Warning: the sub domain " << i << " ref = " << subdomains[i].ref 
-									  << " is previouly defined with "  <<mark[it] << " ref = " << subdomains[mark[it]].ref
-									  << " skip this def " << endl;
-								}
 								break;
 							}
@@ -3106,6 +2721,4 @@
 								tt=tt->link;
 							  } while (tt!=t);
-							if(verbosity>7)
-							 cout << "     Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl;
 							break;
 						}
@@ -3118,6 +2731,5 @@
 
 				if (inew < NbSubDomains) {
-					if (verbosity>5) 
-					 cout << "     Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl;
+					if (verbosity>5) printf("WARNING: %i SubDomains are being removed\n",NbSubDomains-inew);
 					NbSubDomains=inew;}
 
@@ -3133,10 +2745,4 @@
 			for (it=0;it<nbt;it++) 
 			 if(!triangles[it].link)  NbOutT++;
-			if (verbosity> 4)
-			 cout << "    " ;
-			if (verbosity> 2)
-			 cout << " Nb of Sub borned Domain  = " <<  NbSubDomTot << " NbOutTriangles = " << NbOutT <<endl;
-
-
 	}
 	/*}}}1*/
@@ -3236,6 +2842,4 @@
 			while (t0 != (t=t->link));
 		  }
-		if (verbosity>9)
-		 cout << " number of inside triangles " << k << " nbt = " << nbt << endl;
 		// take is same numbering if possible    
 		if(justcompress)
@@ -3300,11 +2904,5 @@
 			while (t0 != (t=t->link));
 		  }
-		//  NbOutT = nbt - k;
-		if (verbosity>5) 
-		 cout << " Nb of Sub Domain =" << NbSubDomains  << " Nb of In Triangles " << k 
-			<< " Nbt = " << nbt << " Out Triangles = " << nbt - k <<  endl;
-
 		return k;   
-
 	}
 	/*}}}1*/
@@ -3375,5 +2973,4 @@
 				Gh[i].to = vertices + nbv;// save Geom -> Th
 				VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
-				// cout << "--------- "  <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
 				nbv++;
 			}
@@ -3400,11 +2997,4 @@
 			Gh.UnMarkEdges();	
 			int bfind=0;
-			/*
-				cout << " nb curves = " << Gh.NbOfCurves << endl;
-				for(int i=0;i<Gh.NbOfCurves ;i++)
-				{
-				cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb 
-				<< "  end e= " << Gh.Number(Gh.curves[i].ee) << " k=" << Gh.curves[i].ke << endl;
-				}*/
 			for (int i=0;i<Gh.NbOfCurves;i++)
 			  {
@@ -3420,7 +3010,4 @@
 					 // a begin of curve 
 					 int nc = ei.on->CurveNumber;
-
-					 //cout << "curve " <<  nc << " v " << Gh.Number((GeometricalVertex *) *ei[je].on) << " "
-					 //     << " e "  << " " << Gh.Number(ei.on) << " vc " << Gh.Number((*Gh.curves[nc].be)[Gh.curves[nc].kb]) << endl;
 					 if(
 								 ei.on==Gh.curves[nc].be    && 
@@ -3428,5 +3015,4 @@
 						)     
 						{ 
-						 // cout << " find " << endl;
 						 bcurve[nc]=iedge*2+je;
 						 bfind++;	
@@ -3449,13 +3035,4 @@
 			Int4 iedge;
 			Gh.UnMarkEdges();	
-			/*   add Curve loop  FH    
-			// find a starting back groud edges to walk 
-			for (iedge=0;iedge<BTh.nbe;iedge++) {      
-			Edge & ei = BTh.edges[iedge];
-			for(int jedge=0;jedge<2;jedge++) // for the 2 extremites
-			if (!ei.on->Mark() && ei[jedge].on->IsRequiredVertex() )
-			{
-			*/  
-			// new code FH 2004 
 			Real8 L=0;
 			for (int icurve=0;icurve<Gh.NbOfCurves;icurve++)
@@ -3472,6 +3049,4 @@
 				Int4 NbCreatePointOnCurve=0;// Nb of new points on curve     (phase==1) 
 
-
-				//    cout.precision(16);
 				for(int phase=0;phase<=step;phase++) 
 				  {
@@ -3564,13 +3139,7 @@
 									 A1->ReferenceNumber = eeequi.ref;
 									 A1->DirOfSearch =NoDirOfSearch;
-									 //cout << icurveequi << " " << i << " " <<  *A1 << endl;
 									 e->on = ongequi;
 									 e->v[0]=  A0;
 									 e->v[1]=  A1;
-									 if(verbosity>99)
-									  cout << i << "+ New P "<< nbv-1 << " "  <<sNew<< " L0=" << L0 
-										 << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= "  
-										 << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl;
-
 									 e->ref = eeequi.ref;
 									 e->adj[0]=PreviousNewEdge;
@@ -3581,6 +3150,4 @@
 									 A0=A1;
 									 sNew += Lstep;
-									 //   cout << " sNew = " << sNew << " L = " << L 
-									 //        << "  ------" <<NbCreatePointOnCurve << " " << i <<  endl;
 									 if (++i== NbCreatePointOnCurve) break;
 								 }
@@ -3590,5 +3157,4 @@
 								 throw ErrorException(__FUNCT__,exprintf("ee.on->CurveNumber!=ei.on->CurveNumber"));
 							 }
-							 if(verbosity>98) cout <<  BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() <<  endl;
 							 if ( ee[k1].on->IsRequiredVertex()) {
 								 if (!eeequi[k1equi].on->IsRequiredVertex()){
@@ -3614,7 +3180,4 @@
 						  {
 							Edge *e = edges + nbe++;
-							if (verbosity>10) 
-							 cout << " Fin curve A1" << *A1 << " " << icurve << " <=> " << icurveequi <<"-----" <<
-								NbCreatePointOnCurve << " == " <<i<<endl;
 							e->on  = ongequi;
 							e->v[0]=  A0;
@@ -3645,7 +3208,4 @@
 							NbOfNewPoints += NbCreatePointOnCurve;
 						  }
-						if(verbosity>5)
-						 cout << icurve << " NbSegOnCurve = " <<  NbSegOnCurve << " Lstep=" 
-							<< Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve <<endl;
 						// do'nt 
 						//  if(NbCreatePointOnCurve<1) break;
@@ -3719,5 +3279,4 @@
 			 Gh[i].to = vertices + nbv;// save Geom -> Th
 			 VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);
-			 //  cout << "--------- "  <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
 			 nbv++;
 		 }
@@ -3738,5 +3297,4 @@
 			nbe = 0;
 			Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
-			//  cout <<  "  -------------- step =" << step << endl;
 			Gh.UnMarkEdges();	
 			NbOfCurves = 0;
@@ -3747,5 +3305,4 @@
 				  if (!ei.Mark() && ei[j].Required()) { 
 					  // warning ei.Mark() can be change in loop for(j=0;j<2;j++) 
-					  //  cout << " New curve = " << NbOfCurves << endl;
 					  Int4 nbvend  =0;
 
@@ -3794,5 +3351,4 @@
 							  va = a->to;
 							  e->SetMark();
-							  //  cout << " New curve " ;
 
 							  // if SameGeo  We have go in the background geometry 
@@ -3818,6 +3374,4 @@
 									  A= a->r;
 									  Metric MAs =MA,MBs;
-									  // cout << " lSubEdge old=" << ledge 
-									  //      << " new " << A << MA << endl;
 									  ledge = 0; 
 									  Real8 x =0, xstep= 1. /  NbSubEdge;
@@ -3828,8 +3382,5 @@
 										  AB = A-B;
 										  lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);
-										  // cout << "     " << lSubEdge[kk] << " x " << x  
-										  //      << " " << A << B << MA << MB<< endl ;
 									  }
-									  //  cout << endl;
 								  }
 
@@ -3884,7 +3435,4 @@
 								  lcurve = lcurveb;
 								  e->SetMark();
-								  // cout << e-Gh.edges << ", " << k << " " 
-								  //      <<(*e)[k].r <<" " <<(*e)[1-k].r <<" " 
-								  //      << lcurve<< ";; " ;                          
 								  a=b;
 								  if (b->Required() ) break;
@@ -3897,5 +3445,4 @@
 								 }// for(;;)
 							  vb = b->to;
-							  //            cout << endl;
 							  NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1);
 							  NbNewPoints = NbEdgeCurve-1;
@@ -3907,6 +3454,4 @@
 
 								  lstep = lcurve / NbEdgeCurve;
-								  //   cout <<"lstep " << lstep << " lcurve " 
-								  //    << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " <<NbVerticesOnGeomEdge0<<" " <<NbVerticesOnGeomEdge<<" step =" <<step<<  endl;
 							 } 
 						  // end of curve --
@@ -3979,12 +3524,6 @@
 			 int j= ii;
 			 while (!(*on)[j].Required()) { 
-				 //	cout << i << " " << ii << " j= " << j << " curve = " 
-				 //           <<  on->CurveNumber << "  " << Gh.Number(on) << " on " << j 
-				 //   << " s0 " << Gh.Number( (*on)[0]) << " s1  " << Gh.Number( (*on)[1]) 
-				 //   << " ->  " ;
 				 Adj(on,j); // next geom edge
 				 j=1-j;
-				 //       cout << Gh.Number(on) << "  " << j  << " e[ON] =  " <<  e[Gh.Number(on)] 
-				 //    << " s0 " << Gh.Number( (*on)[0]) << " s1  " << Gh.Number( (*on)[1]) << endl; 
 				 if (e[Gh.Number(on)])  break; // optimisation     
 				 e[Gh.Number(on)] = ei; 
@@ -3993,10 +3532,10 @@
 
 		int kk=0;
-		for ( i=0;i<Gh.nbe ; i++)
-		 if (!e[i]) 
-		  if(kk++<10) {
-			  cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber 
-				 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1  " << Gh.Number( Gh.edges[i][1]) << endl; 
-		  }
+		for ( i=0;i<Gh.nbe ; i++){
+			if (!e[i]){
+				kk++;
+				if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);
+			}
+		}
 		if(kk) throw ErrorException(__FUNCT__,exprintf("See above"));
 
@@ -4217,5 +3756,4 @@
 				 for (Int4 icount=2; icount<nbvb; icount++) {
 					 Vertex *vi  = ordre[icount];
-					 //	  cout << " Add vertex " <<  Number(vi) << endl;
 					 Icoor2 dete[3];
 					 Triangle *tcvi = FindTriangleContening(vi->i,dete);
@@ -4322,16 +3860,8 @@
 				 triangles = savetriangles;
 				 subdomains = savesubdomains;
-				 //	 cout <<  triangles << " <> " << OutSidesTriangles << endl; 
-				 /*	 k=0;
-						 for (i=0;i<nbt;i++)
-						 for (int j=0;j<3;j++)
-						 if (!triangles[i].TriangleAdj(j))
-						 k++;
-						 */
 				 if (k) {
 					 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));
 				 }
 				 FindSubDomain();
-				 // cout << " NbTOld = " << NbTold << " ==  " << nbt - NbOutT << " " << nbt << endl;
 
 				 delete edge4;
@@ -4347,13 +3877,5 @@
 					 if (!edges[i].adj[j])
 					  if(!edges[i][j].on->IsRequiredVertex()) {
-						  cerr << " Erreur adj et sommet requis edges [" << i <<  "][ " << j << "]= "
-							 <<  Number(edges[i][j]) << " : "  << " on = " << Gh.Number(edges[i].on) ;
-						  if (edges[i][j].on->OnGeomVertex())
-							cerr << " vertex " << Gh.Number(edges[i][j].on->gv);
-						  else if (edges[i][j].on->OnGeomEdge())
-							cerr << "Edges " << Gh.Number(edges[i][j].on->ge);
-						  else
-							cerr << " = " << edges[i][j].on ;
-						  cerr << endl;
+						  throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)"));
 					  }
 		  }
@@ -4372,5 +3894,4 @@
 		jp = aa[jp]&3;
 		do {
-			//    cout << *t << " " <<  j  << "\n\t try swap " ;
 			while (t->swap(j,koption))
 			  {
@@ -4382,5 +3903,4 @@
 				t=  tp->at[jp];      // set unchange t qnd j for previous triangles
 				j=  NextEdge[tp->aa[jp]&3];
-				//   cout << "\n\t s  " <<  *t << " " << j << endl;
 			  }
 			// end on this  Triangle 
@@ -4416,6 +3936,5 @@
 		for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 
 		 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
-		if(verbosity>2) 
-		 cout << "  -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl;
+		if(verbosity>2) printf("   SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);
 		for (k=0;k<nbiter;k++)
 		  {
@@ -4429,7 +3948,5 @@
 			  if (tstart[i] != &vide) // not a boundary vertex 
 				NbSwap += vertices[i].Optim(1);
-			if (verbosity>3)
-			 cout << "    Move max = " <<  sqrt(delta) << " iteration = " 
-				<< k << " Nb of Swap = " << NbSwap << endl;
+			if (verbosity>3) printf("      move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);
 		  }
 
@@ -4441,6 +3958,4 @@
 	void Triangles::MakeQuadTree() {  
 		long int verbosity=0;
-		if(verbosity>8)
-		 cout << "      MakeQuadTree" << endl;
 		if (  !quadtree )  quadtree = new QuadTree(this);
 
@@ -4454,10 +3969,4 @@
 		D2xD2 Br(D2xD2(Beq,Heq).t());
 		D2xD2 B1r(Br.inv());
-		/*   D2xD2 BB = Br.t()*Br;
-			  cout << " BB = " << BB << " " << Br*B1r <<  endl; 
-			  MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y);
-			  MatVVP2x2 VMM(MMM);
-			  cout << " " << VMM.lambda1 << " " << VMM.lambda2 <<  endl; 
-			  */
 		double gammamn=1e100,hmin=1e100;
 		double gammamx=0,hmax=0;
@@ -4484,5 +3993,4 @@
 			 MatVVP2x2 VMK(MK);
 			 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
-			 // cout << B_K << " * " << B1r << " == " << BK << " " << B_K*B_K.inv() << endl;
 			 Real8 betaK=0;
 
@@ -4499,5 +4007,4 @@
 				 MetricAnIso M1(BMB.x.x,BMB.x.y,BMB.y.y);
 				 MatVVP2x2 VM1(M1);
-				 //cout << B_K <<" " <<  M << " " <<  he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<<   endl; 
 				 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
 				 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);		
@@ -4505,5 +4012,4 @@
 			 betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
 			 Marea+= betaK;
-			 // cout << betaK << " " << area3 << " " << beta << " " << beta0 << " " << area3*3*3*3 <<endl;
 			 beta=min(beta,betaK);
 			 beta0=max(beta0,betaK);
@@ -4512,9 +4018,10 @@
 		gammamn=sqrt(gammamn);
 		gammamx=sqrt(gammamx);    
-		cout << "  -- adaptmesh Regulary:  Nb triangles " << nt <<  " , h  min " << hmin  << " , h max " << hmax << endl;  
-		cout << "     area =  " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl; 
-		cout << "     infiny-regulaty:  min " << gammamn << "  max " << gammamx << endl;  
-		cout << "     anisomax  "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh) 
-		  << " min  "<<  1./sqrt(beta0/aireKh)  << endl;
+		printf("   Adaptmesh info:\n");
+		printf("      number of triangles = %i\n",nt);
+		printf("      hmin = %g, hmax=%g\n",hmin,hmax);
+		printf("      area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt));
+		printf("      infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx);
+		printf("      anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5));
 	}
 	/*}}}1*/
@@ -4549,29 +4056,16 @@
 				}
 			}  
-		cout << "  -- Histogram of the unit mesh,  nb of edges" << nbedges << endl <<endl;
-
-		cout << "        length of edge in   | % of edge  | Nb of edges " << endl;
-		cout << "        ------------------- | ---------- | ----------- " << endl;
-		for   (i=0;i<=kmax;i++)
-		  { 
-			cout << "    " ;
-			cout.width(10);
-			if (i==0) cout  << " 0 " ;
-			else cout  << exp(lmin+i/delta) ;
-			cout.width(); cout << "," ;
-			cout.width(10);
-			if (i==kmax) cout << " +infty " ;
-			else cout  << exp(lmin+(i+1)/delta) ;
-			cout.width();cout << "   |   " ;
-
-			cout.precision(4);
-			cout.width(6);
-			cout <<  ((long)  ((10000.0 * histo[i])/ nbedges))/100.0 ;
-			cout.width();
-			cout.precision();
-			cout <<  "   |   " << histo[i] <<endl;
+		printf(" --- Histogram of the unit mesh,  nb of edges = %i\n",nbedges);
+		printf("      length of edge in   | %% of edge  | Nb of edges \n"); 
+		printf("      --------------------+-------------+-------------\n"); 
+		for   (i=0;i<=kmax;i++){ 
+			if (i==0) printf("      %10i",0);
+			else      printf("      %10g",exp(lmin+i/delta));
+			if (i==kmax) printf("          +inf   ");
+			else printf("      %10g",exp(lmin+(i+1)/delta));
+			printf("|  %10g |\n",((long)  ((10000.0 * histo[i])/ nbedges))/100.0);
+			printf("  %g\n",histo[i]);
 		  }
-		cout << "        ------------------- | ---------- | ----------- " << endl <<endl;
-
+		printf("      --------------------+-------------+-------------\n"); 
 	}
 	/*}}}1*/
@@ -4605,5 +4099,5 @@
 		if( k==0) return 0;
 		CurrentTh = this;
-		cout << " Nb of Cracked Edges = " << k << endl;
+		printf("      number of Cracked Edges = %i\n",k);
 		NbCrackedEdges =k;
 		CrackedEdges = new  CrackedEdge[k];
@@ -4689,13 +4183,10 @@
 
 		//  set the ref 
-		cout << " set the ref " <<  endl ;
 		NbCrackedVertices =   nbcrakev;
 		// int nbvo = nbv;
 		nbv = vlast - vertices;
 		int nbnewv =  nbv - nbv; // nb of new vrtices 
-		if (nbcrakev && verbosity > 1 )
-		 cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " <<   nbnewv<< endl;
+		if (nbcrakev && verbosity > 1 ) printf("      number of Cracked vertices = %i, number of created vertices = %i\n",nbcrakev,nbnewv);
 		// all the new vertices are on geometry 
-		//  BOFBO--  A VOIR
 		if (nbnewv)
 		  { // 
@@ -4738,5 +4229,4 @@
 			throw ErrorException(__FUNCT__,exprintf("Not enougth vertices: to crack the mesh we need %i vertices",nbv));
 		  }
-		cout << "  NbCrackedVertices " <<  NbCrackedVertices << endl;
 		CurrentTh = CurrentThOld;
 		return  NbCrackedVertices;
@@ -4835,7 +4325,4 @@
 	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);
@@ -4854,5 +4341,4 @@
 			 Vertex V;
 			 VertexOnGeom GV;
-			 // cerr << Number(edges[i]) << " " << ss[j] << endl;
 			 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
 				{
@@ -4870,9 +4356,6 @@
 				 }
 				 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
-				 //cerr << " : " ;
 				 Metric MVp(Vp);
-				 // cerr << " : "  << MVp  << endl;
 				 edges[i][j].m.IntersectWith(MVp);
-				 //cerr << " . " << endl;
 				}
 
@@ -4886,6 +4369,5 @@
 
 		double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
-		if (verbosity > 1) 
-		 cout << "  -- BoundAnisotropy by  " << anisomax << endl; 
+		if (verbosity > 1)  printf("   BoundAnisotropy by %g\n",anisomax);
 		Real8 h1=1.e30,h2=1e-30,rx=0;
 		Real8 coef = 1./(anisomax*anisomax);
@@ -4913,10 +4395,7 @@
 		  }
 
-		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;
+		if (verbosity>2){
+			printf("      input:  Hmin = %g, Hmax = %g, factor of anisotropy max  = %g\n",pow(h2,-0.5),pow(h1,-0.5),pow(rx,0.5));
+			printf("      output: Hmin = %g, Hmax = %g, factor of anisotropy max  = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5));
 		  }
 	}
@@ -4960,17 +4439,13 @@
 		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;
+		if(verbosity>1) {
+			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv);
+			printf("      coef = %g\n",coef); 
+			printf("      hmin = %g hmax = %g\n",hmin,hmax); 
+			printf("      anisomax = %g nb Jacobi = %i, power = %i\n",anisomax,NbJacobi,power); 
+			if (RelativeMetric) printf("      RelativeErr with CutOff= %g\n",CutOff);
+			else printf("      Absolute error\n");
 		  }
-		double *ss=(double*)s;//, *ssiii = ss;
+		double *ss=(double*)s;
 
 		double sA,sB,sC;
@@ -5070,14 +4545,8 @@
 			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;
+			if(verbosity>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield);
+
+			if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) {
+				if (verbosity>2) printf("      Solution %i is constant, skipping...\n");
 				continue;
 			  }
@@ -5144,10 +4613,4 @@
 								 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] ;
 
@@ -5184,9 +4647,7 @@
 						 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)
 							{
@@ -5245,5 +4706,4 @@
 					 Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
 					 MatVVP2x2 Vp(M);
-					 //cout <<iv <<  "  M  = " <<  M <<  " aniso= " << Vp.Aniso() ;
 					 Vp.Abs();
 					 M = Vp;
@@ -5251,5 +4711,4 @@
 					 dxdy[iv] = M.a21;
 					 dydy[iv] = M.a22;
-					 //  cout << " (abs)  iv M  = " <<  M <<  " aniso= " << Vp.Aniso() <<endl;
 					}
 				 else kk++;
@@ -5311,11 +4770,5 @@
 				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)
@@ -5329,7 +4782,4 @@
 					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);
@@ -5358,11 +4808,8 @@
 					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;
+				if (verbosity>2) { 
+					printf("      Field %i of solution %i\n",nufield,nusol);
+					printf("         before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5));
+					printf("         after  bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5));
 				  }
 			  } //  end of for all field
@@ -5386,6 +4833,5 @@
 
 		const  Real8 maxsubdiv2 = maxsubdiv*maxsubdiv;
-		if(verbosity>1)
-		 cout << "  -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv <<   endl  ;
+		if(verbosity>1) printf("   Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);
 		// for all the edges 
 		// if the len of the edge is to long 
@@ -5434,8 +4880,7 @@
 			  }
 		  }
-		if(verbosity>3)
-		 cout << "    Nb of metric change = " << nbchange 
-			<< " Max  of the subdivision of a edges before change  = " << sqrt(lmax) << endl;
-
+		if(verbosity>3){
+			printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
+		}
 	}
 	/*}}}1*/
@@ -5445,6 +4890,5 @@
 
 		if(raisonmax<1.1) return;
-		if(verbosity > 1)
-		 cout << "  -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl;
+		if(verbosity > 1) printf("   Triangles::SmoothMetric raisonmax = %g\n",raisonmax);
 		ReMakeTriangleContainingTheVertex();
 		Int4 i,j,kch,kk,ip;
@@ -5464,5 +4908,4 @@
 			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;
 				if (!t){
@@ -5491,5 +4934,4 @@
 							  {
 								Real8 dh=(hj-hi)/ll;
-								//cout << " dh = " << dh << endl;
 								if (dh>logseuil) {
 									vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
@@ -5502,9 +4944,5 @@
 						  {
 							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;
@@ -5517,9 +4955,6 @@
 			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;
+		}
+		if(verbosity>2) printf("      number of iterations = %i\n",kch); 
 		delete [] first_np_or_next_t0;
 		delete [] first_np_or_next_t1;
@@ -5599,9 +5034,7 @@
 			strcat(strcat(identity,";"),countbuf);
 			strcat(identity,buf);
-			// cout << "New MAILLAGE "<< identity << endl;
 		  } 
 
 		quadtree=0;
-		//  edgescomponante=0;
 		edges=0;
 		VerticesOnGeomVertex=0;
@@ -5609,14 +5042,6 @@
 		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*/
@@ -5724,5 +5149,4 @@
 			I2 I=vI, J=vJ, IJ= J-I;
 			IJ_IA = (IJ ,(A-I));
-			//   cout << A << vI.i << vJ.i << edge << " " <<  IJ_IA << " dir " << dir <<endl;
 			if (IJ_IA<0) {
 				if (dir>0) {a=1;b=0;return edge;}// change of signe => I
@@ -5740,5 +5164,4 @@
 							a= IJ_AJ/IJ2;
 							b= IJ_IA/IJ2;
-							//    cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl;
 							return edge;
 		  } 
@@ -5786,5 +5209,4 @@
 				if ( (d2=(double) ACAC)  <  dd2) 
 				  {
-					//  cout << " A "  << d2  << " " <<  dd2;
 					er = ta;
 					l0 = ACAC;
@@ -5798,5 +5220,4 @@
 				if ( (d2=(double) BCBC)  <  dd2) 
 				  {
-					// cout << " B "  << d2  << " " <<  dd2;
 					dd2 = d2;
 					er = Adj(ta); // other direction
@@ -5822,6 +5243,4 @@
 					s = 0;
 					cas = -1;
-					//	 cout << " ABAC " <<  ABAC << " ABAC " << ABAC
-					//	      << " AB2 " << AB2 << endl;
 					b = ((double) ABAC/(double) AB2);
 					a = 1 - b;
@@ -5911,5 +5330,4 @@
 			if ( o < oo) 
 			 o=oo,pi=BigPrimeNumber[i];}
-			//  cout << " AGoodNumberPrimeWith " << n << " " <<pi << " "<< o << endl;
 			return pi; 
 	}
@@ -6050,5 +5468,4 @@
 			}
 			det2 = det(*v2,a,b);
-			//   cout << " No Change try the next" << endl;
 		}
 
@@ -6068,5 +5485,4 @@
 				}
 				Icoor2 detss = 0,l=0,ks;
-				// cout << "Real ForcingEdge " << *va << *vb << detss << endl;
 				while ((ks=SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
 				 if(l++ > 10000000) {