Index: /issm/trunk/src/c/Bamgx/QuadTree.h
===================================================================
--- /issm/trunk/src/c/Bamgx/QuadTree.h	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/QuadTree.h	(revision 2854)
@@ -1,2 +1,6 @@
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+
 namespace bamg {
 
@@ -58,12 +62,11 @@
   void  Add( Vertex & w);
 
-  QuadTreeBox* NewQuadTreeBox()
-  {
-    ///cout << "NewQuadTreeBox " << sb << " " << sb->bc << " " 
-    //<< sb->be << " " <<lenStorageQuadTreeBox <<endl;
+  QuadTreeBox* NewQuadTreeBox(){
     if(! (sb->bc<sb->be)) 
-	sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
+	  sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
 
-    assert(sb && (sb->bc->n == 0));
+	 if (!sb || (sb->bc->n != 0)){
+		 throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));
+	 }
     NbQuadTreeBox++;
     return sb->bc++;
Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp	(revision 2854)
@@ -8,4 +8,8 @@
 #include "../QuadTree.h"
 #include "../SetOfE4.h"
+
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
@@ -27,6 +31,10 @@
 		// Real8 t1t1 = t1*t1;
 		Real8 tt = theta*theta;
-		assert( theta >=0);
-		assert( theta <=1);
+		if ( theta<0){
+			throw ErrorException(__FUNCT__,exprintf("theta<0"));
+		}
+		if ( theta>1){
+			throw ErrorException(__FUNCT__,exprintf("theta>1"));
+		}
 		if (TgA()) 
 		 if (TgB()) // interpolation d'hermite
@@ -92,6 +100,10 @@
 	  { R2 A=v[0]->r,B=v[1]->r;
 		Real8 ca,cb,cta,ctb;
-		assert( theta >=-1e-12);
-		assert( theta <=1+1e-12);
+		if ( theta<-1e-12){
+			throw ErrorException(__FUNCT__,exprintf("theta<-1e-12"));
+		}
+		if ( theta>1+1e-12){
+			throw ErrorException(__FUNCT__,exprintf("theta>1+1e-12"));
+		}
 		if (TgA()) 
 		 if (TgB()) // interpolation d'hermite
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2854)
@@ -42,8 +42,7 @@
 		for (i=0;i<NbSubDomains;i++)
 		 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
-
-		//    for (i=0;i<nbt;i++)
-		//      triangles[i].Set(Gh.triangles[i],Gh,*this);
-		assert(!nbt);   
+		if (nbt);  {
+			throw ErrorException(__FUNCT__,exprintf("nbt"));
+		}
 	}
 	/*}}}1*/
@@ -52,5 +51,7 @@
 		long int verbosity=0;
 
-		assert(NbRef<=0);
+		if (NbRef>0){
+			throw ErrorException(__FUNCT__,exprintf("NbRef>0"));
+		}
 		if(verbosity>9)
 		 cout << "DELETE      ~Geometry "<< this  << endl;
@@ -512,8 +513,10 @@
 		//  cout << P ;
 		int k=0;
-		while(pon != on)
-		  {  
+		while(pon != on){  
+			k++;
 			pon = on;
-			assert(k++<100);
+			if (k>=100){
+				throw ErrorException(__FUNCT__,exprintf("k>=100"));
+			}
 			R2 A= (*on)[0];
 			R2 B= (*on)[1];
@@ -538,7 +541,11 @@
 retry:    
 		s=save_s;
-		GeometricalEdge * on = e.on;
-		assert(on);
-		assert( e[0].on &&  e[1].on);
+		GeometricalEdge* on = e.on;
+		if (!on){
+			throw ErrorException(__FUNCT__,exprintf("!on"));
+		}
+		if (!e[0].on ||  !e[1].on){
+			throw ErrorException(__FUNCT__,exprintf("!e[0].on ||  !e[1].on"));
+		}
 		const Vertex &v0=e[0],&v1=e[1];
 		V.m = Metric(1.0-s, v0,s, v1);
@@ -586,5 +593,7 @@
 					Number(eg0->Adj[1]) <<"," ;
 				ge[--bge] =eg0 = eg0->Adj[sens0];
-				assert(bge>=0 && bge <= mxe);
+				if (bge<0 || bge>mxe){
+					throw ErrorException(__FUNCT__,exprintf("bge<0 || bge>mxe"));
+				}
 				sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]);
 				if(NbTry)
@@ -615,5 +624,7 @@
 			ge[++tge] =eg1 = eg1->Adj[sens1];
 			sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1];
-			assert(tge>=0 && tge <= mxe);
+			if (tge<0 || tge>mxe){
+				throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)"));
+			}
 			if(NbTry)
 			 cout << "  Edge "  <<  Number(eg1) << " " << sens1 << " S "
@@ -648,7 +659,8 @@
 			int i=bge;
 			Real8 ll=0;
-			for(i=bge;i<tge;i++) 
-			  {
-				assert( i>=0 && i <= mxe);
+			for(i=bge;i<tge;i++){
+				if ( i<0 || i>mxe){
+					throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe"));
+				}
 				BB =  (*ge[i])[sensge[i]];
 				lge[i]=ll += Norme2(AA-BB);
@@ -660,5 +672,7 @@
 				//	     <<" on = " << Number(ge[tge]) <<  " sens= " << sensge[tge] << endl;
 				// search the geometrical edge
-				assert(s <= 1.0);
+				if (s>1.0){
+					throw ErrorException(__FUNCT__,exprintf("s>1.0"));
+				}
 				Real8 ls= s*ll;
 				on =0;
@@ -668,5 +682,7 @@
 				i=bge;
 				while (  (l1=lge[i]) < ls ) {
-					assert(i >= 0 && i <= mxe);
+					if (i<0 || i>mxe){
+						throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe"));
+					}
 					i++,s0=1-(s1=sensge[i]),l0=l1;}
 					on=ge[i];
@@ -679,6 +695,7 @@
 					sg =  s0 * (1.0-s) +  s * s1;    
 		} 
-		assert(on);
-		// assert(sg && sg-1);
+		if (!on){
+			throw ErrorException(__FUNCT__,exprintf("!on"));
+		}
 		V.r= on->F(sg);
 		//  if (eg0 != eg1) 
@@ -716,5 +733,7 @@
 				GeometricalVertex * vg = (GeometricalVertex  *) (void *) v;
 				int j = vg-v0g;
-				assert( v ==  & (Vertex &) vertices[j]);
+				if ( v !=  &(Vertex &) vertices[j]){
+					throw ErrorException(__FUNCT__,exprintf(" v !=  &(Vertex &) vertices[j]"));
+				}
 				vertices[i].link = vertices + j;
 				k++;	      
@@ -949,6 +968,9 @@
 					  }
 
-				  }} 
-				 assert(nbgem && nbe);
+				  }
+			 } 
+			if (nbgem==0 || nbe==0){
+				throw ErrorException(__FUNCT__,exprintf("nbgem==0 || nbe==0"));
+			}
 
 				 if(step==0) {
@@ -961,9 +983,12 @@
 			//GeometricalEdge * ee=curves[i].ee, *eqee=be->link;
 			curves[i].master=true;
-			if(be->Equi() || be->ReverseEqui() ) 
-			  {
-				assert(eqbe);
+			if(be->Equi() || be->ReverseEqui() ){
+				if (!eqbe){
+					throw ErrorException(__FUNCT__,exprintf("!eqbe"));
+				}
 				int nc = eqbe->CurveNumber;
-				assert(i!=nc);
+				if (i==nc){
+					throw ErrorException(__FUNCT__,exprintf("i==nc"));
+				}
 				curves[i].next=curves[nc].next;
 				curves[i].master=false;
Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp	(revision 2854)
@@ -67,5 +67,7 @@
 			while (t->det <0) { // intersection boundary edge and a,b,
 				k=(*t)(0) ?  ((  (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1  )) : 0;
-				assert(k>=0);
+				if (k<0){
+					throw ErrorException(__FUNCT__,exprintf("k<0"));
+				}
 				ocut = OppositeEdge[k];
 				i=VerticesOfTriangularEdge[ocut][0];
@@ -114,5 +116,7 @@
 				j=VerticesOfTriangularEdge[iedge][1];
 				Real8 dij = detj-deti;
-				assert(i+j+k == 0 + 1 +2);
+				if (i+j+k != 0 + 1 +2){
+					throw ErrorException(__FUNCT__,exprintf("i+j+k != 0 + 1 +2"));
+				}
 				ba[j] =  detj/dij;
 				ba[i] = -deti/dij;
@@ -268,5 +272,7 @@
 	Real8  ListofIntersectionTriangles::Length(){
 		//  cout << " n= " << Size << ":" ;
-		assert(Size>0);
+		if (Size<=0){
+			throw ErrorException(__FUNCT__,exprintf("Size<=0"));
+		}
 		// computation of the length      
 		R2 C;
Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2854)
@@ -8,4 +8,8 @@
 #include "../QuadTree.h"
 #include "../SetOfE4.h"
+
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
@@ -206,5 +210,7 @@
 		}
 		// warning for optimisation S is in [0:0.5] not in [0:1]
-		assert(i<512);
+		if (i>=512){
+			throw ErrorException(__FUNCT__,exprintf("i>=512"));
+		}
 		LastMetricInterpole.lab=l;
 		LastMetricInterpole.opt=i;
@@ -243,5 +249,7 @@
 			 r =  2*(S[i]*(L[j]-l)+ S[j]*(l-L[i]))/(L[j]-L[i]);
 		  }
-		assert(r<=1 && r>=0);
+		if (r>1 || r<0){
+			throw ErrorException(__FUNCT__,exprintf("r>1 || r<0"));
+		}
 		return r ;
 	}
Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp	(revision 2854)
@@ -5,4 +5,8 @@
 #include "../Mesh2.h"
 #include "../QuadTree.h"
+
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../../toolkits/toolkits.h"
 
 #undef __FUNCT__ 
@@ -33,5 +37,7 @@
 	 sb =new StorageQuadTreeBox(lenStorageQuadTreeBox);
 	 root=NewQuadTreeBox();
-	 assert( MaxISize > MaxICoor);
+	 if ( MaxISize <= MaxICoor){
+		 throw ErrorException(__FUNCT__,exprintf("MaxISize <= MaxICoor"));
+	 }
 	 for (Int4 i=0;i<nbv;i++) 
 	  Add(t->vertices[i]);
@@ -373,5 +379,7 @@
 			if (b->n > 0 &&  b->v[0] == &w) return;
 		}
-		assert(l);
+		if (l==0){
+			throw ErrorException(__FUNCT__,exprintf("l==0"));
+		}
 		while ((b= *pb) && (b->n == 4)) // the QuadTreeBox is full
 		  { 
@@ -412,5 +420,7 @@
 		bc =b;
 		be = b +ll;
-		assert(b);
+		if (!b){
+			throw ErrorException(__FUNCT__,exprintf("!b"));
+		}
 	}
 	/*}}}1*/
Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 2854)
@@ -37,5 +37,7 @@
 			t = ttc;
 			j = NextEdge[jc];
-			assert(k<2000);
+			if (k>=2000){
+				throw ErrorException(__FUNCT__,exprintf("k>=2000"));
+			}
 		  } while ( (this!= t)); 
 		return TriangleAdjacent(0,0);
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2854)
@@ -2454,7 +2454,4 @@
 		Int4 BeginNewPoint[3];
 		Int4 EndNewPoint[3];
-#ifdef TRACETRIANGLE
-		Int4 trace=0;
-#endif
 		int step[3];
 		Int4 *first_np_or_next_t = new Int4[nbtx];
@@ -2480,11 +2477,9 @@
 			i=Headt;
 			next_t=-first_np_or_next_t[i];
-			for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) 
-			  { // for each triangle  t
+			for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) { // for each triangle  t
 				// we can change first_np_or_next_t[i]
-#ifdef TRACETRIANGLE
-				trace =  TRACETRIANGLE <0 ? 1 : i == TRACETRIANGLE;
-#endif
-				assert(i>=0 && i < nbt );
+				if (i<0 || i >= nbt ){
+					throw ErrorException(__FUNCT__,exprintf("i<0 || i >= nbt"));
+				}
 				first_np_or_next_t[i] = nbv; // to save the fist new point of triangle
 				for(j=0;j<3;j++)
@@ -2500,10 +2495,4 @@
 					Vertex & vB = * tj.EdgeVertex(1);
 
-#ifdef TRACETRIANGLE
-					if(trace) {
-						cout << " i " << Number(vA) <<" j "<<  Number(vB) 
-						  << " "  << t->Locked(j) ;
-					}
-#endif
 					if (!t->link) continue;// boundary
 					if (t->det <0) continue;
@@ -2523,21 +2512,13 @@
 					const Vertex & vC2 = *tadjj.OppositeVertex();
 
-#ifdef TRACETRIANGLE
-					trace = trace || k == TRACETRIANGLE;
-					if(trace) {
-						cout << "Test Arete " << i << " AB = " << A << B 
-						  << "i "  <<Number(vA)<< "j" <<Number(vB); 
-						cout << " link" <<(int)t->link << " ta=" << Number( ta) 
-						  << " det " <<ta->det ;
-						cout << " hA = " <<vA.m.h << " hB = " <<vB.m.h ;
-						cout << " loc " << ta->Locked(j) << endl;
-					}
-#endif
-
 					if(first_np_or_next_t[k]>0) { // this edge is done before 
 						// find the color of the edge and begin , end of newpoint
 						register int kk = t->NuEdgeTriangleAdj(j);
-						assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1]));
-						assert ((*t)(VerticesOfTriangularEdge[j][1]) == (*ta)(VerticesOfTriangularEdge[kk][0]));
+						if ((*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])){
+							throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])"));
+						}
+						if ((*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])){
+							throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])"));
+						}
 						register Int4 kolor =3*k + kk;
 						ColorEdge[j]=kolor;
@@ -2578,17 +2559,4 @@
 					EndNewPoint[j] = nbv-1;
 				  } // end loop for each edge 
-
-#ifdef TRACETRIANGLE
-				if(trace) {
-					// verification des point cree 
-					cout << "\n ------------ " << t->link << " " << t->det 
-					  << " b " << BeginNewPoint[0] << " " << BeginNewPoint[1]
-					  << " " << BeginNewPoint[2] << " " 
-					  << " e " << EndNewPoint[0] << " " << EndNewPoint[1] 
-					  << " " << EndNewPoint[2] << " " 
-					  << " s " << step[0] << " " << step[1] << " " << step[2] << " " 
-					  <<  endl;
-				}
-#endif
 
 				if (!t->link) continue;// boundary
@@ -2664,21 +2632,19 @@
 					 Int4 i1;
 					 kstack =0; 
-					 while( (i1=stack[kstack++]) >= 0) 
-						{ // the two parameter is i0 and i1 
-						 assert(i1 < nbv && i1 >= 0);
-						 assert(i0 < nbv && i0 >= 0);
-						 assert(i1 != i0);
+					 while( (i1=stack[kstack++]) >= 0){ // the two parameter is i0 and i1 
+						 if (i1 >= nbv || i1 < 0){
+							 throw ErrorException(__FUNCT__,exprintf("i1 >= nbv || i1 < 0"));
+						 }
+						 if (i0 >= nbv || i0 < 0){
+							 throw ErrorException(__FUNCT__,exprintf("i0 >= nbv || i0 < 0"));
+						 }
+						 if (i1 == i0){
+							 throw ErrorException(__FUNCT__,exprintf("i1 == i0"));
+						 }
 						 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0];
 						 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01));
-
-
-#ifdef TRACETRIANGLE
-						 if(trace) {
-							 cout << "\n test j" << j <<" "  << i0 
-								<< " " << i1 << " d01=" << d01 <<endl;}
-#endif
-						 assert (i0 >= nbvold);
-						 assert (i1 >= nbvold);
-						 assert(i0 != i1);
+						 if (i0< nbvold || i1<nbvold){
+							 throw ErrorException(__FUNCT__,exprintf("i0< nbvold || i1<nbvold"));
+						 }
 						 if (d01 == 0) 
 						  break; 
@@ -2708,9 +2674,4 @@
 							  =  vertices[i0].r *w0 + vertices[i1].r* w1;
 
-#ifdef TRACETRIANGLE
-							  if(trace) {
-								  cout << "\n ref = "<< vertices[i0].ref << " " <<vertices[i1].ref <<endl;
-							  }
-#endif    
 							  // update the new point points of the list 
 							  for  (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0)
@@ -2767,5 +2728,7 @@
 
 					quadtree->Add(*vi); //
-					assert (tcvi->det >= 0) ;// internal 
+					if (tcvi->det<0){// internal
+						throw ErrorException(__FUNCT__,exprintf("tcvi->det<0"));
+					}
 					Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);          
 				  }  
@@ -2791,5 +2754,7 @@
 					if (first_np_or_next_t[kt]>0) 
 					 first_np_or_next_t[kt]=-Headt,Headt=kt;
-					assert( ta.EdgeVertex(0) == s);
+					if (ta.EdgeVertex(0) != s){
+						throw ErrorException(__FUNCT__,exprintf("ta.EdgeVertex(0) != s"));
+					}
 					ta = Next(Adj(ta));
 				} while ( (tbegin != (Triangle*) ta)); 
@@ -3098,5 +3063,7 @@
 						//    else if(mark[it] == -2 ) triangles[it].Draw(999);
 						it++;} // end white (it<nbt)
-						assert(k== NbSubDomains);
+						if (k!=NbSubDomains){
+							throw ErrorException(__FUNCT__,exprintf("k!=NbSubDomains"));
+						}
 						if(OutSide) 
 						  {
@@ -3179,5 +3146,7 @@
 					// see routine MakeGeometricalEdgeToEdge
 					Edge &e = *GeometricalEdgetoEdge[Gh.Number(eg)];
-					assert(&e);
+					if (&e==NULL){
+						throw ErrorException(__FUNCT__,exprintf("&e==NULL"));
+					}
 					Vertex * v0 =  e(0),*v1 = e(1);
 					Triangle *t  = v0->t;
@@ -3186,16 +3155,17 @@
 					//	cout << " geom edge = " <<  Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref 
 					//     << " ref edge =" << eg.ref << " sens " << sens ;
-					if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0)
-					 sens = -sens ;
+					if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) sens = -sens ;
 					subdomains[i].sens = sens;
 					subdomains[i].edge = &e;
-					//	cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh  " << e[0].r << e[1].r << endl;
-					//	cout << "  v0 , v1 = " << Number(v0) << " "  << Number(v1) << endl;
-					assert(t && sens);
+					if (!t || !sens){
+						throw ErrorException(__FUNCT__,exprintf("!t || !sens"));
+					}
 
 					TriangleAdjacent  ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges
 
 					while (1) {
-						assert( v0 == ta.EdgeVertex(1) );
+						if ( v0 != ta.EdgeVertex(1) ){
+							throw ErrorException(__FUNCT__,exprintf("v0 != ta.EdgeVertex(1)"));
+						}
 						//	 cout << " recherche " << Number( ta.EdgeVertex(0)) << endl;
 						if (ta.EdgeVertex(0) == v1) { // ok we find the edge
@@ -3225,5 +3195,7 @@
 							  {
 								kkk++;
-								assert(mark[Number(tt)]<0);
+								if (mark[Number(tt)]>=0){
+									throw ErrorException(__FUNCT__,exprintf("mark[Number(tt)]>=0"));
+								}
 								mark[Number(tt)]=i;
 								tt=tt->link;
@@ -3344,9 +3316,15 @@
 		  { 
 			t=t0=subdomains[i].head;
-			assert(t0); // no empty sub domain
+			if (!t0){ // not empty sub domain
+				throw ErrorException(__FUNCT__,exprintf("!t0"));
+			}
 			do { 
 				Int4 kt = Number(t);
-				assert(kt>=0 && kt < nbt );
-				assert(renu[kt]==-1);
+				if (kt<0 || kt >= nbt ){
+					throw ErrorException(__FUNCT__,exprintf("kt<0 || kt >= nbt"));
+				}
+				if (renu[kt]!=-1){
+					throw ErrorException(__FUNCT__,exprintf("renu[kt]!=-1"));
+				}
 				renu[kt]=k++;
 			}
@@ -3362,9 +3340,10 @@
 
 		// put the outside triangles at the end
-		for ( it=0;it<nbt;it++) 
-		 if (renu[it]==-1) 
-		  renu[it]=k++;
-
-		assert(k == nbt);
+		for ( it=0;it<nbt;it++){
+			if (renu[it]==-1) renu[it]=k++;
+		}
+		if (k != nbt){
+			throw ErrorException(__FUNCT__,exprintf("k != nbt"));
+		}
 		// do the change on all the pointeur 
 		for ( it=0;it<nbt;it++)
@@ -3396,21 +3375,21 @@
 	/*}}}1*/
 	/*FUNCTION Triangles::ConsRefTriangle{{{1*/
-	Int4  Triangles::ConsRefTriangle(Int4 *reft) const {
+	Int4  Triangles::ConsRefTriangle(Int4* reft) const {
 		long int verbosity=0;
-		assert(reft);
 		register Triangle *t0,*t;
 		register Int4 k=0, num;   
-		for (Int4 it=0;it<nbt;it++) 
-		 reft[it]=-1; // outside triangle 
-		for (Int4 i=0;i<NbSubDomains;i++)
-		  { 
+		for (Int4 it=0;it<nbt;it++) reft[it]=-1; // outside triangle 
+		for (Int4 i=0;i<NbSubDomains;i++){ 
 			t=t0=subdomains[i].head;
-			assert(t0); // no empty sub domain
+			if (!t0){ // no empty sub domai{
+				throw ErrorException(__FUNCT__,exprintf("!t0"));
+			}
 			// register Int4 color=i+1;// because the color 0 is outside triangle
 			do { k++;
 				num = Number(t);
-				assert(num>=0 &&num < nbt);
+				if (num<0 || num>=nbt){
+					throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt"));
+				}
 				reft[num]=i;
-				//  cout << Number(t0) << " " <<Number(t)<< " "  << i << endl;
 			}
 			while (t0 != (t=t->link));
@@ -3446,7 +3425,8 @@
 1 on GeometricalEdge + abcisse
 2 internal 
-
 		 *************************************************************************/
-		assert(&BTh.Gh == &Gh);
+		if (&BTh.Gh != &Gh){
+			throw ErrorException(__FUNCT__,exprintf("&BTh.Gh != &Gh"));
+		}
 
 		long int verbosity=0;
@@ -3481,30 +3461,36 @@
 			throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx));
 		}
-		assert(vertices);
-		for (i=0;i<Gh.nbv;i++)
-		 if (Gh[i].Required()) {//Gh  vertices Required
-			 vertices[nbv] =  Gh[i];
-			 vertices[nbv].i = I2(0,0);
-			 Gh[i].to = vertices + nbv;// save Geom -> Th
-			 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
-			 // cout << "--------- "  <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
-			 nbv++;}
-		 else Gh[i].to=0;
-		// 
-		for (i=0;i<BTh.NbVerticesOnGeomVertex;i++)
-		  { 
+		if (vertices==NULL){
+			throw ErrorException(__FUNCT__,exprintf("vertices==NULL"));
+		}
+		for (i=0;i<Gh.nbv;i++){
+			if (Gh[i].Required()) {//Gh  vertices Required
+				vertices[nbv] =  Gh[i];
+				vertices[nbv].i = I2(0,0);
+				Gh[i].to = vertices + nbv;// save Geom -> Th
+				VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
+				// cout << "--------- "  <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
+				nbv++;
+			}
+			else Gh[i].to=0;
+		}
+		for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){ 
 			VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i];
 			if (vog.IsRequiredVertex()) {
 				GeometricalVertex * gv = vog;
 				Vertex *bv = vog;
-				assert(gv->to);// use of Geom -> Th
+				if (!gv->to){// use of Geom -> Th
+					throw ErrorException(__FUNCT__,exprintf("!gv->to"));
+				}
 				VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv);
 				gv->to->m = bv->m; // for taking the metrix of the background mesh
-				;}
+				;
+			}
 		  }
-		assert(NbVertexOnBThVertex == NbVerticesOnGeomVertex);
+		if (NbVertexOnBThVertex != NbVerticesOnGeomVertex){
+			throw ErrorException(__FUNCT__,exprintf("NbVertexOnBThVertex != NbVerticesOnGeomVertex"));
+		}
 		// new stuff FH with curve
 		//  find the begin of the curve in BTh
-		  {
 			Gh.UnMarkEdges();	
 			int bfind=0;
@@ -3543,6 +3529,8 @@
 					}
 			  } 
-			assert( bfind==Gh.NbOfCurves);
-		  }          
+			if ( bfind!=Gh.NbOfCurves){
+				throw ErrorException(__FUNCT__,exprintf("bfind!=Gh.NbOfCurves"));
+			}
+
 		// method in 2 + 1 step 
 		//  0.0) compute the length and the number of vertex to do allocation
@@ -3610,7 +3598,8 @@
 						VertexOnGeom *GA1;
 						Edge * PreviousNewEdge = 0;
-						//  cout << "  --------------New Curve phase " << phase 
-						//       << "---------- A0=" << *A0 << ei[k0]  <<endl;
-						assert (A0-vertices>=0 && A0-vertices <nbv);
+						// New Curve phase 
+						if (A0-vertices<0 || A0-vertices>=nbv){
+							throw ErrorException(__FUNCT__,exprintf("A0-vertices<0 || A0-vertices>=nbv"));
+						}
 						if(ongequi->Required() ) 
 						  {
@@ -3619,7 +3608,5 @@
 						  }       
 						else 
-						 for(;;) 
-							{
-							 //   assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe);
+						 for(;;){
 							 Edge &ee=*pe; 
 							 Edge &eeequi=*peequi; 
@@ -3627,5 +3614,7 @@
 							 k1equi= 1 - k0equi;
 
-							 assert(pe  && ee.on);
+							 if (!pe  || !ee.on){
+								 throw ErrorException(__FUNCT__,exprintf("!pe  || !ee.on"));
+							 }
 							 ee.on->SetMark();
 							 Vertex & v0=ee[0], & v1=ee[1];
@@ -3636,13 +3625,19 @@
 							 if (phase) {// computation of the new points
 								 while ((i!=NbCreatePointOnCurve) && sNew <= L) { 
-									 //    cout  << " L0= " << L0 << " L " << L << " sN=" 
-									 //         << sNew << " LAB=" << LAB << " NBPC =" <<NbCreatePointOnCurve<< " i " << i  << endl;
-									 assert (sNew >= L0);
-									 assert(LAB);
-
-
-									 assert(vertices && nbv<nbvx);
-									 assert(edges && nbe < nbex);
-									 assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge < NbVerticesOnGeomEdgex);
+									 if (sNew<L0){
+										 throw ErrorException(__FUNCT__,exprintf("sNew<L0"));
+									 }
+									 if (!LAB){
+										 throw ErrorException(__FUNCT__,exprintf("!LAB"));
+									 }
+									 if (!vertices || nbv>=nbvx){
+										 throw ErrorException(__FUNCT__,exprintf("!vertices || nbv>=nbvx"));
+									 }
+									 if (!edges || nbe>=nbex){
+										 throw ErrorException(__FUNCT__,exprintf("!edges || nbe>=nbex"));
+									 }
+									 if (!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex){
+										 throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex"));
+									 }
 									 // new vertex on edge
 									 A1=vertices+nbv++;
@@ -3650,7 +3645,11 @@
 									 Edge *e = edges + nbe++;
 									 Real8 se= (sNew-L0)/LAB;
-									 assert(se>=0 && se < 1.000000001);
+									 if (se<0 || se>=1.000000001){
+										 throw ErrorException(__FUNCT__,exprintf("se<0 || se>=1.000000001"));
+									 }
 									 se =  abscisseInterpole(v0.m,v1.m,AB,se,1);
-									 assert(se>=0 && se <= 1);
+									 if (se<0 || se>1){
+										 throw ErrorException(__FUNCT__,exprintf("se<0 || se>1"));
+									 }
 									 //((k1==1) != (k1==k1equi))
 									 se = k1 ? se : 1. - se;
@@ -3683,11 +3682,17 @@
 
 							 }               
-							 assert(ee.on->CurveNumber==ei.on->CurveNumber);
+							 if (ee.on->CurveNumber!=ei.on->CurveNumber){
+								 throw ErrorException(__FUNCT__,exprintf("ee.on->CurveNumber!=ei.on->CurveNumber"));
+							 }
 							 if(verbosity>98) cout <<  BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() <<  endl;
 							 if ( ee[k1].on->IsRequiredVertex()) {
-								 assert(eeequi[k1equi].on->IsRequiredVertex());
+								 if (!eeequi[k1equi].on->IsRequiredVertex()){
+									 throw ErrorException(__FUNCT__,exprintf("!eeequi[k1equi].on->IsRequiredVertex()"));
+								 }
 								 register GeometricalVertex * GA1 = *eeequi[k1equi].on;
 								 A1=GA1->to;// the vertex in new mesh
-								 assert (A1-vertices>=0 && A1-vertices <nbv);
+								 if (A1-vertices<0 || A1-vertices>=nbv){
+									 throw ErrorException(__FUNCT__,exprintf("A1-vertices<0 || A1-vertices>=nbv"));
+								 }
 								 break;}
 								 if (!ee.adj[k1]) {
@@ -3717,5 +3722,7 @@
 							PreviousNewEdge = e;
 
-							assert(i==NbCreatePointOnCurve);
+							if (i!=NbCreatePointOnCurve){
+								throw ErrorException(__FUNCT__,exprintf("i!=NbCreatePointOnCurve"));
+							}
 
 						  }
@@ -3763,6 +3770,7 @@
 		  }
 		  } // for(step;;)
-		assert(nbe);
-
+		if (nbe==0){
+			throw ErrorException(__FUNCT__,exprintf("nbe==0"));
+		}
 		delete [] bcurve;
 
@@ -3809,5 +3817,4 @@
 			 nbv++;
 		 }
-		//  assert( Gh.nbv < nbvx);
 
 		// Method in 2 step:  0 and 1 
@@ -3815,5 +3822,7 @@
 		// 2) construct the edge    
 		// generation of the curves
-		assert(! edges);
+		if (edges){
+			throw ErrorException(__FUNCT__,exprintf("edges"));
+		}
 		// 2 step 
 		// --step=0 to compute the number of edges + alloc at end
@@ -3849,5 +3858,7 @@
 								a=ei(0)->The();
 								b=ei(1)->The();
-								assert(edges);
+								if (!edges){
+									throw ErrorException(__FUNCT__,exprintf("!edges"));
+								}
 								edges[nbe].v[0]=a->to;
 								edges[nbe].v[1]=b->to;;
@@ -3868,5 +3879,7 @@
 							  NbNewPoints=0;
 							  NbEdgeCurve=0;
-							  assert(nbvend < nbvx); 
+							  if (nbvend>=nbvx){
+								  throw ErrorException(__FUNCT__,exprintf("nbvend>=nbvx"));
+							  }
 							  lcurve =0;
 							  s = lstep;
@@ -3934,5 +3947,7 @@
 										  else
 											kk0=kkk,ll0=llk;}
-										  assert(kk1 != kk0);
+										  if (kk1 == kk0){
+											  throw ErrorException(__FUNCT__,exprintf("kk1 == kk0"));
+										  }
 
 										  Real8 sbb = (ss-ll0  )/(ll1-ll0);
@@ -3972,5 +3987,7 @@
 								  k = e->SensAdj[kprev];// next vertices
 								  e = e->Adj[kprev];
-								  assert(e);
+								  if (!e){
+									  throw ErrorException(__FUNCT__,exprintf("!e"));
+								  }
 								 }// for(;;)
 							  vb = b->to;
@@ -4006,20 +4023,22 @@
 			} // for (i=0;i<nbe;i++)
 			if(!step) {
-				// cout << "edges " << edges << " VerticesOnGeomEdge " <<VerticesOnGeomEdge << endl;
-				assert(!edges);
-				assert(!VerticesOnGeomEdge);
+				if (edges){
+					throw ErrorException(__FUNCT__,exprintf("edges"));
+				}
+				if (VerticesOnGeomEdge){
+					throw ErrorException(__FUNCT__,exprintf("VerticesOnGeomEdge"));
+				}
 				edges = new Edge[nbex=nbe];
 				if(NbVerticesOnGeomEdge0)
 				 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
-				assert(edges);
-				assert(VerticesOnGeomEdge || NbVerticesOnGeomEdge0 ==0);
+				if (!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0){
+					throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0"));
+				}
 				// do the vertex on a geometrical vertex
 				NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;       
 			}
-			else 
-			 assert(NbVerticesOnGeomEdge == NbVerticesOnGeomEdge0);
-			//     cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe 
-			//	  << "== " << nbex << "  nbv = " << nbv <<  endl;
-			assert(nbex=nbe);
+			else if (NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0){
+				throw ErrorException(__FUNCT__,exprintf("NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0"));
+			}
 		  } // for (step=0;step<2;step++)
 
@@ -4035,5 +4054,7 @@
 	/*FUNCTION Triangles::MakeGeometricalEdgeToEdge{{{1*/
 	Edge** Triangles::MakeGeometricalEdgeToEdge() {
-		assert(Gh.nbe);
+		if (!Gh.nbe){
+			throw ErrorException(__FUNCT__,exprintf("!Gh.nbe"));
+		}
 		Edge **e= new (Edge* [Gh.nbe]);
 
@@ -4072,5 +4093,4 @@
 			  cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber 
 				 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1  " << Gh.Number( Gh.edges[i][1]) << endl; 
-			  //	 assert( e[i]);
 		  }
 		if(kk) throw ErrorException(__FUNCT__,exprintf("See above"));
@@ -4096,5 +4116,7 @@
 		pmax = pmax+DD; 
 		coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
-		assert(coefIcoor >0);
+		if (coefIcoor<=0){
+			throw ErrorException(__FUNCT__,exprintf("coefIcoor<=0"));
+		}
 
 		// generation of integer coord  
@@ -4142,5 +4164,7 @@
 
 			//initialize ordre
-			assert(ordre);
+			if (!ordre){
+				throw ErrorException(__FUNCT__,exprintf("!ordre"));
+			}
 			for (i=0;i<nbv;i++) ordre[i]=0;
 
@@ -4174,5 +4198,7 @@
 					}
 					else if(st[k]>=0) {
-						assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)));
+						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)))"));
+						}
 
 						triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
@@ -4329,5 +4355,7 @@
 								 Vertex *v1= ta.EdgeVertex(1);
 								 Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
-								 assert(st[k] >=0); 
+								 if (st[k]<0){
+									 throw ErrorException(__FUNCT__,exprintf("st[k]<0"));
+								 }
 								 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
 								 ta.SetLock();
@@ -4348,6 +4376,7 @@
 						}
 				 }
-				 // cout <<      savenbt+NbTfillHoll << " " <<  savenbtx  << endl;
-				 assert(savenbt+NbTfillHoll <= savenbtx );
+				 if (savenbt+NbTfillHoll>savenbtx ){
+					 throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx"));
+				 }
 				 // copy of the outside triangles in saveTriangles 
 				 for (i=0;i<nbt;i++){
@@ -4442,5 +4471,8 @@
 			  {
 				NbSwap++;
-				assert(k++<20000);
+				k++;
+				if (k>=20000){
+					throw ErrorException(__FUNCT__,exprintf("k>=20000"));
+				}
 				t=  tp->at[jp];      // set unchange t qnd j for previous triangles
 				j=  NextEdge[tp->aa[jp]&3];
@@ -4641,7 +4673,8 @@
 	/*FUNCTION Triangles::Crack{{{1*/
 	int  Triangles::Crack() { 
-		assert(NbCrackedEdges ==0 || NbCrackedVertices >0); 
-		for (int i=0;i<NbCrackedEdges;i++)
-		 CrackedEdges[i].Crack();
+		if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{
+			throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges!=0 && NbCrackedVertices<=0"));
+		}
+		for (int i=0;i<NbCrackedEdges;i++) CrackedEdges[i].Crack();
 		return NbCrackedEdges;
 	}
@@ -4649,5 +4682,7 @@
 	/*FUNCTION Triangles::UnCrack{{{1*/
 	int Triangles::UnCrack() { 
-		assert(NbCrackedEdges ==0 || NbCrackedVertices >0); 
+		if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{
+			throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges ==0 || NbCrackedVertices >0"));
+		}
 		for (int i=0;i<NbCrackedEdges;i++)
 		 CrackedEdges[i].UnCrack();
@@ -4703,5 +4738,7 @@
 			Triangle * tbegin = v.t;
 			int i  = v.vint;       
-			assert(tbegin && (i >= 0 ) && (i <3));
+			if (!tbegin || (i<0) || (i>=3)){
+				throw ErrorException(__FUNCT__,exprintf("!tbegin || (i<0) || (i>=3)"));
+			}
 			// turn around the vertex v
 			TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge
@@ -4714,5 +4751,7 @@
 				  {   
 					TriangleAdjacent tta=(ta.Adj());
-					assert(tta.Cracked());
+					if (!tta.Cracked()){
+						throw ErrorException(__FUNCT__,exprintf("!tta.Cracked()"));
+					}
 					if ( kk == 0) tbegin=ta,kkk=0;  //  begin by a cracked edge  => restart                
 					if (  kkk ) { kc =1;vv = vlast++;  kkk = 0; } // new vertex if use 
@@ -4721,5 +4760,7 @@
 				if ( tt->link ) { // if good triangles store the value 
 					int it = Number(tt);
-					assert(it < nt);
+					if (it>=nt){
+						throw ErrorException(__FUNCT__,exprintf("(it>=nt)"));
+					}
 					(*tt)(kv)= vv; //   Change the vertex of triangle 
 					if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref 
@@ -4732,5 +4773,7 @@
 				ta = Next(ta).Adj(); 
 			} while ( (tbegin != ta)); 
-			assert(k);
+			if (!k){
+				throw ErrorException(__FUNCT__,exprintf("!k"));
+			}
 			if (kc)  nbcrakev++;
 		  }
@@ -4799,9 +4842,9 @@
 		Triangle * t=0;	
 		int j,jp,jn,jj;
-		if (tstart) 
-		 t=tstart;
-		else 
-		  {
-			assert(quadtree);
+		if (tstart) t=tstart;
+		else {
+			if (!quadtree){
+				throw ErrorException(__FUNCT__,exprintf("!quadtree"));
+			}
 			Vertex *a = quadtree->NearestVertex(B.x,B.y) ;
 
@@ -4812,17 +4855,20 @@
 				throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContening"));
 			}
-			assert(a>= vertices && a < vertices+nbv);
-			//  int k=0;
+			if (a<vertices || a>=vertices+nbv){
+				throw ErrorException(__FUNCT__,exprintf("a<vertices || a>=vertices+nbv"));
+			}
 			t = a->t;
-			assert(t>= triangles && t < triangles+nbt);
-
-		  }
+			if (t<triangles || t>=triangles+nbt){
+				throw ErrorException(__FUNCT__,exprintf("t<triangles || t>=triangles+nbt"));
+			}
+		}
 		Icoor2  detop ;
 		int kkkk =0; // number of test triangle 
 
-		while ( t->det < 0) 
-		  { // the initial triangles is outside  
+		while ( t->det < 0){ // the initial triangles is outside  
 			int k0=(*t)(0) ?  ((  (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1  )) : 0;
-			assert(k0>=0); // k0 the NULL  vertex 
+			if (k0<0){ // k0 the NULL  vertex
+				throw ErrorException(__FUNCT__,exprintf("k0<0"));
+			}
 			int k1=NextVertex[k0],k2=PreviousVertex[k0];
 			dete[k0]=det(B,(*t)[k1],(*t)[k2]);
@@ -4831,5 +4877,8 @@
 			 return t; 
 			t = t->TriangleAdj(OppositeEdge[k0]);
-			assert(kkkk++ < 2);
+			kkkk++;
+			if (kkkk>=2){
+				throw ErrorException(__FUNCT__,exprintf("kkkk>=2"));
+			}
 		  }
 
@@ -4837,7 +4886,9 @@
 		detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
 
-		while(t->det  > 0 ) 
-		  { 
-			assert( kkkk++ < 2000 ); 
+		while(t->det  > 0 ) { 
+			kkkk++;
+			if (kkkk>=2000){
+				throw ErrorException(__FUNCT__,exprintf("kkkk>=2000"));
+			}
 			j= OppositeVertex[jj];
 			dete[j] = detop;  //det(*b,*s1,*s2);
@@ -4856,9 +4907,9 @@
 			// 2 => two way go in way 1 or 2 randomly
 
-			if (k==0) 
-			 break;
-			if (k == 2 && BinaryRand())
-			 Exchange(ii[0],ii[1]);
-			assert ( k  < 3);
+			if (k==0) break;
+			if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]);
+			if ( k>=3){
+				throw ErrorException(__FUNCT__,exprintf("k>=3"));
+			}
 			TriangleAdjacent t1 = t->Adj(jj=ii[0]);
 			if ((t1.det() < 0 ) && (k == 2))
@@ -4887,5 +4938,7 @@
 		Real8 hmin = Gh.MinimalHmin();
 		Real8 maxaniso = 1e6;
-		assert(hmax>0);
+		if (hmax<=0){
+			throw ErrorException(__FUNCT__,exprintf("hmax<=0"));
+		}
 		SetVertexFieldOn();
 		if (errC > 1) errC = 1;
@@ -4902,17 +4955,13 @@
 				 Real8 s = GV;
 				 R2 tg;
-				 //	   cerr << i << " " << j << " " << Number(V) << " on = " 
-				 //	<< Gh.Number(eg) << " at s = " << s << " " << endl;
 				 Real8  R1= eg->R1tg(s,tg);
-				 // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x " 
-				 //    << V.r << errC/ Max(R1,1e-20) <<  " hold=" <<V.m(tg) << " "  << endl;
 				 Real8 ht = hmax;
-				 if (R1>1.0e-20) 
-					{  // err relative to the length of the edge
+				 if (R1>1.0e-20) {  // err relative to the length of the edge
 					 ht = Min(Max(errC/R1,hmin),hmax);
 					}
 				 Real8 hn = iso? ht : Min(hmax,ht*maxaniso);
-				 //cerr << ht << " " << hn << "m=" << edges[i][j].m <<  endl;
-				 assert(ht>0 && hn>0);
+				 if (ht<=0 || hn<=0){
+					 throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
+				 }
 				 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
 				 //cerr << " : " ;
@@ -5217,13 +5266,7 @@
 						 Real8 cAB   =  det3x3(taa[0],taa[1],bb);
 
-						 assert(det33);
-						 //	det33=1;
-						 // verif
-						 //	cout << " " << (taa[0][0]*cBC +  taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ;
-						 //	cout << " " << (taa[0][1]*cBC +  taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1];
-						 //	cout << " " << (taa[0][2]*cBC +  taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2] 
-						 //	     << "  -- " ;
-						 //cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB +  llb*llc* cBC + llc*lla*cCA)/det33 
-						 //   << " == " << llf <<  endl;
+						 if (!det33){
+							 throw ErrorException(__FUNCT__,exprintf("!det33"));
+						 }
 						 // computation of the gradient in the element 
 
@@ -5514,22 +5557,25 @@
 		while (Head0>=0&& kk++<100) {
 			kch=0;
-			for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1)
-			  {  //  pour tous les triangles autour du sommet s
+			for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
+				//  pour tous les triangles autour du sommet s
 				// 	cout << kk << " i = " << i << " " << ip << endl;
-				register Triangle * t= vertices[i].t;
-				assert(t);
+				register Triangle* t= vertices[i].t;
+				if (!t){
+					throw ErrorException(__FUNCT__,exprintf("!t"));
+				}
 				Vertex & vi = vertices[i];
 				TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
 				Vertex *pvj0 = ta.EdgeVertex(0);
 				while (1) {
-					//	  cout << i << " " <<  Number(ta.EdgeVertex(0)) << " "
-					//      << Number(ta.EdgeVertex(1)) << "  ---> " ;
 					ta=Previous(Adj(ta));
-					// cout <<  Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl;
-					assert(vertices+i == ta.EdgeVertex(1));
+					if (vertices+i != ta.EdgeVertex(1)){
+						throw ErrorException(__FUNCT__,exprintf("vertices+i != ta.EdgeVertex(1)"));
+					}
 					Vertex & vj = *(ta.EdgeVertex(0));
 					if ( &vj ) {
 						j= &vj-vertices;
-						assert(j>=0 && j < nbv);
+						if (j<0 || j >= nbv){
+							throw ErrorException(__FUNCT__,exprintf("j<0 || j >= nbv"));
+						}
 						R2 Aij = (R2) vj - (R2) vi;
 						Real8 ll =  Norme2(Aij);
@@ -5611,9 +5657,16 @@
 		if (inbvx) {
 			vertices=new Vertex[nbvx];
-			assert(vertices);
+			if (!vertices){
+				throw ErrorException(__FUNCT__,exprintf("!vertices"));
+			}
 			ordre=new (Vertex* [nbvx]);
-			assert(ordre);
+			if (!ordre){
+				throw ErrorException(__FUNCT__,exprintf("!ordre"));
+			}
 			triangles=new Triangle[nbtx];
-			assert(triangles);}
+			if (!triangles){
+				throw ErrorException(__FUNCT__,exprintf("!triangles"));
+			}
+		}
 		else {
 			vertices=0;
@@ -5622,6 +5675,5 @@
 			nbtx=0;
 		}
-		if ( name || inbvx)
-		  {
+		if ( name || inbvx) {
 			time_t timer =time(0);
 			char buf[70];     
@@ -5713,5 +5765,4 @@
 	Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside){
 		CurrentTh=&Th;
-		assert(&Th);
 		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];
@@ -5738,5 +5789,7 @@
 				k = ta;
 				Exchange(aa,bb);
-				assert(tc->link);
+				if (!tc->link){
+					throw ErrorException(__FUNCT__,exprintf("!tc->link"));
+				}
 			  }
 			a[VerticesOfTriangularEdge[k][0]] = aa;
@@ -5751,12 +5804,15 @@
 		int k=(*t)(0) ?  ((  (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1  )) : 0;
 		int dir=0;
-		assert(k>=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)))) 
-		  {  
-
-			assert(kkk++<1000);      
+		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);
@@ -5774,5 +5830,7 @@
 							continue;}}// go in direction j
 							double IJ2 = IJ_IA + IJ_AJ;
-							assert(IJ2);
+							if (IJ2==0){
+								throw ErrorException(__FUNCT__,exprintf("IJ2==0"));
+							}
 							a= IJ_AJ/IJ2;
 							b= IJ_IA/IJ2;
@@ -5790,5 +5848,7 @@
 		// restart:
 		//   int dir=0;
-		assert(t->link == 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 
@@ -5864,5 +5924,7 @@
 			  }
 		  }
-		assert(cas !=-2);
+		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
@@ -5877,7 +5939,9 @@
 
 			Triangle * tt=t=edge=Adj(Previous(edge));
-			do  {  // loop around vertex s
-
-				assert(edge.EdgeVertex(0)==s && kkk++<10000);
+			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;
@@ -5909,5 +5973,7 @@
 			} while (t!=tt);
 
-			assert((Triangle *) er);
+			if (!(Triangle *) er){
+				throw ErrorException(__FUNCT__,exprintf("!(Triangle *) er"));
+			}
 			I2 A((I2)*er.EdgeVertex(0));
 			I2 B((I2)*er.EdgeVertex(1));
@@ -5980,5 +6046,7 @@
 		Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
 		Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
-		assert ( a1 >= 0 && a1 < 3 );
+		if ( a1<0 || a1>=3 ){
+			throw ErrorException(__FUNCT__,exprintf("a1<0 || a1>=3"));
+		}
 
 		Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
@@ -5990,6 +6058,10 @@
 		Icoor2 det1=t1->det , det2=t2->det ;
 		Icoor2 detT = det1+det2;
-		assert((det1>0 ) && (det2 > 0));
-		assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb
+		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;
@@ -6052,5 +6124,7 @@
 	int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)  { 
 		int NbSwap =0;
-		assert(a.t && b.t); // the 2 vertex is in a mesh 
+		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 
@@ -6067,5 +6141,7 @@
 			v2 = tta.EdgeVertex(0);
 			vbegin =v2;
-			assert(v2);
+			if (!v2){
+				throw ErrorException(__FUNCT__,exprintf("!v2"));
+			}
 			det2 = det(*v2,a,b);
 			//   cout << " No Change try the next" << endl;
@@ -6083,5 +6159,7 @@
 				Vertex * va = &a, *vb = &b;
 				tc = Previous(tc);
-				assert ( v1 && v2);
+				if (!v1 || !v2){
+					throw ErrorException(__FUNCT__,exprintf("!v1 || !v2"));
+				}
 				Icoor2 detss = 0,l=0,ks;
 				// cout << "Real ForcingEdge " << *va << *vb << detss << endl;
@@ -6105,5 +6183,8 @@
 			}
 			tta = tc;
-			assert(k++<2000);
+			k++;
+			if (k>=2000){
+				throw ErrorException(__FUNCT__,exprintf("k>=2000"));
+			}
 			if ( vbegin == v2 ) return -1;// error 
 		}
Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2853)
+++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp	(revision 2854)
@@ -28,6 +28,4 @@
 		register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc;
 		R2 P(s->r),PNew(0,0);
-		//  cout << BTh.quadtree << " " <<  BTh.quadtree->root << endl;
-		// assert(BTh.quadtree && BTh.quadtree->root);
 		do {
 			k++; 
@@ -46,5 +44,7 @@
 			tria = ttc;
 			j = NextEdge[jc];
-			assert(k<2000);
+			if (k>=2000){
+				throw ErrorException(__FUNCT__,exprintf("k>=2000"));
+			}
 		} while ( tbegin != tria); 
 		if (kk<4) return 0;
@@ -117,5 +117,7 @@
 				tria = ttc;
 				j = NextEdge[jc];
-				assert(k<2000);
+				if (k>=2000){
+					throw ErrorException(__FUNCT__,exprintf("k>=2000"));
+				}
 			} while ( tbegin != tria); 
 			if (ok && loop) vP=vPsave; // no move