Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgVertex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/BamgVertex.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/BamgVertex.cpp	(revision 12509)
@@ -13,8 +13,8 @@
 	void BamgVertex::Echo(void){
 
-		printf("Vertex:\n");
-		printf("  integer   coordinates i.x: %i, i.y: %i\n",i.x,i.y);
-		printf("  Euclidean coordinates r.x: %g, r.y: %g\n",r.x,r.y);
-		printf("  ReferenceNumber = %i\n",ReferenceNumber);
+		_printLine_("Vertex:");
+		_printLine_("  integer   coordinates i.x: " << i.x << ", i.y: " << i.y);
+		_printLine_("  Euclidean coordinates r.x: " << r.x << ", r.y: " << r.y);
+		_printLine_("  ReferenceNumber = " << ReferenceNumber);
 		m.Echo();
 
Index: /issm/trunk-jpl/src/c/objects/Bamg/Edge.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/Edge.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/Edge.cpp	(revision 12509)
@@ -26,9 +26,9 @@
 	/*FUNCTION Edge::Echo {{{*/
 	void Edge::Echo(void){ 
-		printf("Edge:\n");
-		printf("   pointers towards two vertices: %p %p\n",v[0],v[1]);
-		printf("   ReferenceNumber = %i\n",ReferenceNumber);
-		printf("   GeomEdgeHook = %p\n",GeomEdgeHook);
-		printf("   two adjacent edges on the same curve: %p %p\n",adj[0],adj[1]);
+		_printLine_("Edge:");
+		_printLine_("   pointers towards two vertices: " << v[0] << " " << v[1]);
+		_printLine_("   ReferenceNumber = " << ReferenceNumber);
+		_printLine_("   GeomEdgeHook = " << GeomEdgeHook);
+		_printLine_("   two adjacent edges on the same curve: " << adj[0] << " " << adj[1]);
 	}
 	/*}}}*/
Index: /issm/trunk-jpl/src/c/objects/Bamg/EigenMetric.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/EigenMetric.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/EigenMetric.cpp	(revision 12509)
@@ -114,9 +114,9 @@
 	void EigenMetric::Echo(void){
 
-		printf("EigenMetric:\n");
-		printf("   lambda1: %g\n",lambda1);
-		printf("   lambda2: %g\n",lambda2);
-		printf("   v.x: %g\n",v.x);
-		printf("   v.y: %g\n",v.y);
+		_printLine_("EigenMetric:");
+		_printLine_("   lambda1: " << lambda1);
+		_printLine_("   lambda2: " << lambda2);
+		_printLine_("   v.x: " << v.x);
+		_printLine_("   v.y: " << v.y);
 
 		return;
Index: /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp	(revision 12509)
@@ -48,5 +48,5 @@
 	Geometry::~Geometry() {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
-		if(NbRef>0){   printf("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
+		if(NbRef>0){   _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
 		if(vertices)   delete [] vertices;  vertices=0;
 		if(edges)      delete [] edges;     edges=0;
@@ -81,5 +81,5 @@
 		//Vertices
 		if (bamggeom->Vertices){
-			if(verbose>5) printf("      processing Vertices\n");
+			if(verbose>5) _printLine_("      processing Vertices");
 			if (bamggeom->VerticesSize[1]!=3) _error2_("Vertices should have 3 columns");
 			vertices = new GeomVertex[nbv];
@@ -126,5 +126,5 @@
 			double* verticeslength=NULL;
 
-			if(verbose>5) printf("      processing Edges\n");
+			if(verbose>5) _printLine_("      processing Edges");
 			if (bamggeom->EdgesSize[1]!=3) _error2_("Edges should have 3 columns");
 			edges = new GeomEdge[nbe];
@@ -181,5 +181,5 @@
 		//hVertices
 		if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){
-			if(verbose>5) printf("      processing hVertices\n");
+			if(verbose>5) _printLine_("      processing hVertices");
 			for (i=0;i< nbv;i++){
 				if (!isnan(bamgopts->hVertices[i])){
@@ -191,5 +191,5 @@
 		//MetricVertices
 		if(bamgopts->metric && bamgopts->metric[0]==nbv){
-			if(verbose>5) printf("      processing MetricVertices\n");
+			if(verbose>5) _printLine_("      processing MetricVertices");
 			for (i=0;i< nbv;i++) {
 				vertices[i].m = Metric((double)bamgopts->metric[i*3+0],(double)bamgopts->metric[i*3+1],(double)bamgopts->metric[i*3+2]);
@@ -199,5 +199,5 @@
 		//MaxCornerAngle
 		if (bamgopts->MaxCornerAngle){
-			if(verbose>5) printf("      processing MaxCornerAngle\n");
+			if(verbose>5) _printLine_("      processing MaxCornerAngle");
 			MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;
 		}
@@ -205,5 +205,5 @@
 		//TangentAtEdges
 		if (bamggeom->TangentAtEdges){
-			if(verbose>5) printf("      processing TangentAtEdges");
+			if(verbose>5) _printString_("      processing TangentAtEdges");
 			if (bamggeom->TangentAtEdgesSize[1]!=4) _error2_("TangentAtEdges should have 4 columns");
 			int n,i,j,k;
@@ -224,5 +224,5 @@
 		//Corners
 		if(bamggeom->Corners){
-			if(verbose>5) printf("      processing Corners");
+			if(verbose>5) _printString_("      processing Corners");
 			if (bamggeom->CornersSize[1]!=1) _error2_("Corners should have 1 column");
 			n=bamggeom->CornersSize[0];
@@ -238,5 +238,5 @@
 		//RequiredVertices
 		if(bamggeom->RequiredVertices){
-			if(verbose>5) printf("      processing RequiredVertices\n");
+			if(verbose>5) _printLine_("      processing RequiredVertices");
 			if (bamggeom->RequiredVerticesSize[1]!=1) _error2_("RequiredVertices should have 1 column");
 			n=bamggeom->RequiredVerticesSize[0];
@@ -250,5 +250,5 @@
 		//RequiredEdges
 		if(bamggeom->RequiredEdges){
-			if(verbose>5) printf("      processing RequiredEdges\n");
+			if(verbose>5) _printLine_("      processing RequiredEdges");
 			if (bamggeom->RequiredEdgesSize[1]!=1) _error2_("RequiredEdges should have 1 column");
 			n=bamggeom->RequiredEdgesSize[0];
@@ -262,5 +262,5 @@
 		//SubDomain
 		if(bamggeom->SubDomains){
-			if(verbose>5) printf("      processing SubDomains\n");
+			if(verbose>5) _printLine_("      processing SubDomains");
 			if (bamggeom->SubDomainsSize[1]!=4) _error2_("SubDomains should have 4 columns");
 			nbsubdomains=bamggeom->SubDomainsSize[0];
@@ -294,5 +294,5 @@
 
 		/*Vertices*/
-		if(verbose>5) printf("      writing Vertices\n");
+		if(verbose>5) _printLine_("      writing Vertices");
 		bamggeom->VerticesSize[0]=nbv;
 		bamggeom->VerticesSize[1]=3;
@@ -310,5 +310,5 @@
 
 		/*Edges*/
-		if(verbose>5) printf("      writing Edges\n");
+		if(verbose>5) _printLine_("      writing Edges");
 		bamggeom->EdgesSize[0]=nbe;
 		bamggeom->EdgesSize[1]=3;
@@ -328,5 +328,5 @@
 
 		/*RequiredEdges*/
-		if(verbose>5) printf("      writing %i RequiredEdges\n",nbreq);
+		if(verbose>5) _printLine_("      writing " << nbreq << " RequiredEdges");
 		bamggeom->RequiredEdgesSize[0]=nbreq;
 		bamggeom->RequiredEdgesSize[1]=1;
@@ -345,5 +345,5 @@
 
 		/*RequiredVertices*/
-		if(verbose>5) printf("      writing %i RequiredVertices\n",nbreqv);
+		if(verbose>5) _printLine_("      writing " << nbreqv << " RequiredVertices");
 		bamggeom->RequiredVerticesSize[0]=nbreqv;
 		bamggeom->RequiredVerticesSize[1]=1;
@@ -360,5 +360,5 @@
 
 		/*SubDomains*/
-		if(verbose>5) printf("      writing SubDomains\n");
+		if(verbose>5) _printLine_("      writing SubDomains");
 		bamggeom->SubDomainsSize[0]=nbsubdomains;
 		bamggeom->SubDomainsSize[1]=4;
@@ -374,5 +374,5 @@
 
 		/*TangentAtEdges*/
-		if(verbose>5) printf("      writing TangentAtEdges\n");
+		if(verbose>5) _printLine_("      writing TangentAtEdges");
 		bamggeom->TangentAtEdgesSize[0]=nbtan;
 		bamggeom->TangentAtEdgesSize[1]=4;
@@ -403,18 +403,18 @@
 	void Geometry::Echo(void){
 
-		printf("Geometry:\n");
-		printf("   nbv  (number of vertices) : %i\n",nbv);
-		printf("   nbe  (number of edges)    : %i\n",nbe);
-		printf("   nbsubdomains: %i\n",nbsubdomains);
-		printf("   nbcurves: %i\n",nbcurves);
-		printf("   vertices: %p\n",vertices);
-		printf("   edges: %p\n",edges);
-		printf("   quadtree: %p\n",quadtree);
-		printf("   subdomains: %p\n",subdomains);
-		printf("   curves: %p\n",curves);
-		printf("   pmin (x,y): (%g %g)\n",pmin.x,pmin.y);
-		printf("   pmax (x,y): (%g %g)\n",pmax.x,pmax.y);
-		printf("   coefIcoor: %g\n",coefIcoor);
-		printf("   MaxCornerAngle: %g\n",MaxCornerAngle);
+		_printLine_("Geometry:");
+		_printLine_("   nbv  (number of vertices) : " << nbv);
+		_printLine_("   nbe  (number of edges)    : " << nbe);
+		_printLine_("   nbsubdomains: " << nbsubdomains);
+		_printLine_("   nbcurves: " << nbcurves);
+		_printLine_("   vertices: " << vertices);
+		_printLine_("   edges: " << edges);
+		_printLine_("   quadtree: " << quadtree);
+		_printLine_("   subdomains: " << subdomains);
+		_printLine_("   curves: " << curves);
+		_printLine_("   pmin (x,y): (" << pmin.x << " " << pmin.y << ")");
+		_printLine_("   pmax (x,y): (" << pmax.x << " " << pmax.y << ")");
+		_printLine_("   coefIcoor: " << coefIcoor);
+		_printLine_("   MaxCornerAngle: " << MaxCornerAngle);
 
 		return;
@@ -527,6 +527,6 @@
 			/*if there is a vertex found that is to close to vertices[i] -> error*/
 			if( v && Norme1(v->r - vertices[i].r) < eps ){
-				printf("reference numbers: %i %i\n",v->ReferenceNumber,vertices[i].ReferenceNumber);
-				printf("Id: %i\n",i+1);
+				_printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
+				_printLine_("Id: " << i+1);
 				delete [] next_p;
 				delete [] head_v;
@@ -870,9 +870,9 @@
 			if (bge<=0) {
 				if(NbTry) {
-					printf("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve\n");
-					printf("That bug might come from:\n");
-					printf(" 1)  a mesh edge  containing more than %i geometrical edges\n",mxe/2);
-					printf(" 2)  code bug : be sure that we call   Mesh::SetVertexFieldOn() before\n");
-					printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n");
+					_printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
+					_printLine_("That bug might come from:");
+					_printLine_(" 1)  a mesh edge  containing more than " << mxe/2 << " geometrical edges");
+					_printLine_(" 2)  code bug : be sure that we call   Mesh::SetVertexFieldOn() before");
+					_printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
 					_error2_("see above");
 				}
@@ -887,12 +887,12 @@
 		while (eg1 != (GeomEdge*) vg1  &&  (*eg1)(direction1) != (GeomVertex*) vg1) { 
 			if(tge>=mxe ) { 
-				printf("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n");
+				_printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
 				NbTry++;
 				if (NbTry<2) goto retry;
-				printf("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve\n");
-				printf("That bug might come from:\n");
-				printf(" 1)  a mesh edge  contening more than %i geometrical edges\n",mxe/2);
-				printf(" 2)  code bug : be sure that we call   Mesh::SetVertexFieldOn() before\n");
-				printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n");
+				_printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
+				_printLine_("That bug might come from:");
+				_printLine_(" 1)  a mesh edge  contening more than " << mxe/2 << " geometrical edges");
+				_printLine_(" 2)  code bug : be sure that we call   Mesh::SetVertexFieldOn() before");
+				_printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
 				_error2_("see above");
 			}
Index: /issm/trunk-jpl/src/c/objects/Bamg/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 12509)
@@ -167,5 +167,5 @@
 			//x.Echo();
 			//y.Echo();
-			//printf("cx = %g, cy=%g\n",cx,cy);
+			//_printLine_("cx = " << cx << ", cy=" << cy);
 
 			si += sint;
@@ -186,5 +186,5 @@
 			MaxNbSeg *= 2;
 			if (verbosity>3){
-				printf("   reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
+				_printLine_("   reshape lSegsI from " << mneo << " to " << MaxNbSeg);
 			}
 			_assert_(lSegsI && NbSeg<MaxNbSeg);
@@ -211,5 +211,5 @@
 		for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];       
 		long int verbosity=0;
-		if(verbosity>3) printf("   ListofIntersectionTriangles  ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg);
+		if(verbosity>3) _printLine_("   ListofIntersectionTriangles  ReShape Maxsize " << MaxSize << " -> " << MaxNbSeg);
 		MaxSize = newsize; 
 		delete [] lIntTria;// remove old
Index: /issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp	(revision 12509)
@@ -29,5 +29,5 @@
 		if(bamggeom->Edges==NULL) {
 			/*Recreate geometry if needed*/
-			printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
+			_printLine_("WARNING: mesh present but no geometry found. Reconstructing...");
 			BuildGeometryFromMesh(bamgopts);
 			Gh.PostRead();
@@ -101,7 +101,7 @@
 			  if (kk[i]>=0) kk[i]=k++;
 			}
-		  printf("   number of vertices %i, remove = %i\n",k,Tho.nbv - k);
-		  printf("   number of triangles %i, remove = %i\n",kt,nbInT-kt);
-		  printf("   number of New boundary edge %i\n",nbNewBedge);
+		  _printLine_("   number of vertices " << k << ", remove = " << Tho.nbv - k);
+		  _printLine_("   number of triangles " << kt << ", remove = " << nbInT-kt);
+		  _printLine_("   number of New boundary edge " << nbNewBedge);
 		  long imaxnbv =k;
 		  Init(imaxnbv);
@@ -280,5 +280,5 @@
 
 		//Vertices
-		if (verbose) printf("Reading vertices (%i)\n",nbv);
+		if (verbose) _printLine_("Reading vertices (" << nbv << ")");
 		vertices=xNew<BamgVertex>(nbv);
 		orderedvertices=xNew<BamgVertex*>(nbv);
@@ -294,5 +294,5 @@
 
 		//Triangles
-		if (verbose) printf("Reading triangles (%i)\n",nbt);
+		if (verbose) _printLine_("Reading triangles (" << nbt << ")");
 		triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since 
 		nodeflags=xNew<bool>(nbv);
@@ -310,7 +310,7 @@
 
 		/*Recreate geometry: */
-		if (verbose) printf("Building Geometry\n");
+		if (verbose) _printLine_("Building Geometry");
 		BuildGeometryFromMesh();
-		if (verbose) printf("Completing geometry\n");
+		if (verbose) _printLine_("Completing geometry");
 		Gh.PostRead();
 
@@ -319,5 +319,5 @@
 		for(i=0;i<nbv;i++){
 			if(!nodeflags[i]){
-				printf("Vertex %i does not belong to any element\n",i+1);
+				_printLine_("Vertex " << i+1 << " does not belong to any element");
 				isorphan=true;
 			}
@@ -346,5 +346,5 @@
 		//Vertices
 		if(bamgmesh->Vertices){
-			if(verbose>5) printf("      processing Vertices\n");
+			if(verbose>5) _printLine_("      processing Vertices");
 
 			vertices=xNew<BamgVertex>(nbv);
@@ -367,5 +367,5 @@
 		//Triangles
 		if(bamgmesh->Triangles){
-			if(verbose>5) printf("      processing Triangles\n");
+			if(verbose>5) _printLine_("      processing Triangles");
 			triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since 
 			//other triangles will be added for each edge
@@ -385,5 +385,5 @@
 		//Quadrilaterals
 		if(bamgmesh->Quadrilaterals){
-			if(verbose>5) printf("      processing Quadrilaterals\n");
+			if(verbose>5) _printLine_("      processing Quadrilaterals");
 			long i1,i2,i3,i4,iref;
 			triangles =new Triangle[nbt];
@@ -407,5 +407,5 @@
 		//VerticesOnGeomEdge
 		if(bamgmesh->VerticesOnGeomEdge){
-			if(verbose>5) printf("      processing VerticesOnGeomEdge\n");
+			if(verbose>5) _printLine_("      processing VerticesOnGeomEdge");
 			NbVerticesOnGeomEdge=bamgmesh->VerticesOnGeomEdgeSize[0];
 			VerticesOnGeomEdge= new  VertexOnGeom[NbVerticesOnGeomEdge] ;
@@ -422,5 +422,5 @@
 		//VerticesOnGeomVertex
 		if(bamgmesh->VerticesOnGeomVertexSize[0]){
-			if(verbose>5) printf("      processing VerticesOnGeomVertex\n");
+			if(verbose>5) _printLine_("      processing VerticesOnGeomVertex");
 			NbVerticesOnGeomVertex=bamgmesh->VerticesOnGeomVertexSize[0];
 			VerticesOnGeomVertex  = new  VertexOnGeom[NbVerticesOnGeomVertex] ;
@@ -438,5 +438,5 @@
 			double* len=NULL;
 
-			if(verbose>5) printf("      processing Edges\n");
+			if(verbose>5) _printLine_("      processing Edges");
 			nbe=bamgmesh->EdgesSize[0];
 			edges= new Edge[nbe];
@@ -498,5 +498,5 @@
 		//EdgeOnGeomEdge
 		if(bamgmesh->EdgesOnGeomEdge){
-			if(verbose>5) printf("      processing EdgesOnGeomEdge\n");
+			if(verbose>5) _printLine_("      processing EdgesOnGeomEdge");
 			int i1,i2,i,j;
 			i2=bamgmesh->EdgesOnGeomEdgeSize[0];
@@ -515,5 +515,5 @@
 		if(bamgmesh->SubDomains){
 			long i3,head,direction;
-			if(verbose>5) printf("      processing SubDomains\n");
+			if(verbose>5) _printLine_("      processing SubDomains");
 			nbsubdomains=bamgmesh->SubDomainsSize[0];
 			subdomains = new SubDomain [ nbsubdomains ];
@@ -584,5 +584,5 @@
 
 		/*Vertices*/
-		if(verbose>5) printf("      writing Vertices\n");
+		if(verbose>5) _printLine_("      writing Vertices");
 		bamgmesh->VerticesSize[0]=nbv;
 		bamgmesh->VerticesSize[1]=3;
@@ -597,5 +597,5 @@
 
 		/*Edges*/
-		if(verbose>5) printf("      writing Edges\n");
+		if(verbose>5) _printLine_("      writing Edges");
 		bamgmesh->EdgesSize[0]=nbe;
 		bamgmesh->EdgesSize[1]=3;
@@ -614,5 +614,5 @@
 
 		/*Element edges*/
-		if(verbose>5) printf("      writing element edges\n");
+		if(verbose>5) _printLine_("      writing element edges");
 		SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv);
 		double* elemedge=NULL;
@@ -670,5 +670,5 @@
 
 		/*IssmSegments*/
-		if(verbose>5) printf("      writing IssmSegments\n");
+		if(verbose>5) _printLine_("      writing IssmSegments");
 		bamgmesh->IssmSegmentsSize[0]=NumIssmSegments;
 		bamgmesh->IssmSegmentsSize[1]=4;
@@ -713,5 +713,5 @@
 
 		/*Triangles*/
-		if(verbose>5) printf("      writing Triangles\n");
+		if(verbose>5) _printLine_("      writing Triangles");
 		k=nbInT-nbq*2;
 		num=0;
@@ -734,5 +734,5 @@
 
 		/*Quadrilaterals*/
-		if(verbose>5) printf("      writing Quadrilaterals\n");
+		if(verbose>5) _printLine_("      writing Quadrilaterals");
 		bamgmesh->QuadrilateralsSize[0]=nbq;
 		bamgmesh->QuadrilateralsSize[1]=5;
@@ -755,5 +755,5 @@
 
 		/*SubDomains*/
-		if(verbose>5) printf("      writing SubDomains\n");
+		if(verbose>5) _printLine_("      writing SubDomains");
 		bamgmesh->SubDomainsSize[0]=nbsubdomains;
 		bamgmesh->SubDomainsSize[1]=4;
@@ -769,5 +769,5 @@
 
 		/*SubDomainsFromGeom*/
-		if(verbose>5) printf("      writing SubDomainsFromGeom\n");
+		if(verbose>5) _printLine_("      writing SubDomainsFromGeom");
 		bamgmesh->SubDomainsFromGeomSize[0]=Gh.nbsubdomains;
 		bamgmesh->SubDomainsFromGeomSize[1]=4;
@@ -783,5 +783,5 @@
 
 		/*VerticesOnGeomVertex*/
-		if(verbose>5) printf("      writing VerticesOnGeomVertex\n");
+		if(verbose>5) _printLine_("      writing VerticesOnGeomVertex");
 		bamgmesh->VerticesOnGeomVertexSize[0]=NbVerticesOnGeomVertex;
 		bamgmesh->VerticesOnGeomVertexSize[1]=2;
@@ -797,5 +797,5 @@
 
 		/*VertexOnGeomEdge*/
-		if(verbose>5) printf("      writing VerticesOnGeomEdge\n");
+		if(verbose>5) _printLine_("      writing VerticesOnGeomEdge");
 		bamgmesh->VerticesOnGeomEdgeSize[0]=NbVerticesOnGeomEdge;
 		bamgmesh->VerticesOnGeomEdgeSize[1]=3;
@@ -814,5 +814,5 @@
 
 		/*EdgesOnGeomEdge*/
-		if(verbose>5) printf("      writing EdgesOnGeomEdge\n");
+		if(verbose>5) _printLine_("      writing EdgesOnGeomEdge");
 		k=0;
 		for (i=0;i<nbe;i++){
@@ -834,5 +834,5 @@
 
 		/*Element Connectivity*/
-		if(verbose>5) printf("      writing Element connectivity\n");
+		if(verbose>5) _printLine_("      writing Element connectivity");
 		bamgmesh->ElementConnectivitySize[0]=nbt-nbtout;
 		bamgmesh->ElementConnectivitySize[1]=3;
@@ -854,5 +854,5 @@
 
 		/*ElementNodal Connectivity*/
-		if(verbose>5) printf("      writing Nodal element connectivity\n");
+		if(verbose>5) _printLine_("      writing Nodal element connectivity");
 		bamgmesh->NodalElementConnectivitySize[0]=nbv;
 		bamgmesh->NodalElementConnectivitySize[1]=connectivitymax_1;
@@ -869,5 +869,5 @@
 
 		/*Nodal Connectivity*/
-		if(verbose>5) printf("      writing Nodal connectivity\n");
+		if(verbose>5) _printLine_("      writing Nodal connectivity");
 		//chaining algorithm (again...)
 		int* head_2=NULL;
@@ -921,5 +921,5 @@
 
 		/*Cracked vertices*/
-		if(verbose>5) printf("      writing Cracked vertices\n");
+		if(verbose>5) _printLine_("      writing Cracked vertices");
 		bamgmesh->CrackedVerticesSize[0]=NbCrackedVertices;
 		bamgmesh->CrackedVerticesSize[1]=2;
@@ -933,5 +933,5 @@
 
 		/*Cracked vertices*/
-		if(verbose>5) printf("      writing Cracked vertices\n");
+		if(verbose>5) _printLine_("      writing Cracked vertices");
 		bamgmesh->CrackedEdgesSize[0]=NbCrackedEdges;
 		bamgmesh->CrackedEdgesSize[1]=4;
@@ -961,5 +961,5 @@
 		int  i,j;
 
-		if(bamgopts->verbose>3) printf("      processing metric\n");
+		if(bamgopts->verbose>3) _printLine_("      processing metric");
 		double hmin = Max(bamgopts->hmin,MinimalHmin());
 		double hmax = Min(bamgopts->hmax,MaximalHmax());
@@ -1264,5 +1264,5 @@
 
 		//display info
-		if (verbose > 1)  printf("   BoundAnisotropy by %g\n",anisomax);
+		if (verbose > 1)  _printLine_("   BoundAnisotropy by " << anisomax);
 
 		double h1=1.e30,h2=1e-30;
@@ -1291,6 +1291,6 @@
 		//display info
 		if (verbose>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));
+			_printLine_("      input:  Hmin = " << pow(h2,-0.5)  << ", Hmax = " << pow(h1,-0.5) << ", factor of anisotropy max  = " << pow(rx,0.5));
+			_printLine_("      output: Hmin = " << pow(hn2,-0.5) << ", Hmax = " << pow(hn1,-0.5)<< ", factor of anisotropy max  = " <<pow(rnx,0.5));
 		}
 	}
@@ -1314,5 +1314,5 @@
 
 		//display info
-		if (verbose>1) printf("   construction of the geometry from the 2d mesh\n");
+		if (verbose>1) _printLine_("   construction of the geometry from the 2d mesh");
 
 		//check that the mesh is not empty
@@ -1374,8 +1374,8 @@
 				//else (see 3 lines above), the edge has been called more than twice: return error
 				else {
-					printf("The edge (%i,%i) belongs to more than 2 triangles (%i)\n",GetId(triangles[i][VerticesOfTriangularEdge[j][0]]),GetId(triangles[i][VerticesOfTriangularEdge[j][1]]),k);
-					printf("Edge %i of triangle %i\n",j,i);
-					printf("Edge %i of triangle %i\n",(-st[k]+2)%3,(-st[k]+2)/3);
-					printf("Edge %i of triangle %i\n",triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3)),GetId(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))));
+					_printLine_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << "," << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles (" << k << ")");
+					_printLine_("Edge " << j << " of triangle " << i);
+					_printLine_("Edge " << (-st[k]+2)%3 << " of triangle " << (-st[k]+2)/3);
+					_printLine_("Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3)) << " of triangle " << GetId(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))));
 					_error2_("An edge belongs to more than 2 triangles");
 				}	
@@ -1389,10 +1389,10 @@
 		//display info
 		if(verbose>5) {
-			printf("         info on Mesh:\n");
-			printf("            - number of vertices    = %i \n",nbv); 
-			printf("            - number of triangles   = %i \n",nbt); 
-			printf("            - number of given edges = %i \n",nbe); 
-			printf("            - number of all edges   = %i \n",nbedges); 
-			printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-nbedges+nbv); 
+			_printLine_("         info on Mesh:");
+			_printLine_("            - number of vertices    = " << nbv); 
+			_printLine_("            - number of triangles   = " << nbt); 
+			_printLine_("            - number of given edges = " << nbe); 
+			_printLine_("            - number of all edges   = " << nbedges); 
+			_printLine_("            - Euler number 1 - nb of holes = " << nbt-nbedges+nbv); 
 		}
 
@@ -1425,5 +1425,5 @@
 
 			//display info
-			if(verbose>4) printf("   Construction of the edges %i\n",nbe);
+			if(verbose>4) _printLine_("   Construction of the edges " << nbe);
 
 			for (i=0;i<nbedges;i++){ 
@@ -1451,5 +1451,5 @@
 						edges[add].ReferenceNumber=edgessave[i].ReferenceNumber; 		      
 						edges[add].GeomEdgeHook=edgessave[i].GeomEdgeHook; //  HACK to get required edges
-						printf("oh no...\n");
+						_printLine_("oh no...");
 					}
 					else
@@ -1566,5 +1566,5 @@
 			}
 		}
-		if (verbose> 3) printf("      The Number of sub domain = %i\n",nbsubdomains); 
+		if (verbose> 3) _printLine_("      The Number of sub domain = " << nbsubdomains); 
 
 		//build subdomains
@@ -1618,5 +1618,5 @@
 		Gh.nbsubdomains = nbsubdomains;
 		Gh.subdomains = new GeomSubDomain[nbsubdomains];
-		if (verbose>3) printf("   number of vertices = %i\n   number of edges = %i\n",Gh.nbv,Gh.nbe);
+		if (verbose>3) _printLine_("   number of vertices = " << Gh.nbv << "\n   number of edges = " << Gh.nbe);
 		NbVerticesOnGeomVertex = Gh.nbv;
 		VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
@@ -1790,5 +1790,5 @@
 		//display infos
 		if(verbose>1) {
-			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv);
+			_printLine_("   Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")");
 		}
 
@@ -1862,9 +1862,9 @@
 
 			//display info
-			if(verbose>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g\n",nusol,smin,smax,sdelta);
+			if(verbose>2) _printLine_("      Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta);
 
 			//skip constant field
 			if (sdelta < 1.0e-10*Max(absmax,1e-20)){
-				printf("      Solution %i is constant, skipping...\n",nusol);
+				_printLine_("      Solution " << nusol << " is constant, skipping...");
 				continue;
 			}
@@ -1992,5 +1992,5 @@
 		//display infos
 		if(verbose>1) {
-			printf("   Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",nbsol,nbt,nbv);
+			_printLine_("   Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")");
 		}
 
@@ -2075,9 +2075,9 @@
 
 			//display info
-			if(verbose>2) printf("      Solution %i, Min = %g, Max = %g, Delta = %g, number of fields = %i\n",nusol,smin,smax,sdelta,nbsol);
+			if(verbose>2) _printLine_("      Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << ", number of fields = " << nbsol);
 
 			//skip constant field
 			if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
-				if (verbose>2) printf("      Solution %i is constant, skipping...\n",nusol);
+				if (verbose>2) _printLine_("      Solution " << nusol << " is constant, skipping...");
 				continue;
 			}
@@ -2264,5 +2264,5 @@
 		//Return if no edge is cracked
 		if(k==0) return;
-		if (verbose>4) printf("      number of Cracked Edges = %i\n",k);
+		if (verbose>4) _printLine_("      number of Cracked Edges = " << k);
 
 		//Initialize Cracked edge
@@ -2305,5 +2305,5 @@
 
 		//Add new vertices
-		if (verbose>4) printf("      number of Cracked Vertices = %i\n",NbCrackedVertices);
+		if (verbose>4) _printLine_("      number of Cracked Vertices = " << NbCrackedVertices);
 		if (NbCrackedVertices){
 			CrackedVertices=xNew<long>(2*NbCrackedVertices);
@@ -2374,5 +2374,5 @@
 					}
 				}
-				//printf("%i -> %i %i %i, edge [%i->%i %i->%i]\n",element_renu[GetId(ta.t)],GetId((*ta.t)[0])+1,GetId((*ta.t)[1])+1,GetId((*ta.t)[2])+1,i1,j1,i2,j2);
+				//_printLine_("" << element_renu[GetId(ta.t)] << " -> " << GetId((*ta.t)[0])+1 << " " << GetId((*ta.t)[1])+1 << " " << GetId((*ta.t)[2])+1 << ", edge [" << i1 << "->" << j1 << " " << i2 << "->" << j2 << "]");
 				ta = Next(ta).Adj(); 
 				if (count++>50) _error2_("Maximum number of iteration exceeded");
@@ -2398,10 +2398,10 @@
 		int i;
 
-		printf("Mesh Echo:\n");
-		printf("   nbv = %i\n",nbv);
-		printf("   nbt = %i\n",nbt);
-		printf("   nbe = %i\n",nbe);
-		printf("   nbq = %i\n",nbq);
-		printf("   index:\n");
+		_printLine_("Mesh Echo:");
+		_printLine_("   nbv = " << nbv);
+		_printLine_("   nbt = " << nbt);
+		_printLine_("   nbe = " << nbe);
+		_printLine_("   nbq = " << nbq);
+		_printLine_("   index:");
 		for (i=0;i<nbt;i++){
 			_printLine_("   " << setw(4) << i+1 << ": [" 
@@ -2410,5 +2410,5 @@
 						<< setw(4) << (((BamgVertex *)triangles[i](0))?GetId(triangles[i][2])+1:0) << "]");
 		}
-		printf("   coordinates:\n");
+		_printLine_("   coordinates:");
 		for (i=0;i<nbv;i++){
 			_printLine_("   " << setw(4) << i+1 << ": [" << vertices[i].r.x << " " << vertices[i].r.y << "]");
@@ -2426,5 +2426,5 @@
 
 			//display
-			if (verbose > 2) printf("   ForceBoundary  nb of edge: %i\n",nbe);
+			if (verbose > 2) _printLine_("   ForceBoundary  nb of edge: " << nbe);
 
 			//check that there is no triangle with 0 determinant
@@ -2457,5 +2457,5 @@
 				Nbswap +=  vertices[j].Optim(1,0);
 			}
-			if (verbose > 3) printf("      number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 
+			if (verbose > 3) _printLine_("      number of inforced edge = " << nbfe << ", number of swap= " << Nbswap); 
 		}
 	/*}}}*/
@@ -2467,6 +2467,6 @@
 
 		if (verbose >2){
-			if (OutSide) printf("   Find all external sub-domain\n"); 
-			else printf("   Find all internal sub-domain\n");
+			if (OutSide) _printLine_("   Find all external sub-domain"); 
+			else _printLine_("   Find all internal sub-domain");
 		  }
 		short * HeapArete = new short[nbt];
@@ -2613,5 +2613,5 @@
 								 }//while (t)
 								}
-							if(verbose>4) printf("      Number of removes subdomains (OutSideMesh) = %i\n",nbsubdomains-j);
+							if(verbose>4) _printLine_("      Number of removes subdomains (OutSideMesh) = " << nbsubdomains-j);
 							nbsubdomains=j;
 						  }
@@ -2692,5 +2692,5 @@
 
 				if (inew < nbsubdomains) {
-					if (verbose>5) printf("WARNING: %i SubDomains are being removed\n",nbsubdomains-inew);
+					if (verbose>5) _printLine_("WARNING: " << nbsubdomains-inew << " SubDomains are being removed");
 					nbsubdomains=inew;}
 
@@ -2800,5 +2800,5 @@
 
 		//Display info
-		if (verbose>2) printf("   Insert initial %i vertices\n",nbv);
+		if (verbose>2) _printLine_("   Insert initial " << nbv << " vertices");
 
 		//Compute integer coordinates for the existing vertices
@@ -2891,5 +2891,5 @@
 		/*Now, add the vertices One by One*/
 		long NbSwap=0;
-		if (verbose>3) printf("   Begining of insertion process...\n");
+		if (verbose>3) _printLine_("   Begining of insertion process...");
 
 		for (int icount=2; icount<nbv; icount++) {
@@ -2914,6 +2914,6 @@
 		//Display info
 		if (verbose>3) {
-			printf("      NbSwap of insertion: %i\n",NbSwap);
-			printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
+			_printLine_("      NbSwap of insertion: " << NbSwap);
+			_printLine_("      NbSwap/nbv:          " << NbSwap/nbv);
 		}
 
@@ -2929,6 +2929,6 @@
 			 NbSwap += orderedvertices[is1]->Optim(0,0);
 			if (verbose>3) {
-				printf("      Optim Loop: %i\n",Nbloop);
-				printf("      NbSwap/nbv:          %i\n",NbSwap/nbv);
+				_printLine_("      Optim Loop: " << Nbloop);
+				_printLine_("      NbSwap/nbv:          " << NbSwap/nbv);
 			}
 			if(!NbSwap) break;
@@ -2953,5 +2953,5 @@
 
 		//display info if required
-		if (verbose>5) printf("      Try to Insert %i new points\n",nbvnew);
+		if (verbose>5) _printLine_("      Try to Insert " << nbvnew << " new points");
 
 		//return if no new points
@@ -3003,12 +3003,12 @@
 		} 
 		if (verbose>3) {
-			printf("         number of new points: %i\n",iv);
-			printf("         number of to close (?) points: %i\n",nbv-iv);
-			printf("         number of swap after: %i\n",NbSwap);
+			_printLine_("         number of new points: " << iv);
+			_printLine_("         number of to close (?) points: " << nbv-iv);
+			_printLine_("         number of swap after: " << NbSwap);
 		}
 		nbv = iv;
 
 		for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);  
-		if (verbose>3) printf("   NbSwap=%i\n",NbSwap);
+		if (verbose>3) _printLine_("   NbSwap=" << NbSwap);
 
 		NbTSwap +=  NbSwap ;
@@ -3056,5 +3056,5 @@
 			if (!e[i]){
 				kk++;
-				if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);
+				if(kk<10) _printLine_("BUG: the geometrical edge " << i << " is on no edge curve");
 			}
 		}
@@ -3070,8 +3070,8 @@
 		long int verbose=0;
 
-		if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta);
+		if (verbose>2) _printLine_("MakeQuadrangles costheta = " << costheta);
 
 		if (costheta >1) {
-			if (verbose>5) printf("   do nothing: costheta > 1\n");
+			if (verbose>5) _printLine_("   do nothing: costheta > 1");
 		}
 
@@ -3100,7 +3100,7 @@
 			nbq = kk;
 			if (verbose>2){
-				printf("   number of quadrilaterals    = %i\n",nbq);
-				printf("   number of triangles         = %i\n",nbt-nbtout- nbq*2);
-				printf("   number of outside triangles = %i\n",nbtout);
+				_printLine_("   number of quadrilaterals    = " << nbq);
+				_printLine_("   number of triangles         = " << nbt-nbtout- nbq*2);
+				_printLine_("   number of outside triangles = " << nbtout);
 			}
 			delete [] qq;
@@ -3128,5 +3128,5 @@
 
 		const  double maxsubdiv2 = maxsubdiv*maxsubdiv;
-		if(verbose>1) printf("   Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);
+		if(verbose>1) _printLine_("   Limit the subdivision of a edges in the new mesh by " << maxsubdiv);
 		// for all the edges 
 		// if the len of the edge is to long 
@@ -3171,5 +3171,5 @@
 		}
 		if(verbose>3){
-			printf("      number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));
+			_printLine_("      number of metric changes = " << nbchange << ", maximum number of subdivision of a edges before change = " << pow(lmax,0.5));
 		}
 	}
@@ -3225,5 +3225,5 @@
 		/*First, insert old points if requested*/
 		if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){
-			if (verbose>5) printf("         Inserting initial mesh points\n");
+			if (verbose>5) _printLine_("         Inserting initial mesh points");
 			for (i=0;i<Bh.nbv;i++){ 
 				BamgVertex &bv=Bh[i];
@@ -3244,5 +3244,5 @@
 		// Big loop (most time consuming)
 		int iter=0;
-		if (verbose>5) printf("         Big loop\n");
+		if (verbose>5) _printLine_("         Big loop");
 		do {
 			/*Update variables*/
@@ -3475,5 +3475,5 @@
 	// find extrema coordinates of vertices pmin,pmax
 	long i;
-	if(verbose>2) printf("      Reconstruct mesh of %i vertices\n",nbv); 
+	if(verbose>2) _printLine_("      Reconstruct mesh of " << nbv << " vertices"); 
 
 	//initialize orderedvertices
@@ -3532,10 +3532,10 @@
 	//Display info if required
 	if(verbose>5) {
-		printf("         info of Mesh:\n");
-		printf("            - number of vertices    = %i \n",nbv); 
-		printf("            - number of triangles   = %i \n",nbt); 
-		printf("            - number of given edges = %i \n",nbe); 
-		printf("            - number of all edges   = %i \n"  ,edge4->nb()); 
-		printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv); 
+		_printLine_("         info of Mesh:");
+		_printLine_("            - number of vertices    = " << nbv << " "); 
+		_printLine_("            - number of triangles   = " << nbt << " "); 
+		_printLine_("            - number of given edges = " << nbe << " "); 
+		_printLine_("            - number of all edges   = " << edge4->nb()); 
+		_printLine_("            - Euler number 1 - nb of holes = " << nbt-edge4->nb()+nbv); 
 	}
 
@@ -3554,8 +3554,8 @@
 				if (k<10) {
 					//print only 10 edges
-					printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
+					_printLine_("Lost boundary edges " << i << " : " << edge4->i(i) << " " << edge4->j(i));
 				}
 				else if (k==10){
-					printf("Other lost boundary edges not shown...\n");
+					_printLine_("Other lost boundary edges not shown...");
 				}
 			}
@@ -3745,14 +3745,14 @@
 					/*Check that the 2 vertices are on geometry AND required*/
 					if(!edges[i][j].GeomEdgeHook->IsRequiredVertex()){
-						printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,GetId(edges[i][0])+1,GetId(edges[i][1])+1);
-						printf("This edge is on geometrical edge number %i\n",Gh.GetId(edges[i].GeomEdgeHook)+1);
+						_printLine_("ReconstructExistingMesh error message: problem with the edge number " << i+1 << ": [" << GetId(edges[i][0])+1 << " " << GetId(edges[i][1])+1 << "]");
+						_printLine_("This edge is on geometrical edge number " << Gh.GetId(edges[i].GeomEdgeHook)+1);
 						if (edges[i][j].GeomEdgeHook->OnGeomVertex())
-						 printf("the vertex number %i of this edge is a geometric BamgVertex number %i\n",GetId(edges[i][j])+1,Gh.GetId(edges[i][j].GeomEdgeHook->gv)+1);
+						 _printLine_("the vertex number " << GetId(edges[i][j])+1 << " of this edge is a geometric BamgVertex number " << Gh.GetId(edges[i][j].GeomEdgeHook->gv)+1);
 						else if (edges[i][j].GeomEdgeHook->OnGeomEdge())
-						 printf("the vertex number %i of this edge is a geometric Edge number %i\n",GetId(edges[i][j])+1,Gh.GetId(edges[i][j].GeomEdgeHook->ge)+1);
+						 _printLine_("the vertex number " << GetId(edges[i][j])+1 << " of this edge is a geometric Edge number " << Gh.GetId(edges[i][j].GeomEdgeHook->ge)+1);
 						else
-						 printf("Its pointer is %p\n",edges[i][j].GeomEdgeHook);
-
-						printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
+						 _printLine_("Its pointer is " << edges[i][j].GeomEdgeHook);
+
+						_printLine_("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required");
 						_error2_("See above (might be cryptic...)");
 					}
@@ -3842,13 +3842,13 @@
 		long it,ie,i;
 
-		printf("renumbering triangles\n");
+		_printLine_("renumbering triangles");
 		for ( it=0;it<nbt;it++) 
 		 triangles[it].Renumbering(vertices,ve,renu);
 
-		printf("renumbering edges\n");
+		_printLine_("renumbering edges");
 		for ( ie=0;ie<nbe;ie++) 
 		 edges[ie].Renumbering(vertices,ve,renu);
 
-		printf("renumbering vertices on geom\n");
+		_printLine_("renumbering vertices on geom");
 		for (i=0;i< NbVerticesOnGeomVertex;i++)
 		  {
@@ -3858,5 +3858,5 @@
 		  }
 
-		printf("renumbering vertices on edge\n");
+		_printLine_("renumbering vertices on edge");
 		for (i=0;i< NbVerticesOnGeomEdge;i++)
 		  {
@@ -3866,5 +3866,5 @@
 		  }
 
-		printf("renumbering vertices on Bth vertex\n");
+		_printLine_("renumbering vertices on Bth vertex");
 		for (i=0;i< NbVertexOnBThVertex;i++)
 		  {
@@ -3956,5 +3956,5 @@
 				number_of_errors++;
 				if (number_of_errors<20){
-					printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det);
+					_printLine_("Area of Triangle " << i+1 << " < 0 (det=" << triangles[i].det << ")");
 				}
 			}
@@ -4024,10 +4024,10 @@
 	gammamn=sqrt(gammamn);
 	gammamx=sqrt(gammamx);    
-	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));
+	_printLine_("   Adaptmesh info:");
+	_printLine_("      number of triangles = " << nt);
+	_printLine_("      hmin = " << hmin << ", hmax=" << hmax);
+	_printLine_("      area = " << area << ", M area = " << Marea << ", M area/( |Khat| nt) = " <<  Marea/(aireKh*nt));
+	_printLine_("      infinite-regularity(?): min = " << gammamn << ", max = " << gammamx);
+	_printLine_("      anisomax = " << pow(alpha2,0.5) << ", beta max = " << 1./pow(beta/aireKh,0.5) << ", min = " << 1./pow(beta0/aireKh,0.5));
 }
 /*}}}*/
@@ -4063,16 +4063,16 @@
 			}
 		}  
-	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"); 
+	_printLine_(" --- Histogram of the unit mesh,  nb of edges = " << nbedges);
+	_printLine_("      length of edge in   | %% of edge  | Nb of edges "); 
+	_printLine_("      --------------------+-------------+-------------"); 
 	for   (i=0;i<=kmax;i++){ 
 		if (i==0) _printString_( "      " << setw(10) << 0.);
 		else      _printString_( "      " << setw(10) << exp(lmin+i/delta));
-		if (i==kmax) printf("          +inf   ");
+		if (i==kmax) _printString_("          +inf   ");
 		else      _printString_( "      " << setw(10) << exp(lmin+(i+1)/delta));
 		_printLine_("|  " << setw(10) << (long((10000. * histo[i])/ nbedges)/100.) << " |");
-		printf("  %i\n",histo[i]);
-	}
-	printf("      --------------------+-------------+-------------\n"); 
+		_printLine_("  " << histo[i]);
+	}
+	_printLine_("      --------------------+-------------+-------------"); 
 }
 /*}}}*/
@@ -4100,5 +4100,5 @@
 	for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 
 	 tstart[ GetId(VerticesOnGeomEdge[k].meshvertex)]=&vide;
-	if(verbose>2) printf("   SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);
+	if(verbose>2) _printLine_("   SmoothingVertex: nb Iteration = " << nbiter << ", Omega=" << omega);
 	for (k=0;k<nbiter;k++)
 	  {
@@ -4112,5 +4112,5 @@
 		  if (tstart[i] != &vide) // not a boundary vertex 
 			NbSwap += vertices[i].Optim(1);
-		if (verbose>3) printf("      move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);
+		if (verbose>3) _printLine_("      move max = " << pow(delta << ", iteration = " << 0.5) << ", nb of swap = " << k);
 	  }
 
@@ -4126,5 +4126,5 @@
 
 	if(raisonmax<1.1) return;
-	if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
+	if(verbose > 1) _printLine_("   Mesh::SmoothMetric raisonmax = " << raisonmax);
 	CreateSingleVertexToTriangleConnectivity();
 	long i,j,kch,kk,ip;
@@ -4192,5 +4192,5 @@
 		Exchange(first_np_or_next_t0,first_np_or_next_t1);
 	}
-	if(verbose>2) printf("      number of iterations = %i\n",kch); 
+	if(verbose>2) _printLine_("      number of iterations = " << kch); 
 	delete [] first_np_or_next_t0;
 	delete [] first_np_or_next_t1;
@@ -4382,5 +4382,5 @@
 										||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0  
 										||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
-							printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd);
+							_printLine_("" << ke + nbvold << " not in triangle " << i << " In= " << !!t.link << " " << aa << " " << bb << " " << cc << " " << dd);
 							_error2_("Number of triangles with P2 interpolation Problem");
 						}
@@ -4390,5 +4390,5 @@
 										||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0 
 										||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
-							printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd);
+							_printLine_("" << ke + nbvold << " not in triangle " << ii << " In= " << !!tt.link << " " << aa << " " << bb << " " << cc << " " << dd);
 							_error2_("Number of triangles with P2 interpolation Problem");
 						}
@@ -4685,7 +4685,7 @@
 
 		if (verbose>2){
-			printf("   number of quadrilaterals    = %i\n",nbq);
-			printf("   number of triangles         = %i\n",nbt-nbtout- nbq*2);
-			printf("   number of outside triangles = %i\n",nbtout);
+			_printLine_("   number of quadrilaterals    = " << nbq);
+			_printLine_("   number of triangles         = " << nbt-nbtout- nbq*2);
+			_printLine_("   number of outside triangles = " << nbtout);
 		}
 
@@ -4754,5 +4754,5 @@
 			Triangle *tcvi = TriangleFindFromCoord(vi.i,det3);
 			if (tcvi && !tcvi->link) {
-				printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
+				_printLine_("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)");
 			}
 
@@ -4767,11 +4767,11 @@
 	}
 	if (verbose>3) {
-		printf("   number of points: %i\n",iv);
-		printf("   number of swap to  split internal edges with border vertices: %i\n",NbSwap);
+		_printLine_("   number of points: " << iv);
+		_printLine_("   number of swap to  split internal edges with border vertices: " << NbSwap);
 		nbv = iv;
 	}
 }
-if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices  to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold));
-if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge);
+if (NbSplitEdge>nbv-nbvold) _printLine_("WARNING: not enough vertices  to split all internal edges, we lost " << NbSplitEdge - ( nbv-nbvold) << " edges...");
+if (verbose>2) _printLine_("SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge);
 
 return  NbSplitEdge;
@@ -4943,5 +4943,5 @@
 
 	//Vertices
-	if(verbose) printf("Reading vertices (%i)\n",nbv);
+	if(verbose) _printLine_("Reading vertices (" << nbv << ")");
 	for (i=0;i<nbv;i++){
 		vertices[i].r.x=x[i];
@@ -5255,20 +5255,20 @@
 
 		//Insert points inside existing triangles
-		if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
-		if (verbose>3) printf("      Creating initial Constrained Delaunay Triangulation...\n");
-		if (verbose>3) printf("         Inserting boundary points\n");
+		if (verbose>4) _printLine_("      -- current number of vertices = " << nbv);
+		if (verbose>3) _printLine_("      Creating initial Constrained Delaunay Triangulation...");
+		if (verbose>3) _printLine_("         Inserting boundary points");
 		Insert();
 
 		//Force the boundary
-		if (verbose>3) printf("         Forcing boundaries\n");
+		if (verbose>3) _printLine_("         Forcing boundaries");
 		ForceBoundary();
 
 		//Extract SubDomains
-		if (verbose>3) printf("         Extracting subdomains\n");
+		if (verbose>3) _printLine_("         Extracting subdomains");
 		FindSubDomain();
 
-		if (verbose>3) printf("      Inserting internal points\n");
+		if (verbose>3) _printLine_("      Inserting internal points");
 		NewPoints(*this,bamgopts,0) ;
-		if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
+		if (verbose>4) _printLine_("      -- current number of vertices = " << nbv);
 	}
 	/*}}}*/
@@ -5373,9 +5373,9 @@
 					int nc=ei.GeomEdgeHook->CurveNumber;
 					
-					//printf("Dealing with curve number %i\n",nc);
-					//printf("edge on geometry is same as GhCurve? %s\n",(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no");
+					//_printLine_("Dealing with curve number " << nc);
+					//_printLine_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no");
 					//if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
-					//	printf("Do we have the right extremity? curve first vertex -> %s\n",((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no");
-					//	printf("Do we have the right extremity? curve last  vertex -> %s\n",((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no");
+					//	_printLine_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no");
+					//	_printLine_("Do we have the right extremity? curve last  vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no");
 					//}
 					//BUG FIX from original bamg
@@ -5585,20 +5585,20 @@
 
 		//Insert points inside existing triangles
-		if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
-		if (verbose>3) printf("      Creating initial Constrained Delaunay Triangulation...\n");
-		if (verbose>3) printf("         Inserting boundary points\n");
+		if (verbose>4) _printLine_("      -- current number of vertices = " << nbv);
+		if (verbose>3) _printLine_("      Creating initial Constrained Delaunay Triangulation...");
+		if (verbose>3) _printLine_("         Inserting boundary points");
 		Insert();
 
 		//Force the boundary
-		if (verbose>3) printf("         Forcing boundaries\n");
+		if (verbose>3) _printLine_("         Forcing boundaries");
 		ForceBoundary();
 
 		//Extract SubDomains
-		if (verbose>3) printf("         Extracting subdomains\n");
+		if (verbose>3) _printLine_("         Extracting subdomains");
 		FindSubDomain();
 
-		if (verbose>3) printf("      Inserting internal points\n");
+		if (verbose>3) _printLine_("      Inserting internal points");
 		NewPoints(BTh,bamgopts,KeepVertices) ;
-		if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
+		if (verbose>4) _printLine_("      -- current number of vertices = " << nbv);
 	}
 	/*}}}*/
Index: /issm/trunk-jpl/src/c/objects/Bamg/Metric.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/Metric.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/Metric.cpp	(revision 12509)
@@ -72,6 +72,6 @@
 	void Metric::Echo(void){
 
-		printf("Metric:\n");
-		printf("   [a11 a21 a22]: [%g %g %g]\n",a11,a21,a22);
+		_printLine_("Metric:");
+		_printLine_("   [a11 a21 a22]: [" << a11 << " " << a21 << " " << a22 << "]");
 
 		return;
@@ -201,5 +201,5 @@
 		LastMetricInterpole.lab=l;
 		LastMetricInterpole.opt=i;
-		if (i>200 && kkk++<10) printf("WARNING: LengthInterpole: ( i=%i l=%i sss=%g ) %g\n",i,l,sss,sstop); 
+		if (i>200 && kkk++<10) _printLine_("WARNING: LengthInterpole: ( i=" << i << " l=" << l << " sss=" << sss << " ) " << sstop); 
 		return l;
 	}
Index: /issm/trunk-jpl/src/c/objects/Bamg/Triangle.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Bamg/Triangle.cpp	(revision 12508)
+++ /issm/trunk-jpl/src/c/objects/Bamg/Triangle.cpp	(revision 12509)
@@ -107,18 +107,18 @@
 		int i;
 
-		printf("Triangle:\n");
-		printf("   vertices pointer towards three vertices\n");
-		printf("      vertices[0] vertices[1] vertices[2] = %p %p %p\n",vertices[0],vertices[1],vertices[2]);
-		printf("   adj pointer towards three adjacent triangles\n");
-		printf("      adj[0] adj[1] adj[2] = %p %p %p\n",adj[0],adj[1],adj[2]);
-		printf("   det (integer triangle determinant) = %i\n",det);
+		_printLine_("Triangle:");
+		_printLine_("   vertices pointer towards three vertices");
+		_printLine_("      vertices[0] vertices[1] vertices[2] = " << vertices[0] << " " << vertices[1] << " " << vertices[2]);
+		_printLine_("   adj pointer towards three adjacent triangles");
+		_printLine_("      adj[0] adj[1] adj[2] = " << adj[0] << " " << adj[1] << " " << adj[2]);
+		_printLine_("   det (integer triangle determinant) = " << det);
 		if (link){
-			printf("   link (pointer toward duplicate triangle)= %p\n",link);
+			_printLine_("   link (pointer toward duplicate triangle)= " << link);
 		}
 		else{
-			printf("   color = %i\n",color);
-		}
-
-		printf("\nThree vertices:\n");
+			_printLine_("   color = " << color);
+		}
+
+		_printLine_("\nThree vertices:");
 		for(i=0;i<3;i++){
 			if (vertices[i]){
@@ -126,5 +126,5 @@
 			}
 			else{
-				printf("   vertex %i does not exist\n",i+1);
+				_printLine_("   vertex " << i+1 << " does not exist");
 			}
 		}