Index: /issm/trunk/src/c/Bamgx/BamgObjects.h
===================================================================
--- /issm/trunk/src/c/Bamgx/BamgObjects.h	(revision 3237)
+++ /issm/trunk/src/c/Bamgx/BamgObjects.h	(revision 3238)
@@ -34,68 +34,4 @@
 	
 	/*INLINE functions{{{1*/
-	// to sort in descending order
-	template<class T>  inline void  HeapSort(T *c,long n){
-		/*Intermediary*/
-		int l,j,r,i;
-		T   crit;
-		c--;                    //the array must starts at 1 and not 0 
-		if(n<=1) return;        //return if size <=1
-		l=n/2+1;                //initialize l and r
-		r=n;
-		for(;;){
-			if(l<=1){
-				crit  =c[r];
-				c[r--]=c[1];
-				if (r==1){c[1]=crit; return;}
-			}
-			else  crit = c[--l]; 
-			j=l;
-			for(;;){
-				i=j;
-				j=2*j;
-				if  (j>r) {c[i]=crit;break;}
-				if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
-				if (crit < c[j]) c[i]=c[j];       //c[j]  > crit -> stock this large value in i(<j)
-				else{c[i]=crit;break;}            //c[j]  < crit -> stock crit in i (<j)
-			}
-		}
-	}
-	// to sort in descending order and return ordering
-	template<class T>  inline void  HeapSort(int** porder,T* c,int n){
-		/*Intermediary*/
-		int  l,j,r,i;
-		T    crit;
-		int  pos;
-		int* order=NULL;
-		order=(int*)xmalloc(n*sizeof(int));
-		for(i=0;i<n;i++) order[i]=i+1;
-		c--;                    //the array must starts at 1 and not 0 
-		order--;
-		if(n<=1) return;        //return if size <=1
-		l=n/2+1;                //initialize l and r
-		r=n;
-		for(;;){
-			if(l<=1){
-				crit  =c[r]; pos=order[r];
-				c[r--]=c[1]; order[r+1]=order[1];
-				if (r==1){
-					c[1]=crit; order[1]=pos;
-					order++;
-					*porder=order;
-					return;
-				}
-			}
-			else  {crit=c[--l]; pos=order[l];}
-			j=l;
-			for(;;){
-				i=j;
-				j=2*j;
-				if  (j>r) {c[i]=crit;order[i]=pos;break;}
-				if ((j<r) && (c[j] < c[j+1]))j++;
-				if (crit < c[j]) {c[i]=c[j];order[i]=order[j];} 
-				else{c[i]=crit;order[i]=pos;break;}
-			}
-		}
-	}
 	inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){
 		return A[0]*( B[1]*C[2]-B[2]*C[1])
Index: /issm/trunk/src/c/Bamgx/meshtype.h
===================================================================
--- /issm/trunk/src/c/Bamgx/meshtype.h	(revision 3237)
+++ /issm/trunk/src/c/Bamgx/meshtype.h	(revision 3238)
@@ -58,4 +58,66 @@
 	template<class T> inline T Max3 (const T &a,const T & b,const T & c){return Max(Max(a,b),c);}
 	template<class T> inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);}
+	template<class T> inline void  HeapSort(T *c,long n){
+		/*Intermediary*/
+		int l,j,r,i;
+		T   crit;
+		c--;                    //the array must starts at 1 and not 0 
+		if(n<=1) return;        //return if size <=1
+		l=n/2+1;                //initialize l and r
+		r=n;
+		for(;;){
+			if(l<=1){
+				crit  =c[r];
+				c[r--]=c[1];
+				if (r==1){c[1]=crit; return;}
+			}
+			else  crit = c[--l]; 
+			j=l;
+			for(;;){
+				i=j;
+				j=2*j;
+				if  (j>r) {c[i]=crit;break;}
+				if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
+				if (crit < c[j]) c[i]=c[j];       //c[j]  > crit -> stock this large value in i(<j)
+				else{c[i]=crit;break;}            //c[j]  < crit -> stock crit in i (<j)
+			}
+		}
+	}
+	template<class T> inline void  HeapSort(int** porder,T* c,int n){
+		/*Intermediary*/
+		int  l,j,r,i;
+		T    crit;
+		int  pos;
+		int* order=NULL;
+		order=(int*)xmalloc(n*sizeof(int));
+		for(i=0;i<n;i++) order[i]=i+1;
+		c--;                    //the array must starts at 1 and not 0 
+		order--;
+		if(n<=1) return;        //return if size <=1
+		l=n/2+1;                //initialize l and r
+		r=n;
+		for(;;){
+			if(l<=1){
+				crit  =c[r]; pos=order[r];
+				c[r--]=c[1]; order[r+1]=order[1];
+				if (r==1){
+					c[1]=crit; order[1]=pos;
+					order++;
+					*porder=order;
+					return;
+				}
+			}
+			else  {crit=c[--l]; pos=order[l];}
+			j=l;
+			for(;;){
+				i=j;
+				j=2*j;
+				if  (j>r) {c[i]=crit;order[i]=pos;break;}
+				if ((j<r) && (c[j] < c[j+1]))j++;
+				if (crit < c[j]) {c[i]=c[j];order[i]=order[j];} 
+				else{c[i]=crit;order[i]=pos;break;}
+			}
+		}
+	}
 
 	//Inline functions
@@ -77,5 +139,4 @@
 	}
 
-
 }
 #endif
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3237)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3238)
@@ -5359,463 +5359,3 @@
 /*}}}1*/
 
-	/*Intermediary*/
-	/*FUNCTION swap{{{1*/
-	void  swap(Triangle *t1,Int1 a1, Triangle *t2,Int1 a2, Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2){ 
-		// --------------------------------------------------------------
-		// Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles
-		//                               
-		//               sb                     sb    
-		//             / | \                   /   \                      !
-		//         as1/  |  \                 /a2   \                     !
-		//           /   |   \               /    t2 \                    !
-		//       s1 /t1  | t2 \s2  -->   s1 /___as2___\s2                 !
-		//          \  a1|a2  /             \   as1   /  
-		//           \   |   /               \ t1    /   
-		//            \  |  / as2             \   a1/    
-		//             \ | /                   \   /     
-		//              sa                       sa   
-		//  -------------------------------------------------------------
-		int as1 = NextEdge[a1];
-		int as2 = NextEdge[a2];
-		int ap1 = PreviousEdge[a1];
-		int ap2 = PreviousEdge[a2];
-		(*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb
-		(*t2)(VerticesOfTriangularEdge[a2][1]) = s1  ; // avant sa
-		// mise a jour des 2 adjacences externes 
-		TriangleAdjacent taas1 = t1->Adj(as1),
-							  taas2 = t2->Adj(as2),
-							  tas1(t1,as1), tas2(t2,as2),
-							  ta1(t1,a1),ta2(t2,a2);
-		// externe haut gauche
-		taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap());
-		// externe bas droite
-		taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap());
-		// remove the Mark  UnMarkSwap 
-		t1->SetUnMarkUnSwap(ap1);
-		t2->SetUnMarkUnSwap(ap2);
-		// interne 
-		tas1.SetAdj2(tas2);
-
-		t1->det = det1;
-		t2->det = det2;
-
-		t1->SetTriangleContainingTheVertex();
-		t2->SetTriangleContainingTheVertex();
-	} // end swap 
-	/*}}}1*/
-	/*FUNCTION FindTriangle{{{1*/
-	Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside){
-		I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); 
-		Icoor2 dete[3];
-		Triangle & tb = *Th.FindTriangleContaining(I,dete);
-
-		if  (tb.link) 
-		  { // internal point in a true triangles
-			a[0]= (Real8) dete[0]/ tb.det;
-			a[1]= (Real8) dete[1] / tb.det;
-			a[2] = (Real8) dete[2] / tb.det;
-			inside = 1;	 
-			return Th.Number(tb);
-		  } 
-		else 
-		  {
-			inside = 0; 
-			double aa,bb;
-			TriangleAdjacent  ta=CloseBoundaryEdgeV2(I,&tb,aa,bb);	 
-			int k = ta;
-			Triangle * tc = ta;
-			if (!tc->link) 
-			  { ta = ta.Adj();
-				tc=ta;
-				k = ta;
-				Exchange(aa,bb);
-				if (!tc->link){
-					throw ErrorException(__FUNCT__,exprintf("!tc->link"));
-				}
-			  }
-			a[VerticesOfTriangularEdge[k][0]] = aa;
-			a[VerticesOfTriangularEdge[k][1]] = bb;
-			a[OppositeVertex[k]] = 1- aa -bb;
-			return Th.Number(tc);
-		  }
-	}
-	/*}}}1*/
-	/*FUNCTION CloseBoundaryEdge{{{1*/
-	TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) {
-		int k=(*t)(0) ?  ((  (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1  )) : 0;
-		int dir=0;
-		if (k<0){
-			throw ErrorException(__FUNCT__,exprintf("k<0"));
-		}
-		int kkk=0;  
-		Icoor2 IJ_IA,IJ_AJ;
-		TriangleAdjacent edge(t,OppositeEdge[k]);          
-		for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) {  
-			kkk++;
-			if (kkk>=1000){
-				throw ErrorException(__FUNCT__,exprintf("kkk>=1000"));
-			}
-			Vertex  &vI =  *edge.EdgeVertex(0);
-			Vertex  &vJ =  *edge.EdgeVertex(1);
-			I2 I=vI, J=vJ, IJ= J-I;
-			IJ_IA = (IJ ,(A-I));
-			if (IJ_IA<0) {
-				if (dir>0) {a=1;b=0;return edge;}// change of signe => I
-				else {dir=-1;
-					continue;}};// go in direction i 
-					IJ_AJ = (IJ ,(J-A));
-					if (IJ_AJ<0) {
-						if(dir<0)  {a=0;b=1;return edge;}            
-						else {dir = 1;
-							continue;}}// go in direction j
-							double IJ2 = IJ_IA + IJ_AJ;
-							if (IJ2==0){
-								throw ErrorException(__FUNCT__,exprintf("IJ2==0"));
-							}
-							a= IJ_AJ/IJ2;
-							b= IJ_IA/IJ2;
-							return edge;
-		  } 
-	}
-	/*}}}1*/
-	/*FUNCTION CloseBoundaryEdgeV2{{{1*/
-	TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) { 
-		// walk around the vertex 
-		// version 2 for remove the probleme if we fill the hole
-		//int bug=1;
-		//  Triangle *torigine = t;
-		// restart:
-		//   int dir=0;
-		if (t->link != 0){
-			throw ErrorException(__FUNCT__,exprintf("t->link != 0"));
-		}
-		// to have a starting edges 
-		// try the 3 edge bourna-- in case of internal hole 
-		// and choice  the best 
-		// 
-		// the probleme is in case of  the fine and long internal hole
-		// for exemple neart the training edge of a wing
-		Vertex * s=0,*s1=0, *s0=0;
-		Icoor2 imax = MaxICoor22;
-		Icoor2 l0 = imax,l1 = imax;
-		double dd2 =  imax;// infinity
-		TriangleAdjacent er; 
-		int  cas=-2;
-		for (int j=0;j<3;j++)
-		  { 
-			TriangleAdjacent ta=t->FindBoundaryEdge(j);
-			if  (! (Triangle *) ta) continue;
-			s0 = ta.EdgeVertex(0);
-			s1 = ta.EdgeVertex(1);
-			I2 A = * s0;
-			I2 B = *ta.EdgeVertex(1);
-			I2 AB = B-A,AC=C-A,BC=B-C;
-			Icoor2  ACAC = (AC,AC), BCBC = (BC,BC);
-			Icoor2  AB2  =   Norme2_2(AB); //  ||AB||^2
-			Icoor2  ABAC  =   (AB,AC);         //  AB.AC|
-
-			double d2;
-			if ( ABAC < 0 )   // DIST A
-			  {
-				if ( (d2=(double) ACAC)  <  dd2) 
-				  {
-					er = ta;
-					l0 = ACAC;
-					l1 = BCBC;
-					cas = 0;
-					s = s0;
-				  }
-			  }
-			else if (ABAC > AB2)  // DIST B
-			  {
-				if ( (d2=(double) BCBC)  <  dd2) 
-				  {
-					dd2 = d2;
-					er = Adj(ta); // other direction
-					l0 = BCBC;
-					l1 = ACAC;
-					cas = 1;
-					s = s1;
-				  }
-			  }
-			else  // DIST AB
-			  { 
-
-				double det_2 =  (double) Det(AB,AC); 
-				det_2 *= det_2; // square of area*2 of triangle ABC
-				d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC      
-
-				if (d2 < dd2) 
-				  {
-					dd2 = d2;
-					er = ta;
-					l0 = (AC,AC);
-					l1 = (BC,BC);
-					s = 0;
-					cas = -1;
-					b = ((double) ABAC/(double) AB2);
-					a = 1 - b;
-				  }
-			  }
-		  }
-		if (cas ==-2){
-			throw ErrorException(__FUNCT__,exprintf("cas==-2"));
-		}
-		// l1 = ||C s1||  , l0 = ||C s0||
-		// where s0,s1 are the vertex of the edge er
-
-		if ( s) 
-		  { 
-			t=er;
-			TriangleAdjacent edge(er); 
-
-			int kkk=0;  
-			int linkp = t->link == 0;
-
-			Triangle * tt=t=edge=Adj(Previous(edge));
-			do  {  // loop over vertex s
-				kkk++;
-				if (edge.EdgeVertex(0)!=s && kkk>=10000){
-					throw ErrorException(__FUNCT__,exprintf("edge.EdgeVertex(0)!=s && kkk>=10000"));
-				}
-
-				int link = tt->link == 0;
-				if ((link + linkp) == 1) 
-				  { // a boundary edge 
-					Vertex * st = edge.EdgeVertex(1);
-					I2 I=*st;
-					Icoor2  ll = Norme2_2 (C-I);
-					if (ll < l1) {  // the other vertex is neart 
-						s1=st;
-						l1=ll;
-						er = edge;
-						if(ll<l0) { // change of direction --
-							s1=s;
-							l1=l0;
-							s=st;
-							l0=ll;
-							t=tt;
-							edge=Adj(edge);
-							link=linkp;
-							er = edge;
-						}
-					}
-				  }
-
-				linkp=link;
-				edge=Adj(Previous(edge));
-				tt = edge;
-			} while (t!=tt);
-
-			if (!(Triangle *) er){
-				throw ErrorException(__FUNCT__,exprintf("!(Triangle *) er"));
-			}
-			I2 A((I2)*er.EdgeVertex(0));
-			I2 B((I2)*er.EdgeVertex(1));
-			I2 AB=B-A,AC=C-A,CB=B-C;
-			double aa =  (double) (AB,AC);
-			double bb =  (double) (AB,CB);
-			if (aa<0)       a=1,b=0;
-			else if(bb<0)   a=0,b=1;
-			else  
-			  {
-				a  = bb/(aa+bb);
-				b  = aa/(aa+bb);
-			  }
-		  }
-		return er;
-	} 
-	/*}}}1*/
-	/*FUNCTION AGoodNumberPrimeWith{{{1*/
-	Int4 AGoodNumberPrimeWith(Int4 n){
-
-		//list of big prime numbers
-		const Int4 BigPrimeNumber[] ={ 567890359L,
-			567890431L,  567890437L,  567890461L,  567890471L,
-			567890483L,  567890489L,  567890497L,  567890507L,
-			567890591L,  567890599L,  567890621L,  567890629L , 0};
-
-		//initialize o and pi
-		Int4 o =0;
-		Int4 pi=BigPrimeNumber[1];
-
-		//loop until BigPrimeNumber[i]==0 (end of BigPrimeNumber)
-		for (int i=0; BigPrimeNumber[i]; i++){
-
-			//compute r, rest of the remainder of the division of BigPrimeNumber[i] by n
-			Int4 r = BigPrimeNumber[i] % n;
-
-			/*compute oo = min ( r , n-r , |n - 2r|, |n-3r|)*/
-			Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));
-			if ( o < oo){
-				o=oo;
-				pi=BigPrimeNumber[i];
-			}
-		}
-		return pi; 
-	}
-	/*}}}1*/
-	/*FUNCTION SwapForForcingEdge{{{1*/
-	int SwapForForcingEdge(Vertex   *  & pva ,Vertex  * &   pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {
-		// l'arete ta coupe l'arete pva pvb
-		// de cas apres le swap sa coupe toujours
-		// on cherche l'arete suivante 
-		// on suppose que detsa >0 et detsb <0
-		// attention la routine echange pva et pvb 
-
-		if(tt1.Locked()) return 0; // frontiere croise 
-
-		TriangleAdjacent tt2 = Adj(tt1);
-		Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
-		Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
-		if ( a1<0 || a1>=3 ){
-			throw ErrorException(__FUNCT__,exprintf("a1<0 || a1>=3"));
-		}
-
-		Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
-		Vertex & s1= (*t1)[OppositeVertex[a1]];
-		Vertex & s2= (*t2)[OppositeVertex[a2]];
-
-
-		Icoor2 dets2 = det(*pva,*pvb,s2);
-		Icoor2 det1=t1->det , det2=t2->det ;
-		Icoor2 detT = det1+det2;
-		if ((det1<=0 ) || (det2<=0)){
-			throw ErrorException(__FUNCT__,exprintf("(det1<=0 ) || (det2<=0)"));
-		}
-		if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb
-			throw ErrorException(__FUNCT__,exprintf("(detsa>=0) || (detsb<=0)"));
-		}
-		Icoor2 ndet1 = bamg::det(s1,sa,s2);
-		Icoor2 ndet2 = detT - ndet1;
-
-		int ToSwap =0; //pas de swap
-		if ((ndet1 >0) && (ndet2 >0)) 
-		  { // on peut swaper  
-			if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
-			 ToSwap =1; 
-			else // swap alleatoire 
-			 if (BinaryRand()) 
-			  ToSwap =2; 
-		  }
-		if (ToSwap) NbSwap++,
-		 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
-
-		int ret=1;
-
-		if (dets2 < 0) {// haut
-			dets1 = ToSwap ? dets1 : detsa ;
-			detsa = dets2; 
-			tt1 =  Previous(tt2) ;}
-		else if (dets2 > 0){// bas 
-			dets1 = ToSwap ? dets1 : detsb ;
-			detsb = dets2;
-			//xxxx tt1 = ToSwap ? tt1 : Next(tt2);
-			if(!ToSwap) tt1 =  Next(tt2);
-		}
-		else { // changement de sens 
-			ret = -1;
-			Exchange(pva,pvb);
-			Exchange(detsa,detsb);
-			Exchange(dets1,dets2);
-			Exchange(tt1,tt2);
-			dets1=-dets1;
-			dets2=-dets2;
-			detsa=-detsa;
-			detsb=-detsb;
-
-			if (ToSwap) 
-			 if (dets2 < 0) {// haut
-				 dets1 = (ToSwap ? dets1 : detsa) ;
-				 detsa = dets2; 
-				 tt1 =  Previous(tt2) ;}
-			 else if (dets2 > 0){// bas 
-				 dets1 = (ToSwap ? dets1 : detsb) ;
-				 detsb =  dets2;
-				 if(!ToSwap) tt1 =  Next(tt2);
-			 }
-			 else {// on a fin ???
-				 tt1 = Next(tt2);
-				 ret =0;}
-
-		}
-		return ret;
-	}
-	/*}}}1*/
-	/*FUNCTION ForceEdge{{{1*/
-
-	int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)  { 
-		int NbSwap =0;
-		if (!a.t || !b.t){ // the 2 vertex is in a mesh
-			throw ErrorException(__FUNCT__,exprintf("!a.t || !b.t"));
-		}
-		int k=0;
-		taret=TriangleAdjacent(0,0); // erreur 
-
-		TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]);
-		Vertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
-		// we turn around a in the  direct sens  
-
-		Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1;
-		if(v2) // normal case 
-		 det2 = det(*v2,a,b);
-		else { // no chance infini vertex try the next
-			tta= Previous(Adj(tta));
-			v2 = tta.EdgeVertex(0);
-			vbegin =v2;
-			if (!v2){
-				throw ErrorException(__FUNCT__,exprintf("!v2"));
-			}
-			det2 = det(*v2,a,b);
-		}
-
-		while (v2 != &b) {
-			TriangleAdjacent tc = Previous(Adj(tta));    
-			v1 = v2; 
-			v2 = tc.EdgeVertex(0);
-			det1 = det2;
-			det2 =  v2 ? det(*v2,a,b): det2; 
-
-			if((det1 < 0) && (det2 >0)) { 
-				// try to force the edge 
-				Vertex * va = &a, *vb = &b;
-				tc = Previous(tc);
-				if (!v1 || !v2){
-					throw ErrorException(__FUNCT__,exprintf("!v1 || !v2"));
-				}
-				Icoor2 detss = 0,l=0,ks;
-				while ((ks=SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
-				 if(l++ > 10000000) {
-					 throw ErrorException(__FUNCT__,exprintf("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i) too big",NbSwap,l));
-				 }
-				Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
-				if (( aa == &a ) && (bb == &b) ||  (bb ==  &a ) && (aa == &b)) {
-					tc.SetLock();
-					a.Optim(1,0);
-					b.Optim(1,0);
-					taret = tc;
-					return NbSwap;
-				}
-				else 
-				  {
-					taret = tc;
-					return -2; // error  boundary is crossing
-				  }
-			}
-			tta = tc;
-			k++;
-			if (k>=2000){
-				throw ErrorException(__FUNCT__,exprintf("k>=2000"));
-			}
-			if ( vbegin == v2 ) return -1;// error 
-		}
-
-		tta.SetLock();
-		taret=tta;
-		a.Optim(1,0);
-		b.Optim(1,0);
-		return NbSwap; 
-	}
-	/*}}}1*/
-
 }