Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp	(revision 16966)
+++ /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp	(revision 16967)
@@ -20,4 +20,5 @@
 	this->omega=0;
 	this->power=0;
+	this->random=0;
 	this->verbose=0;
 
Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/BamgOpts.h	(revision 16966)
+++ /issm/trunk-jpl/src/c/bamg/BamgOpts.h	(revision 16967)
@@ -25,4 +25,5 @@
 		double  omega;
 		double  power;
+		bool    random;
 		int     verbose;
 
Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16966)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16967)
@@ -2775,5 +2775,5 @@
 	/*}}}*/
 	/*FUNCTION Mesh::Insert{{{*/
-	void Mesh::Insert() {
+	void Mesh::Insert(bool random) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
 
@@ -2822,5 +2822,6 @@
 		//Get Prime number
 		const long PrimeNumber= BigPrimeNumber(nbv);
-		int   k0=rand()%nbv; 
+		int temp = rand(); if(random) temp = 756804110;
+		int  k0=temp%nbv; 
 
 		//Build orderedvertices
@@ -2904,28 +2905,8 @@
 			_printf_("      NbSwap/nbv:          " << NbSwap/nbv << "\n");
 		}
-
-#ifdef NBLOOPOPTIM
-
-		k0 = rand()%nbv ; 
-		for (int is4=0; is4<nbv; is4++) 
-		 orderedvertices[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];
-
-		for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){
-			long  NbSwap = 0;
-			for (int is1=0; is1<nbv; is1++) 
-			 NbSwap += orderedvertices[is1]->Optim(0,0);
-			if (verbose>3) {
-				_printf_("      Optim Loop: " << Nbloop << "\n");
-				_printf_("      NbSwap/nbv:          " << NbSwap/nbv << "\n");
-			}
-			if(!NbSwap) break;
-		}
-		CreateSingleVertexToTriangleConnectivity(); 
-		// because we break the TriangleContainingTheVertex
-#endif
 	}
 	/*}}}*/
 	/*FUNCTION Mesh::InsertNewPoints{{{*/
-	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {
+	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap,bool random) {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
 
@@ -2948,5 +2929,6 @@
 		const long PrimeNumber= BigPrimeNumber(nbv)  ;
 		//remainder of the division of rand() by nbvnew
-		long k3 = rand()%nbvnew;
+		int  temp = rand(); if(!random) temp = 1098566905;
+		long k3 = temp%nbvnew;
 		//loop over the new points
 		for (int is3=0; is3<nbvnew; is3++){
@@ -3212,5 +3194,5 @@
 			}
 			Bh.CreateSingleVertexToTriangleConnectivity();     
-			InsertNewPoints(nbvold,NbTSwap)   ;            
+			InsertNewPoints(nbvold,NbTSwap,bamgopts->random);
 		}  
 		else Bh.CreateSingleVertexToTriangleConnectivity();     
@@ -3270,5 +3252,5 @@
 			  }// for triangle   
 
-			if (!InsertNewPoints(nbvold,NbTSwap)) break;
+			if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break;
 			for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
 			Headt = nbt; // empty list 
@@ -4231,5 +4213,5 @@
 
 	/*Insert Vertices*/
-	Insert();
+	Insert(true);
 }
 /*}}}*/
@@ -4534,5 +4516,5 @@
 		if (verbose>3) _printf_("      Creating initial Constrained Delaunay Triangulation...\n");
 		if (verbose>3) _printf_("         Inserting boundary points\n");
-		Insert();
+		Insert(bamgopts->random);
 
 		//Force the boundary
@@ -4862,5 +4844,5 @@
 		if (verbose>3) _printf_("      Creating initial Constrained Delaunay Triangulation...\n");
 		if (verbose>3) _printf_("         Inserting boundary points\n");
-		Insert();
+		Insert(bamgopts->random);
 
 		//Force the boundary
Index: /issm/trunk-jpl/src/c/bamg/Mesh.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.h	(revision 16966)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.h	(revision 16967)
@@ -77,5 +77,5 @@
 			R2 I2ToR2(const I2 & P) const;
 			void AddVertex(BamgVertex & s,Triangle * t,Icoor2 *  =0) ;
-			void Insert();
+			void Insert(bool random);
 			void Echo(void);
 			void ForceBoundary();
@@ -92,5 +92,5 @@
 			void MakeBamgQuadtree();
 			void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1);
-			long InsertNewPoints(long nbvold,long & NbTSwap) ; 
+			long InsertNewPoints(long nbvold,long & NbTSwap,bool random); 
 			void TrianglesRenumberBySubDomain(bool justcompress=false);
 			void SmoothingVertex(int =3,double=0.3);
Index: /issm/trunk-jpl/src/m/mesh/bamg.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 16966)
+++ /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 16967)
@@ -325,4 +325,5 @@
 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
 bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0);
+bamg_options.random=getfieldvalue(options,'rand',true);
 bamg_options.verbose=getfieldvalue(options,'verbose',1);
 %}}}
Index: /issm/trunk-jpl/src/m/mesh/bamg.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 16966)
+++ /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 16967)
@@ -309,4 +309,5 @@
 	bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1)
 	bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0)
+	bamg_options['random']=options.getfieldvalue('rand',True)
 	bamg_options['verbose']=options.getfieldvalue('verbose',1)
 	#}}}