Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3232)
@@ -18,5 +18,5 @@
 #include <string.h>
 #include "Mesh2.h"
-#include "QuadTree.h"
+#include "objects/QuadTree.h"
 
 using namespace bamg;
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 3232)
@@ -8,816 +8,26 @@
 
 #include "meshtype.h"
-#include "Metric.h"
+
+#include "objects/Metric.h"
+#include "objects/DoubleAndInt4.h"
+#include "objects/Direction.h"
+#include "objects/Vertex.h"
+#include "objects/TriangleAdjacent.h"
+#include "objects/Edge.h"
+#include "objects/GeometricalVertex.h"
+#include "objects/GeometricalEdge.h"
+#include "objects/Curve.h"
+#include "objects/Triangle.h"
+#include "objects/ListofIntersectionTriangles.h"
+#include "objects/GeometricalSubDomain.h"
+#include "objects/SubDomain.h"
+#include "objects/VertexOnGeom.h"
+#include "objects/VertexOnVertex.h"
+#include "objects/VertexOnEdge.h"
+#include "objects/CrackedEdge.h"
+#include "objects/Triangles.h"
+#include "objects/Geometry.h"
 
 namespace bamg {
-
-	//classes
-	class Geometry;
-	class Triangles;
-	class Triangle;
-	class QuadTree;
-	class GeometricalEdge;
-	class VertexOnGeom;
-	class VertexOnEdge;
-
-	/*CLASS: DoubleAndInt4 {{{1*/
-	class DoubleAndInt4 {
-		//class used by Triangles::MakeQuadrangles
-
-		public:
-			double q;
-			Int4 i3j;
-
-			//Operators
-			int operator<(DoubleAndInt4 a){return q > a.q;}
-	};
-	/*}}}1*/
-	/*CLASS: Direction{{{1*/
-	class Direction { //   
-		private:
-			Icoor1 dir;
-
-		public:
-			//Methods
-			Direction(): dir(MaxICoor){}; //  no direction set
-			Direction(Icoor1 i,Icoor1 j) {
-				Icoor2 n2 = 2*(Abs(i)+Abs(j));  
-				Icoor2 r  = MaxICoor* (Icoor2) i;
-				Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number 
-				dir = (j>0) ? r1 : r1+1;          // odd-> j>0 even-> j<0
-			}
-			int sens(Icoor1 i,Icoor1 j) {
-				int r =1; 
-				if (dir!= MaxICoor) {
-					Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1);
-					r = (x*i + y*j) >=0;
-				}
-				return r;
-			}
-	};
-	/*}}}1*/
-	/*CLASS: Vertex{{{1*/
-	class Vertex {
-
-		public:
-			I2 i;  // integer coordinates
-			R2 r;  // real coordinates
-			Metric m;
-			Int4 ReferenceNumber;
-			Direction DirOfSearch;
-			Int1 vint;  // the vertex number in triangle; varies between 0 and 2 in t
-			union {
-				Triangle* t; // one triangle which is containing the vertex
-				Int4 color;
-				Vertex* to;  // use in geometry Vertex to now the Mesh Vertex associed 
-				VertexOnGeom* onGeometry;       // if vint == 8; // set with Triangles::SetVertexFieldOn()
-				Vertex* onBackgroundVertex;     // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
-				VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge
-			};
-
-			//Operators
-			operator  I2() const {return i;}             // Cast operator
-			operator  const R2 & () const {return r;}    // Cast operator
-			operator Metric () const {return m;}         // Cast operator
-			Real8 operator()(R2 x) const { return m(x);} // Get x in the metric m
-
-			//methods (No constructor and no destructors...)
-			Int4  Optim(int =1,int =0); 
-			Real8 Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,Real8 =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);
-			void  Echo();
-			int   ref() const { return ReferenceNumber;}
-
-			//inline functions
-			inline void Set(const Vertex & rec,const Triangles &,Triangles &);
-	};
-	/*}}}1*/
-	inline Vertex* TheVertex(Vertex * a); // for remove crak in mesh 
-	double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &);
-	/*CLASS: TriangleAdjacent{{{1*/
-	class TriangleAdjacent {
-
-		public:
-			Triangle* t; //pointer toward the triangle
-			int  a;      //Edge number
-
-			//Constructors
-			TriangleAdjacent() {};
-			TriangleAdjacent(Triangle  * tt,int  aa): t(tt),a(aa &3) {};
-
-			//Operators
-			operator Triangle * () const {return t;}
-			operator Triangle & () const {return *t;}
-			operator int() const {return a;}
-			TriangleAdjacent & operator++(){ a= NextEdge[a]; return *this; }
-			TriangleAdjacent operator--(){ a= PreviousEdge[a]; return *this; }
-
-			//Methods
-			int swap();
-
-			//Inline methods
-			inline  TriangleAdjacent  Adj() const ;
-			inline void SetAdj2(const TriangleAdjacent& , int =0);
-			inline Vertex *  EdgeVertex(const int &) const ;
-			inline Vertex *  OppositeVertex() const ;
-			inline Icoor2 & det() const;
-			inline int Locked() const  ;
-			inline int GetAllFlag_UnSwap() const ;
-			inline void SetLock();
-			inline int MarkUnSwap()  const;
-			inline void SetMarkUnSwap();
-			inline void SetCracked();
-			inline int Cracked() const ;
-	};
-	/*}}}1*/
-	/*CLASS: Edge{{{1*/
-	class Edge {
-
-		public:
-			Vertex* v[2];
-			Int4 ref;
-			GeometricalEdge* onGeometry;
-			Edge* adj[2]; // the 2 adj edges if on the same curve 
-
-			//Operators
-			Vertex & operator[](int i){return *v[i];};
-			Vertex * operator()(int i){return v[i];};
-			R2       operator()(double t) const; // return the point 
-			const Vertex & operator[](int i) const { return *v[i];};
-
-			//Methods
-			void ReNumbering(Vertex *vb,Vertex *ve, Int4 *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)){
-					throw ErrorException(__FUNCT__,exprintf("Intersection bug"));
-				}
-				if (adj[0]!=&e && adj[1]!=&e){
-					throw ErrorException(__FUNCT__,exprintf("adj[0]!=&e && adj[1]!=&e"));
-				}
-				return adj[0]==&e ? 0 : 1;
-			}
-
-			//Inline methods
-			inline void Set(const Triangles &,Int4,Triangles &);
-
-	};
-	/*}}}1*/
-	inline TriangleAdjacent FindTriangleAdjacent(Edge &E);
-	/*CLASS: GeometricalVertex{{{1*/
-	class GeometricalVertex : public Vertex { 
-
-		public:
-			friend class Geometry;
-
-			int cas;
-			GeometricalVertex* link; //  link all the same GeometricalVertex circular (Crack) 
-
-			//Methods
-			int  Corner() const {return cas&4;}
-			int  Required()const {return cas&6;}// a corner is required
-			int  IsThe() const { return link == this;}  
-			void SetCorner(){ cas |= 4;}
-			void SetRequired(){ cas |= 2;}
-			void Set(){cas=0;}
-			GeometricalVertex() :cas(0), link(this) {};
-			GeometricalVertex* The() {
-				if (!link){ throw ErrorException(__FUNCT__,exprintf("!link"));}
-				return link;
-			}// return a unique vertex
-
-			//Inline methods
-			inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew);
-	};
-	/*}}}1*/
-	/*CLASS: GeometricalEdge{{{1*/
-	class GeometricalEdge {
-
-		public:
-			GeometricalVertex* v[2];
-			Int4 ref;
-			Int4 CurveNumber;
-			R2   tg[2]; // the 2 tangentes (tg[0] =0 => no continuity)
-			GeometricalEdge* Adj[2]; 
-			int DirAdj[2];
-			int flag ;
-			GeometricalEdge* link; // if   Cracked() or Equi()
-
-			//Operators
-			GeometricalVertex       & operator[](int i){return *v[i];};
-			const GeometricalVertex & operator[](int i) const { return *v[i];};
-			GeometricalVertex       * operator()(int i){return v[i];};  
-
-			//Methods
-			R2 F(Real8 theta) const ; // parametrization of the curve edge
-			Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente
-			int  Cracked() const {return flag & 1;}
-			int  Dup() const { return flag & 32;}
-			int  Equi()const {return flag & 2;}
-			int  ReverseEqui()const {return flag & 128;}
-			int  TgA()const {return flag &4;}
-			int  TgB()const {return flag &8;}
-			int  Tg(int i) const{return i==0 ? TgA() : TgB();}
-			int  Mark()const {return flag &16;}
-			int  Required() { return flag & 64;}
-			void SetCracked() { flag |= 1;}
-			void SetDup()     { flag |= 32;} // not a real edge 
-			void SetEqui()    { flag |= 2;}
-			void SetTgA()     { flag|=4;}
-			void SetTgB()     { flag|=8;}
-			void SetMark()    { flag|=16;}
-			void SetUnMark()  { flag &= 1007 /* 1023-16*/;}
-			void SetRequired() { flag |= 64;}
-			void SetReverseEqui() {flag |= 128;}
-
-			//Inline methods
-			inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
-	};
-	/*}}}1*/
-	/*CLASS: Curve{{{1*/
-	class Curve {
-		public:
-			GeometricalEdge * be,*ee; // begin et end edge
-			int kb,ke;  //  begin vetex and end vertex
-			Curve *next; // next curve equi to this
-			bool master; // true => of equi curve point on this curve  
-
-			//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 
-
-			//Inline methods
-			inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
-	};
-	/*}}}1*/
-	/*CLASS: Triangle{{{1*/
-	class Triangle {
-
-		friend class TriangleAdjacent;
-
-		private:
-			Vertex*   TriaVertices[3];      // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
-			Triangle* TriaAdjTriangles[3];  // 3 pointers toward the adjacent triangles
-			Int1      TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1
-
-		public: 
-			Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
-			union { 
-				Triangle * link ;
-				Int4 color;
-			};
-
-			//Constructors/Destructors
-			Triangle() {}
-			Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k);
-			Triangle(Vertex *v0,Vertex *v1,Vertex *v2);
-
-			//Operators
-			const Vertex & operator[](int i) const {return *TriaVertices[i];};
-			Vertex & operator[](int i)  {return *TriaVertices[i];};
-			const Vertex * operator()(int i) const {return TriaVertices[i];};
-			Vertex * & operator()(int i)  {return TriaVertices[i];};
-
-			//Methods
-			void   Echo();
-			int    swap(Int2 a1,int=0);
-			Int4   Optim(Int2 a,int =0);
-			int    Locked(int a)const { return TriaAdjSharedEdge[a]&4;} 
-			int    Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} 
-			int    Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;}
-			int    GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;}
-			void   SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);}
-			double QualityQuad(int a,int option=1) const;
-			Int1   NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the  adjacent edge in adj tria  
-			TriangleAdjacent FindBoundaryEdge(int i) const;
-			TriangleAdjacent Adj(int i)  const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
-			Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
-			Triangle* Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
-			void  ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
-				if (link  >=tb && link  <te) link  = tb + renu[link -tb];
-				if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
-				if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb];
-				if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];    
-			}
-			void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
-				if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
-				if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
-				if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb];    
-			}
-			void SetAdjAdj(Int1 a){
-				a &= 3;
-				register Triangle *tt=TriaAdjTriangles[a];
-				TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap
-				register Int1 aatt = TriaAdjSharedEdge[a] & 3;
-				if(tt){ 
-					tt->TriaAdjTriangles[aatt]=this;
-					tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark 
-			  }
-			void SetAdj2(Int1 a,Triangle *t,Int1 aat){
-				TriaAdjTriangles[a]=t;   //the adjacent triangle to the edge a is t
-				TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle
-				if(t) { //if t!=NULL add adjacent triangle to t (this)
-					t->TriaAdjTriangles[aat]=this;
-					t->TriaAdjSharedEdge[aat]=a;
-				}
-			}
-			void SetTriangleContainingTheVertex() { 
-				if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
-				if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);
-				if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2);
-			}
-			void SetHidden(int a){
-				//Get Adjacent Triangle number a
-				register Triangle* t = TriaAdjTriangles[a];
-				//if it exist
-				//C|=D -> C=(C|D) bitwise inclusive OR
-				if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16;
-				TriaAdjSharedEdge[a] |= 16;
-			}
-			void SetCracked(int a){
-				register Triangle * t = TriaAdjTriangles[a];
-				if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32;
-				TriaAdjSharedEdge[a] |= 32;
-			}
-
-			void SetLocked(int a){
-				//mark the edge as on Boundary
-				register Triangle * t = TriaAdjTriangles[a];
-				t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4;
-				TriaAdjSharedEdge[a] |= 4;
-			}
-			void SetMarkUnSwap(int a){
-				register Triangle * t = TriaAdjTriangles[a];
-				t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8;
-				TriaAdjSharedEdge[a] |=8 ;
-			}
-			void SetUnMarkUnSwap(int a){ 
-				register Triangle * t = TriaAdjTriangles[a];
-				t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 
-				TriaAdjSharedEdge[a] &=55 ;
-			}
-			void SetDet() {
-				if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2])    det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
-				else det = -1; }
-
-			//Inline methods
-			inline Real4 qualite() ;
-			inline void  Set(const Triangle &,const Triangles &,Triangles &);
-			inline int   In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
-
-
-	};  // end of Triangle class 
-	/*}}}1*/
-	/*CLASS: ListofIntersectionTriangles{{{1*/
-	class ListofIntersectionTriangles {
-
-		class IntersectionTriangles {
-			public: 
-				Triangle* t;
-				Real8  bary[3];  // use if t != 0
-				R2 x;
-				Metric m;
-				Real8 s; // curvilinear coordinate 
-				Real8 sp;// length of the previous segment in m
-				Real8 sn;// length of the next segment in m
-		};
-
-		class SegInterpolation {
-			public:
-				GeometricalEdge * e;
-				Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
-				Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
-				int last;// last index  in ListofIntersectionTriangles for this Sub seg of edge
-
-				//Methods
-				R2 F(Real8 s){ 
-					Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01;
-					if (lBegin>s || s>lEnd){
-						throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd"));
-					}
-					return e->F(sBegin*c0+sEnd*c1);
-				}
-		};
-
-		public:
-
-			int MaxSize;
-			int Size;
-			Real8 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;
-			} 
-
-			//Operators
-			IntersectionTriangles & operator[](int i) {return lIntTria[i];}
-			operator int&() {return Size;}
-
-			//Methods
-			void  init(){state=0;len=0;Size=0;}
-			int   NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
-			int   NewItem(R2,const Metric & );
-			void  SplitEdge(const Triangles & ,const R2 &,const R2  &,int nbegin=0); 
-			Real8 Length(); 
-			Int4  NewPoints(Vertex *,Int4 & nbv,Int4 nbvx);
-			void  NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 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){
-						throw ErrorException(__FUNCT__,exprintf("!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){ throw ErrorException(__FUNCT__,exprintf("!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
-			}
-
-	};
-	/*}}}1*/
-	/*CLASS: GeometricalSubDomain{{{1*/
-	class GeometricalSubDomain {
-		public:
-			GeometricalEdge *edge;
-			int sens; // -1 or 1
-			Int4 ref;
-
-			//Inline methods
-			inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &);
-
-	};
-	/*}}}1*/
-	/*CLASS: SubDomain{{{1*/
-	class SubDomain {
-		public:
-			Triangle * head;
-			Int4  ref;  
-			int sens; // -1 or 1
-			Edge * edge; // to  geometrical 	
-
-			//Inline methods
-			inline void Set(const Triangles &,Int4,Triangles &);
-	};
-	/*}}}1*/
-	/*CLASS: VertexOnGeom{{{1*/
-	class VertexOnGeom{
-
-		public:
-
-			Vertex* mv;
-			Real8 abscisse;  
-			union{ 
-				GeometricalVertex * gv; // if abscisse <0; 
-				GeometricalEdge * ge;  // if abscisse in [0..1]
-			};
-
-			//Constructors/Destructors
-			VertexOnGeom(): mv(0),abscisse(0){gv=0;} 
-			VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
-			VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;}
-
-			//Operators
-			operator Vertex * () const  {return mv;}
-			operator GeometricalVertex * () const  {return gv;}
-			operator GeometricalEdge * () const  {return ge;}
-			operator const Real8 & () const {return abscisse;}
-
-			//Methods
-			int  OnGeomVertex()const {return this? abscisse <0 :0;}
-			int  OnGeomEdge() const {return this? abscisse >=0 :0;}
-			int  IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
-			void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;}
-
-			//Inline methods
-			inline void Set(const Triangles &,Int4,Triangles &);
-			inline void Set(const VertexOnGeom&,const Triangles &,Triangles &);  
-
-	};
-	/*}}}1*/
-	/*CLASS: VertexOnVertex{{{1*/
-	class VertexOnVertex {
-
-		public:
-			Vertex* v;
-			Vertex* bv;
-
-			//Constructors
-			VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){}
-			VertexOnVertex() {};
-
-			//Methods
-			void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}
-
-			//Inline methods
-			inline void Set(const Triangles &,Int4,Triangles &);
-	};
-	/*}}}1*/
-	/*CLASS: VertexOnEdge{{{1*/
-	class VertexOnEdge {
-
-		public:
-			Vertex* v;
-			Edge*   be;
-			Real8 abcisse;
-
-			//Constructors
-			VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {}
-			VertexOnEdge(){}
-
-			//Operators
-			operator Real8 () const { return abcisse;}
-			operator Vertex* () const { return v;}  
-			operator Edge* () const { return be;}  
-			Vertex & operator[](int i) const { return (*be)[i];}
-
-			//Methods
-			void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;}  
-
-			//Inline methods
-			inline void Set(const Triangles &,Int4,Triangles &);  
-	};
-	/*}}}1*/
-	/*CLASS: CrackedEdge{{{1*/
-	class CrackedEdge {
-
-		friend class Triangles;
-
-		class CrackedTriangle {
-			friend class Triangles;
-			friend class CrackedEdge;
-			Triangle* t; // edge of triangle t
-			int i;       //  edge number of in triangle
-			Edge *edge;  // the  2 edge 
-			Vertex *New[2]; // new vertex number 
-
-			//Constructors
-			CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;} 
-			CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;} 
-
-			//Methods
-			void Crack(){ 
-				Triangle & T(*t); 
-				int i0=VerticesOfTriangularEdge[i][0];
-				int i1=VerticesOfTriangularEdge[i][0];
-				if (!New[0] && !New[1]){
-					throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
-				}
-				T(i0) = New[0];
-				T(i1) = New[1];
-			}    
-			void UnCrack(){ 
-				Triangle & T(*t); 
-				int i0=VerticesOfTriangularEdge[i][0];
-				int i1=VerticesOfTriangularEdge[i][0];
-				if (!New[0] && !New[1]){
-					throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
-				}
-				T(i0) = TheVertex(T(i0));
-				T(i1) = TheVertex(T(i1));
-			} 
-			void Set() {
-				TriangleAdjacent ta ( FindTriangleAdjacent(*edge));
-				t = ta;
-				i = ta;
-
-				New[0] = ta.EdgeVertex(0);
-				New[1] = ta.EdgeVertex(1);
-				// warning the ref 
-			}    
-		};
-
-		public:  
-			CrackedTriangle a,b; 
-
-			//Constructors
-			CrackedEdge() :a(),b() {}
-			CrackedEdge(Edge * start, Int4  i,Int4 j) : a(start+i),b(start+j) {};
-			CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {};
-
-			//Methods
-			void Crack() { a.Crack(); b.Crack();}
-			void UnCrack() { a.UnCrack(); b.UnCrack();}
-			void Set() { a.Set(), b.Set();}
-	};
-	/*}}}1*/
-	/*CLASS: Triangles{{{1*/
-	class Triangles {
-		public:
-
-			Geometry & Gh;   // Geometry
-			Triangles & BTh; // Background Mesh Bth==*this =>no  background 
-			Int4 NbRef;      // counter of ref on the this class if 0 we can delete
-			Int4 nbvx,nbtx;  // nombre max  de sommets , de  triangles
-			Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,
-			Int4 NbOfQuad;  // nb of quadrangle 
-			Int4 NbSubDomains;
-			Int4 NbOutT;    // Nb of oudeside triangle
-			Int4 NbOfTriangleSearchFind;
-			Int4 NbOfSwapTriangle;
-			Vertex* vertices;
-			Int4 NbVerticesOnGeomVertex;
-			VertexOnGeom * VerticesOnGeomVertex;
-			Int4 NbVerticesOnGeomEdge;
-			VertexOnGeom * VerticesOnGeomEdge;
-			Int4 NbVertexOnBThVertex;
-			VertexOnVertex *VertexOnBThVertex;
-			Int4 NbVertexOnBThEdge;
-			VertexOnEdge *VertexOnBThEdge;
-			Int4 NbCrackedVertices;
-			Int4 NbCrackedEdges;
-			CrackedEdge *CrackedEdges;
-			R2 pmin,pmax;    // extrema
-			Real8 coefIcoor; // coef to integer Icoor1;
-			Triangle* triangles;
-			Edge* edges; 
-			QuadTree *quadtree;
-			Vertex** ordre;
-			SubDomain* subdomains;
-			ListofIntersectionTriangles  lIntTria;
-
-			//Constructors/Destructors
-			Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
-			Triangles(double* index,double* x,double* y,int nods,int nels);
-			Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); //copy operator
-			Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
-			Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
-				try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
-				catch(...) { this->~Triangles(); throw; }
-			}
-			Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
-				try { GeomToTriangles0(nbvx,bamgopts);}
-				catch(...) { this->~Triangles(); throw; }
-			}
-			~Triangles(); 
-
-			//Operators
-			const Vertex & operator[]  (Int4 i) const { return vertices[i];};
-			Vertex & operator[](Int4 i) { return vertices[i];};
-			const Triangle & operator()  (Int4 i) const { return triangles[i];};
-			Triangle & operator()(Int4 i) { return triangles[i];};
-
-			//Methods
-			void SetIntCoor(const char * from =0);
-			Real8 MinimalHmin() {return 2.0/coefIcoor;}
-			Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
-			I2 toI2(const R2 & P) const {
-				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
-							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
-			R2 toR2(const I2 & P) const {
-				return  R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
-			}
-			void AddVertex(Vertex & s,Triangle * t,Icoor2 *  =0) ;
-			void Insert();
-			void ForceBoundary();
-			void Renumber(BamgOpts* bamgopts);
-			void FindSubDomain(int OutSide=0);
-			Int4 TriangleReferenceList(Int4 *) const;
-			void ShowHistogram() const;
-			void ShowRegulaty() const;
-			void ReMakeTriangleContainingTheVertex();
-			void UnMarkUnSwapTriangle();
-			void SmoothMetric(Real8 raisonmax) ;
-			void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ;
-			void MaxSubDivision(Real8 maxsubdiv);
-			Edge** MakeGeometricalEdgeToEdge();
-			void SetVertexFieldOn();  
-			void SetVertexFieldOnBTh();
-			Int4 SplitInternalEdgeWithBorderVertices();
-			void MakeQuadrangles(double costheta);
-			int  SplitElement(int choice);
-			void MakeQuadTree();
-			void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1);
-			Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; 
-			void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
-			void ReNumberingVertex(Int4 * renu);
-			void SmoothingVertex(int =3,Real8=0.3);
-			Metric MetricAt (const R2 &) const;
-			GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex &  A, Vertex & B,Real8 theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
-			Int4 Number(const Triangle & t) const  { return &t - triangles;}
-			Int4 Number(const Triangle * t) const  { return t - triangles;}
-			Int4 Number(const Vertex & t) const  { return &t - vertices;}
-			Int4 Number(const Vertex * t) const  { return t - vertices;}
-			Int4 Number(const Edge & t) const  { return &t - edges;}
-			Int4 Number(const Edge * t) const  { return t - edges;}
-			Int4 Number2(const Triangle * t) const  {
-				return t - triangles;
-			}
-			Vertex* NearestVertex(Icoor1 i,Icoor1 j) ;
-			Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
-			void ReadMesh(double* index,double* x,double* y,int nods,int nels);
-			void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
-			void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
-			void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
-			void WriteMetric(BamgOpts* bamgopts);
-			void AddMetric(BamgOpts* bamgopts);
-			void BuildMetric0(BamgOpts* bamgopts);
-			void BuildMetric1(BamgOpts* bamgopts);
-			void AddGeometryMetric(BamgOpts* bamgopts);
-			int  isCracked() const {return NbCrackedVertices != 0;}
-			int  Crack();
-			int  UnCrack();
-			void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
-			void GenerateMeshProperties() ;
-			int  CrackMesh();
-
-		private:
-			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);
-
-	};
-	/*}}}1*/
-   /*CLASS: Geometry{{{1*/
-	class Geometry { 
-
-		public:
-			Int4 NbRef;     // counter of ref on the this class if 0 we can delete
-			Int4 nbvx,nbtx; // maximum number of vertices
-			Int4 nbv,nbt,nbiv,nbe; // number of vertices
-			Int4 NbSubDomains; // 
-			Int4 NbOfCurves;
-			GeometricalVertex* vertices;
-			Triangle* triangles; 
-			GeometricalEdge* edges;
-			QuadTree* quadtree;
-			GeometricalSubDomain* subdomains;
-			Curve* curves;
-			R2 pmin,pmax; // extrema
-			Real8 coefIcoor;  // coef to integer Icoor1;
-			Real8 MaximalAngleOfCorner;
-
-			//Constructor/Destructors
-			~Geometry(); 
-			Geometry(const Geometry & Gh); //Copy  Operator 
-			Geometry(int nbg,const Geometry** ag); // intersection operator 
-
-			//Operators
-			const GeometricalVertex & operator[]  (Int4 i) const { return vertices[i];};
-			GeometricalVertex & operator[](Int4 i) { return vertices[i];};
-			const  GeometricalEdge & operator()  (Int4 i) const { return edges[i];};
-			GeometricalEdge & operator()(Int4 i) { return edges[i];}; 
-
-			//Methods
-			int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
-			void Echo();
-			I2 toI2(const R2 & P) const {
-				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
-							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
-			Real8 MinimalHmin() {return 2.0/coefIcoor;}
-			Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
-			void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
-			void EmptyGeometry();
-			Geometry() {EmptyGeometry();}// empty Geometry
-			void AfterRead();
-			Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}
-			Int4 Number(const GeometricalVertex & t) const  { return &t - vertices;}
-			Int4 Number(const GeometricalVertex * t) const  { return t - vertices;}
-			Int4 Number(const GeometricalEdge & t) const  { return &t - edges;}
-			Int4 Number(const GeometricalEdge * t) const  { return t - edges;}
-			Int4 Number(const Curve * c) const  { return c - curves;}
-			void UnMarkEdges() {for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
-			GeometricalEdge *  ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
-			GeometricalEdge *  Contening(const R2 P,  GeometricalEdge * start) const;
-			void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
-	};
-	/*}}}1*/
 	
 	/*INLINE functions{{{1*/
@@ -980,21 +190,4 @@
 	/*}}}1*/
 
-   /*INLINE functions of CLASS Vertex{{{1*/
-	inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){ 
-		*this  = rec;
-	  }
-	Int4 inline  Vertex::Optim(int i,int koption){ 
-		Int4 ret=0;
-		if ( t && (vint >= 0 ) && (vint <3) )
-		  {
-			ret = t->Optim(vint,koption);
-			if(!i) 
-			  {
-				t =0; // for no future optime 
-				vint= 0; }
-		  }
-		return ret;
-	  }
-	/*}}}1*/
 	/*INLINE functions of CLASS GeometricalVertex{{{1*/
 	inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){ 
@@ -1189,4 +382,20 @@
 	  }
 	/*}}}1*/
+	/*INLINE functions of CLASS Vertex{{{1*/
+	Int4 inline Vertex::Optim(int i,int koption){ 
+		Int4 ret=0;
+		if ( t && (vint >= 0 ) && (vint <3) ){
+			ret = t->Optim(vint,koption);
+			if(!i){
+				t =0; // for no future optime 
+				vint= 0;
+			}
+		}
+		return ret;
+	}
+	inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){ 
+		*this  = rec;
+	}
+	/*}}}1*/
 	/*INLINE functions of CLASS VertexOnEdge{{{1*/
 	inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
Index: sm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 3231)
+++ 	(revision )
@@ -1,124 +1,0 @@
-#ifndef _METRIC_H
-#define _METRIC_H
-
-#include "R2.h"
-
-namespace bamg {
-
-	typedef P2<double,double>    D2;
-	typedef P2xP2<double,double> D2xD2;
-
-	class Metric;
-	class MatVVP2x2;
-
-
-	class Metric{
-
-		public:
-			//fields
-			Real8 a11,a21,a22;
-			//friends
-			friend class MatVVP2x2;
-			//functions
-			Metric(){};
-			Metric(const MatVVP2x2);
-			Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}
-			Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){}
-			Metric( Real8  a,const  Metric ma, Real8  b,const  Metric mb);
-			Metric(const Real8  a[3],const  Metric m0,const  Metric m1,const  Metric m2 );
-			void  Echo();
-			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);}
-			int   IntersectWith(const Metric M2);
-			inline void Box(Real4 &hx,Real4 &hy) const ;  
-			//operators
-			Metric operator*(Real8 c) const {Real8 c2=c*c;return  Metric(a11*c2,a21*c2,a22*c2);} 
-			Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return  Metric(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()(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;};
-
-	};
-
-	class MatVVP2x2{
-		public:
-			//fields
-			double lambda1,lambda2;
-			D2     v;
-			//friends
-			friend  class Metric;
-			//functions
-			MatVVP2x2(const Metric );
-			MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
-			void  Echo();
-			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) {Min(1.0/(h*h));}
-			void Maxh(double h) {Max(1.0/(h*h));}
-			void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
-			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();}
-			Real8 Aniso() const  { return sqrt( Aniso2());}
-			void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
-			inline void BoundAniso2(const Real8 coef);
-			//operators
-			void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
-	};
-
-	class SaveMetricInterpole {
-		friend Real8 LengthInterpole(const Metric ,const  Metric , R2 );
-		friend Real8 abscisseInterpole(const Metric ,const  Metric , R2 ,Real8 ,int );
-		int opt;
-		Real8 lab;
-		Real8 L[1024],S[1024];
-	};
-
-	extern SaveMetricInterpole  LastMetricInterpole; // for optimization 
-	//Functions
-	void  SimultaneousMatrixReduction( Metric M1,  Metric M2,D2xD2 &V);
-	Real8 LengthInterpole(const Metric Ma,const  Metric Mb, R2 AB);
-	Real8 abscisseInterpole(const Metric Ma,const  Metric Mb, R2 AB,Real8 s,int optim=0);
-
-	//inlines
-	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{  //TO BE CHECKED
-			if (lambda1 > lambda2)
-			 lambda1 = bamg::Min(lambda1,lambda2*coef);
-			else
-			 lambda2 = bamg::Min(lambda2,lambda1*coef);
-		}
-	}
-	inline Metric::Metric(const MatVVP2x2 M) {
-		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  Metric::Box(Real4 &hx,Real4 &hy) const {
-		Real8 d=  a11*a22-a21*a21;
-		hx = sqrt(a22/d);
-		hy = sqrt(a11/d);
-	}
-	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);
-	}
-}
-#endif
Index: sm/trunk/src/c/Bamgx/QuadTree.h
===================================================================
--- /issm/trunk/src/c/Bamgx/QuadTree.h	(revision 3231)
+++ 	(revision )
@@ -1,68 +1,0 @@
-#ifndef _QUADTREE_H
-#define _QUADTREE_H
-
-#include "../shared/shared.h"
-#include "../include/macros.h"
-#include "../toolkits/toolkits.h"
-
-namespace bamg {
-
-	typedef long IntQuad;
-
-	const int MaxDeep = 30;
-	const IntQuad MaxISize = ( 1L << MaxDeep); 
-
-	class Triangles;
-	class Vertex;
-
-	class QuadTree{
-		private:
-			class QuadTreeBox { 
-				public:
-					long n; 
-					//contains only one object form the list (either Vertex or QuadTreeBox)
-					// if n < 4 => Vertex else =>  QuadTreeBox;
-					union{
-						QuadTreeBox* b[4];
-						Vertex* v[4];
-					};
-			};
-			class StorageQuadTreeBox {
-				public:
-					QuadTreeBox *b,*bc,*be;
-					long len;
-					StorageQuadTreeBox* n; // next StorageQuadTreeBox
-					StorageQuadTreeBox(long ,StorageQuadTreeBox* =NULL);
-					~StorageQuadTreeBox() {
-						if(n) delete n;
-						delete [] b;
-					}
-					long  SizeOf() const {return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);}
-			};
-			StorageQuadTreeBox* sb;
-			long                lenStorageQuadTreeBox;
-
-		public:
-			//fields
-			QuadTreeBox* root;
-			Triangles*   th;
-			//functions
-			~QuadTree();
-			QuadTree(Triangles * t,long nbv=-1);
-			QuadTree();
-			long    NbQuadTreeBox,NbVertices;
-			long    NbQuadTreeBoxSearch,NbVerticesSearch;
-			Vertex* NearestVertex(Icoor1 i,Icoor1 j);
-			Vertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);
-			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++;
-			}
-	};
-}
-#endif
Index: sm/trunk/src/c/Bamgx/R2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/R2.h	(revision 3231)
+++ 	(revision )
@@ -1,101 +1,0 @@
-#ifndef _R2_H
-#define _R2_H
-
-namespace bamg {
-
-	template <class R,class RR> class P2xP2;
-	template <class R,class RR>
-	  class P2 {
-		  public:  
-			  //fields
-			  R x,y;
-			  //functions
-			  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) {}
-			  void Echo(){
-				  printf("Member of P2:\n");
-				  printf("   x: %g\n",x);
-				  printf("   y: %g\n",y);
-			  }
-			  //operators
-			  RR       operator,(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} //scalar product
-			  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);}
-			  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 {
-		  private:
-			  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:
-			  //fields
-			  P2<R,RR> x,y; 
-			  //functions
-			  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) {}
-			  void Echo(){
-				  printf("Member of P2xP2:\n");
-				  printf("   x.x: %g   x.y: %g\n",x.x,x.y);
-				  printf("   y.x: %g   y.x: %g\n",y.x,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);} 
-			  //Operators
-			  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);
-			  }
-	  };  
-
-	//inline functions
-	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 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 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 P2<R,RR> Orthogonal (const P2<R,RR> x) {
-		  return  P2<R,RR>(-x.y,x.x);
-	  } 
-}
-#endif
Index: sm/trunk/src/c/Bamgx/SetOfE4.h
===================================================================
--- /issm/trunk/src/c/Bamgx/SetOfE4.h	(revision 3231)
+++ 	(revision )
@@ -1,38 +1,0 @@
-#ifndef _SetOfEdge4_h
-#define _SetOfEdge4_h
-
-namespace bamg {
-
-	class SetOfEdges4;
-
-	class Int4Edge{
-		friend class SetOfEdges4;
-		public:
-		Int4 i,j;
-		Int4 next; 
-	};
-
-	class SetOfEdges4 {
-
-		private:
-			Int4 nx,nbax,NbOfEdges;
-			Int4* head; 
-			Int4Edge* Edges;
-
-		public:
-			SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 
-			~SetOfEdges4() {delete [] head; delete [] Edges;}
-			Int4 add (Int4 ii,Int4 jj);
-			Int4 SortAndAdd (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
-			Int4  nb(){return NbOfEdges;}
-			Int4 find (Int4 ii,Int4 jj);
-			Int4 SortAndFind (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;}
-
-			//operators
-			Int4Edge & operator[](Int4 k){return  Edges[k];}
-	};
-}
-#endif 
Index: /issm/trunk/src/c/Bamgx/meshtype.h
===================================================================
--- /issm/trunk/src/c/Bamgx/meshtype.h	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/meshtype.h	(revision 3232)
@@ -2,5 +2,10 @@
 #define MESHTYPE_H
 
-#include "R2.h"
+#include "objects/R2.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <ctime>
 
 namespace bamg {
Index: /issm/trunk/src/c/Bamgx/objects/CrackedEdge.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/CrackedEdge.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/CrackedEdge.h	(revision 3232)
@@ -0,0 +1,83 @@
+#ifndef _CRACKEDEDGE_H_
+#define _CRACKEDEDGE_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+#include "TriangleAdjacent.h"
+#include "Edge.h"
+#include "Triangle.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+
+	class CrackedEdge {
+
+		friend class Triangles;
+
+		class CrackedTriangle {
+			friend class Triangles;
+			friend class CrackedEdge;
+			Triangle* t; // edge of triangle t
+			int i;       //  edge number of in triangle
+			Edge *edge;  // the  2 edge 
+			Vertex *New[2]; // new vertex number 
+
+			//Constructors
+			CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;} 
+			CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;} 
+
+			//Methods
+			void Crack(){ 
+				Triangle & T(*t); 
+				int i0=VerticesOfTriangularEdge[i][0];
+				int i1=VerticesOfTriangularEdge[i][0];
+				if (!New[0] && !New[1]){
+					throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
+				}
+				T(i0) = New[0];
+				T(i1) = New[1];
+			}    
+			void UnCrack(){ 
+				Triangle & T(*t); 
+				int i0=VerticesOfTriangularEdge[i][0];
+				int i1=VerticesOfTriangularEdge[i][0];
+				if (!New[0] && !New[1]){
+					throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
+				}
+				T(i0) = TheVertex(T(i0));
+				T(i1) = TheVertex(T(i1));
+			} 
+			void Set() {
+				TriangleAdjacent ta ( FindTriangleAdjacent(*edge));
+				t = ta;
+				i = ta;
+
+				New[0] = ta.EdgeVertex(0);
+				New[1] = ta.EdgeVertex(1);
+				// warning the ref 
+			}    
+		};
+
+		public:  
+			CrackedTriangle a,b; 
+
+			//Constructors
+			CrackedEdge() :a(),b() {}
+			CrackedEdge(Edge * start, Int4  i,Int4 j) : a(start+i),b(start+j) {};
+			CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {};
+
+			//Methods
+			void Crack() { a.Crack(); b.Crack();}
+			void UnCrack() { a.UnCrack(); b.UnCrack();}
+			void Set() { a.Set(), b.Set();}
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/Curve.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Curve.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Curve.h	(revision 3232)
@@ -0,0 +1,33 @@
+#ifndef _CURVE_H_
+#define _CURVE_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "GeometricalEdge.h"
+
+namespace bamg {
+
+	//classes
+	class Geometry;
+
+	class Curve {
+		public:
+			GeometricalEdge * be,*ee; // begin et end edge
+			int kb,ke;  //  begin vetex and end vertex
+			Curve* next; // next curve equi to this
+			bool master; // true => of equi curve point on this curve  
+
+			//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 
+
+			//Inline methods
+			inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/Direction.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Direction.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Direction.h	(revision 3232)
@@ -0,0 +1,32 @@
+#ifndef _DIRECTION_H_
+#define _DIRECTION_H_
+
+#include "../meshtype.h"
+
+namespace bamg {
+
+	class Direction {
+		private:
+			Icoor1 dir;
+
+		public:
+			//Methods
+			Direction(): dir(MaxICoor){}; //  no direction set
+			Direction(Icoor1 i,Icoor1 j) {
+				Icoor2 n2 = 2*(Abs(i)+Abs(j));  
+				Icoor2 r  = MaxICoor* (Icoor2) i;
+				Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number 
+				dir = (j>0) ? r1 : r1+1;          // odd-> j>0 even-> j<0
+			}
+			int sens(Icoor1 i,Icoor1 j) {
+				int r =1; 
+				if (dir!= MaxICoor) {
+					Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1);
+					r = (x*i + y*j) >=0;
+				}
+				return r;
+			}
+	};
+}
+#endif
+
Index: /issm/trunk/src/c/Bamgx/objects/DoubleAndInt4.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/DoubleAndInt4.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/DoubleAndInt4.h	(revision 3232)
@@ -0,0 +1,19 @@
+#ifndef _DOUBLEANDINT4_H_
+#define _DOUBLEANDINT4_H_
+
+#include "../meshtype.h"
+
+namespace bamg {
+
+	class DoubleAndInt4 {
+		//class used by Triangles::MakeQuadrangles
+
+		public:
+			double q;
+			Int4 i3j;
+
+			//Operators
+			int operator<(DoubleAndInt4 a){return q > a.q;}
+	};
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/Edge.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Edge.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Edge.h	(revision 3232)
@@ -0,0 +1,58 @@
+#ifndef _EDGE_H_
+#define _EDGE_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+#include "TriangleAdjacent.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+	class GeometricalEdge;
+
+	class Edge {
+
+		public:
+			Vertex* v[2];
+			Int4 ref;
+			GeometricalEdge* onGeometry;
+			Edge* adj[2]; // the 2 adj edges if on the same curve 
+
+			//Operators
+			Vertex & operator[](int i){return *v[i];};
+			Vertex * operator()(int i){return v[i];};
+			R2       operator()(double t) const; // return the point 
+			const Vertex & operator[](int i) const { return *v[i];};
+
+			//Methods
+			void ReNumbering(Vertex *vb,Vertex *ve, Int4 *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)){
+					throw ErrorException(__FUNCT__,exprintf("Intersection bug"));
+				}
+				if (adj[0]!=&e && adj[1]!=&e){
+					throw ErrorException(__FUNCT__,exprintf("adj[0]!=&e && adj[1]!=&e"));
+				}
+				return adj[0]==&e ? 0 : 1;
+			}
+
+			//Inline methods
+			inline void Set(const Triangles &,Int4,Triangles &);
+
+	};
+
+	//FOR NOW
+	inline TriangleAdjacent FindTriangleAdjacent(Edge &E);
+
+}
+#endif
+
Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp	(revision 3232)
@@ -4,11 +4,5 @@
 #include <time.h>
 
-#include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
+#include "GeometricalEdge.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h	(revision 3232)
@@ -0,0 +1,63 @@
+#ifndef _GEOMETRICALEDGE_H_
+#define _GEOMETRICALEDGE_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Edge.h"
+#include "GeometricalVertex.h"
+#include "GeometricalEdge.h"
+
+namespace bamg {
+
+	//classes
+	class Geometry;
+
+	class GeometricalEdge {
+
+		public:
+			GeometricalVertex* v[2];
+			Int4 ref;
+			Int4 CurveNumber;
+			R2   tg[2]; // the 2 tangentes (tg[0] =0 => no continuity)
+			GeometricalEdge* Adj[2]; 
+			int DirAdj[2];
+			int flag ;
+			GeometricalEdge* link; // if   Cracked() or Equi()
+
+			//Operators
+			GeometricalVertex       & operator[](int i){return *v[i];};
+			const GeometricalVertex & operator[](int i) const { return *v[i];};
+			GeometricalVertex       * operator()(int i){return v[i];};  
+
+			//Methods
+			R2 F(Real8 theta) const ; // parametrization of the curve edge
+			Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente
+			int  Cracked() const {return flag & 1;}
+			int  Dup() const { return flag & 32;}
+			int  Equi()const {return flag & 2;}
+			int  ReverseEqui()const {return flag & 128;}
+			int  TgA()const {return flag &4;}
+			int  TgB()const {return flag &8;}
+			int  Tg(int i) const{return i==0 ? TgA() : TgB();}
+			int  Mark()const {return flag &16;}
+			int  Required() { return flag & 64;}
+			void SetCracked() { flag |= 1;}
+			void SetDup()     { flag |= 32;} // not a real edge 
+			void SetEqui()    { flag |= 2;}
+			void SetTgA()     { flag|=4;}
+			void SetTgB()     { flag|=8;}
+			void SetMark()    { flag|=16;}
+			void SetUnMark()  { flag &= 1007 /* 1023-16*/;}
+			void SetRequired() { flag |= 64;}
+			void SetReverseEqui() {flag |= 128;}
+
+			//Inline methods
+			inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/GeometricalSubDomain.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/GeometricalSubDomain.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/GeometricalSubDomain.h	(revision 3232)
@@ -0,0 +1,28 @@
+#ifndef _GEOMETRICALSUBDOMAIN_H_
+#define _GEOMETRICALSUBDOMAIN_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "GeometricalEdge.h"
+
+namespace bamg {
+
+	//classes
+	class Geometry;
+
+	class GeometricalSubDomain {
+		public:
+			GeometricalEdge *edge;
+			int sens; // -1 or 1
+			Int4 ref;
+
+			//Inline methods
+			inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &);
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/GeometricalVertex.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/GeometricalVertex.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/GeometricalVertex.h	(revision 3232)
@@ -0,0 +1,42 @@
+#ifndef _GEOMETRICALVERTEX_H_
+#define _GEOMETRICALVERTEX_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+
+namespace bamg {
+
+	//classes
+	class Geometry;
+
+	class GeometricalVertex : public Vertex { 
+
+		public:
+			friend class Geometry;
+
+			int cas;
+			GeometricalVertex* link; //  link all the same GeometricalVertex circular (Crack) 
+
+			//Methods
+			int  Corner() const {return cas&4;}
+			int  Required()const {return cas&6;}// a corner is required
+			int  IsThe() const { return link == this;}  
+			void SetCorner(){ cas |= 4;}
+			void SetRequired(){ cas |= 2;}
+			void Set(){cas=0;}
+			GeometricalVertex() :cas(0), link(this) {};
+			GeometricalVertex* The() {
+				if (!link){ throw ErrorException(__FUNCT__,exprintf("!link"));}
+				return link;
+			}// return a unique vertex
+
+			//Inline methods
+			inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew);
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3232)
@@ -5,10 +5,4 @@
 
 #include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.h	(revision 3232)
@@ -0,0 +1,79 @@
+#ifndef _GEOMETRY_H_
+#define _GEOMETRY_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+#include "Edge.h"
+#include "GeometricalVertex.h"
+#include "GeometricalEdge.h"
+#include "Curve.h"
+#include "Triangle.h"
+#include "GeometricalSubDomain.h"
+#include "SubDomain.h"
+#include "VertexOnGeom.h"
+#include "QuadTree.h"
+
+namespace bamg {
+
+   /*CLASS: Geometry{{{1*/
+	class Geometry { 
+
+		public:
+			Int4 NbRef;     // counter of ref on the this class if 0 we can delete
+			Int4 nbvx,nbtx; // maximum number of vertices
+			Int4 nbv,nbt,nbiv,nbe; // number of vertices
+			Int4 NbSubDomains; // 
+			Int4 NbOfCurves;
+			GeometricalVertex* vertices;
+			Triangle* triangles; 
+			GeometricalEdge* edges;
+			QuadTree* quadtree;
+			GeometricalSubDomain* subdomains;
+			Curve* curves;
+			R2 pmin,pmax; // extrema
+			Real8 coefIcoor;  // coef to integer Icoor1;
+			Real8 MaximalAngleOfCorner;
+
+			//Constructor/Destructors
+			~Geometry(); 
+			Geometry(const Geometry & Gh); //Copy  Operator 
+			Geometry(int nbg,const Geometry** ag); // intersection operator 
+
+			//Operators
+			const GeometricalVertex & operator[]  (Int4 i) const { return vertices[i];};
+			GeometricalVertex & operator[](Int4 i) { return vertices[i];};
+			const  GeometricalEdge & operator()  (Int4 i) const { return edges[i];};
+			GeometricalEdge & operator()(Int4 i) { return edges[i];}; 
+
+			//Methods
+			int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
+			void Echo();
+			I2 toI2(const R2 & P) const {
+				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
+							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
+			Real8 MinimalHmin() {return 2.0/coefIcoor;}
+			Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
+			void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
+			void EmptyGeometry();
+			Geometry() {EmptyGeometry();}// empty Geometry
+			void AfterRead();
+			Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}
+			Int4 Number(const GeometricalVertex & t) const  { return &t - vertices;}
+			Int4 Number(const GeometricalVertex * t) const  { return t - vertices;}
+			Int4 Number(const GeometricalEdge & t) const  { return &t - edges;}
+			Int4 Number(const GeometricalEdge * t) const  { return t - edges;}
+			Int4 Number(const Curve * c) const  { return c - curves;}
+			void UnMarkEdges() {for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
+			GeometricalEdge *  ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
+			GeometricalEdge *  Contening(const R2 P,  GeometricalEdge * start) const;
+			void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
+	};
+	/*}}}1*/
+	
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 3232)
@@ -5,10 +5,4 @@
 
 #include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.h	(revision 3232)
@@ -0,0 +1,122 @@
+#ifndef _LISTOFINTERSECTIONTRIANGLES_H_
+#define _LISTOFINTERSECTIONTRIANGLES_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Metric.h"
+#include "Vertex.h"
+#include "Edge.h"
+#include "GeometricalEdge.h"
+#include "Triangle.h"
+
+namespace bamg {
+
+	class Triangles;
+
+	class ListofIntersectionTriangles {
+
+		class IntersectionTriangles {
+			public: 
+				Triangle* t;
+				Real8  bary[3];  // use if t != 0
+				R2 x;
+				Metric m;
+				Real8 s; // curvilinear coordinate 
+				Real8 sp;// length of the previous segment in m
+				Real8 sn;// length of the next segment in m
+		};
+
+		class SegInterpolation {
+			public:
+				GeometricalEdge * e;
+				Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
+				Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
+				int last;// last index  in ListofIntersectionTriangles for this Sub seg of edge
+
+				//Methods
+				R2 F(Real8 s){ 
+					Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01;
+					if (lBegin>s || s>lEnd){
+						throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd"));
+					}
+					return e->F(sBegin*c0+sEnd*c1);
+				}
+		};
+
+		public:
+
+			int MaxSize;
+			int Size;
+			Real8 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;
+			} 
+
+			//Operators
+			IntersectionTriangles & operator[](int i) {return lIntTria[i];}
+			operator int&() {return Size;}
+
+			//Methods
+			void  init(){state=0;len=0;Size=0;}
+			int   NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
+			int   NewItem(R2,const Metric & );
+			void  SplitEdge(const Triangles & ,const R2 &,const R2  &,int nbegin=0); 
+			Real8 Length(); 
+			Int4  NewPoints(Vertex *,Int4 & nbv,Int4 nbvx);
+			void  NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 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){
+						throw ErrorException(__FUNCT__,exprintf("!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){ throw ErrorException(__FUNCT__,exprintf("!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
+			}
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp	(revision 3232)
@@ -4,11 +4,5 @@
 #include <ctime>
 
-#include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
+#include "Metric.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/Metric.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Metric.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/Metric.cpp	(revision 3232)
@@ -4,11 +4,5 @@
 #include <time.h>
 
-#include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
+#include "Metric.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Metric.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Metric.h	(revision 3232)
@@ -0,0 +1,130 @@
+#ifndef _METRIC_H
+#define _METRIC_H
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "R2.h"
+
+namespace bamg {
+
+	typedef P2<double,double>    D2;
+	typedef P2xP2<double,double> D2xD2;
+
+	class Metric;
+	class MatVVP2x2;
+
+
+	class Metric{
+
+		public:
+			//fields
+			Real8 a11,a21,a22;
+			//friends
+			friend class MatVVP2x2;
+			//functions
+			Metric(){};
+			Metric(const MatVVP2x2);
+			Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}
+			Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){}
+			Metric( Real8  a,const  Metric ma, Real8  b,const  Metric mb);
+			Metric(const Real8  a[3],const  Metric m0,const  Metric m1,const  Metric m2 );
+			void  Echo();
+			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);}
+			int   IntersectWith(const Metric M2);
+			inline void Box(Real4 &hx,Real4 &hy) const ;  
+			//operators
+			Metric operator*(Real8 c) const {Real8 c2=c*c;return  Metric(a11*c2,a21*c2,a22*c2);} 
+			Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return  Metric(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()(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;};
+
+	};
+
+	class MatVVP2x2{
+		public:
+			//fields
+			double lambda1,lambda2;
+			D2     v;
+			//friends
+			friend  class Metric;
+			//functions
+			MatVVP2x2(const Metric );
+			MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
+			void  Echo();
+			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) {Min(1.0/(h*h));}
+			void Maxh(double h) {Max(1.0/(h*h));}
+			void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
+			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();}
+			Real8 Aniso() const  { return sqrt( Aniso2());}
+			void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
+			inline void BoundAniso2(const Real8 coef);
+			//operators
+			void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
+	};
+
+	class SaveMetricInterpole {
+		friend Real8 LengthInterpole(const Metric ,const  Metric , R2 );
+		friend Real8 abscisseInterpole(const Metric ,const  Metric , R2 ,Real8 ,int );
+		int opt;
+		Real8 lab;
+		Real8 L[1024],S[1024];
+	};
+
+	extern SaveMetricInterpole  LastMetricInterpole; // for optimization 
+	//Functions
+	void  SimultaneousMatrixReduction( Metric M1,  Metric M2,D2xD2 &V);
+	Real8 LengthInterpole(const Metric Ma,const  Metric Mb, R2 AB);
+	Real8 abscisseInterpole(const Metric Ma,const  Metric Mb, R2 AB,Real8 s,int optim=0);
+
+	//inlines
+	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{  //TO BE CHECKED
+			if (lambda1 > lambda2)
+			 lambda1 = bamg::Min(lambda1,lambda2*coef);
+			else
+			 lambda2 = bamg::Min(lambda2,lambda1*coef);
+		}
+	}
+	inline Metric::Metric(const MatVVP2x2 M) {
+		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  Metric::Box(Real4 &hx,Real4 &hy) const {
+		Real8 d=  a11*a22-a21*a21;
+		hx = sqrt(a22/d);
+		hy = sqrt(a11/d);
+	}
+	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);
+	}
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp	(revision 3232)
@@ -3,10 +3,6 @@
 #include <stdlib.h>
 
+//#include "QuadTree.h"
 #include "../Mesh2.h"
-#include "../QuadTree.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/QuadTree.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/QuadTree.h	(revision 3232)
@@ -0,0 +1,69 @@
+#ifndef _QUADTREE_H
+#define _QUADTREE_H
+
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+
+namespace bamg {
+
+	typedef long IntQuad;
+
+	const int MaxDeep = 30;
+	const IntQuad MaxISize = ( 1L << MaxDeep); 
+
+	class Triangles;
+	class Vertex;
+
+	class QuadTree{
+		private:
+			class QuadTreeBox { 
+				public:
+					long n; 
+					//contains only one object form the list (either Vertex or QuadTreeBox)
+					// if n < 4 => Vertex else =>  QuadTreeBox;
+					union{
+						QuadTreeBox* b[4];
+						Vertex* v[4];
+					};
+			};
+			class StorageQuadTreeBox {
+				public:
+					QuadTreeBox *b,*bc,*be;
+					long len;
+					StorageQuadTreeBox* n; // next StorageQuadTreeBox
+					StorageQuadTreeBox(long ,StorageQuadTreeBox* =NULL);
+					~StorageQuadTreeBox() {
+						if(n) delete n;
+						delete [] b;
+					}
+					long  SizeOf() const {return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);}
+			};
+			StorageQuadTreeBox* sb;
+			long                lenStorageQuadTreeBox;
+
+		public:
+			//fields
+			QuadTreeBox* root;
+			Triangles*   th;
+			//functions
+			~QuadTree();
+			QuadTree(Triangles * t,long nbv=-1);
+			QuadTree();
+			long    NbQuadTreeBox,NbVertices;
+			long    NbQuadTreeBoxSearch,NbVerticesSearch;
+			Vertex* NearestVertex(Icoor1 i,Icoor1 j);
+			Vertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);
+			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++;
+			}
+	};
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/R2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/R2.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/R2.h	(revision 3232)
@@ -0,0 +1,101 @@
+#ifndef _R2_H
+#define _R2_H
+
+namespace bamg {
+
+	template <class R,class RR> class P2xP2;
+	template <class R,class RR>
+	  class P2 {
+		  public:  
+			  //fields
+			  R x,y;
+			  //functions
+			  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) {}
+			  void Echo(){
+				  printf("Member of P2:\n");
+				  printf("   x: %g\n",x);
+				  printf("   y: %g\n",y);
+			  }
+			  //operators
+			  RR       operator,(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} //scalar product
+			  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);}
+			  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 {
+		  private:
+			  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:
+			  //fields
+			  P2<R,RR> x,y; 
+			  //functions
+			  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) {}
+			  void Echo(){
+				  printf("Member of P2xP2:\n");
+				  printf("   x.x: %g   x.y: %g\n",x.x,x.y);
+				  printf("   y.x: %g   y.x: %g\n",y.x,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);} 
+			  //Operators
+			  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);
+			  }
+	  };  
+
+	//inline functions
+	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 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 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 P2<R,RR> Orthogonal (const P2<R,RR> x) {
+		  return  P2<R,RR>(-x.y,x.x);
+	  } 
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp	(revision 3232)
@@ -1,7 +1,4 @@
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
-#include "../meshtype.h"
-#include "../SetOfE4.h"
+
+#include "../Mesh2.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/SetOfE4.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/SetOfE4.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/SetOfE4.h	(revision 3232)
@@ -0,0 +1,38 @@
+#ifndef _SetOfEdge4_h
+#define _SetOfEdge4_h
+
+namespace bamg {
+
+	class SetOfEdges4;
+
+	class Int4Edge{
+		friend class SetOfEdges4;
+		public:
+		Int4 i,j;
+		Int4 next; 
+	};
+
+	class SetOfEdges4 {
+
+		private:
+			Int4 nx,nbax,NbOfEdges;
+			Int4* head; 
+			Int4Edge* Edges;
+
+		public:
+			SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 
+			~SetOfEdges4() {delete [] head; delete [] Edges;}
+			Int4 add (Int4 ii,Int4 jj);
+			Int4 SortAndAdd (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
+			Int4  nb(){return NbOfEdges;}
+			Int4 find (Int4 ii,Int4 jj);
+			Int4 SortAndFind (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;}
+
+			//operators
+			Int4Edge & operator[](Int4 k){return  Edges[k];}
+	};
+}
+#endif 
Index: /issm/trunk/src/c/Bamgx/objects/SubDomain.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/SubDomain.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/SubDomain.h	(revision 3232)
@@ -0,0 +1,30 @@
+#ifndef _SUBDOMAIN_H_
+#define _SUBDOMAIN_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Edge.h"
+#include "Triangle.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+
+	class SubDomain {
+		public:
+			Triangle * head;
+			Int4  ref;  
+			int sens; // -1 or 1
+			Edge* edge; // to  geometrical 	
+
+			//Inline methods
+			inline void Set(const Triangles &,Int4,Triangles &);
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 3232)
@@ -5,10 +5,4 @@
 
 #include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/Triangle.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangle.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Triangle.h	(revision 3232)
@@ -0,0 +1,135 @@
+#ifndef _TRIANGLE_H_
+#define _TRIANGLE_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+#include "TriangleAdjacent.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+
+	class Triangle {
+
+		friend class TriangleAdjacent;
+
+		private:
+			Vertex*   TriaVertices[3];      // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
+			Triangle* TriaAdjTriangles[3];  // 3 pointers toward the adjacent triangles
+			Int1      TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1
+
+		public: 
+			Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
+			union { 
+				Triangle * link ;
+				Int4 color;
+			};
+
+			//Constructors/Destructors
+			Triangle() {}
+			Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k);
+			Triangle(Vertex *v0,Vertex *v1,Vertex *v2);
+
+			//Operators
+			const Vertex & operator[](int i) const {return *TriaVertices[i];};
+			Vertex & operator[](int i)  {return *TriaVertices[i];};
+			const Vertex * operator()(int i) const {return TriaVertices[i];};
+			Vertex * & operator()(int i)  {return TriaVertices[i];};
+
+			//Methods
+			void   Echo();
+			int    swap(Int2 a1,int=0);
+			Int4   Optim(Int2 a,int =0);
+			int    Locked(int a)const { return TriaAdjSharedEdge[a]&4;} 
+			int    Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} 
+			int    Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;}
+			int    GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;}
+			void   SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);}
+			double QualityQuad(int a,int option=1) const;
+			Int1   NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the  adjacent edge in adj tria  
+			TriangleAdjacent FindBoundaryEdge(int i) const;
+			TriangleAdjacent Adj(int i)  const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
+			Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
+			Triangle* Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
+			void  ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
+				if (link  >=tb && link  <te) link  = tb + renu[link -tb];
+				if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
+				if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb];
+				if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];    
+			}
+			void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
+				if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
+				if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
+				if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb];    
+			}
+			void SetAdjAdj(Int1 a){
+				a &= 3;
+				register Triangle *tt=TriaAdjTriangles[a];
+				TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap
+				register Int1 aatt = TriaAdjSharedEdge[a] & 3;
+				if(tt){ 
+					tt->TriaAdjTriangles[aatt]=this;
+					tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark 
+			  }
+			void SetAdj2(Int1 a,Triangle *t,Int1 aat){
+				TriaAdjTriangles[a]=t;   //the adjacent triangle to the edge a is t
+				TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle
+				if(t) { //if t!=NULL add adjacent triangle to t (this)
+					t->TriaAdjTriangles[aat]=this;
+					t->TriaAdjSharedEdge[aat]=a;
+				}
+			}
+			void SetTriangleContainingTheVertex() { 
+				if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
+				if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);
+				if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2);
+			}
+			void SetHidden(int a){
+				//Get Adjacent Triangle number a
+				register Triangle* t = TriaAdjTriangles[a];
+				//if it exist
+				//C|=D -> C=(C|D) bitwise inclusive OR
+				if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16;
+				TriaAdjSharedEdge[a] |= 16;
+			}
+			void SetCracked(int a){
+				register Triangle * t = TriaAdjTriangles[a];
+				if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32;
+				TriaAdjSharedEdge[a] |= 32;
+			}
+
+			void SetLocked(int a){
+				//mark the edge as on Boundary
+				register Triangle * t = TriaAdjTriangles[a];
+				t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4;
+				TriaAdjSharedEdge[a] |= 4;
+			}
+			void SetMarkUnSwap(int a){
+				register Triangle * t = TriaAdjTriangles[a];
+				t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8;
+				TriaAdjSharedEdge[a] |=8 ;
+			}
+			void SetUnMarkUnSwap(int a){ 
+				register Triangle * t = TriaAdjTriangles[a];
+				t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 
+				TriaAdjSharedEdge[a] &=55 ;
+			}
+			void SetDet() {
+				if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2])    det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
+				else det = -1; }
+
+			//Inline methods
+			inline Real4 qualite() ;
+			inline void  Set(const Triangle &,const Triangles &,Triangles &);
+			inline int   In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
+
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/TriangleAdjacent.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/TriangleAdjacent.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/TriangleAdjacent.h	(revision 3232)
@@ -0,0 +1,53 @@
+#ifndef _TRIANGLEADJACENT_H_
+#define _TRIANGLEADJACENT_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+
+namespace bamg {
+
+	//classes
+	class Triangle;
+
+	class TriangleAdjacent {
+
+		public:
+			Triangle* t; //pointer toward the triangle
+			int  a;      //Edge number
+
+			//Constructors
+			TriangleAdjacent() {};
+			TriangleAdjacent(Triangle* tt,int  aa): t(tt),a(aa &3) {};
+
+			//Operators
+			operator Triangle * () const {return t;}
+			operator Triangle & () const {return *t;}
+			operator int() const {return a;}
+			TriangleAdjacent & operator++(){ a= NextEdge[a]; return *this; }
+			TriangleAdjacent operator--(){ a= PreviousEdge[a]; return *this; }
+
+			//Methods
+			int swap();
+
+			//Inline methods
+			inline TriangleAdjacent  Adj() const ;
+			inline void SetAdj2(const TriangleAdjacent& , int =0);
+			inline Vertex *  EdgeVertex(const int &) const ;
+			inline Vertex *  OppositeVertex() const ;
+			inline Icoor2 & det() const;
+			inline int Locked() const  ;
+			inline int GetAllFlag_UnSwap() const ;
+			inline void SetLock();
+			inline int MarkUnSwap()  const;
+			inline void SetMarkUnSwap();
+			inline void SetCracked();
+			inline int Cracked() const ;
+	};
+}
+#endif
+
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3232)
@@ -5,10 +5,6 @@
 
 #include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
+#include "QuadTree.h"
+#include "SetOfE4.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3232)
@@ -0,0 +1,162 @@
+#ifndef _TRIANGLES_H_
+#define _TRIANGLES_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Metric.h"
+#include "SetOfE4.h"
+#include "DoubleAndInt4.h"
+#include "Direction.h"
+#include "Vertex.h"
+#include "TriangleAdjacent.h"
+#include "Edge.h"
+#include "GeometricalVertex.h"
+#include "GeometricalEdge.h"
+#include "Curve.h"
+#include "Triangle.h"
+#include "ListofIntersectionTriangles.h"
+#include "GeometricalSubDomain.h"
+#include "SubDomain.h"
+#include "VertexOnGeom.h"
+#include "VertexOnVertex.h"
+#include "VertexOnEdge.h"
+#include "CrackedEdge.h"
+#include "QuadTree.h"
+
+namespace bamg {
+
+	//classes
+	class Geometry;
+
+	class Triangles {
+		public:
+
+			Geometry & Gh;   // Geometry
+			Triangles & BTh; // Background Mesh Bth==*this =>no  background 
+			Int4 NbRef;      // counter of ref on the this class if 0 we can delete
+			Int4 nbvx,nbtx;  // nombre max  de sommets , de  triangles
+			Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,
+			Int4 NbOfQuad;  // nb of quadrangle 
+			Int4 NbSubDomains;
+			Int4 NbOutT;    // Nb of oudeside triangle
+			Int4 NbOfTriangleSearchFind;
+			Int4 NbOfSwapTriangle;
+			Vertex* vertices;
+			Int4 NbVerticesOnGeomVertex;
+			VertexOnGeom * VerticesOnGeomVertex;
+			Int4 NbVerticesOnGeomEdge;
+			VertexOnGeom * VerticesOnGeomEdge;
+			Int4 NbVertexOnBThVertex;
+			VertexOnVertex *VertexOnBThVertex;
+			Int4 NbVertexOnBThEdge;
+			VertexOnEdge *VertexOnBThEdge;
+			Int4 NbCrackedVertices;
+			Int4 NbCrackedEdges;
+			CrackedEdge *CrackedEdges;
+			R2 pmin,pmax;    // extrema
+			Real8 coefIcoor; // coef to integer Icoor1;
+			Triangle* triangles;
+			Edge* edges; 
+			QuadTree* quadtree;
+			Vertex** ordre;
+			SubDomain* subdomains;
+			ListofIntersectionTriangles lIntTria;
+
+			//Constructors/Destructors
+			Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
+			Triangles(double* index,double* x,double* y,int nods,int nels);
+			Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); //copy operator
+			Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
+			Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
+				try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
+				catch(...) { this->~Triangles(); throw; }
+			}
+			Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
+				try { GeomToTriangles0(nbvx,bamgopts);}
+				catch(...) { this->~Triangles(); throw; }
+			}
+			~Triangles(); 
+
+			//Operators
+			const Vertex & operator[]  (Int4 i) const { return vertices[i];};
+			Vertex & operator[](Int4 i) { return vertices[i];};
+			const Triangle & operator()  (Int4 i) const { return triangles[i];};
+			Triangle & operator()(Int4 i) { return triangles[i];};
+
+			//Methods
+			void SetIntCoor(const char * from =0);
+			Real8 MinimalHmin() {return 2.0/coefIcoor;}
+			Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
+			I2 toI2(const R2 & P) const {
+				return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
+							,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
+			R2 toR2(const I2 & P) const {
+				return  R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
+			}
+			void AddVertex(Vertex & s,Triangle * t,Icoor2 *  =0) ;
+			void Insert();
+			void ForceBoundary();
+			void Renumber(BamgOpts* bamgopts);
+			void FindSubDomain(int OutSide=0);
+			Int4 TriangleReferenceList(Int4 *) const;
+			void ShowHistogram() const;
+			void ShowRegulaty() const;
+			void ReMakeTriangleContainingTheVertex();
+			void UnMarkUnSwapTriangle();
+			void SmoothMetric(Real8 raisonmax) ;
+			void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ;
+			void MaxSubDivision(Real8 maxsubdiv);
+			Edge** MakeGeometricalEdgeToEdge();
+			void SetVertexFieldOn();  
+			void SetVertexFieldOnBTh();
+			Int4 SplitInternalEdgeWithBorderVertices();
+			void MakeQuadrangles(double costheta);
+			int  SplitElement(int choice);
+			void MakeQuadTree();
+			void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1);
+			Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; 
+			void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
+			void ReNumberingVertex(Int4 * renu);
+			void SmoothingVertex(int =3,Real8=0.3);
+			Metric MetricAt (const R2 &) const;
+			GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex &  A, Vertex & B,Real8 theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
+			Int4 Number(const Triangle & t) const  { return &t - triangles;}
+			Int4 Number(const Triangle * t) const  { return t - triangles;}
+			Int4 Number(const Vertex & t) const  { return &t - vertices;}
+			Int4 Number(const Vertex * t) const  { return t - vertices;}
+			Int4 Number(const Edge & t) const  { return &t - edges;}
+			Int4 Number(const Edge * t) const  { return t - edges;}
+			Int4 Number2(const Triangle * t) const  {
+				return t - triangles;
+			}
+			Vertex* NearestVertex(Icoor1 i,Icoor1 j) ;
+			Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
+			void ReadMesh(double* index,double* x,double* y,int nods,int nels);
+			void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
+			void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
+			void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
+			void WriteMetric(BamgOpts* bamgopts);
+			void AddMetric(BamgOpts* bamgopts);
+			void BuildMetric0(BamgOpts* bamgopts);
+			void BuildMetric1(BamgOpts* bamgopts);
+			void AddGeometryMetric(BamgOpts* bamgopts);
+			int  isCracked() const {return NbCrackedVertices != 0;}
+			int  Crack();
+			int  UnCrack();
+			void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
+			void GenerateMeshProperties() ;
+			int  CrackMesh();
+
+		private:
+			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);
+
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 3231)
+++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 3232)
@@ -5,10 +5,4 @@
 
 #include "../Mesh2.h"
-#include "../QuadTree.h"
-#include "../SetOfE4.h"
-
-#include "../../shared/shared.h"
-#include "../../include/macros.h"
-#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/Bamgx/objects/Vertex.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Vertex.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/Vertex.h	(revision 3232)
@@ -0,0 +1,62 @@
+#ifndef _VERTEX_H_
+#define _VERTEX_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Direction.h"
+#include "Metric.h"
+
+namespace bamg {
+
+	//classes
+	class Triangle;
+	class Triangles;
+	class VertexOnGeom;
+	class VertexOnEdge;
+
+	class Vertex {
+
+		public:
+			I2 i;  // integer coordinates
+			R2 r;  // real coordinates
+			Metric m;
+			Int4 ReferenceNumber;
+			Direction DirOfSearch;
+			Int1 vint;  // the vertex number in triangle; varies between 0 and 2 in t
+			union {
+				Triangle* t; // one triangle which is containing the vertex
+				Int4 color;
+				Vertex* to;  // use in geometry Vertex to now the Mesh Vertex associed 
+				VertexOnGeom* onGeometry;       // if vint == 8; // set with Triangles::SetVertexFieldOn()
+				Vertex* onBackgroundVertex;     // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
+				VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge
+			};
+
+			//Operators
+			operator  I2() const {return i;}             // Cast operator
+			operator  const R2 & () const {return r;}    // Cast operator
+			operator Metric () const {return m;}         // Cast operator
+			Real8 operator()(R2 x) const { return m(x);} // Get x in the metric m
+
+			//methods (No constructor and no destructors...)
+			Real8 Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,Real8 =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);
+			void  Echo();
+			int   ref() const { return ReferenceNumber;}
+
+			//inline functions
+			inline Int4 Optim(int =1,int =0); 
+			inline void Set(const Vertex & rec,const Triangles &,Triangles &);
+
+
+	};
+
+	//FOR NOW (WARNING)
+	inline Vertex* TheVertex(Vertex * a); // for remove crak in mesh 
+	double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &);
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/VertexOnEdge.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/VertexOnEdge.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/VertexOnEdge.h	(revision 3232)
@@ -0,0 +1,43 @@
+#ifndef _VERTEXONEDGE_H_
+#define _VERTEXONEDGE_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+#include "Edge.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+
+	class VertexOnEdge {
+
+		public:
+			Vertex* v;
+			Edge*   be;
+			Real8 abcisse;
+
+			//Constructors
+			VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {}
+			VertexOnEdge(){}
+
+			//Operators
+			operator Real8 () const { return abcisse;}
+			operator Vertex* () const { return v;}  
+			operator Edge* () const { return be;}  
+			Vertex & operator[](int i) const { return (*be)[i];}
+
+			//Methods
+			void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;}  
+
+			//Inline methods
+			inline void Set(const Triangles &,Int4,Triangles &);  
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h	(revision 3232)
@@ -0,0 +1,54 @@
+#ifndef _VERTEXONGEOM_H_
+#define _VERTEXONGEOM_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+#include "GeometricalVertex.h"
+#include "GeometricalEdge.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+
+	class VertexOnGeom{
+
+		public:
+
+			Vertex* mv;
+			Real8 abscisse;  
+			union{ 
+				GeometricalVertex * gv; // if abscisse <0; 
+				GeometricalEdge * ge;  // if abscisse in [0..1]
+			};
+
+			//Constructors/Destructors
+			VertexOnGeom(): mv(0),abscisse(0){gv=0;} 
+			VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
+			VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;}
+
+			//Operators
+			operator Vertex * () const  {return mv;}
+			operator GeometricalVertex * () const  {return gv;}
+			operator GeometricalEdge * () const  {return ge;}
+			operator const Real8 & () const {return abscisse;}
+
+			//Methods
+			int  OnGeomVertex()const {return this? abscisse <0 :0;}
+			int  OnGeomEdge() const {return this? abscisse >=0 :0;}
+			int  IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
+			void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;}
+
+			//Inline methods
+			inline void Set(const Triangles &,Int4,Triangles &);
+			inline void Set(const VertexOnGeom&,const Triangles &,Triangles &);  
+
+	};
+
+}
+#endif
Index: /issm/trunk/src/c/Bamgx/objects/VertexOnVertex.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/VertexOnVertex.h	(revision 3232)
+++ /issm/trunk/src/c/Bamgx/objects/VertexOnVertex.h	(revision 3232)
@@ -0,0 +1,35 @@
+#ifndef _VERTEXONVERTEX_H_
+#define _VERTEXONVERTEX_H_
+
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
+
+#include "../meshtype.h"
+#include "Vertex.h"
+
+namespace bamg {
+
+	//classes
+	class Triangles;
+
+	class VertexOnVertex {
+
+		public:
+			Vertex* v;
+			Vertex* bv;
+
+			//Constructors
+			VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){}
+			VertexOnVertex() {};
+
+			//Methods
+			void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}
+
+			//Inline methods
+			inline void Set(const Triangles &,Int4,Triangles &);
+	};
+
+}
+#endif