Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3147)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3148)
@@ -95,6 +95,4 @@
 		if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
 		Th.ReNumberingTheTriangleBySubDomain();
-		//if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
-		Th.SmoothingVertex(3,omega);
 
 		//Build output
@@ -159,5 +157,5 @@
 		if (verbosity>1) printf("   Generating Mesh...\n");
 		Thr=&BTh,Thb=0;
-		Triangles & Th( *(0 ?  new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) :  new Triangles(maxnbv,BTh,bamgopts->KeepVertices)));
+		Triangles & Th( *(0 ?  new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) :  new Triangles(maxnbv,BTh,bamgopts,bamgopts->KeepVertices)));
 		if (Thr != &BTh) delete Thr;
 		if(costheta<=1.0) Th.MakeQuadrangles(costheta);
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3147)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3148)
@@ -367,11 +367,11 @@
 		class IntersectionTriangles {
 			public: 
-				Triangle *t;
+				Triangle* t;
 				Real8  bary[3];  // use if t != 0
 				R2 x;
 				Metric m;
-				Real8 s;// abscisse curviline
-				Real8 sp; // len of the previous seg in m
-				Real8 sn;// len of the  next seg in m
+				Real8 s; // curvilinear coordinate 
+				Real8 sp;// length of the previous segment in m
+				Real8 sn;// length of the next segment in m
 		};
 
@@ -380,5 +380,5 @@
 				GeometricalEdge * e;
 				Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
-				Real8 lBegin,lEnd; // length abscisse  set in ListofIntersectionTriangles::Length
+				Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
 				int last;// last index  in ListofIntersectionTriangles for this Sub seg of edge
 				R2 F(Real8 s){ 
@@ -390,7 +390,7 @@
 		};
 
-		int MaxSize; // 
-		int Size; //
-		Real8 len; //
+		int MaxSize;
+		int Size;
+		Real8 len;
 		int state;
 		IntersectionTriangles * lIntTria;
@@ -403,9 +403,8 @@
 
 		ListofIntersectionTriangles(int n=256,int m=16)
-		  :   MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
-		  NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m])  
-		{
-		 long int verbosity=0;
-		 if (verbosity>9) printf("   construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
+		  : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
+		  NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
+			  long int verbosity=0;
+			  if (verbosity>9) printf("   construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
 		};
 
@@ -417,6 +416,5 @@
 			int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
 			int NewItem(R2,const Metric & );
-			void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1) 
-			  { 
+			void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){ 
 				long int verbosity=0;
 
@@ -436,6 +434,5 @@
 					lSegsI = lEn;        
 				}
-				if (NbSeg) 
-				 lSegsI[NbSeg-1].last=Size;
+				if (NbSeg) lSegsI[NbSeg-1].last=Size;
 				lSegsI[NbSeg].e=e;
 				lSegsI[NbSeg].sBegin=s0;
@@ -443,7 +440,4 @@
 				NbSeg++;           
 			  }
-
-			//  void CopyMetric(int i,int j){ lIntTria[j].m=lIntTria[i].m;}
-			//  void CopyMetric(const Metric & mm,int j){ lIntTria[j].m=mm;}
 
 			void ReShape() { 
@@ -453,6 +447,6 @@
 					throw ErrorException(__FUNCT__,exprintf("!nw"));
 				}
-				for (int i=0;i<MaxSize;i++) // recopy
-				 nw[i] = lIntTria[i];       
+				// recopy
+				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);
@@ -597,5 +591,4 @@
 		public:
 
-			int static counter; // to kown the number of mesh in memory 
 			Geometry & Gh;   // Geometry
 			Triangles & BTh; // Background Mesh Bth==*this =>no  background 
@@ -650,6 +643,6 @@
 			Triangles(double* index,double* x,double* y,int nods,int nels);
 
-			Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
-				try {GeomToTriangles1(nbvx,keepBackVertices);}
+			Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
+				try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
 				catch(...) { this->~Triangles(); throw; }
 			}
@@ -700,7 +693,6 @@
 			int SplitElement(int choice);
 			void MakeQuadTree();
-			void NewPoints( Triangles &,int KeepVertices =1 );
+			void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1);
 			Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; 
-			void NewPoints(int KeepVertices=1){ NewPoints(*this,KeepVertices);}
 			void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
 			void ReNumberingVertex(Int4 * renu);
@@ -744,5 +736,5 @@
 			int  CrackMesh();
 		private:
-			void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption
+			void GeomToTriangles1(Int4 nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
 			void GeomToTriangles0(Int4 nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
 			void PreInit(Int4);
Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 3147)
+++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 3148)
@@ -263,17 +263,18 @@
 	/*FUNCTION ListofIntersectionTriangles::Length{{{1*/
 	Real8  ListofIntersectionTriangles::Length(){
+		// computation of the length
+
+		// check Size
 		if (Size<=0){
 			throw ErrorException(__FUNCT__,exprintf("Size<=0"));
 		}
-		// computation of the length      
-		R2 C;
+
 		Metric Mx,My;
 		int ii,jj;
 		R2 x,y,xy;
 
-		SegInterpolation *SegI=lSegsI;
+		SegInterpolation* SegI=lSegsI;
 		SegI=lSegsI;
-		lSegsI[NbSeg].last=Size;// improvement 
-
+		lSegsI[NbSeg].last=Size;
 		int EndSeg=Size;
 
@@ -285,19 +286,19 @@
 		for (jj=0,ii=1;ii<Size;jj=ii++) {  
 			// seg jj,ii
-			x=y;
-			y = lIntTria[ii].x;
+			x  = y;
+			y  = lIntTria[ii].x;
 			xy = y-x;
 			Mx = lIntTria[ii].m;
 			My = lIntTria[jj].m;
-			sxy =  LengthInterpole(Mx,My,xy);
+			sxy= LengthInterpole(Mx,My,xy);
 			s += sxy;
 			lIntTria[ii].s = s;
-			if (ii == EndSeg) 
-			 SegI->lEnd=s,
-				SegI++,
-				EndSeg=SegI->last,
+			if (ii == EndSeg){
+				SegI->lEnd=s;
+				SegI++;
+				EndSeg=SegI->last;
 				SegI->lBegin=s;
-
-		  }
+			}
+		}
 		len = s;
 		SegI->lEnd=s;
@@ -307,30 +308,30 @@
 	/*}}}1*/
 	/*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
-	Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4  nbvx){
-		long int verbosity=0;
-
-
-		const Int4 nbvold = nbv;
-		Real8 s = Length();
-		if (s <  1.5 ) return 0;
-		//////////////////////   
+	Int4 ListofIntersectionTriangles::NewPoints(Vertex* vertices,Int4 &nbv,Int4 nbvx){
+
+		//If length<1.5, do nothing
+		Real8 s=Length();
+		if (s<1.5) return 0;
+
+		const Int4 nbvold=nbv;
 		int ii = 1 ;
 		R2 y,x;
 		Metric My,Mx ;
 		Real8 sx =0,sy;
-		int nbi = Max(2,(int) (s+0.5));
-		Real8 sint = s/nbi;
-		Real8 si = sint;
+		int nbi=Max(2,(int) (s+0.5));
+		Real8 sint=s/nbi;
+		Real8 si  =sint;
 
 		int EndSeg=Size;
-		SegInterpolation *SegI=0;
-		if (NbSeg) 
-		 SegI=lSegsI,EndSeg=SegI->last;
-
-		for (int k=1;k<nbi;k++)
-		  {
-			while ((ii < Size) && ( lIntTria[ii].s <= si )) 
-			 if (ii++ == EndSeg) 
-			  SegI++,EndSeg=SegI->last;
+		SegInterpolation* SegI=NULL;
+		if (NbSeg) SegI=lSegsI,EndSeg=SegI->last;
+
+		for (int k=1;k<nbi;k++){
+			while ((ii < Size) && ( lIntTria[ii].s <= si )){
+				if (ii++ == EndSeg){
+					SegI++;
+					EndSeg=SegI->last;
+				}
+			}
 
 			int ii1=ii-1;
@@ -347,4 +348,8 @@
 			Real8 cx = 1-cy;
 			C = SegI ? SegI->F(si): x * cx + y *cy;
+			//C.Echo();
+			//x.Echo();
+			//y.Echo();
+			//printf("cx = %g, cy=%g\n",cx,cy);
 
 			si += sint;
Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 3147)
+++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 3148)
@@ -98,23 +98,22 @@
 	/*FUNCTION Triangle::Optim{{{1*/
 	Int4  Triangle::Optim(Int2 i,int koption) {
-		// turne around in positif sens
+		// turn around (positive direction)
+		Triangle *t=this;
 		Int4 NbSwap =0;
-		Triangle  *t = this;
-		int k=0,j =OppositeEdge[i];
-		int jp = PreviousEdge[j];
-		// initialise   tp, jp the previous triangle & edge
-		Triangle *tp= TriaAdjTriangles[jp];
+		int  k = 0;
+		int  j = OppositeEdge[i];
+		int  jp= PreviousEdge[j];
+
+		// initialize tp, jp the previous triangle & edge
+		Triangle *tp=TriaAdjTriangles[jp];
 		jp = TriaAdjSharedEdge[jp]&3;
 		do {
-			while (t->swap(j,koption))
-			  {
+			while (t->swap(j,koption)){
+				if (k>=20000) throw ErrorException(__FUNCT__,exprintf("k>=20000"));
 				NbSwap++;
 				k++;
-				if (k>=20000){
-					throw ErrorException(__FUNCT__,exprintf("k>=20000"));
-				}
 				t=  tp->TriaAdjTriangles[jp];      // set unchange t qnd j for previous triangles
 				j=  NextEdge[tp->TriaAdjSharedEdge[jp]&3];
-			  }
+			}
 			// end on this  Triangle 
 			tp = t;
@@ -187,6 +186,5 @@
 						 break;
 					 }
-					 else 
-						{	
+					 else {	
 						 // critere de Delaunay anisotrope 
 						 Real8 som;
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3147)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3148)
@@ -19,5 +19,4 @@
 	static const  Direction NoDirOfSearch=Direction();
 	long NbUnSwap =0;
-	int  Triangles::counter = 0;
 
 	/*Constructors/Destructors*/
@@ -3151,5 +3150,5 @@
 		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 points\n");
+		if (verbose>3) printf("         Inserting boundary points\n");
 		Insert();
 
@@ -3164,10 +3163,14 @@
 		// NewPointsOld(*this) ;
 		if (verbose>3) printf("      Inserting internal points\n");
-		NewPoints(*this,0) ;
+		NewPoints(*this,bamgopts,0) ;
 		if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
 	}
 	/*}}}1*/
 	/*FUNCTION Triangles::GeomToTriangles1{{{1*/
-	void Triangles::GeomToTriangles1(Int4 inbvx,int KeepVertices) { 
+	void Triangles::GeomToTriangles1(Int4 inbvx,BamgOpts* bamgopts,int KeepVertices){ 
+
+		/*Get options*/
+		int verbosity=bamgopts->verbose;
+
 		Gh.NbRef++;// add a ref to Gh
 
@@ -3194,5 +3197,4 @@
 		}
 
-		long int verbosity=0;
 		BTh.NbRef++; // add a ref to BackGround Mesh
 		PreInit(inbvx);
@@ -3501,5 +3503,5 @@
 		FindSubDomain();
 
-		NewPoints(BTh,KeepVertices) ;
+		NewPoints(BTh,bamgopts,KeepVertices) ;
 	}
 	/*}}}1*/
@@ -3973,30 +3975,28 @@
 /*}}}1*/
 	/*FUNCTION Triangles::NewPoints{{{1*/
-	void  Triangles::NewPoints(Triangles & Bh,int KeepVertices) {
-
-		long int verbosity=2;
-		Int4 nbtold(nbt),nbvold(nbv);
-		Int4 i,k;
-		int j;
-		Int4 *first_np_or_next_t = new Int4[nbtx];
-		Int4 NbTSwap =0;
-
-		//display info
-		if (verbosity>2)  printf("   Triangles::NewPoints\n");
-
-		/*First, insert old points*/
-
-		nbtold = nbt;
-
+	void  Triangles::NewPoints(Triangles & Bh,BamgOpts* bamgopts,int KeepVertices){
+
+		int i,j,k;
+		Int4 NbTSwap=0;
+		Int4 nbtold=nbt;
+		Int4 nbvold=nbv;
+		Int4 Headt=0;
+		Int4 next_t;
+		Int4* first_np_or_next_t=new Int4[nbtx];
+		Triangle* t=NULL;
+
+		/*Recover options*/
+		int verbosity=bamgopts->verbose;
+
+		/*First, insert old points if requested*/
 		if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){
-			//   Bh.SetVertexFieldOn();
+			if (verbosity>5) printf("         Inserting initial mesh points\n");
 			for (i=0;i<Bh.nbv;i++){ 
 				Vertex &bv=Bh[i];
 				if (!bv.onGeometry){
-					vertices[nbv].r = bv.r;
+					vertices[nbv].r   = bv.r;
 					vertices[nbv++].m = bv.m;
 				}
 			}
-			int nbv1=nbv;
 			Bh.ReMakeTriangleContainingTheVertex();     
 			InsertNewPoints(nbvold,NbTSwap)   ;            
@@ -4004,69 +4004,63 @@
 		else Bh.ReMakeTriangleContainingTheVertex();     
 
-		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 traingle on i is  -first_np_or_next_t[i]
+		for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
+		// the next traingle of i is -first_np_or_next_t[i]
+
+		// Big loop (most time consuming)
 		int iter=0;
-
-		// Big loop (most time consuming)
+		if (verbosity>5) printf("         Big loop\n");
 		do {
+			/*Update variables*/
 			iter++;
-			nbtold = nbt;
-			nbvold = nbv;
-
-			// default size of  IntersectionTriangle
+			nbtold=nbt;
+			nbvold=nbv;
+
+			/*We test all triangles*/
 			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]
+
+				//check i
 				if (i<0 || i>=nbt ){
-					throw ErrorException(__FUNCT__,exprintf("i<0 || i>=nbt"));
-				}
+					throw ErrorException(__FUNCT__,exprintf("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1));
+				}
+				//change first_np_or_next_t[i]
 				first_np_or_next_t[i] = iter; 
+
+				//Loop over the edges of t
 				for(j=0;j<3;j++){
-					// for each edge 
 					TriangleAdjacent tj(t,j);
-					Vertex & vA = * tj.EdgeVertex(0);
-					Vertex & vB = * tj.EdgeVertex(1);
-
-					if (!t->link) continue;// boundary
-					if (t->det <0) continue;
+					Vertex &vA = *tj.EdgeVertex(0);
+					Vertex &vB = *tj.EdgeVertex(1);
+
+					//if t is a boundary triangle, or tj locked, continue
+					if (!t->link)     continue;
+					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;
-
+					Triangle* ta=tadjj;
+
+					//if the adjacent triangle is a boundary triangle, continur
+					if (ta->det<0) continue;	  
+
+					R2 A=vA;
+					R2 B=vB;
 					k=Number(ta);
 
-					if(first_np_or_next_t[k]==iter)  // this edge is done before 
-					 continue; // next edge of the triangle 
-
-					//const Int4 NbvOld = nbv;
+					//if this edge has already been done, go to next edge of triangle
+					if(first_np_or_next_t[k]==iter) continue;
+
 					lIntTria.SplitEdge(Bh,A,B);
 					lIntTria.NewPoints(vertices,nbv,nbvx);
 				  } // end loop for each edge 
-
 			  }// for triangle   
 
-			if (!InsertNewPoints(nbvold,NbTSwap)) 
-			 break;
-
+			if (!InsertNewPoints(nbvold,NbTSwap)) break;
 			for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
-
 			Headt = nbt; // empty list 
 
-			// for all the triangle contening the vertex i
+			// for all the triangle containing the vertex i
 			for (i=nbvold;i<nbv;i++){ 
 				Vertex*          s  = vertices + i;
@@ -4085,5 +4079,5 @@
 					ta = Next(Adj(ta));
 				} while ( (tbegin != (Triangle*) ta)); 
-			}   
+			}
 
 		} while (nbv!=nbvold);