Index: /issm/trunk/src/c/objects/Bamg/BamgVertex.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/BamgVertex.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/BamgVertex.h	(revision 5130)
@@ -19,4 +19,5 @@
 		public:
 
+			/*Fields*/
 			I2        i;                 // integer coordinates
 			R2        r;                 // real coordinates
@@ -35,5 +36,5 @@
 			};
 
-			//Operators
+			/*Operators*/
 			operator I2() const {return i;}             // Cast operator
 			operator const R2 & () const {return r;}    // Cast operator
@@ -41,5 +42,5 @@
 			double operator()(R2 x) const { return m(x);} // Get x in the metric m
 
-			//methods (No constructor and no destructors...)
+			/*methods (No constructor and no destructors...)*/
 			double Smoothing(Mesh & ,const Mesh & ,Triangle  * & ,double =1);
 			void   MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
Index: /issm/trunk/src/c/objects/Bamg/Curve.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Curve.cpp	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Curve.cpp	(revision 5130)
@@ -10,6 +10,18 @@
 
 	/*Constructors/Destructors*/
+	/*FUNCTION Curve::Curve(){{{1*/
+	Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0),master(true){
+
+	} 
+	/*}}}*/
 
 	/*Methods*/
+	/*FUNCTION Curve::Reverse {{{1*/
+	void Curve::Reverse() {
+		/*reverse the direction of the curve */
+		Exchange(be,ee);
+		Exchange(kb,ke);
+	}
+	/*}}}*/
 	/*FUNCTION Curve::Set {{{1*/
 	void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew){
Index: /issm/trunk/src/c/objects/Bamg/Curve.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Curve.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Curve.h	(revision 5130)
@@ -19,6 +19,6 @@
 
 			//Methods
-			Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {} 
-			void Reverse() { Exchange(be,ee); Exchange(kb,ke);} //  revese the sens of the curse 
+			Curve();
+			void Reverse(void);
 			void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
 	};
Index: /issm/trunk/src/c/objects/Bamg/Direction.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Direction.cpp	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Direction.cpp	(revision 5130)
@@ -9,4 +9,9 @@
 
 	/*Constructors/Destructors*/
+	/*FUNCTION Direction() {{{1*/
+	Direction::Direction():
+		dir(MaxICoor){
+
+	}/*}}}*/
 	/*FUNCTION Direction(Icoor1 i,Icoor1 j) {{{1*/
 	Direction::Direction(Icoor1 i,Icoor1 j) {
Index: /issm/trunk/src/c/objects/Bamg/Direction.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Direction.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Direction.h	(revision 5130)
@@ -13,5 +13,5 @@
 		public:
 			//Methods
-			Direction(): dir(MaxICoor){}; //  no direction set
+			Direction();
 			Direction(Icoor1 i,Icoor1 j);
 			int sens(Icoor1 i,Icoor1 j);
Index: /issm/trunk/src/c/objects/Bamg/Edge.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Edge.cpp	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Edge.cpp	(revision 5130)
@@ -33,4 +33,22 @@
 	}
 	/*}}}*/
+	/*FUNCTION Edge::Renumbering{{{1*/
+	void Edge::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
+
+		if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
+		if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
+
+	}
+	/*}}}*/
+	/*FUNCTION Edge::Intersection{{{1*/
+	int Edge::Intersection(const  Edge & e){ 
+
+		/*some shecks*/
+		if (!(adj[0]==&e || adj[1]==&e)){ ISSMERROR("Intersection bug"); }
+		ISSMASSERT(adj[0]==&e || adj[1]==&e);
+
+		return adj[0]==&e?0:1;
+	}
+	/*}}}*/
 
 } 
Index: /issm/trunk/src/c/objects/Bamg/Edge.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Edge.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Edge.h	(revision 5130)
@@ -28,17 +28,6 @@
 
 			//Methods
-			void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
-				if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
-				if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
-			}
-			int Intersection(const  Edge & e) const { 
-				if (!(adj[0]==&e || adj[1]==&e)){
-					ISSMERROR("Intersection bug");
-				}
-				if (adj[0]!=&e && adj[1]!=&e){
-					ISSMERROR("adj[0]!=&e && adj[1]!=&e");
-				}
-				return adj[0]==&e ? 0 : 1;
-			}
+			void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu);
+			int  Intersection(const  Edge & e);
 			void Set(const Mesh &,long,Mesh &);
 			void Echo(void);
Index: /issm/trunk/src/c/objects/Bamg/GeometricalEdge.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/GeometricalEdge.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/GeometricalEdge.h	(revision 5130)
@@ -30,13 +30,13 @@
 			double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
 			int    Tg(int i) const{return i==0 ? TgA() : TgB();}
-			int    Cracked() const    {return flag & 1;  }
-			int    TgA()const         {return flag & 4;  }
-			int    TgB()const         {return flag & 8;  }
-			int    Mark()const        {return flag & 16; }
-			int    Required()         {return flag & 64; }
-			void   SetCracked()     { flag |= 1;  }
-			void   SetTgA()         { flag |=4;   }
-			void   SetTgB()         { flag |=8;   }
-			void   SetMark()        { flag |=16;  }
+			int    Cracked() const{return flag & 1;}
+			int    TgA()const{return flag & 4; }
+			int    TgB()const {return flag & 8;}
+			int    Mark()const {return flag & 16;}
+			int    Required() {return flag & 64;}
+			void   SetCracked()     { flag |= 1;}
+			void   SetTgA()         { flag |=4; }
+			void   SetTgB()         { flag |=8; }
+			void   SetMark()        { flag |=16;}
 			void   SetUnMark()      { flag &= 1007 /* 1023-16*/;}
 			void   SetRequired()    { flag |= 64; }
Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5130)
@@ -13,7 +13,10 @@
 
 	/*Constructors/Destructors*/
+	/*FUNCTION Geometry::Geometry(){{{1*/
 	Geometry::Geometry(){
 		Init();
 	}
+	/*}}}*/
+	/*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{1*/
 	Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
 		Init();
@@ -21,5 +24,6 @@
 		AfterRead();
 	}
-	/*FUNCTION  Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/
+	/*}}}*/
+	/*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/
 	Geometry::Geometry(const Geometry & Gh) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/
@@ -405,5 +409,5 @@
 
 	/*Methods*/
-	/*FUNCTION  Geometry::AfterRead(){{{1*/
+	/*FUNCTION Geometry::AfterRead(){{{1*/
 	void Geometry::AfterRead(){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
@@ -719,5 +723,5 @@
 	}
 	/*}}}1*/
-	/*FUNCTION  Geometry::Containing{{{1*/
+	/*FUNCTION Geometry::Containing{{{1*/
 	GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
@@ -745,5 +749,5 @@
 	}
 	/*}}}1*/
-	/*FUNCTION  Geometry::Echo {{{1*/
+	/*FUNCTION Geometry::Echo {{{1*/
 
 	void Geometry::Echo(void){
@@ -768,5 +772,5 @@
 	}
 	/*}}}*/
-	/*FUNCTION  Geometry::Init{{{1*/
+	/*FUNCTION Geometry::Init{{{1*/
 	void Geometry::Init(void){
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
@@ -786,5 +790,33 @@
 	}
 	/*}}}1*/
-	/*FUNCTION  Geometry::ProjectOnCurve {{{1*/
+	/*FUNCTION Geometry::MinimalHmin{{{1*/
+	double Geometry::MinimalHmin() {
+		return 2.0/coefIcoor;
+	}/*}}}*/
+	/*FUNCTION Geometry::MaximalHmax{{{1*/
+	double Geometry::MaximalHmax() {
+		return Max(pmax.x-pmin.x,pmax.y-pmin.y);
+	}/*}}}*/
+	/*FUNCTION Geometry::Number(const GeometricalVertex &t){{{1*/
+	long Geometry::Number(const GeometricalVertex & t) const  {
+		return &t - vertices;
+	}/*}}}*/
+	/*FUNCTION Geometry::Number(const GeometricalVertex * t){{{1*/
+	long Geometry::Number(const GeometricalVertex * t) const  {
+		return t - vertices;
+	}/*}}}*/
+	/*FUNCTION Geometry::Number(const GeometricalEdge & t){{{1*/
+	long Geometry::Number(const GeometricalEdge & t) const  {
+		return &t - edges;
+	}/*}}}*/
+	/*FUNCTION Geometry::Number(const GeometricalEdge * t){{{1*/
+	long Geometry::Number(const GeometricalEdge * t) const  {
+		return t - edges;
+	}/*}}}*/
+	/*FUNCTION Geometry::Number(const Curve * c){{{1*/
+	long Geometry::Number(const Curve * c) const  {
+		return c - curves;
+	}/*}}}*/
+	/*FUNCTION Geometry::ProjectOnCurve {{{1*/
 	GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
@@ -844,5 +876,5 @@
 					printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n");
 					ISSMERROR("see above");
-				  }
+				}
 				NbTry++;
 				goto retry;
@@ -920,4 +952,12 @@
 	}
 	/*}}}1*/
-
+	/*FUNCTION Geometry::toI2{{{1*/
+	I2 Geometry::toI2(const R2 & P) const {
+		return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
+					,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
+	}/*}}}*/
+	/*FUNCTION Geometry::UnMarkEdges{{{1*/
+	void Geometry::UnMarkEdges() {
+		for (int i=0;i<nbe;i++) edges[i].SetUnMark();
+	}/*}}}*/
 } 
Index: /issm/trunk/src/c/objects/Bamg/Geometry.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.h	(revision 5130)
@@ -42,29 +42,26 @@
 
 			//Operators
-			const GeometricalVertex & operator[]  (long i) const { return vertices[i];};
-			GeometricalVertex & operator[](long i) { return vertices[i];};
-			const  GeometricalEdge & operator()  (long i) const { return edges[i];};
-			GeometricalEdge & operator()(long i) { return edges[i];}; 
+			const GeometricalVertex &operator[](long i) const { return vertices[i]; };
+			GeometricalVertex       &operator[](long i) { return vertices[i];       };
+			const GeometricalEdge   &operator()(long i) const { return edges[i];    };
+			GeometricalEdge         &operator()(long  i) { return edges[i];                };
 
 			//Methods
-			void Echo();
-			I2 toI2(const R2 & P) const {
-				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
-							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
-			}
-			double MinimalHmin() {return 2.0/coefIcoor;}
-			double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
-			void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
-			void Init(void);
-			void AfterRead();
-			long Number(const GeometricalVertex & t) const  { return &t - vertices;}
-			long Number(const GeometricalVertex * t) const  { return t - vertices;}
-			long Number(const GeometricalEdge & t) const  { return &t - edges;}
-			long Number(const GeometricalEdge * t) const  { return t - edges;}
-			long Number(const Curve * c) const  { return c - curves;}
-			void UnMarkEdges() {for (int i=0;i<nbe;i++) edges[i].SetUnMark();}
-			GeometricalEdge *  ProjectOnCurve(const Edge & ,double,BamgVertex &,VertexOnGeom &) const ;
-			GeometricalEdge *  Containing(const R2 P,  GeometricalEdge * start) const;
-			void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
+			void             Echo();
+			I2               toI2(const R2 &P) const;
+			double           MinimalHmin();
+			double           MaximalHmax();
+			void             ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
+			void             Init(void);
+			void             AfterRead();
+			long             Number(const GeometricalVertex &t) const;
+			long             Number(const GeometricalVertex *t) const;
+			long             Number(const GeometricalEdge &t) const;
+			long             Number(const GeometricalEdge *t) const;
+			long             Number(const Curve *c) const;
+			void             UnMarkEdges();
+			GeometricalEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const;
+			GeometricalEdge *Containing(const R2 P, GeometricalEdge *start) const;
+			void             WriteGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
 	};
 	
Index: /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp	(revision 5130)
@@ -8,5 +8,216 @@
 namespace bamg {
 
+	/*Constructors Destructors*/
+	/*FUNCTION ListofIntersectionTriangles::ListofIntersectionTriangles{{{1*/
+	ListofIntersectionTriangles::ListofIntersectionTriangles(int n,int m)
+	  : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
+	  NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
+	  }
+	/*}}}*/
+	/*FUNCTION ListofIntersectionTriangles::~ListofIntersectionTriangles{{{1*/
+	ListofIntersectionTriangles::~ListofIntersectionTriangles(){
+		if (lIntTria) delete [] lIntTria,lIntTria=0;
+		if (lSegsI) delete [] lSegsI,lSegsI=0;
+	} 
+	/*}}}*/
+
 	/*Methods*/
+	/*FUNCTION ListofIntersectionTriangles::Init{{{1*/
+	void ListofIntersectionTriangles::Init(void){
+		state=0;
+		len=0;
+		Size=0;
+	}
+	/*}}}*/
+	/*FUNCTION ListofIntersectionTriangles::Length{{{1*/
+	double  ListofIntersectionTriangles::Length(){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Length)*/
+
+		// computation of the length
+
+		// check Size
+		if (Size<=0){
+			ISSMERROR("Size<=0");
+		}
+
+		Metric Mx,My;
+		int ii,jj;
+		R2 x,y,xy;
+
+		SegInterpolation* SegI=lSegsI;
+		SegI=lSegsI;
+		lSegsI[NbSeg].last=Size;
+		int EndSeg=Size;
+
+		y = lIntTria[0].x;
+		double sxy, s = 0;
+		lIntTria[0].s =0;
+		SegI->lBegin=s;
+
+		for (jj=0,ii=1;ii<Size;jj=ii++) {  
+			// seg jj,ii
+			x  = y;
+			y  = lIntTria[ii].x;
+			xy = y-x;
+			Mx = lIntTria[ii].m;
+			My = lIntTria[jj].m;
+			sxy= LengthInterpole(Mx,My,xy);
+			s += sxy;
+			lIntTria[ii].s = s;
+			if (ii == EndSeg){
+				SegI->lEnd=s;
+				SegI++;
+				EndSeg=SegI->last;
+				SegI->lBegin=s;
+			}
+		}
+		len = s;
+		SegI->lEnd=s;
+
+		return s;
+	}
+	/*}}}1*/
+	/*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {{{1*/
+	int  ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) { 
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
+
+		register int n;
+		R2 x(0,0);
+		if ( d0) x =      (*tt)[0].r * d0;
+		if ( d1) x = x +  (*tt)[1].r * d1;
+		if ( d2) x = x +  (*tt)[2].r * d2;
+		// newer add same point 
+		if(!Size ||  Norme2_2(lIntTria[Size-1].x-x)) {
+			if (Size==MaxSize) ReShape();
+			lIntTria[Size].t=tt;
+			lIntTria[Size].bary[0]=d0;
+			lIntTria[Size].bary[1]=d1;
+			lIntTria[Size].bary[2]=d2;
+			lIntTria[Size].x = x;
+			Metric m0,m1,m2;
+			register BamgVertex * v;
+			if ((v=(*tt)(0))) m0    = v->m;
+			if ((v=(*tt)(1))) m1    = v->m;
+			if ((v=(*tt)(2))) m2    = v->m;
+			lIntTria[Size].m =  Metric(lIntTria[Size].bary,m0,m1,m2);
+			n=Size++;}
+		else n=Size-1;
+		return n;
+	}
+	/*}}}1*/
+	/*FUNCTION ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm){{{1*/
+	int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
+
+		register int n;
+		if(!Size ||  Norme2_2(lIntTria[Size-1].x-A)) {
+			if (Size==MaxSize) ReShape();
+			lIntTria[Size].t=0;
+			lIntTria[Size].x=A;
+			lIntTria[Size].m=mm;
+			n=Size++;
+		}
+		else  n=Size-1;
+		return  n; 
+	}
+	/*}}}1*/
+	/*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
+	long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
+
+		//If length<1.5, do nothing
+		double s=Length();
+		if (s<1.5) return 0;
+
+		const long nbvold=nbv;
+		int ii = 1 ;
+		R2 y,x;
+		Metric My,Mx ;
+		double sx =0,sy;
+		int nbi=Max(2,(int) (s+0.5));
+		double sint=s/nbi;
+		double si  =sint;
+
+		int EndSeg=Size;
+		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;
+			x  =lIntTria[ii1].x;
+			sx =lIntTria[ii1].s;
+			Metric Mx=lIntTria[ii1].m;
+			y  =lIntTria[ii].x;
+			sy =lIntTria[ii].s;
+			Metric My=lIntTria[ii].m;
+			double lxy = sy-sx;
+			double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
+
+			R2 C;
+			double 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;
+			if ( nbv<maxnbv) {
+				vertices[nbv].r = C;
+				vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
+			}
+			else return nbv-nbvold;
+		  }
+		return nbv-nbvold;
+	}
+	/*}}}1*/
+	/*FUNCTION ListofIntersectionTriangles::NewSubSeg{{{1*/
+	void  ListofIntersectionTriangles::NewSubSeg(GeometricalEdge *e,double s0,double s1){ 
+		long int verbosity=0;
+		if (NbSeg>=MaxNbSeg) {
+			int mneo= MaxNbSeg;
+			MaxNbSeg *= 2;
+			if (verbosity>3){
+				printf("   reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
+			}
+			SegInterpolation * lEn =  new SegInterpolation[MaxNbSeg];
+			if (!lSegsI || NbSeg>=MaxNbSeg){
+				ISSMERROR("!lSegsI || NbSeg>=MaxNbSeg");
+			}
+			for (int i=0;i< NbSeg;i++) 
+			 lEn[i] = lSegsI[MaxNbSeg]; // copy old to new            
+			delete []  lSegsI; // remove old
+			lSegsI = lEn;        
+		}
+		if (NbSeg) lSegsI[NbSeg-1].last=Size;
+		lSegsI[NbSeg].e=e;
+		lSegsI[NbSeg].sBegin=s0;
+		lSegsI[NbSeg].sEnd=s1;     
+		NbSeg++;           
+	}
+	/*}}}*/
+	/*FUNCTION ListofIntersectionTriangles::ReShape{{{1*/
+	void ListofIntersectionTriangles::ReShape(){ 
+
+		register int newsize = MaxSize*2;
+		IntersectionTriangles* nw = new IntersectionTriangles[newsize];
+		ISSMASSERT(nw);
+
+		// 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);
+		MaxSize = newsize; 
+		delete [] lIntTria;// remove old
+		lIntTria = nw; // copy pointer
+	}
+	/*}}}*/
 	/*FUNCTION ListofIntersectionTriangles::SplitEdge{{{1*/
 	void ListofIntersectionTriangles::SplitEdge(const Mesh & Bh, const R2 &A,const R2  &B,int nbegin) {
@@ -33,5 +244,5 @@
 			 ifirst = ilast;}  
 		else {// not optimisation 
-			init();
+			Init();
 			t=tbegin = Bh.TriangleFindFromCoord(a,deta);
 			if( t->det>=0)
@@ -180,187 +391,38 @@
 				}
 
-						k = OppositeVertex[ocut];
-
-						Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
-
-
-						if (detbij >= 0) { //we find the triangle contening b
-							dt[0]=bamg::det((*t)[1],(*t)[2],b);
-							dt[1]=bamg::det((*t)[2],(*t)[0],b);
-							dt[2]=bamg::det((*t)[0],(*t)[1],b);
-							double dd = t->det;
-							NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);      
-							return ;}
-						else { // next triangle by  adjacent by edge ocut 
-							deti = dt[i];
-							detj = dt[j];
-							double dij = detj-deti;
-							ba[i] =  detj/dij;
-							ba[j] = -deti/dij;
-							ba[3-i-j ] = 0;
-							ilast=NewItem(t, ba[0],ba[1],ba[2]);      
-
-							TriangleAdjacent ta =t->Adj(ocut);
-							t = ta;
-							iedge= ta; 
-							if (t->det <= 0)  {
-								double ba,bb;
-								TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
-								NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
-								return;
-							}
-						}// we  go outside of omega 
+				k = OppositeVertex[ocut];
+
+				Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
+
+
+				if (detbij >= 0) { //we find the triangle contening b
+					dt[0]=bamg::det((*t)[1],(*t)[2],b);
+					dt[1]=bamg::det((*t)[2],(*t)[0],b);
+					dt[2]=bamg::det((*t)[0],(*t)[1],b);
+					double dd = t->det;
+					NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);      
+					return ;}
+				else { // next triangle by  adjacent by edge ocut 
+					deti = dt[i];
+					detj = dt[j];
+					double dij = detj-deti;
+					ba[i] =  detj/dij;
+					ba[j] = -deti/dij;
+					ba[3-i-j ] = 0;
+					ilast=NewItem(t, ba[0],ba[1],ba[2]);      
+
+					TriangleAdjacent ta =t->Adj(ocut);
+					t = ta;
+					iedge= ta; 
+					if (t->det <= 0)  {
+						double ba,bb;
+						TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
+						NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
+						return;
+					}
+				}// we  go outside of omega 
 		} // for(;;)
 	}
 	/*}}}1*/
-	/*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {{{1*/
-	int  ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) { 
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
-
-		register int n;
-		R2 x(0,0);
-		if ( d0) x =      (*tt)[0].r * d0;
-		if ( d1) x = x +  (*tt)[1].r * d1;
-		if ( d2) x = x +  (*tt)[2].r * d2;
-		// newer add same point 
-		if(!Size ||  Norme2_2(lIntTria[Size-1].x-x)) {
-			if (Size==MaxSize) ReShape();
-			lIntTria[Size].t=tt;
-			lIntTria[Size].bary[0]=d0;
-			lIntTria[Size].bary[1]=d1;
-			lIntTria[Size].bary[2]=d2;
-			lIntTria[Size].x = x;
-			Metric m0,m1,m2;
-			register BamgVertex * v;
-			if ((v=(*tt)(0))) m0    = v->m;
-			if ((v=(*tt)(1))) m1    = v->m;
-			if ((v=(*tt)(2))) m2    = v->m;
-			lIntTria[Size].m =  Metric(lIntTria[Size].bary,m0,m1,m2);
-			n=Size++;}
-		else n=Size-1;
-		return n;
-	}
-	/*}}}1*/
-	/*FUNCTION ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm){{{1*/
-	int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
-
-		register int n;
-		if(!Size ||  Norme2_2(lIntTria[Size-1].x-A)) {
-			if (Size==MaxSize) ReShape();
-			lIntTria[Size].t=0;
-			lIntTria[Size].x=A;
-			lIntTria[Size].m=mm;
-			n=Size++;
-		}
-		else  n=Size-1;
-		return  n; 
-	}
-	/*}}}1*/
-	/*FUNCTION ListofIntersectionTriangles::Length{{{1*/
-	double  ListofIntersectionTriangles::Length(){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Length)*/
-
-		// computation of the length
-
-		// check Size
-		if (Size<=0){
-			ISSMERROR("Size<=0");
-		}
-
-		Metric Mx,My;
-		int ii,jj;
-		R2 x,y,xy;
-
-		SegInterpolation* SegI=lSegsI;
-		SegI=lSegsI;
-		lSegsI[NbSeg].last=Size;
-		int EndSeg=Size;
-
-		y = lIntTria[0].x;
-		double sxy, s = 0;
-		lIntTria[0].s =0;
-		SegI->lBegin=s;
-
-		for (jj=0,ii=1;ii<Size;jj=ii++) {  
-			// seg jj,ii
-			x  = y;
-			y  = lIntTria[ii].x;
-			xy = y-x;
-			Mx = lIntTria[ii].m;
-			My = lIntTria[jj].m;
-			sxy= LengthInterpole(Mx,My,xy);
-			s += sxy;
-			lIntTria[ii].s = s;
-			if (ii == EndSeg){
-				SegI->lEnd=s;
-				SegI++;
-				EndSeg=SegI->last;
-				SegI->lBegin=s;
-			}
-		}
-		len = s;
-		SegI->lEnd=s;
-
-		return s;
-	}
-	/*}}}1*/
-	/*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
-	long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
-
-		//If length<1.5, do nothing
-		double s=Length();
-		if (s<1.5) return 0;
-
-		const long nbvold=nbv;
-		int ii = 1 ;
-		R2 y,x;
-		Metric My,Mx ;
-		double sx =0,sy;
-		int nbi=Max(2,(int) (s+0.5));
-		double sint=s/nbi;
-		double si  =sint;
-
-		int EndSeg=Size;
-		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;
-			x  =lIntTria[ii1].x;
-			sx =lIntTria[ii1].s;
-			Metric Mx=lIntTria[ii1].m;
-			y  =lIntTria[ii].x;
-			sy =lIntTria[ii].s;
-			Metric My=lIntTria[ii].m;
-			double lxy = sy-sx;
-			double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
-
-			R2 C;
-			double 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;
-			if ( nbv<maxnbv) {
-				vertices[nbv].r = C;
-				vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
-			}
-			else return nbv-nbvold;
-		  }
-		return nbv-nbvold;
-	}
-	/*}}}1*/
 
 }
Index: /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h	(revision 5130)
@@ -11,20 +11,22 @@
 
 		class IntersectionTriangles {
+
 			public: 
-				Triangle* t;
-				double  bary[3];  // use if t != 0
-				R2 x;
-				Metric m;
-				double s; // curvilinear coordinate 
-				double sp;// length of the previous segment in m
-				double sn;// length of the next segment in m
+				Triangle *t;
+				double    bary[3];   // use if t != 0
+				R2        x;
+				Metric    m;
+				double    s;         // curvilinear coordinate
+				double    sp;        // length of the previous segment in m
+				double    sn;        // length of the next segment in m
 		};
 
 		class SegInterpolation {
+
 			public:
-				GeometricalEdge * e;
-				double sBegin,sEnd; // abscisse of the seg on edge parameter
-				double lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
-				int last;// last index  in ListofIntersectionTriangles for this Sub seg of edge
+				GeometricalEdge *e;
+				double           sBegin  ,sEnd; // abscisse of the seg on edge parameter
+				double           lBegin  ,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
+				int              last;          // last index in ListofIntersectionTriangles for this Sub seg of edge
 
 				//Methods
@@ -40,24 +42,16 @@
 		public:
 
-			int MaxSize;
-			int Size;
-			double len;
-			int state;
-			IntersectionTriangles * lIntTria;
-			int NbSeg;
-			int MaxNbSeg;
-			SegInterpolation * lSegsI;
+			int                    MaxSize;
+			int                    Size;
+			double                 len;
+			int                    state;
+			IntersectionTriangles *lIntTria;
+			int                    NbSeg;
+			int                    MaxNbSeg;
+			SegInterpolation      *lSegsI;
 
 			//Constructors/Destructors
-			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);
-			}
-			~ListofIntersectionTriangles(){
-				if (lIntTria) delete [] lIntTria,lIntTria=0;
-				if (lSegsI) delete [] lSegsI,lSegsI=0;
-			} 
+			ListofIntersectionTriangles(int n=256,int m=16);
+			~ListofIntersectionTriangles();
 
 			//Operators
@@ -66,45 +60,12 @@
 
 			//Methods
-			void  init(){state=0;len=0;Size=0;}
-			int   NewItem(Triangle * tt,double d0,double d1,double d2);
-			int   NewItem(R2,const Metric & );
-			void  SplitEdge(const Mesh & ,const R2 &,const R2  &,int nbegin=0); 
-			double Length(); 
-			long  NewPoints(BamgVertex *,long & nbv,long maxnbv);
-			void  NewSubSeg(GeometricalEdge *e,double s0,double s1){ 
-				long int verbosity=0;
-				if (NbSeg>=MaxNbSeg) {
-					int mneo= MaxNbSeg;
-					MaxNbSeg *= 2;
-					if (verbosity>3){
-						printf("   reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
-					}
-					SegInterpolation * lEn =  new SegInterpolation[MaxNbSeg];
-					if (!lSegsI || NbSeg>=MaxNbSeg){
-						ISSMERROR("!lSegsI || NbSeg>=MaxNbSeg");
-					}
-					for (int i=0;i< NbSeg;i++) 
-					 lEn[i] = lSegsI[MaxNbSeg]; // copy old to new            
-					delete []  lSegsI; // remove old
-					lSegsI = lEn;        
-				}
-				if (NbSeg) lSegsI[NbSeg-1].last=Size;
-				lSegsI[NbSeg].e=e;
-				lSegsI[NbSeg].sBegin=s0;
-				lSegsI[NbSeg].sEnd=s1;     
-				NbSeg++;           
-			}
-			void ReShape(){ 
-				register int newsize = MaxSize*2;
-				IntersectionTriangles* nw = new IntersectionTriangles[newsize];
-				if (!nw) ISSMERROR("!nw");
-				// 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);
-				MaxSize = newsize; 
-				delete [] lIntTria;// remove old
-				lIntTria = nw; // copy pointer
-			}
+			void   Init();
+			int    NewItem(Triangle *tt,double d0,double d1,double d2);
+			int    NewItem(R2 ,const Metric &);
+			void   SplitEdge(const Mesh &,const R2 &,const R2 &,int nbegin=0);
+			double Length();
+			long   NewPoints(BamgVertex *,long &nbv,long maxnbv);
+			void   NewSubSeg(GeometricalEdge *e,double s0,double s1);
+			void   ReShape();
 	};
 
Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5130)
@@ -219,4 +219,14 @@
 
 	  }
+	/*}}}1*/
+	/*FUNCTION Mesh::Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices){{{1*/
+	Mesh::Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices) :Gh(BT.Gh),BTh(BT) {
+		TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);
+	}
+	/*}}}1*/
+	/*FUNCTION Mesh::Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts){{{1*/
+	Mesh::Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
+		TriangulateFromGeom0(maxnbv,bamgopts);
+	}
 	/*}}}1*/
 	/*FUNCTION Mesh::~Mesh(){{{1*/
Index: /issm/trunk/src/c/objects/Bamg/Mesh.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.h	(revision 5130)
@@ -62,19 +62,13 @@
 			Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator
 			Mesh(const Mesh &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
-			Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
-				try {TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);}
-				catch(...) { this->~Mesh(); throw; }
-			}
-			Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
-				try { TriangulateFromGeom0(maxnbv,bamgopts);}
-				catch(...) { this->~Mesh(); throw; }
-			}
+			Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1);
+			Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts);
 			~Mesh(); 
 
 			//Operators
-			const BamgVertex & operator[]  (long i) const { return vertices[i];};
-			BamgVertex & operator[](long i) { return vertices[i];};
-			const Triangle & operator()  (long i) const { return triangles[i];};
-			Triangle & operator()(long i) { return triangles[i];};
+			const BamgVertex &operator[](long i) const { return vertices[i];  };
+			BamgVertex       &operator[](long i) { return vertices[i];        };
+			const Triangle   &operator()(long i) const { return triangles[i]; };
+			Triangle         &operator()(long  i) { return triangles[i];             };
 
 			//Methods
@@ -159,5 +153,4 @@
 			void TriangulateFromGeom0(long maxnbv,BamgOpts* bamgopts);// the real constructor mesh generator
 			void Init(long);
-
 	};
 
Index: /issm/trunk/src/c/objects/Bamg/include.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/include.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/include.h	(revision 5130)
@@ -6,9 +6,6 @@
 #define  _INCLUDE2_H_
 
-
 #include "./macros.h"
 #include "./typedefs.h"
 
-
 #endif //ifndef _INCLUDE2_H_
-
Index: /issm/trunk/src/c/objects/Bamg/macros.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/macros.h	(revision 5129)
+++ /issm/trunk/src/c/objects/Bamg/macros.h	(revision 5130)
@@ -19,9 +19,9 @@
 	static const short NextVertex[3] = {1,2,0};
 	static const short PreviousVertex[3] = {2,0,1};
-#if LONG_BIT > 63
+	#if LONG_BIT > 63
 	const  Icoor1 MaxICoor   = 1073741823; // 2^30-1 =111...111 (29 times)
-#else
+	#else
 	const  Icoor1 MaxICoor   = 8388608;    // 2^23
-#endif
+	#endif
 	const  Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ;
 }