Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2957)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2958)
@@ -144,6 +144,6 @@
 
 		public:
-			Triangle* t; // le triangle 
-			int  a;      // le numero de l arete
+			Triangle* t; //pointer toward the triangle
+			int  a;      //Edge number
 
 			TriangleAdjacent(Triangle  * tt,int  aa): t(tt),a(aa &3) {};
@@ -153,25 +153,25 @@
 			operator Triangle & () const {return *t;}
 			operator int() const {return a;}
-			TriangleAdjacent & operator++() 
-			  {
+			TriangleAdjacent & operator++(){
 				a= NextEdge[a];
-				return *this;}
-				TriangleAdjacent operator--()
-				  { 
-					a= PreviousEdge[a];
-					return *this;}
-					inline  TriangleAdjacent  Adj() const ;
-					int swap();
-					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 ;
+				return *this;
+			}
+			TriangleAdjacent operator--(){ 
+				a= PreviousEdge[a];
+				return *this;
+			}
+			inline  TriangleAdjacent  Adj() const ;
+			int swap();
+			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*/
@@ -298,5 +298,5 @@
 		Vertex* ns [3];    // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
 		Triangle* at [3];  // 3 adjacent triangles (nu)
-		Int1  aa[3];       // les nu des arete dans le triangles (mod 4)
+		Int1  aa[3];       // les nu des arete dans le triangles (mod 3)
 		public: 
 		Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
@@ -354,8 +354,15 @@
 		  }
 
-		void SetAdj2(Int1 a,Triangle *t,Int1 aat)
-		  {  at[a]=t;aa[a]=aat;
-			if(t) {t->at[aat]=this;t->aa[aat]=a;}
-		  }
+		void SetAdj2(Int1 a,Triangle *t,Int1 aat){
+			//at = pointer toward the adjacent triangles of this
+			//aa = position of the edge in the triangle (mod 3)
+			at[a]=t;   //the adjacent triangle to the edge a is t
+			aa[a]=aat; //position of the edge in the adjacent triangle
+			//if t!=NULL add adjacent triangle to t (this)
+			if(t) {
+				t->at[aat]=this;
+				t->aa[aat]=a;
+			}
+		}
 
 		void SetTriangleContainingTheVertex()
@@ -377,5 +384,10 @@
 
 		void SetHidden(int a){
-			register Triangle * t = at[a];
+
+			//Get Adjacent Triangle number a
+			register Triangle* t = at[a];
+			
+			//if it exist
+			//C|=D -> C=(C|D) bitwise inclusive OR
 			if(t) t->aa[aa[a] & 3] |=16;
 			aa[a] |= 16;
@@ -391,4 +403,5 @@
 
 		void SetLocked(int a){
+			//mark the edge as on Boundary
 			register Triangle * t = at[a];
 			t->aa[aa[a] & 3] |=4;
@@ -1071,14 +1084,15 @@
 	  }	       
 
-	inline  void  TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l  )
-	  { // set du triangle adjacent 
+	inline  void  TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l  ){
+		//set Adjacent Triangle of a triangle
 		if(t) {
 			t->at[a]=ta.t;
-			t->aa[a]=ta.a|l;}
-			if(ta.t) {
-				ta.t->at[ta.a] = t ;
-				ta.t->aa[ta.a] = a| l ;
-			}
-	  }
+			t->aa[a]=ta.a|l;
+		}
+		if(ta.t) {
+			ta.t->at[ta.a] = t ;
+			ta.t->aa[ta.a] = a| l ;
+		}
+	}
 
 
Index: /issm/trunk/src/c/Bamgx/SetOfE4.h
===================================================================
--- /issm/trunk/src/c/Bamgx/SetOfE4.h	(revision 2957)
+++ /issm/trunk/src/c/Bamgx/SetOfE4.h	(revision 2958)
@@ -14,15 +14,15 @@
 	class SetOfEdges4 {
 		Int4 nx,nbax,NbOfEdges;
-		Int4 * tete; 
-		Int4Edge * Edges;
+		Int4* head; 
+		Int4Edge* Edges;
 
 		public:
 		SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 
-		~SetOfEdges4() {delete [] tete; delete [] Edges;}
+		~SetOfEdges4() {delete [] head; delete [] Edges;}
 		Int4 add (Int4 ii,Int4 jj);
-		Int4 addtrie (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
+		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 findtrie (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj)  : find (jj,ii) ;}
+		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;}
Index: /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp	(revision 2957)
+++ /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp	(revision 2958)
@@ -14,12 +14,18 @@
 	/*FUNCTION  SetOfEdges4::SetOfEdges4(Int4 mmx,Int4 nnx){{{1*/
 	SetOfEdges4::SetOfEdges4(Int4 mmx,Int4 nnx){
-		nx=nnx;
-		nbax=mmx;
-		NbOfEdges = 0;
-		tete= new Int4 [nx];
-		Int4 i=nx;
-		while (i--)
-		 tete[i] = -1;// vide 
-		Edges =new Int4Edge[nbax];
+
+		/*Intermediary*/
+		int i;
+
+		//initialize fields
+		nx   =nnx;   //number of vertices
+		nbax =mmx; // 3 * number of triangles
+		NbOfEdges=0;
+		head = new Int4 [nx];
+		Edges= new Int4Edge[nbax];
+
+		//initialize head (-1 everywhere)
+		i=nx;
+		while (i--) head[i]=-1;
 	}
 	/*}}}1*/
@@ -28,37 +34,67 @@
 	/*FUNCTION  SetOfEdges4::find {{{1*/
 	Int4 SetOfEdges4::find(Int4 ii,Int4 jj) { 
-		if (tete == 0 ) {
-			throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::find no more tete de liste (?)"));
+
+		/*Intermediary*/
+		int n;
+
+		//check that head is not empty
+		if (head==NULL) {
+			throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)"));
 		}
-		Int4 n = tete[ Abs( ii ) % nx ];
 
-		while (n >= 0) 
-		 if (ii == Edges[n].i && jj == Edges[n].j) return n;
-		 else n = Edges[n].next;
-		return -1; //do not exist, return -1
+		//get n from h (usually h=ii)
+		n=head[Abs(ii)%nx];
+
+		//go through the existing edges that holds h (=ii) and return position in Edge
+		while (n >= 0){
+
+			//if the edge ii jj is already in Edges, return n
+			if (ii == Edges[n].i && jj == Edges[n].j) return n;
+
+			//else go to next edge that holds ii
+			else n = Edges[n].next;
+		}
+
+		//if we reach this point, the edge does not exist return -1
+		return -1;
 	}
 	/*}}}1*/
 	/*FUNCTION  SetOfEdges4::add{{{1*/
 	Int4 SetOfEdges4::add(Int4 ii,Int4 jj) {
-		if (tete == 0 ) {
-			throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no more tete de liste (?)"));
+
+		/*Intermediary*/
+		int h,n;
+
+		//check that head is not empty
+		if (head==NULL) {
+			throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)"));
 		}
-		Int4 h;
-		Int4 n = tete[ h = Abs( ii ) % nx ];
-		while (n >= 0) 
-		 if (ii == Edges[n].i && jj == Edges[n].j)
-		  return n;
-		 else n = Edges[n].next;
+
+		//get n from h (usually h=ii)
+		n=head[h=Abs(ii)%nx];
+
+		//go through the existing edges that holds h (=ii) and check that 
+		//the edge ii jj is not already in Edge
+		while (n >= 0){
+
+			//if the edge ii jj is already in Edges, return n
+			if (ii == Edges[n].i && jj == Edges[n].j) return n;
+
+			//else go to next edge that holds ii
+			else n = Edges[n].next;
+		}
+
+		//check that nbax <=NbOfEdges
 		if (nbax <=NbOfEdges ) {
 			throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add overflow: NbOfEdges=%i > nbax=%i",NbOfEdges,nbax));
 		}
 
+		//update chain
 		Edges[NbOfEdges].i=ii;
 		Edges[NbOfEdges].j=jj;
-		Edges[NbOfEdges].next= tete[h];
-		tete[h] = NbOfEdges;
+		Edges[NbOfEdges].next= head[h];
+		head[h] = NbOfEdges;
 		return NbOfEdges ++;
 	}
 	/*}}}1*/
-
 }
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2957)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2958)
@@ -881,5 +881,5 @@
 
 		/*Intermediary*/
-		int i,j,k;
+		int i,j,k,kk,it,jt;
 
 		/*Recover options*/
@@ -899,34 +899,52 @@
 		if (cutoffradian>=0) Gh.MaximalAngleOfCorner = cutoffradian;
 
-		//Construction of the edges 
+		/*Construction of the edges*/
+
+		//initialize st and edge4
 		SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
 		Int4*        st   = new Int4[nbt*3];
-		for (i=0;i<nbt*3;i++)
-		 st[i]=-1;
-		Int4 kk =0;
-
-		Int4 nbeold = nbe;
+
+		//initialize st as -1 (chaining algorithm)
+		for (i=0;i<nbt*3;i++) st[i]=-1;
+
+		//build edge4 (chain)
 		for (i=0;i<nbe;i++){
-			edge4->addtrie(Number(edges[i][0]),Number(edges[i][1]));
-		}
+			edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]));
+		}
+		//check that there is no double edge
 		if (nbe !=  edge4->nb()){ 
 			throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i, nbe4=%i",nbe,edge4->nb())); 
 		}
+		//keep nbe in nbeold
+		Int4 nbeold = nbe;
+
+		//Go through the triangles and ass the edges in edge4 if they are not there yet
 		for (i=0;i<nbt;i++){
+			//3 edges per triangle
 			for  (j=0;j<3;j++) {
-				Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
+				//Add Edge to edge4 (k=numberofedges in edge4)
+				Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
 				Int4 invisible = triangles[i].Hidden(j);
+
+				//if st[k] has not been changed yet, add 3*i+j (= vertex position in the index)
 				if(st[k]==-1) st[k]=3*i+j;
+
+				//else st[k]>=0 -> the edge already exist, check
 				else if(st[k]>=0) {
+					//check that it is not an edge on boundary (should not alrady exist)
 					if (triangles[i].TriangleAdj(j) || triangles[st[k]/3].TriangleAdj((int) (st[k]%3))){
-						throw ErrorException(__FUNCT__,exprintf("triangles[i].TriangleAdj(j) || triangles[st[k]/3].TriangleAdj((int) (st[k]%3))"));
+						throw ErrorException(__FUNCT__,exprintf("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)"));
 					}
+					//OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list
 					triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
 					if (invisible)  triangles[i].SetHidden(j);
+					// if k < nbe mark the adge as on Boundary (Locked)
 					if (k<nbe) {
 						triangles[i].SetLocked(j);
 					}
+					//set st[k] as negative so that the edge should not be called again
 					st[k]=-2-st[k]; 
 				}
+				//else (see 3 lines above), the edge has been called more than twice: return error
 				else {
 					printf("The edge (%i,%i) belongs to more than 2 triangles (%i)\n",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]),k);
@@ -938,10 +956,12 @@
 			}
 		}
+
+		//delete edge4
 		Int4 nbedges = edge4->nb(); // the total number of edges 
-		delete edge4;
-		edge4 =0;
-
+		delete edge4; edge4=NULL;
+
+		//display info
 		if(verbosity>5) {
-			printf("         info on Meshs:\n");
+			printf("         info on Mesh:\n");
 			printf("            - number of vertices    = %i \n",nbv); 
 			printf("            - number of triangles   = %i \n",nbt); 
@@ -950,28 +970,33 @@
 			printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv); 
 		}
-		// check the consistant of edge[].adj and the geometrical required  vertex
-		k=0;
-		kk=0;
-		Int4 it;
-
-		for (i=0;i<nbedges;i++)
-		 if (st[i] <-1) {// edge internal
-			 it =  (-2-st[i])/3;
-			 j  =  (int) ((-2-st[i])%3);
-			 Triangle & tt = * triangles[it].TriangleAdj(j);
-			 if (triangles[it].color != tt.color|| i < nbeold) k++;
-		 }
-		 else if (st[i] >=0) // edge alone 
-		  kk++;
+
+		// check consistency of edge[].adj and geometrical required  vertices
+		k=0; kk=0;
+		for (i=0;i<nbedges;i++){
+			//internal edge
+			if (st[i] <-1) {
+				//get triangle number back
+				it =  (-2-st[i])/3;
+				//get edge position back
+				j  =  (int) ((-2-st[i])%3);
+				Triangle &tt=*triangles[it].TriangleAdj(j);
+				if (triangles[it].color != tt.color|| i < nbeold) k++;
+			}
+			//boundary edge (alone)
+			else if (st[i] >=0) 
+			 kk++;
+		}
+
+		/*Constructions of edges*/
 
 		k += kk;
 		kk=0;
 		if (k) {
-			// construction of the edges 
 			nbe = k;
-			Edge * edgessave = edges;
+			Edge* edgessave=edges;
 			edges = new Edge[nbe];
 			k =0;
-			// construction of the edges 
+
+			//display info
 			if(verbosity>4) printf("   Construction of the edges %i\n",nbe);
 
@@ -979,69 +1004,95 @@
 				Int4  add= -1;
 
-				if (st[i] <-1) // edge internal
-				  { 
+				//internal edge (belongs to two triangles)
+				if (st[i] <-1){ 
 					it =  (-2-st[i])/3;
 					j  =  (int) ((-2-st[i])%3);
 					Triangle & tt = * triangles[it].TriangleAdj(j);
-					if (triangles[it].color !=  tt.color || i < nbeold) // Modif FH 06122055
+					if (triangles[it].color !=  tt.color || i < nbeold)
 					 add=k++;
-				  }
-				else if (st[i] >=0) // edge alone 
-				  {
+				}
+				//boundary edge
+				else if (st[i] >=0){
 					it = st[i]/3;
 					j  = (int) (st[i]%3);
 					add=k++;
-				  }
-
-				if (add>=0 && add < nbe)
-				  {
-
+				}
+				if (add>=0 && add < nbe){
 					edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]];
 					edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]];
 					edges[add].onGeometry=0; 
-					if (i<nbeold) // in file edge // Modif FH 06122055 
-					  {
-						edges[add].ref = edgessave[i].ref; 		      
-						edges[add].onGeometry = edgessave[i].onGeometry; //  HACK pour recuperer les aretes requise midf FH avril 2006 ???? 
-					  }
+					//if already existed
+					if (i<nbeold){
+						edges[add].ref=edgessave[i].ref; 		      
+						edges[add].onGeometry=edgessave[i].onGeometry; //  HACK to get required edges
+					}
 					else
-					 edges[add].ref = Min(edges[add].v[0]->ref(),edges[add].v[1]->ref()); // no a good choice
+					 edges[add].ref=Min(edges[add].v[0]->ref(),edges[add].v[1]->ref());
 				  }
 			}
+
+			//check that we have been through all edges
 			if (k!=nbe){
-				throw ErrorException(__FUNCT__,exprintf("k!=nbe"));
-			}
+				throw ErrorException(__FUNCT__,exprintf("problem in edge construction process: k!=nbe (should not happen)"));
+			}
+			//delete edgessave
 			if (edgessave) delete [] edgessave;
 		}
 
-		// construction of edges[].adj 
-		for (i=0;i<nbv;i++) 
-		 vertices[i].color =0;
-		for (i=0;i<nbe;i++)
-		 for (j=0;j<2;j++) 
-		  edges[i].v[j]->color++;
-
-		for (i=0;i<nbv;i++) 
-		 vertices[i].color = (vertices[i].color ==2) ? -1 : -2;
-		for (i=0;i<nbe;i++)
-		 for (j=0;j<2;j++) 
-			{ 
-			 Vertex *v=edges[i].v[j];
-			 Int4 i0=v->color,j0;
-			 if(i0<0)
-			  edges[i ].adj[ j ]=0;  // Add FH Jan 2008   
-			 if(i0==-1)
-			  v->color=i*2+j;
-			 else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
-				 j0 =  i0%2;
-				 i0 =  i0/2;
-				 if (v!=edges[i0 ].v[j0]){
-					 throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]"));
-				 }
-				 edges[i ].adj[ j ] =edges +i0;
-				 edges[i0].adj[ j0] =edges +i ;
-				 v->color = -3;}
-			}
-		// now reconstruct the sub domain info 
+		/*Color the vertices*/
+
+		//initialize color of all vertices as 0
+		for (i=0;i<nbv;i++) vertices[i].color =0;
+
+		//go through the edges and add a color to corresponding vertices
+		//(A vertex in 4 edges will have a color 4)
+		for (i=0;i<nbe;i++){
+		 for (j=0;j<2;j++) edges[i].v[j]->color++;
+		}
+
+		//change the color: if a vertex belongs to 2 edges -1, else -2
+		for (i=0;i<nbv;i++) {
+			vertices[i].color=(vertices[i].color ==2)? -1 : -2;
+		}
+
+		/*Build edges[].adj*/
+
+		for (i=0;i<nbe;i++){
+			for (j=0;j<2;j++){ 
+				//get current vertex
+				Vertex* v=edges[i].v[j];
+				//get vertex color (i0)
+				Int4 i0=v->color,j0;
+
+				//if color<0 (first time), no adjacent edge
+				if(i0<0) edges[i].adj[j]=NULL;
+
+				//if color=-1 (corner),change the vertex color as 2*i+j (position of the vertex in edges)
+				if(i0==-1) v->color=i*2+j;
+
+				//if color>=0 (i and i0 edge are adjacent by the vertex v)
+				else if (i0>=0) {
+					//get position of v in edges back
+					j0 =  i0%2; //column in edges
+					i0 =  i0/2; //line in edges
+
+					//check that we have the correct vertex
+					if (v!=edges[i0 ].v[j0]){
+						throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]: this should not happen as the vertex belongs to this edge"));
+					}
+
+					//Add adjacence
+					edges[i ].adj[j ]=edges +i0;
+					edges[i0].adj[j0]=edges +i ;
+
+					//change color to -3
+					v->color = -3;
+				}
+			}
+		}
+
+		/*Reconstruct subdomains info*/
+
+		//check that NbSubDomains is empty
 		if (NbSubDomains){
 			throw ErrorException(__FUNCT__,exprintf("NbSubDomains should be 0"));
@@ -1049,84 +1100,89 @@
 		NbSubDomains=0;
 
-		  { 
-			Int4 it;
-			// find all the sub domain
-			Int4 *colorT = new Int4[nbt];
-			Triangle *tt,*t;
-			Int4 k;
-			for ( it=0;it<nbt;it++)
-			 colorT[it]=-1;
-			for (it=0;it<nbt;it++)
-			  {
-				if (colorT[it]<0) 
-				  {
-					colorT[it]=NbSubDomains;
-					Int4 level =1,j,jt,kolor=triangles[it].color;
-					st[0]=it; // stack 
-					st[1]=0;
-					k=1;
-					while (level>0)
-					 if( ( j=st[level]++) <3)
-						{ 
-						 t = &triangles[st[level-1]];
-						 tt=t->TriangleAdj((int)j);
-
-						 if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor))
-							{
-							 colorT[jt]=NbSubDomains;
-							 st[++level]=jt;
-							 st[++level]=0;
-							 k++;
-							}
+		//color the subdomains
+		Int4* colorT= new Int4[nbt];
+		Triangle *tt,*t;
+
+		//initialize the color of each triangle as -1
+		for (it=0;it<nbt;it++) colorT[it]=-1;
+
+		//loop over the triangles
+		for (it=0;it<nbt;it++){
+
+			//if the triangle has not been colored yet:
+			if (colorT[it]<0){
+
+				//color = number of subdomains
+				colorT[it]=NbSubDomains;
+
+				//color all the adjacent triangles of T that share a non marked edge
+				int level =1;
+				int kolor=triangles[it].color;
+				st[0]=it; // stack 
+				st[1]=0;
+				k=1;
+				while (level>0){
+					if( (j=st[level]++)<3 ){ 
+						t = &triangles[st[level-1]];
+						tt=t->TriangleAdj((int)j);
+
+						//color the adjacent triangle
+						if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) {
+							colorT[jt]=NbSubDomains;
+							st[++level]=jt;
+							st[++level]=0;
+							k++;
 						}
-					 else 
-					  level-=2;
-					NbSubDomains++;
-				  }
-			  }
-			if (verbosity> 3) printf("      The Number of sub domain = %i\n",NbSubDomains); 
-
-			Int4 isd;
-			subdomains = new SubDomain[NbSubDomains];
-			for (isd=0;isd<NbSubDomains;isd++)
-			  {
-				subdomains[isd].head =0;
-			  }
-			k=0;
-			for (it=0;it<nbt;it++)
-			 for (int j=0;j<3;j++)
-				{
-				 tt=triangles[it].TriangleAdj(j);
-				 if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head)
-					{
-					 subdomains[isd].head = triangles+it;
-					 subdomains[isd].ref =  triangles[it].color;
-					 subdomains[isd].sens = j; // hack
-					 subdomains[isd].edge = 0;
-					 k++;
 					}
-				}  
-			if (k!= NbSubDomains){
-				throw ErrorException(__FUNCT__,exprintf("k!= NbSubDomains"));
-			}
-
-			delete [] colorT;
-
-
-		  }      
+					else level-=2;
+				}
+				NbSubDomains++;
+			}
+		}
+		if (verbosity> 3) printf("      The Number of sub domain = %i\n",NbSubDomains); 
+
+		//build subdomains
+		Int4 isd;
+		subdomains = new SubDomain[NbSubDomains];
+
+		//initialize subdomains[isd].head as 0
+		for (isd=0;isd<NbSubDomains;isd++) subdomains[isd].head =0;
+		  
+		k=0;
+		for (it=0;it<nbt;it++){
+			for (int j=0;j<3;j++){
+				tt=triangles[it].TriangleAdj(j);
+				if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head){
+					subdomains[isd].head = triangles+it;
+					subdomains[isd].ref =  triangles[it].color;
+					subdomains[isd].sens = j; // hack
+					subdomains[isd].edge = 0;
+					k++;
+				}
+			}
+		}
+		//check that we have been through all subdomains
+		if (k!= NbSubDomains){
+			throw ErrorException(__FUNCT__,exprintf("k!= NbSubDomains"));
+		}
+		//delete colorT and st
+		delete [] colorT;
 		delete [] st;
-		// now make the geometry
-		// 1 compress the vertices 
-		Int4 * colorV = new Int4[nbv];
-		for (i=0;i<nbv;i++) 
-		 colorV[i]=-1;
-		for (i=0;i<nbe;i++)
-		 for ( j=0;j<2;j++)
-		  colorV[Number(edges[i][j])]=0;
+
+		/*Reconstruct Geometry Gh*/
+
+		//build colorV -1 for all vertex and 0 for the vertices belonging to edges
+		Int4* colorV = new Int4[nbv];
+		for (i=0;i<nbv;i++) colorV[i]=-1;
+		for (i=0;i<nbe;i++){
+		 for ( j=0;j<2;j++) colorV[Number(edges[i][j])]=0;
+		}
+		//number the vertices nelonging to edges
 		k=0;
-		for (i=0;i<nbv;i++) 
-		 if(!colorV[i])
-		  colorV[i]=k++;
-
+		for (i=0;i<nbv;i++){
+		 if(!colorV[i]) colorV[i]=k++;
+		}
+
+		//Build Gh
 		Gh.nbv=k;
 		Gh.nbe = nbe;
@@ -1140,23 +1196,16 @@
 		NbVerticesOnGeomEdge =0;
 		VerticesOnGeomEdge =0;
-		  {
-			Int4 j;
-			for (i=0;i<nbv;i++) 
-			 if((j=colorV[i])>=0)
-				{
-
-				 Vertex & v = Gh.vertices[j];
-				 v = vertices[i];
-				 v.color =0;
-				 VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]);
-				}
-
-		  }
-		edge4= new SetOfEdges4(nbe,nbv);  
-
-		Real4 * len = new Real4[Gh.nbv];
-		for(i=0;i<Gh.nbv;i++)
-		 len[i]=0;
-
+
+		//Build VertexOnGeom
+		for (i=0;i<nbv;i++){
+			if((j=colorV[i])>=0){
+				Vertex & v = Gh.vertices[j];
+				v = vertices[i];
+				v.color =0;
+				VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]);
+			}
+		}
+
+		//Buid pmin and pmax of Gh (extrema coordinates)
 		Gh.pmin =  Gh.vertices[0].r;
 		Gh.pmax =  Gh.vertices[0].r;
@@ -1168,36 +1217,48 @@
 			Gh.pmax.y = Max(Gh.pmax.y,Gh.vertices[i].r.y);
 		}
-
 		R2 DD05 = (Gh.pmax-Gh.pmin)*0.05;
 		Gh.pmin -=  DD05;
 		Gh.pmax +=  DD05;
 
+		//Build Gh.coefIcoor
 		Gh.coefIcoor= (MaxICoor)/(Max(Gh.pmax.x-Gh.pmin.x,Gh.pmax.y-Gh.pmin.y));
 		if (Gh.coefIcoor<=0){
-			throw ErrorException(__FUNCT__,exprintf("Gh.coefIcoor<=0"));
-		}
-
+			throw ErrorException(__FUNCT__,exprintf("Gh.coefIcoor<=0 in infered Geometry (this should not happen)"));
+		}
+
+		/*Build Gh.edges*/
+
+		//initialize len as 0
+		Real4 * len = new Real4[Gh.nbv];
+		for(i=0;i<Gh.nbv;i++) len[i]=0;
+
+		//initialize edge4 again
+		edge4= new SetOfEdges4(nbe,nbv);  
 		Real8 hmin = HUGE_VAL;
 		int kreq=0;
-		for (i=0;i<nbe;i++)
-		  {
+		for (i=0;i<nbe;i++){
+
 			Int4 i0 = Number(edges[i][0]);
 			Int4 i1 = Number(edges[i][1]);
-			Int4 j0 =	 colorV[i0];
-			Int4 j1 =  colorV[i1];
+			Int4 j0 = colorV[i0];
+			Int4 j1 = colorV[i1];
 
 			Gh.edges[i].v[0] = Gh.vertices +  j0;
 			Gh.edges[i].v[1] = Gh.vertices +  j1;
+
 			Gh.edges[i].flag = 0;
+
 			Gh.edges[i].tg[0]=R2();
 			Gh.edges[i].tg[1]=R2();
-			bool requis= edges[i].onGeometry; 
-			if(requis) kreq++;
+
+			bool required= edges[i].onGeometry; 
+			if(required) kreq++;
 			edges[i].onGeometry =  Gh.edges + i;
-			if(requis)  {  // correction fevr 2009 JYU ...
+			if(required){
 				Gh.edges[i].v[0]->SetRequired();
 				Gh.edges[i].v[1]->SetRequired();
-				Gh.edges[i].SetRequired(); // fin modif ... 
-			}
+				Gh.edges[i].SetRequired();
+			}
+
 			R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r;
 			Real8 l12=Norme2(x12);        
@@ -1212,27 +1273,28 @@
 			Gh.edges[i].ref  = edges[i].ref;
 
-			k = edge4->addtrie(i0,i1);
+			k = edge4->SortAndAdd(i0,i1);
 			if (k != i){
-				throw ErrorException(__FUNCT__,exprintf("k != i"));
-			}
-
-		  }
-
-
-		for (i=0;i<Gh.nbv;i++) 
+				throw ErrorException(__FUNCT__,exprintf("problem in Edge4 construction: k != i"));
+			}
+		}
+
+		//Build metric for all vertices of Gh
+		for (i=0;i<Gh.nbv;i++){
 		 if (Gh.vertices[i].color > 0) 
 		  Gh.vertices[i].m=  Metric(len[i] /(Real4) Gh.vertices[i].color);
 		 else 
 		  Gh.vertices[i].m=  Metric(hmin);
+		}
+		//delete len
 		delete [] len;
-		for (i=0;i<NbSubDomains;i++)
-		  {
-			Int4 it = Number(subdomains[i].head);
-			int j = subdomains[i].sens;
+
+		//Build Gh.subdomains
+		for (i=0;i<NbSubDomains;i++){
+			it = Number(subdomains[i].head);
+			j = subdomains[i].sens;
 			Int4 i0 = Number(triangles[it][VerticesOfTriangularEdge[j][0]]);
 			Int4 i1 = Number(triangles[it][VerticesOfTriangularEdge[j][1]]);
-			k = edge4->findtrie(i0,i1);
-			if(k>=0)
-			  {
+			k = edge4->SortAndFind(i0,i1);
+			if(k>=0){
 				subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
 				subdomains[i].edge = edges+k;
@@ -1240,5 +1302,5 @@
 				Gh.subdomains[i].sens  =  subdomains[i].sens;
 				Gh.subdomains[i].ref =  subdomains[i].ref;
-			  }
+			}
 			else
 			 throw ErrorException(__FUNCT__,exprintf("%i should be >=0"));
@@ -1247,8 +1309,11 @@
 		delete edge4;
 		delete [] colorV;
-		//  -- unset adj
-		for (i=0;i<nbt;i++)
-		 for ( j=0;j<3;j++)
-		  triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j));
+
+		//unset adj
+		for (i=0;i<nbt;i++){
+			for ( j=0;j<3;j++){
+				triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j));
+			}
+		}
 
 	}
@@ -2006,5 +2071,5 @@
 			Int4 kk =0;
 			for (i=0;i<nbe;i++){
-				kk=kk+(i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
+				kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
 			}
 			if (kk != nbe) { 
@@ -2015,5 +2080,5 @@
 			for (i=0;i<nbt;i++){
 				for (int j=0;j<3;j++) {
-					Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
+					Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
 								Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
 					Int4 invisible = triangles[i].Hidden(j);
@@ -2177,5 +2242,5 @@
 							Vertex *v0= ta.EdgeVertex(0);
 							Vertex *v1= ta.EdgeVertex(1);
-							Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
+							Int4 k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
 							if (st[k]<0){
 								throw ErrorException(__FUNCT__,exprintf("st[k]<0"));
@@ -4504,5 +4569,5 @@
 
 
-			kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
+			kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
 			if (ong) // a geometrical edges 
 			  { 
@@ -4595,5 +4660,5 @@
 				const Vertex & v0 = t[VerticesOfTriangularEdge[j][0]];
 				const Vertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				Int4  ke =edge4->findtrie(Number(v0),Number(v1));
+				Int4  ke =edge4->SortAndFind(Number(v0),Number(v1));
 				if (ke>0) 
 				  {
@@ -4646,5 +4711,5 @@
 				if ( kedge[3*i+j] < 0) 
 				  {
-					Int4  ke =edge4->findtrie(Number(v0),Number(v1));
+					Int4  ke =edge4->SortAndFind(Number(v0),Number(v1));
 					if (ke<0) // new 
 					  {