Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2861)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2862)
@@ -6,176 +6,163 @@
 namespace bamg {
 
-typedef P2<double,double> D2;
-typedef P2xP2<double,double> D2xD2;
+	typedef P2<double,double> D2;
+	typedef P2xP2<double,double> D2xD2;
 
-class  MetricAnIso;
-class MatVVP2x2;
-class MetricIso;
+	class  MetricAnIso;
+	class MatVVP2x2;
+	class MetricIso;
 
-typedef TYPEMETRIX Metric;
+	typedef TYPEMETRIX Metric;
 
 
-class MetricIso{
-  friend class MatVVP2x2;
-   Real4 h;
-public:
-  MetricIso(Real4 a): h(a){}
-  MetricIso(const MetricAnIso M);// {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
-  MetricIso(Real8 a11,Real8 a21,Real8 a22);// {*this=MetricAnIso(a11,a21,a22));}
-  MetricIso(): h(1) {}; // 
-  MetricIso(const Real8  a[3],const  MetricIso m0,
-	   const  MetricIso m1,const  MetricIso m2 )
-    : h(hinterpole 
-	? (a[0]*m0.h+a[1]*m1.h+a[2]*m2.h)
-	: 1/sqrt(a[0]/(m0.h*m0.h)+a[1]/(m1.h*m1.h)+a[2]/(m2.h*m2.h))) {}
-    MetricIso(const Real8  a,const  MetricIso ma,
-	          const Real8  b,const  MetricIso mb)
-    : h(hinterpole
-	? a*ma.h+b*mb.h
-	:1/sqrt(a/(ma.h*ma.h)+b/(mb.h*mb.h))) {}
-  R2 Orthogonal(const R2 A)const {return R2(-h*A.y,h*A.x);}
-  R2 Orthogonal(const I2 A)const {return R2(-h*A.y,h*A.x);}
-//  D2 Orthogonal(const D2 A)const {return D2(-h*A.y,h*A.x);}
-  Real8 operator()(R2 x) const { return sqrt((x,x))/h;};
-  Real8 operator()(R2 x,R2 y) const { return ((x,y))/(h*h);};
-  int  IntersectWith(MetricIso M) {int r=0;if (M.h<h) r=1,h=M.h;return r;}
-  MetricIso operator*(Real8 c) const {return  MetricIso(h/c);} 
-  MetricIso operator/(Real8 c) const {return  MetricIso(h*c);}
-  Real8 det() const {return 1./(h*h*h*h);}    
-  operator D2xD2(){ return D2xD2(1/(h*h),0.,0.,1/(h*h));}
-  void     Box(Real4 & hx,Real4 & hy) { hx=h,hy=h;}
-};
+	class MetricIso{
+		friend class MatVVP2x2;
+		Real4 h;
+		public:
+		MetricIso(Real4 a): h(a){}
+		MetricIso(const MetricAnIso M);// {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
+		MetricIso(Real8 a11,Real8 a21,Real8 a22);// {*this=MetricAnIso(a11,a21,a22));}
+		MetricIso(): h(1) {}; // 
+		MetricIso(const Real8  a[3],const  MetricIso m0,
+					const  MetricIso m1,const  MetricIso m2 )
+		  : h(hinterpole 
+					  ? (a[0]*m0.h+a[1]*m1.h+a[2]*m2.h)
+					  : 1/sqrt(a[0]/(m0.h*m0.h)+a[1]/(m1.h*m1.h)+a[2]/(m2.h*m2.h))) {}
+		MetricIso(const Real8  a,const  MetricIso ma,
+					const Real8  b,const  MetricIso mb)
+		  : h(hinterpole
+					  ? a*ma.h+b*mb.h
+					  :1/sqrt(a/(ma.h*ma.h)+b/(mb.h*mb.h))) {}
+		R2 Orthogonal(const R2 A)const {return R2(-h*A.y,h*A.x);}
+		R2 Orthogonal(const I2 A)const {return R2(-h*A.y,h*A.x);}
+		//  D2 Orthogonal(const D2 A)const {return D2(-h*A.y,h*A.x);}
+		Real8 operator()(R2 x) const { return sqrt((x,x))/h;};
+		Real8 operator()(R2 x,R2 y) const { return ((x,y))/(h*h);};
+		int  IntersectWith(MetricIso M) {int r=0;if (M.h<h) r=1,h=M.h;return r;}
+		MetricIso operator*(Real8 c) const {return  MetricIso(h/c);} 
+		MetricIso operator/(Real8 c) const {return  MetricIso(h*c);}
+		Real8 det() const {return 1./(h*h*h*h);}    
+		operator D2xD2(){ return D2xD2(1/(h*h),0.,0.,1/(h*h));}
+		void     Box(Real4 & hx,Real4 & hy) { hx=h,hy=h;}
+	};
 
 
-class MetricAnIso{ public:
-  friend class MatVVP2x2;
-  Real8 a11,a21,a22;
-  MetricAnIso(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}
-  MetricAnIso(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){}
-  MetricAnIso()  {}; // 
-  MetricAnIso(const Real8  a[3],const  MetricAnIso m0,
-	      const  MetricAnIso m1,const  MetricAnIso m2 );
-  R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);}
-  Real8 det() const {return a11*a22-a21*a21;}  
-  R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
-  R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
-//  D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
-  MetricAnIso( Real8  a,const  MetricAnIso ma,
-	       Real8  b,const  MetricAnIso mb);
-  int  IntersectWith(const MetricAnIso M);
-  MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return  MetricAnIso(a11*c2,a21*c2,a22*c2);} 
-  MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return  MetricAnIso(a11*c2,a21*c2,a22*c2);} 
-  operator D2xD2(){ return D2xD2(a11,a21,a21,a22);}
+	class MetricAnIso{ public:
+		friend class MatVVP2x2;
+		Real8 a11,a21,a22;
+		MetricAnIso(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}
+		MetricAnIso(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){}
+		MetricAnIso()  {}; // 
+		MetricAnIso(const Real8  a[3],const  MetricAnIso m0,
+					const  MetricAnIso m1,const  MetricAnIso m2 );
+		R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);}
+		Real8 det() const {return a11*a22-a21*a21;}  
+		R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
+		R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
+		//  D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
+		MetricAnIso( Real8  a,const  MetricAnIso ma,
+					Real8  b,const  MetricAnIso mb);
+		int  IntersectWith(const MetricAnIso M);
+		MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return  MetricAnIso(a11*c2,a21*c2,a22*c2);} 
+		MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return  MetricAnIso(a11*c2,a21*c2,a22*c2);} 
+		operator D2xD2(){ return D2xD2(a11,a21,a21,a22);}
 
-  Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
-//  Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
-  Real8 operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};
-  inline void  Box(Real4 &hx,Real4 &hy) const ;  
-  MetricAnIso(const MatVVP2x2);
-};
+		Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
+		//  Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
+		Real8 operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};
+		inline void  Box(Real4 &hx,Real4 &hy) const ;  
+		MetricAnIso(const MatVVP2x2);
+	};
 
+	class MatVVP2x2{
+		friend  class MetricAnIso;
+		friend  class MetricIso;
+		public:
+		double lambda1,lambda2;
+		D2 v;
 
-class MatVVP2x2 
-{
-  friend  class MetricAnIso;
-  friend  class MetricIso;
-public:
-  double lambda1,lambda2;
-  D2 v;
+		MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
 
+		void  Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);}
+		void  pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);}
+		void  Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;}
+		void  Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;}
 
-  MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
-  
-  void  Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);}
-  void  pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);}
-  void  Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;}
-  void  Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;}
+		void Minh(double h) {Max(1.0/(h*h));}
+		void Maxh(double h) {Min(1.0/(h*h));}
+		void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
+		MatVVP2x2(const MetricAnIso );
+		MatVVP2x2(const MetricIso M) :  lambda1(1/(M.h*M.h)),lambda2(1/(M.h*M.h)),v(1,0) {}
+		Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));}
+		Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));}
+		Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);}
+		Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);}
+		Real8 Aniso2() const  { return lmax()/lmin();}
+		inline void BoundAniso2(const Real8 coef);
+		Real8 Aniso() const  { return sqrt( Aniso2());}
+		void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
+		void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
+	  };
 
-  void Minh(double h) {Max(1.0/(h*h));}
-  void Maxh(double h) {Min(1.0/(h*h));}
-  void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
-  MatVVP2x2(const MetricAnIso );
-  MatVVP2x2(const MetricIso M) :  lambda1(1/(M.h*M.h)),lambda2(1/(M.h*M.h)),v(1,0) {}
-  Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));}
-  Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));}
-  Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);}
-  Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);}
-  Real8 Aniso2() const  { return lmax()/lmin();}
-  inline void BoundAniso2(const Real8 coef);
-  Real8 Aniso() const  { return sqrt( Aniso2());}
-  void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
-  void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
-};
+	inline void  MatVVP2x2::BoundAniso2(const Real8 coef){
+		if (coef<=1.00000000001) 
+		 if (lambda1 < lambda2)
+		  lambda1 = bamg::Max(lambda1,lambda2*coef);
+		 else
+		  lambda2 = bamg::Max(lambda2,lambda1*coef);
+		else  // a verifier 
+		 if (lambda1 > lambda2)
+		  lambda1 = bamg::Min(lambda1,lambda2*coef);
+		 else
+		  lambda2 = bamg::Min(lambda2,lambda1*coef);
+	  }
 
-inline void  MatVVP2x2::BoundAniso2(const Real8 coef) 
-{
-  if (coef<=1.00000000001) 
-    if (lambda1 < lambda2)
-      lambda1 = bamg::Max(lambda1,lambda2*coef);
-    else
-      lambda2 = bamg::Max(lambda2,lambda1*coef);
-  else  // a verifier 
-    if (lambda1 > lambda2)
-      lambda1 = bamg::Min(lambda1,lambda2*coef);
-    else
-      lambda2 = bamg::Min(lambda2,lambda1*coef);
-}
+	void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
+	MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
+
+	inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) {
+		//     recompose M in: M = V^t lambda V 
+		//     V = ( v,v^\perp)
+		//  where v^\perp = (-v_1,v_0)
+		double v00=M.v.x*M.v.x;
+		double v11=M.v.y*M.v.y;
+		double v01=M.v.x*M.v.y;
+		a11=v00*M.lambda1+v11*M.lambda2;
+		a21=v01*(M.lambda1-M.lambda2);
+		a22=v00*M.lambda2+v11*M.lambda1;
+	  }
+
+	inline   void  MetricAnIso::Box(Real4 &hx,Real4 &hy) const {
+		Real8 d=  a11*a22-a21*a21;
+		hx = sqrt(a22/d);
+		hy = sqrt(a11/d);
+	}
+
+	inline MetricIso::MetricIso(const MetricAnIso M) 
+	  {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
+
+	inline  MetricIso::MetricIso(Real8 a11,Real8 a21,Real8 a22)
+	  {MatVVP2x2 vp(MetricAnIso(a11,a21,a22));h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
 
 
 
+	class SaveMetricInterpole {
+		friend  Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
+		friend Real8 abscisseInterpole(const MetricAnIso ,const  MetricAnIso , R2 ,Real8 ,int );
+		int opt;
+		Real8 lab;
+		Real8 L[1024],S[1024];
+	};
 
-void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
-MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
+	extern SaveMetricInterpole  LastMetricInterpole; // for optimization 
 
-inline MetricAnIso::MetricAnIso(const MatVVP2x2 M)  
-{
- //     recompose M in: M = V^t lambda V 
-  //     V = ( v,v^\perp)
-  //  where v^\perp = (-v_1,v_0)
-  double v00=M.v.x*M.v.x;
-  double v11=M.v.y*M.v.y;
-  double v01=M.v.x*M.v.y;
-  a11=v00*M.lambda1+v11*M.lambda2;
-  a21=v01*(M.lambda1-M.lambda2);
-  a22=v00*M.lambda2+v11*M.lambda1;
-}
+	Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB);
+	Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim=0);
 
+	inline Real8 LengthInterpole(Real8 la,Real8 lb) 
+	  {   return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ?  (la+lb)/2  : la*lb*log(la/lb)/(la-lb);}
 
-inline   void  MetricAnIso::Box(Real4 &hx,Real4 &hy) const {
-  Real8 d=  a11*a22-a21*a21;
-  hx = sqrt(a22/d);
-  hy = sqrt(a11/d);
-}
-
-
-inline MetricIso::MetricIso(const MetricAnIso M) 
-  {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
-  
-inline  MetricIso::MetricIso(Real8 a11,Real8 a21,Real8 a22)
-  {MatVVP2x2 vp(MetricAnIso(a11,a21,a22));h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
-
-
-
-class SaveMetricInterpole {
-  friend  Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
-  friend Real8 abscisseInterpole(const MetricAnIso ,const  MetricAnIso , R2 ,Real8 ,int );
-  int opt;
-  Real8 lab;
-  Real8 L[1024],S[1024];
-};
-
-extern SaveMetricInterpole  LastMetricInterpole; // for optimization 
-
-
- Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB);
- Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim=0);
-
-
-
-inline Real8 LengthInterpole(Real8 la,Real8 lb) 
-{   return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ?  (la+lb)/2  : la*lb*log(la/lb)/(la-lb);}
-
-inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s)
-{ return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20))  ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb);}
+	inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s)
+	  { return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20))  ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb);}
 
 }
Index: /issm/trunk/src/c/Bamgx/QuadTree.h
===================================================================
--- /issm/trunk/src/c/Bamgx/QuadTree.h	(revision 2861)
+++ /issm/trunk/src/c/Bamgx/QuadTree.h	(revision 2862)
@@ -5,74 +5,73 @@
 namespace bamg {
 
-const int MaxDeep = 30;
-typedef  long  IntQuad;
-const IntQuad MaxISize = ( 1L << MaxDeep);
+	const int MaxDeep = 30;
+	typedef  long  IntQuad;
+	const IntQuad MaxISize = ( 1L << MaxDeep);
 
 
-class Triangles;
-class Vertex;
+	class Triangles;
+	class Vertex;
 
-class QuadTree {
- public:
+	class QuadTree {
+		public:
 
-  class QuadTreeBox { 
-  public:
+			class QuadTreeBox { 
+				public:
 
-    long n; // if n < 4 => Vertex else =>  QuadTreeBox;
-    union {
-      QuadTreeBox *b[4];
-      Vertex * v[4];
-    };
-    
-
-  }; // end class QuadTreeBox  /////////////////
-
-  class StorageQuadTreeBox {
-  public:
-    QuadTreeBox *b,*bc,*be;
-    long len;
-    StorageQuadTreeBox *n; // next StorageQuadTreeBox
-    StorageQuadTreeBox(long ,StorageQuadTreeBox * =0);
-    ~StorageQuadTreeBox() {
-      if(n) delete n;
-      delete [] b;
-    }
-    long  SizeOf() const {
-      return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);
-    }
-  }; // end class  StorageQuadTreeBox 
-  
-  StorageQuadTreeBox * sb;
-  
-  
-  long  lenStorageQuadTreeBox;
-
-public:
-  QuadTreeBox * root;
-  Triangles *th;
-  long NbQuadTreeBox,NbVertices;
-  long NbQuadTreeBoxSearch,NbVerticesSearch;
-  Vertex * NearestVertex(Icoor1 i,Icoor1 j);
-  Vertex *  NearestVertexWithNormal(Icoor1 i,Icoor1 j); // new version  
-  Vertex * ToClose(Vertex & ,Real8 ,Icoor1,Icoor1);
-  long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();}
+					long n; // if n < 4 => Vertex else =>  QuadTreeBox;
+					union {
+						QuadTreeBox *b[4];
+						Vertex * v[4];
+					};
 
 
-  void  Add( Vertex & w);
+			}; // end class QuadTreeBox  /////////////////
 
-  QuadTreeBox* NewQuadTreeBox(){
-    if(! (sb->bc<sb->be)) 
-	  sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
+			class StorageQuadTreeBox {
+				public:
+					QuadTreeBox *b,*bc,*be;
+					long len;
+					StorageQuadTreeBox *n; // next StorageQuadTreeBox
+					StorageQuadTreeBox(long ,StorageQuadTreeBox * =0);
+					~StorageQuadTreeBox() {
+						if(n) delete n;
+						delete [] b;
+					}
+					long  SizeOf() const {
+						return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);
+					}
+			}; // end class  StorageQuadTreeBox 
 
-	 if (!sb || (sb->bc->n != 0)){
-		 throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));
-	 }
-    NbQuadTreeBox++;
-    return sb->bc++;
-  }
-  ~QuadTree();
-  QuadTree(Triangles * t,long nbv=-1);
-  QuadTree();
-	
-};
+			StorageQuadTreeBox * sb;
+
+
+			long  lenStorageQuadTreeBox;
+
+		public:
+			QuadTreeBox * root;
+			Triangles *th;
+			long NbQuadTreeBox,NbVertices;
+			long NbQuadTreeBoxSearch,NbVerticesSearch;
+			Vertex * NearestVertex(Icoor1 i,Icoor1 j);
+			Vertex *  NearestVertexWithNormal(Icoor1 i,Icoor1 j); // new version  
+			Vertex * ToClose(Vertex & ,Real8 ,Icoor1,Icoor1);
+			long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();}
+
+
+			void  Add( Vertex & w);
+
+			QuadTreeBox* NewQuadTreeBox(){
+				if(! (sb->bc<sb->be)) 
+				 sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
+
+				if (!sb || (sb->bc->n != 0)){
+					throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));
+				}
+				NbQuadTreeBox++;
+				return sb->bc++;
+			}
+			~QuadTree();
+			QuadTree(Triangles * t,long nbv=-1);
+			QuadTree();
+	};
 }
Index: /issm/trunk/src/c/Bamgx/R2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/R2.h	(revision 2861)
+++ /issm/trunk/src/c/Bamgx/R2.h	(revision 2862)
@@ -2,92 +2,89 @@
 
 namespace bamg {
-template <class R,class RR> class P2xP2;
+	template <class R,class RR> class P2xP2;
 
-template <class R,class RR>
-class P2 {
+	template <class R,class RR>
+	  class P2 {
 
-public:  
-  R x,y;
-  P2 () :x(0),y(0) {};
-  P2 (R a,R b)  :x(a),y(b)  {}
-  P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {}
-  P2<R,RR>   operator+(const P2<R,RR> & cc) const {return P2<R,RR>(x+cc.x,y+cc.y);}
-  P2<R,RR>   operator-(const P2<R,RR> & cc) const {return P2<R,RR>(x-cc.x,y-cc.y);}
-  P2<R,RR>   operator-()  const{return P2<R,RR>(-x,-y);}
-//  RR   operator*(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} // produit scalaire
-  RR   operator,(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} // produit scalaire
-  P2<R,RR>   operator*(R  cc) const {return P2<R,RR>(x*cc,y*cc);}
- // P2<R,RR>   operator*(RR  cc) const {return P2<R,RR>((R)(x*cc),(R)(y*cc));}
-  P2<R,RR>   operator/(R  cc) const {return P2<R,RR>(x/cc,y/cc);}
-  P2<R,RR>  operator+=(const  P2<R,RR> & cc) {x += cc.x;y += cc.y;return *this;}
-  P2<R,RR>  operator/=(const  R r) {x /= r;y /= r;return *this;}
-  P2<R,RR>  operator*=(const  R r) {x *= r;y *= r;return *this;}
-  P2<R,RR>  operator-=(const  P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;}
-//  P2<R,RR> Orthogonal(const   P2<R,RR> r) {return P2<R,RR>(-r.y,r.x);}
- };
+		  public:  
+			  R x,y;
+			  P2 () :x(0),y(0) {};
+			  P2 (R a,R b)  :x(a),y(b)  {}
+			  P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {}
+			  P2<R,RR>   operator+(const P2<R,RR> & cc) const {return P2<R,RR>(x+cc.x,y+cc.y);}
+			  P2<R,RR>   operator-(const P2<R,RR> & cc) const {return P2<R,RR>(x-cc.x,y-cc.y);}
+			  P2<R,RR>   operator-()  const{return P2<R,RR>(-x,-y);}
+			  RR   operator,(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} // produit scalaire
+			  P2<R,RR>   operator*(R  cc) const {return P2<R,RR>(x*cc,y*cc);}
+			  P2<R,RR>   operator/(R  cc) const {return P2<R,RR>(x/cc,y/cc);}
+			  P2<R,RR>  operator+=(const  P2<R,RR> & cc) {x += cc.x;y += cc.y;return *this;}
+			  P2<R,RR>  operator/=(const  R r) {x /= r;y /= r;return *this;}
+			  P2<R,RR>  operator*=(const  R r) {x *= r;y *= r;return *this;}
+			  P2<R,RR>  operator-=(const  P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;}
+	  };
 
-template <class R,class RR>
-class P2xP2 { // x ligne 1 y ligne2 
+	template <class R,class RR>
+	  class P2xP2 { // x ligne 1 y ligne2 
 
-  friend std::ostream& operator <<(std::ostream& f, const P2xP2<R,RR> & c) 
-     { f << '[' << c.x << ',' << c.y << ']' <<std::flush ; return f; }
-     
-  friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc)
-     {return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y);} 
+		  friend std::ostream& operator <<(std::ostream& f, const P2xP2<R,RR> & c) 
+			 { f << '[' << c.x << ',' << c.y << ']' <<std::flush ; return f; }
+
+		  friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc)
+			 {return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y);} 
 
 
- public:
-  P2<R,RR> x,y; 
-  P2xP2 (): x(),y()  {}
-  P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {}
-  P2xP2 (P2<R,RR> a,P2<R,RR> b,P2<R,RR> c ): x(b-a),y(c-a) {}
-  P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {}
-  P2<R,RR> operator*(const P2<R,RR> c) const {return P2<R,RR>(x.x*c.x + x.y*c.y, y.x*c.x + y.y*c.y);}
-  P2xP2<R,RR>  operator*(P2xP2<R,RR> c) const 
-    { return  P2xP2<R,RR>(x.x*c.x.x + x.y*c.y.x,
-			  x.x*c.x.y + x.y*c.y.y,
-			  y.x*c.x.x + y.y*c.y.x,
-			  y.x*c.x.y + y.y*c.y.y);}
-  RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;}
-  P2xP2<R,RR> inv()  const
-     { RR d = (*this).det(); 
-       return P2xP2<R,RR>((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) );
-     };
-   P2xP2<R,RR> t() {return P2xP2<R,RR>(x.x,y.x,x.y,y.y);} //transposer 
-   P2<R,RR>tx() {return P2<R,RR>(x.x,y.x);} 
-   P2<R,RR>ty() {return P2<R,RR>(x.y,y.y);} 
+		  public:
+		  P2<R,RR> x,y; 
+		  P2xP2 (): x(),y()  {}
+		  P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {}
+		  P2xP2 (P2<R,RR> a,P2<R,RR> b,P2<R,RR> c ): x(b-a),y(c-a) {}
+		  P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {}
+		  P2<R,RR> operator*(const P2<R,RR> c) const {return P2<R,RR>(x.x*c.x + x.y*c.y, y.x*c.x + y.y*c.y);}
+		  P2xP2<R,RR>  operator*(P2xP2<R,RR> c) const 
+			 { return  P2xP2<R,RR>(x.x*c.x.x + x.y*c.y.x,
+						 x.x*c.x.y + x.y*c.y.y,
+						 y.x*c.x.x + y.y*c.y.x,
+						 y.x*c.x.y + y.y*c.y.y);}
+		  RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;}
+		  P2xP2<R,RR> inv()  const
+			 { RR d = (*this).det(); 
+			  return P2xP2<R,RR>((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) );
+			 };
+		  P2xP2<R,RR> t() {return P2xP2<R,RR>(x.x,y.x,x.y,y.y);} //transposer 
+		  P2<R,RR>tx() {return P2<R,RR>(x.x,y.x);} 
+		  P2<R,RR>ty() {return P2<R,RR>(x.y,y.y);} 
 
-};  
+	  };  
 
-template  <class R,class RR>  
-inline RR Det(const P2<R,RR> x,const P2<R,RR> y) {
-  return (RR) x.x * (RR) y.y - (RR) x.y * (RR) y.x ;} 
+	template  <class R,class RR>  
+	  inline RR Det(const P2<R,RR> x,const P2<R,RR> y) {
+		  return (RR) x.x * (RR) y.y - (RR) x.y * (RR) y.x ;} 
 
-template  <class R,class RR>  
-inline RR Area2 (const P2<R,RR> a,const P2<R,RR> b,const P2<R,RR> c) {
-  return Det(b-a,c-a) ;} 
+	template  <class R,class RR>  
+	  inline RR Area2 (const P2<R,RR> a,const P2<R,RR> b,const P2<R,RR> c) {
+		  return Det(b-a,c-a) ;} 
 
-template  <class R,class RR>  
-inline R Norme1 (const P2<R,RR> x) {
-  return (Abs(x.x)+Abs(x.y)) ;} 
+	template  <class R,class RR>  
+	  inline R Norme1 (const P2<R,RR> x) {
+		  return (Abs(x.x)+Abs(x.y)) ;} 
 
-template  <class R,class RR>  
-inline R NormeInfini (const P2<R,RR> x) {
-  return Max(Abs(x.x),Abs(x.y)) ;} 
+	template  <class R,class RR>  
+	  inline R NormeInfini (const P2<R,RR> x) {
+		  return Max(Abs(x.x),Abs(x.y)) ;} 
 
-template  <class R,class RR>  
-inline RR Norme2_2 (const P2<R,RR> x) {
-  return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ;} 
+	template  <class R,class RR>  
+	  inline RR Norme2_2 (const P2<R,RR> x) {
+		  return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ;} 
 
-template  <class R,class RR>  
-inline RR Norme2 (const P2<R,RR> x) {
-  return sqrt((RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y) ;} 
+	template  <class R,class RR>  
+	  inline RR Norme2 (const P2<R,RR> x) {
+		  return sqrt((RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y) ;} 
 
-template  <class R,class RR>  
-inline P2<R,RR> Orthogonal (const P2<R,RR> x) {
-  return  P2<R,RR>(-x.y,x.x);} 
+	template  <class R,class RR>  
+	  inline P2<R,RR> Orthogonal (const P2<R,RR> x) {
+		  return  P2<R,RR>(-x.y,x.x);} 
 
-template <class R,class RR>
-inline  std::ostream& operator <<(std::ostream& f, const P2<R,RR> & c)
-  { f << '[' << c.x << ',' << c.y <<']' <<std::flush ; return f; }
+	template <class R,class RR>
+	  inline  std::ostream& operator <<(std::ostream& f, const P2<R,RR> & c)
+		 { f << '[' << c.x << ',' << c.y <<']' <<std::flush ; return f; }
 }
Index: /issm/trunk/src/c/Bamgx/SetOfE4.h
===================================================================
--- /issm/trunk/src/c/Bamgx/SetOfE4.h	(revision 2861)
+++ /issm/trunk/src/c/Bamgx/SetOfE4.h	(revision 2862)
@@ -4,31 +4,30 @@
 namespace bamg {
 
-class SetOfEdges4 ;
-class Int4Edge{
-friend class SetOfEdges4;
-public:
-  Int4 i,j;
-  Int4 next; 
-};
+	class SetOfEdges4 ;
+	class Int4Edge{
+		friend class SetOfEdges4;
+		public:
+		Int4 i,j;
+		Int4 next; 
+	};
 
-class SetOfEdges4 {
-  Int4 nx,nbax,NbOfEdges;
-  Int4 * tete; 
-  Int4Edge * Edges;
+	class SetOfEdges4 {
+		Int4 nx,nbax,NbOfEdges;
+		Int4 * tete; 
+		Int4Edge * Edges;
 
-public:
-  SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 
-  ~SetOfEdges4() {
-  delete [] tete; delete [] Edges;}
-   Int4 add (Int4 ii,Int4 jj);
-  Int4 addtrie (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
-  Int4  nb(){return NbOfEdges;}
-  Int4 find (Int4 ii,Int4 jj);
-  Int4 findtrie (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj)  : find (jj,ii) ;}
-  Int4 i(Int4 k){return Edges[k].i;}
-  Int4 j(Int4 k){return Edges[k].j;}
-  Int4 newarete(Int4 k){return NbOfEdges == k+1;}
-  Int4Edge & operator[](Int4 k){return  Edges[k];}
-};
+		public:
+		SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 
+		~SetOfEdges4() {delete [] tete; delete [] Edges;}
+		Int4 add (Int4 ii,Int4 jj);
+		Int4 addtrie (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
+		Int4  nb(){return NbOfEdges;}
+		Int4 find (Int4 ii,Int4 jj);
+		Int4 findtrie (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj)  : find (jj,ii) ;}
+		Int4 i(Int4 k){return Edges[k].i;}
+		Int4 j(Int4 k){return Edges[k].j;}
+		Int4 newarete(Int4 k){return NbOfEdges == k+1;}
+		Int4Edge & operator[](Int4 k){return  Edges[k];}
+	};
 }