Index: issm/trunk/src/c/bamg/Mesh.cpp
===================================================================
--- issm/trunk/src/c/bamg/Mesh.cpp	(revision 15396)
+++ issm/trunk/src/c/bamg/Mesh.cpp	(revision 16560)
@@ -122,15 +122,15 @@
 		  }
 		  for (i=0;i<Tho.nbt;i++)
-			if(  reft[i] >=0 && flag[i]) 
-			  {
+			if(reft[i] >=0 && flag[i]){
 				const Triangle & t = Tho.triangles[i];
 				int i0 = Tho.GetId(t[0]);
 				int i1 = Tho.GetId(t[1]);
 				int i2 = Tho.GetId(t[2]);
-				if (i0<0 || i1<0 || i2<0){
+				if(i0<0 || i1<0 || i2<0){
 					delete [] refv;
 					_error_("i0<0 || i1<0 || i2< 0");
 				}
-				if (i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv){
+				if(i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv){
+					delete [] refv;
 					_error_("i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv");
 				}
@@ -139,7 +139,4 @@
 				nbt++;           
 			  }
-		  if (kt!=nbt){
-			  _error_("kt!=nbt");
-		  }
 		  if (nbt==0 && nbv==0) {
 			  _error_("All triangles have been removed");
@@ -148,5 +145,4 @@
 		  delete [] reft;
 		  delete [] refv;
-		  //double cutoffradian = 10.0/180.0*Pi;
 		  BuildGeometryFromMesh(bamgopts);
 		  Gh.PostRead(); 
@@ -154,11 +150,8 @@
 		  ReconstructExistingMesh();
 
-		  if (!nbsubdomains){
-			  _error_("nbsubdomains==0");
-		  }
-		  if (!subdomains[0].head || !subdomains[0].head->link){
-			  _error_("!subdomains[0].head || !subdomains[0].head->link");
-		  }
-
+		  /*Final checks*/
+		  _assert_(kt==nbt);
+		  _assert_(nbsubdomains);
+		  _assert_(subdomains[0].head && subdomains[0].head->link); 
 	  }
 	/*}}}*/
@@ -269,5 +262,4 @@
 	void Mesh::ReadMesh(int* index,double* x,double* y,int nods,int nels){
 
-		double Hmin = HUGE_VAL;// the infinie value 
 		long i1,i2,i3;
 		long i;
@@ -523,7 +515,9 @@
 				head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing
 				direction=(int)bamgmesh->SubDomains[i*3+2];
-				if (i3!=23) _error_("Bad Subdomain definition: first number should be 3");
+				if (i3!=3) _error_("Bad Subdomain definition: first number should be 3");
 				if (head<0 || head>=nbt) _error_("Bad Subdomain definition: head should in [1 " << nbt << "] (triangle number)");
 				subdomains[i].head = triangles+head;
+				subdomains[i].direction = direction;
+				subdomains[i].ReferenceNumber = i3;
 			}
 		}
@@ -1030,4 +1024,5 @@
 				Triangle &t=triangles[i];
 				if (t.det>0 && !(t.Hidden(0)||t.Hidden(1) || t.Hidden(2) )){
+					/*Remove triangles that have a bad aspect ratio*/
 					//if(t.Anisotropy()<2 & t.Length()<1.e+5){
 						index[num*3+0]=GetId(t[0])+1; //back to M indexing
@@ -1052,7 +1047,6 @@
 
 		/*Get options*/
-		int    verbose=bamgopts->verbose;
-		double anisomax =bamgopts->anisomax;
-		double errg     =bamgopts->errg;
+		double anisomax = bamgopts->anisomax;
+		double errg     = bamgopts->errg;
 
 		double ss[2]={0.00001,0.99999};
@@ -1062,10 +1056,8 @@
 
 		//check that hmax is positive
-		if (hmax<=0){
-			_error_("hmax<=0");
-		}
+		_assert_(hmax>0); 
 
 		//errC cannot be higher than 1
-		if (errC>1) errC=1;
+		if(errC>1) errC=1;
 
 		//Set all vertices to "on"
@@ -1073,5 +1065,5 @@
 
 		//loop over all the vertices on edges
-		for (int  i=0;i<nbe;i++){
+		for(int  i=0;i<nbe;i++){
 			for (int j=0;j<2;j++){
 
@@ -1160,5 +1152,5 @@
 
 		//some checks
-		if (( infvertexindex <0 ) && (detOld <0) ||  ( infvertexindex >=0  ) && (detOld >0) ){
+		if(((infvertexindex <0 ) && (detOld <0)) ||  ((infvertexindex >=0) && (detOld >0)) ){
 			_error_("inconsistent configuration (Contact ISSM developers)");
 		}
@@ -1727,14 +1719,11 @@
 			long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
 			k = edge4->SortAndFind(i0,i1);
-			if(k>=0){
-				subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
-				subdomains[i].edge = edges+k;
-				Gh.subdomains[i].edge = Gh.edges + k;
-				Gh.subdomains[i].direction  =  subdomains[i].direction;
-				Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
-			}
-			else
-			 _error_("%i should be >=0");
-		  }
+			_assert_(k>=0);
+			subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
+			subdomains[i].edge = edges+k;
+			Gh.subdomains[i].edge = Gh.edges + k;
+			Gh.subdomains[i].direction  =  subdomains[i].direction;
+			Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
+		}
 
 		delete edge4;
@@ -2187,8 +2176,10 @@
 				// correction of second derivative
 				// by a laplacien
-				double* d2[3] = {dxdx, dxdy, dydy};
 				double* dd;
 				for (int xy = 0;xy<3;xy++) {
-					dd = d2[xy];
+					if      (xy==0) dd=dxdx;
+					else if (xy==1) dd=dxdy;
+					else if (xy==2) dd=dydy;
+					else    _error_("not supported yet");
 					// do leat 2 iteration for boundary problem
 					for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
@@ -2620,6 +2611,5 @@
 
 			  }
-			else
-			  { // find the head for all sub domaine
+			else{ // find the head for all subdomains
 				if (Gh.nbsubdomains != nbsubdomains && subdomains)
 				 delete [] subdomains, subdomains=0;
@@ -2627,5 +2617,4 @@
 				 subdomains = new SubDomain[ Gh.nbsubdomains];
 				nbsubdomains =Gh.nbsubdomains;
-				long err=0;
 				CreateSingleVertexToTriangleConnectivity();
 				long * mark = new long[nbt];
@@ -2638,5 +2627,4 @@
 					GeomEdge &eg = *Gh.subdomains[i].edge;
 					subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber;
-					int ssdlab = subdomains[i].ReferenceNumber;
 					// by carefull is not easy to find a edge create from a GeomEdge 
 					// see routine MakeGeomEdgeToEdge
@@ -2773,5 +2761,5 @@
 			vertices=new BamgVertex[maxnbv];
 			_assert_(vertices);
-			orderedvertices=new (BamgVertex* [maxnbv]);
+			orderedvertices=new BamgVertex* [maxnbv];
 			_assert_(orderedvertices);
 			triangles=new Triangle[maxnbt];
@@ -3015,9 +3003,4 @@
 	}
 	/*}}}*/
-	/*FUNCTION Mesh::isCracked{{{*/
-	int Mesh::isCracked() const {
-		return NbCrackedVertices != 0;
-	}
-	/*}}}*/
 	/*FUNCTION Mesh::MakeGeomEdgeToEdge{{{*/
 	Edge** Mesh::MakeGeomEdgeToEdge() {
@@ -3027,5 +3010,5 @@
 			_error_("!Gh.nbe");
 		}
-		Edge **e= new (Edge* [Gh.nbe]);
+		Edge **e= new Edge* [Gh.nbe];
 
 		long i;
@@ -3108,8 +3091,5 @@
 	void Mesh::MakeBamgQuadtree() {  
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/
-
-		long int verbose=0;
-		if (  !quadtree )  quadtree = new BamgQuadtree(this);
-
+		if(!quadtree) quadtree = new BamgQuadtree(this);
 	}
 	/*}}}*/
@@ -3135,5 +3115,5 @@
 			for (int j=0;j<3;j++){
 				Triangle &tt = *t.TriangleAdj(j);
-				if ( ! &tt ||  it < GetId(tt) && ( tt.link || t.link)){
+				if ( (!&tt ||  it < GetId(tt)) && ( tt.link || t.link)){
 					BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
 					BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
@@ -3313,16 +3293,10 @@
 			}
 
-		} while (nbv!=nbvold);
-
-		delete []  first_np_or_next_t;
-
-		long NbSwapf =0,NbSwp;
-
-		NbSwp = NbSwapf;
-		for (i=0;i<nbv;i++)
-		 NbSwapf += vertices[i].Optim(0);
-		NbTSwap +=  NbSwapf ;
-	}
-	/*}}}*/
+		}while(nbv!=nbvold);
+		delete [] first_np_or_next_t;
+
+		long NbSwapf =0;
+		for(i=0;i<nbv;i++) NbSwapf += vertices[i].Optim(0);
+	}/*}}}*/
 	/*FUNCTION Mesh::ProjectOnCurve{{{*/
 	GeomEdge*   Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex &  vA, BamgVertex & vB,
@@ -3765,5 +3739,4 @@
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
 
-		long int verbose=0;
 		long *renu= new long[nbt];
 		register Triangle *t0,*t,*te=triangles+nbt;
@@ -3830,78 +3803,4 @@
 	}
 	/*}}}*/
-	/*FUNCTION Mesh::VerticesRenumber{{{*/
-	void Mesh::VerticesRenumber(long * renu) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
-
-		// warning be carfull because pointer
-		// from on mesh to over mesh 
-		//  --  so do ReNumbering at the beginning
-		BamgVertex * ve = vertices+nbv;
-		long it,ie,i;
-
-		_printf_("renumbering triangles\n");
-		for ( it=0;it<nbt;it++) 
-		 triangles[it].Renumbering(vertices,ve,renu);
-
-		_printf_("renumbering edges\n");
-		for ( ie=0;ie<nbe;ie++) 
-		 edges[ie].Renumbering(vertices,ve,renu);
-
-		_printf_("renumbering vertices on geom\n");
-		for (i=0;i< NbVerticesOnGeomVertex;i++)
-		  {
-			BamgVertex *v = VerticesOnGeomVertex[i].meshvertex;
-			if (v>=vertices && v < ve)
-			 VerticesOnGeomVertex[i].meshvertex=vertices+renu[GetId(v)];
-		  }
-
-		_printf_("renumbering vertices on edge\n");
-		for (i=0;i< NbVerticesOnGeomEdge;i++)
-		  {
-			BamgVertex *v =VerticesOnGeomEdge[i].meshvertex;
-			if (v>=vertices && v < ve)
-			 VerticesOnGeomEdge[i].meshvertex=vertices+renu[GetId(v)];
-		  }
-
-		_printf_("renumbering vertices on Bth vertex\n");
-		for (i=0;i< NbVertexOnBThVertex;i++)
-		  {
-			BamgVertex *v=VertexOnBThVertex[i].v;
-			if (v>=vertices && v < ve)
-			 VertexOnBThVertex[i].v=vertices+renu[GetId(v)];
-		  }
-
-		for (i=0;i< NbVertexOnBThEdge;i++)
-		  {
-			BamgVertex *v=VertexOnBThEdge[i].v;
-			if (v>=vertices && v < ve)
-			 VertexOnBThEdge[i].v=vertices+renu[GetId(v)];
-		  }
-
-		// move the Vertices without a copy of the array 
-		// be carefull not trivial code 
-		long j;
-		for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
-		 if (renu[it] >= 0) // a new sub cycle
-			{ 
-			 i=it;
-			 BamgVertex ti=vertices[i],tj;
-			 while ( (j=renu[i]) >= 0){
-				 // i is old, and j is new 
-				 renu[i] = -1-renu[i]; // mark 
-				 tj = vertices[j];     // save new
-				 vertices[j]= ti;      // new <- old
-				 i=j;     // next 
-				 ti = tj;
-				}  
-			}
-		if (quadtree){
-			delete quadtree;
-			quadtree = new BamgQuadtree(this);
-		}
-
-		for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
-	}
-	/*}}}*/
 /*FUNCTION Mesh::SetIntCoor{{{*/
 void Mesh::SetIntCoor(const char * strfrom) {
@@ -3964,113 +3863,4 @@
 
 	if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above");
-}
-/*}}}*/
-/*FUNCTION Mesh::ShowRegulaty{{{*/
-void  Mesh::ShowRegulaty() const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
-
-	const  double  sqrt32=sqrt(3.)*0.5; 
-	const double  aireKh=sqrt32*0.5;
-	D2  Beq(1,0),Heq(0.5,sqrt32);
-	D2xD2 Br(D2xD2(Beq,Heq).t());
-	D2xD2 B1r(Br.inv());
-	double gammamn=1e100,hmin=1e100;
-	double gammamx=0,hmax=0;
-	double beta=1e100;
-	double beta0=0;
-	double  alpha2=0;
-	double area=0,Marea=0;
-	// double cf= double(coefIcoor);
-	// double cf2= 6.*cf*cf;
-	int nt=0;
-	for (int it=0;it<nbt;it++)
-	 if ( triangles[it].link){
-		 Triangle &K=triangles[it];
-		 double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
-		 area+= area3;
-		 D2xD2 B_Kt(K[0],K[1],K[2]);
-		 D2xD2 B_K(B_Kt.t());
-		 D2xD2 B1K = Br*B_K.inv();
-		 D2xD2 BK =  B_K*B1r;
-		 D2xD2 B1B1 = B1K.t()*B1K;
-		 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
-		 EigenMetric VMK(MK);
-		 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
-		 double betaK=0;
-
-		 for (int j=0;j<3;j++)
-			{
-			 double he= Norme2(R2(K[j],K[(j+1)%3]));
-			 hmin=Min(hmin,he);
-			 hmax=Max(hmax,he);
-			 BamgVertex & v=K[j];
-			 D2xD2 M((Metric)v);
-			 betaK += sqrt(M.det());
-
-			 D2xD2 BMB = BK.t()*M*BK;
-			 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
-			 EigenMetric VM1(M1);
-			 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
-			 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);		
-			}
-		 betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
-		 Marea+= betaK;
-		 beta=min(beta,betaK);
-		 beta0=max(beta0,betaK);
-		}   
-	area*=3; 
-	gammamn=sqrt(gammamn);
-	gammamx=sqrt(gammamx);    
-	_printf_("   Adaptmesh info:\n");
-	_printf_("      number of triangles = " << nt << "\n");
-	_printf_("      hmin = " << hmin << ", hmax=" << hmax << "\n");
-	_printf_("      area = " << area << ", M area = " << Marea << ", M area/( |Khat| nt) = " <<  Marea/(aireKh*nt) << "\n");
-	_printf_("      infinite-regularity(?): min = " << gammamn << ", max = " << gammamx << "\n");
-	_printf_("      anisomax = " << pow(alpha2,0.5) << ", beta max = " << 1./pow(beta/aireKh,0.5) << ", min = " << 1./pow(beta0/aireKh,0.5) << "\n");
-}
-/*}}}*/
-/*FUNCTION Mesh::ShowHistogram{{{*/
-void  Mesh::ShowHistogram() const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
-
-	const long kmax=10;
-	const double llmin = 0.5,llmax=2;
-	const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
-	long histo[kmax+1];
-	long i,it,k, nbedges =0;
-	for (i=0;i<=kmax;i++) histo[i]=0;
-	for (it=0;it<nbt;it++)
-	 if ( triangles[it].link) 
-		{
-
-		 for (int j=0;j<3;j++)
-			{
-			 Triangle *ta = triangles[it].TriangleAdj(j);
-			 if ( !ta || !ta->link || GetId(ta) >= it) 
-				{ 
-				 BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
-				 BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
-				 if ( !& vP || !&vQ) continue;
-				 R2 PQ = vQ.r - vP.r;
-				 double l = log(LengthInterpole(vP,vQ,PQ));
-				 nbedges++;
-				 k = (int) ((l - lmin)*delta);
-				 k = Min(Max(k,0L),kmax);
-				 histo[k]++;
-				}
-			}
-		}  
-	_printf_(" --- Histogram of the unit mesh,  nb of edges = " << nbedges << "\n");
-	_printf_("      length of edge in   | %% of edge  | Nb of edges \n"); 
-	_printf_("      --------------------+-------------+-------------\n"); 
-	for   (i=0;i<=kmax;i++){ 
-		if (i==0) _printf_( "      " << setw(10) << 0.);
-		else      _printf_( "      " << setw(10) << exp(lmin+i/delta));
-		if (i==kmax) _printf_("          +inf   ");
-		else      _printf_( "      " << setw(10) << exp(lmin+(i+1)/delta));
-		_printf_("|  " << setw(10) << (long((10000. * histo[i])/ nbedges)/100.) << " |\n");
-		_printf_("  " << histo[i] << "\n");
-	}
-	_printf_("      --------------------+-------------+-------------\n"); 
 }
 /*}}}*/
@@ -4138,7 +3928,7 @@
 	 first_np_or_next_t1[i]=-1;
 	kk=0;
-	while (Head0>=0&& kk++<100) {
+	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) {
+		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
 			register Triangle* t= vertices[i].t;
@@ -4195,510 +3985,4 @@
 }
 /*}}}*/
-	/*FUNCTION Mesh::SplitElement{{{*/
-	int  Mesh::SplitElement(int choice){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
-
-		long int verbose=0;
-
-		Direction NoDirOfSearch;
-		const  int withBackground = &BTh != this && &BTh;
-
-		TrianglesRenumberBySubDomain();
-		//int nswap =0;
-		const long nfortria( choice ? 4 : 6);
-		if(withBackground) 
-		  {
-			BTh.SetVertexFieldOn();
-			SetVertexFieldOnBTh();
-		  }
-		else
-		 BTh.SetVertexFieldOn();
-
-		long newnbt=0,newnbv=0;
-		long * kedge = 0;
-		long newnbq=nbq;
-		long * ksplit= 0, * ksplitarray=0;
-		long kkk=0;
-		int ret =0;
-		if (maxnbv<nbv+nbe) return 1;//   
-		// 1) create  the new points by spliting the internal edges 
-		// set the 
-		long nbvold = nbv;
-		long nbtold = nbt;
-		long nbtoutold  = nbtout;
-		long  NbEdgeOnGeom=0;
-		long i;
-
-		nbt = nbt - nbtout; // remove all the  the ouside triangles 
-		long nbtsave = nbt;
-		Triangle * lastT = triangles + nbt;
-		for (i=0;i<nbe;i++)
-		 if(edges[i].GeomEdgeHook) NbEdgeOnGeom++;
-		long newnbe=nbe+nbe;
-		//  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
-		long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
-		// long newNbVertexOnBThVertex=NbVertexOnBThVertex;
-		long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
-
-		// do allocation for pointeur to the geometry and background
-		VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
-		VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
-		if (NbVerticesOnGeomEdge)
-		 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
-		if (NbVertexOnBThEdge)
-		 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
-		Edge *newedges = new Edge [newnbe];
-		//  memcpy(newedges,edges,sizeof(Edge)*nbe);
-		SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
-		long k=nbv;
-		long kk=0;
-		long kvb = NbVertexOnBThEdge;
-		long kvg = NbVerticesOnGeomEdge;
-		long ie =0;
-		Edge ** edgesGtoB=0;
-		if (withBackground)
-		 edgesGtoB= BTh.MakeGeomEdgeToEdge();
-		long ferr=0;
-		for (i=0;i<nbe;i++)
-		 newedges[ie].GeomEdgeHook=0;
-
-		for (i=0;i<nbe;i++)
-		  {
-			GeomEdge *ong =  edges[i].GeomEdgeHook;
-
-			newedges[ie]=edges[i];
-			newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
-			newedges[ie].adj[1]=newedges + ie +1;
-			R2 A = edges[i][0],B = edges[i][1];
-
-			kk += (i == edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1])));
-			if (ong) // a geometrical edges 
-			  { 
-				if (withBackground){
-					// walk on back ground mesh 
-					//  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 
-					// a faire -- difficile 
-					// the first PB is to now a background edge between the 2 vertices
-					if (!edgesGtoB){
-						_error_("!edgesGtoB");
-					}
-					ong= ProjectOnCurve(*edgesGtoB[Gh.GetId(edges[i].GeomEdgeHook)],
-								edges[i][0],edges[i][1],0.5,vertices[k],
-								newVertexOnBThEdge[kvb],
-								newVerticesOnGeomEdge[kvg++]);
-					vertices[k].ReferenceNumber= edges[i].ReferenceNumber;
-					vertices[k].DirOfSearch =   NoDirOfSearch;        
-					;
-					// get the Info on background mesh 
-					double s =        newVertexOnBThEdge[kvb];
-					BamgVertex &  bv0  = newVertexOnBThEdge[kvb][0];
-					BamgVertex &  bv1  = newVertexOnBThEdge[kvb][1];
-					// compute the metrix of the new points 
-					vertices[k].m =  Metric(1-s,bv0,s,bv1); 
-					kvb++;
-				  }
-				else 
-				  {
-					ong=Gh.ProjectOnCurve(edges[i],
-								0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
-					// vertices[k].i = R2ToI2( vertices[k].r);
-					vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
-					vertices[k].DirOfSearch = NoDirOfSearch;
-					vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);	      
-				  }  
-			  }
-			else // straigth line edge ---
-			  { 
-				vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
-				vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
-				vertices[k].GeomEdgeHook = 0;
-			  }
-			//vertices[k].i = R2ToI2( vertices[k].r);
-			R2 AB =  vertices[k].r;
-			R2 AA = (A+AB)*0.5;
-			R2 BB = (AB+B)*0.5;
-			vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
-			vertices[k].DirOfSearch = NoDirOfSearch;
-
-			newedges[ie].GeomEdgeHook = Gh.Containing(AA,ong);
-			newedges[ie++].v[1]=vertices+k;
-
-			newedges[ie]=edges[i];
-			newedges[ie].adj[0]=newedges + ie -1;
-			newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
-			newedges[ie].GeomEdgeHook =  Gh.Containing(BB,ong);
-			newedges[ie++].v[0]=vertices+k;
-			k++;
-		  }
-		if (edgesGtoB) delete [] edgesGtoB;
-		edgesGtoB=0;
-
-		newnbv=k;
-		newNbVerticesOnGeomEdge=kvg;
-		if (newnbv> maxnbv) goto Error;// bug 
-
-		nbv = k;
-
-		kedge = new long[3*nbt+1];
-		ksplitarray = new long[nbt+1];
-		ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
-
-		for (i=0;i<3*nbt;i++)
-		 kedge[i]=-1;
-
-		//  
-
-		for (i=0;i<nbt;i++) {
-			Triangle & t = triangles[i];
-			if (!t.link){
-				_error_("!t.link");
-			}
-			for(int j=0;j<3;j++)
-			  {
-				const AdjacentTriangle ta = t.Adj(j);
-				const Triangle & tt = ta;
-				if (&tt >= lastT)
-				 t.SetAdj2(j,0,0);// unset adj
-				const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
-				const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				long  ke =edge4->SortAndFind(GetId(v0),GetId(v1));
-				if (ke>0) 
-				  {
-					long ii = GetId(tt);
-					int  jj = ta;
-					long ks = ke + nbvold;
-					kedge[3*i+j] = ks;
-					if (ii<nbt) // good triangle
-					 kedge[3*ii+jj] = ks;
-					BamgVertex &A=vertices[ks];
-					double aa,bb,cc,dd;
-					if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
-						// warning PB roundoff error 
-						if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0 
-										||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0  
-										||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
-							_printf_(ke + nbvold << " not in triangle " << i << " In= " << !!t.link << " " << aa << " " << bb << " " << cc << " " << dd << "\n");
-							_error_("Number of triangles with P2 interpolation Problem");
-						}
-					}
-					else {
-						if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0 
-										||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0 
-										||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
-							_printf_(ke + nbvold << " not in triangle " << ii << " In= " << !!tt.link << " " << aa << " " << bb << " " << cc << " " << dd << "\n");
-							_error_("Number of triangles with P2 interpolation Problem");
-						}
-					} 
-				  }
-			  }
-		}
-
-		for (i=0;i<nbt;i++){
-			ksplit[i]=1; // no split by default
-			const Triangle & t = triangles[ i];
-			int nbsplitedge =0;
-			int nbinvisible =0;
-			int invisibleedge=0;
-			int kkk[3];      
-			for (int j=0;j<3;j++)
-			  {
-				if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
-
-				const AdjacentTriangle ta = t.Adj(j);
-				const Triangle & tt = ta;
-
-				const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
-				const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				if ( kedge[3*i+j] < 0) 
-				  {
-					long  ke =edge4->SortAndFind(GetId(v0),GetId(v1));
-					if (ke<0) // new 
-					  {
-						if (&tt) // internal triangles all the boundary 
-						  { // new internal edges 
-							long ii = GetId(tt);
-							int  jj = ta;
-
-							kedge[3*i+j]=k;// save the vertex number 
-							kedge[3*ii+jj]=k;
-							if (k<maxnbv) 
-							  {
-								vertices[k].r = ((R2) v0+(R2) v1 )/2;
-								//vertices[k].i = R2ToI2( vertices[k].r);
-								vertices[k].ReferenceNumber=0;
-								vertices[k].DirOfSearch =NoDirOfSearch;
-								vertices[k].m =  Metric(0.5,v0,0.5,v1);
-							  }
-							k++;
-							kkk[nbsplitedge++]=j;		      
-						  } // tt 
-						else
-						 _error_("Bug...");
-					  } // ke<0	       
-					else
-					  { // ke >=0
-						kedge[3*i+j]=nbvold+ke;
-						kkk[nbsplitedge++]=j;// previously splited
-					  }
-				  }
-				else 
-				 kkk[nbsplitedge++]=j;// previously splited
-
-			  } 
-			if (nbinvisible>=2){
-				_error_("nbinvisible>=2");
-			}
-			switch (nbsplitedge) {
-				case 0: ksplit[i]=10; newnbt++; break;   // nosplit
-				case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 
-				case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 
-				case 3:
-						  if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
-						  else   ksplit[i]=10*nfortria,newnbt+=nfortria;
-						  break;
-			} 
-			if (ksplit[i]<40){
-				_error_("ksplit[i]<40");
-			}
-		  }
-		//  now do the element split
-		newnbq = 4*nbq;
-		nbv = k;
-		kkk = nbt;
-		ksplit[-1] = nbt;
-		// look on  old true  triangles 
-
-		for (i=0;i<nbtsave;i++){
-			int  nbmkadj=0;
-			long mkadj [100];
-			mkadj[0]=i;
-			long kk=ksplit[i]/10;
-			int  ke=(int) (ksplit[i]%10);
-			if (kk>=7 || kk<=0){
-				_error_("kk>=7 || kk<=0");
-			}
-
-			// def the numbering   k (edge) i vertex 
-			int k0 = ke;
-			int k1 = NextEdge[k0];
-			int k2 = PreviousEdge[k0];
-			int i0 = OppositeVertex[k0];
-			int i1 = OppositeVertex[k1];
-			int i2 = OppositeVertex[k2];
-
-			Triangle &t0=triangles[i];
-			BamgVertex * v0=t0(i0);           
-			BamgVertex * v1=t0(i1);           
-			BamgVertex * v2=t0(i2);
-
-			if (nbmkadj>=10){
-				_error_("nbmkadj>=10");
-			}
-			// --------------------------
-			AdjacentTriangle ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
-			// save the flag Hidden
-			int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
-			// un set all adj -- save Hidden flag --
-			t0.SetAdj2(0,0,hid[0]);
-			t0.SetAdj2(1,0,hid[1]);
-			t0.SetAdj2(2,0,hid[2]);
-			// --  remake 
-			switch  (kk) {
-				case 1: break;// nothing 
-				case 2: // 
-						  {
-							Triangle &t1=triangles[kkk++];
-							t1=t0;
-							if (kedge[3*i+i0]<0){
-								_error_("kedge[3*i+i0]<0");
-							}
-							BamgVertex * v3 = vertices + kedge[3*i+k0];
-
-							t0(i2) = v3;
-							t1(i1) = v3;
-							t0.SetAllFlag(k2,0);
-							t1.SetAllFlag(k1,0);
-						  } 
-						break; 
-				case 3: //
-						  {
-							Triangle &t1=triangles[kkk++];
-							Triangle &t2=triangles[kkk++];
-							t2=t1=t0;
-							if (kedge[3*i+k1]<0){
-								_error_("kedge[3*i+k1]<0");
-							}
-							if (kedge[3*i+k2]<0){
-								_error_("kedge[3*i+k2]<0");
-							}
-
-							BamgVertex * v01 = vertices + kedge[3*i+k2];
-							BamgVertex * v02 = vertices + kedge[3*i+k1]; 
-							t0(i1) = v01; 
-							t0(i2) = v02; 
-							t1(i2) = v02;
-							t1(i0) = v01; 
-							t2(i0) = v02; 
-							t0.SetAllFlag(k0,0);
-							t1.SetAllFlag(k1,0);
-							t1.SetAllFlag(k0,0);
-							t2.SetAllFlag(k2,0);
-						  } 
-						break;
-				case 4: // 
-				case 6: // split in 4 
-						  {
-							Triangle &t1=triangles[kkk++];
-							Triangle &t2=triangles[kkk++];
-							Triangle &t3=triangles[kkk++];
-							t3=t2=t1=t0;
-							if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
-								_error_("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
-							}
-							BamgVertex * v12 = vertices + kedge[3*i+k0];
-							BamgVertex * v02 = vertices + kedge[3*i+k1]; 
-							BamgVertex * v01 = vertices + kedge[3*i+k2];
-							t0(i1) = v01;
-							t0(i2) = v02;
-							t0.SetAllFlag(k0,hid[k0]);
-
-							t1(i0) = v01;
-							t1(i2) = v12;
-							t0.SetAllFlag(k1,hid[k1]);
-
-							t2(i0) = v02;
-							t2(i1) = v12;
-							t2.SetAllFlag(k2,hid[k2]);
-
-							t3(i0) = v12;
-							t3(i1) = v02;
-							t3(i2) = v01;
-
-							t3.SetAllFlag(0,hid[0]);	   
-							t3.SetAllFlag(1,hid[1]);	   
-							t3.SetAllFlag(2,hid[2]);
-
-							if ( kk == 6)
-							  {
-
-								Triangle &t4=triangles[kkk++];
-								Triangle &t5=triangles[kkk++];
-
-								t4 = t3;
-								t5 = t3;
-
-								t0.SetHidden(k0);
-								t1.SetHidden(k1);
-								t2.SetHidden(k2);
-								t3.SetHidden(0);
-								t4.SetHidden(1);
-								t5.SetHidden(2);
-
-								if (nbv < maxnbv ) 
-								  {
-									vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
-									vertices[nbv].ReferenceNumber =0;
-									vertices[nbv].DirOfSearch =NoDirOfSearch;
-									//vertices[nbv].i = R2ToI2(vertices[nbv].r);
-									double a3[]={1./3.,1./3.,1./3.};
-									vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
-									BamgVertex * vc =  vertices +nbv++;
-									t3(i0) = vc;
-									t4(i1) = vc;
-									t5(i2) = vc;
-
-								  }
-								else
-								 goto Error; 
-							  }
-
-						  } 
-						break;         
-			}
-
-			// save all the new triangles
-			mkadj[nbmkadj++]=i;
-			long jj;
-			if (t0.link) 
-			 for (jj=nbt;jj<kkk;jj++)
-				{
-				 triangles[jj].link=t0.link;
-				 t0.link= triangles+jj;
-				 mkadj[nbmkadj++]=jj;
-				}
-			if (nbmkadj>13){// 13 = 6 + 4 +
-				_error_("nbmkadj>13");
-			}
-
-			if (kk==6)  newnbq+=3;
-			for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
-			ksplit[i]= nbt; // save last adresse of the new triangles
-			kkk = nbt;
-		  }
-
-		for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
-
-		if(withBackground)
-		 for (i=0;i<BTh.nbv;i++)
-		  BTh.vertices[i].m =  BTh.vertices[i].m*2.;
-
-		ret = 2;
-		if (nbt>= maxnbt) goto Error; // bug 
-		if (nbv>= maxnbv) goto Error; // bug 
-		// generation of the new triangles 
-
-		SetIntCoor("In SplitElement"); 
-
-		CreateSingleVertexToTriangleConnectivity();
-		if(withBackground)
-		 BTh.CreateSingleVertexToTriangleConnectivity();
-
-		delete [] edges;
-		edges = newedges;
-		nbe = newnbe;
-		nbq = newnbq;
-
-		for (i=0;i<nbsubdomains;i++)
-		  { 
-			long k = subdomains[i].edge- edges;
-			subdomains[i].edge =  edges+2*k; // spilt all edge in 2 
-		  }
-
-		if (ksplitarray) delete [] ksplitarray;
-		if (kedge) delete [] kedge;
-		if (edge4) delete edge4;
-		if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
-		VerticesOnGeomEdge= newVerticesOnGeomEdge;
-		if(VertexOnBThEdge) delete []  VertexOnBThEdge;
-		VertexOnBThEdge = newVertexOnBThEdge;
-		NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
-		NbVertexOnBThEdge=newNbVertexOnBThEdge;
-		//  CreateSingleVertexToTriangleConnectivity();
-
-		ReconstructExistingMesh();
-
-		if (verbose>2){
-			_printf_("   number of quadrilaterals    = " << nbq << "\n");
-			_printf_("   number of triangles         = " << nbt-nbtout- nbq*2 << "\n");
-			_printf_("   number of outside triangles = " << nbtout << "\n");
-		}
-
-		return 0; //ok
-
-Error:
-		nbv = nbvold;
-		nbt = nbtold;
-		nbtout = nbtoutold;
-		// cleaning memory ---
-		delete [] newedges;
-		if (ksplitarray) delete [] ksplitarray;
-		if (kedge) delete [] kedge;
-		if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
-		if (edge4) delete edge4;
-		if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
-
-		return ret; // ok 
-	}
-	/*}}}*/
 /*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{*/
 long  Mesh::SplitInternalEdgeWithBorderVertices(){
@@ -4888,5 +4172,4 @@
 	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
 
-	long int verbose=0;
 	register Triangle *t0,*t;
 	register long k=0, num;   
@@ -4937,5 +4220,5 @@
 	//Vertices
 	if(verbose) _printf_("Reading vertices (" << nbv << ")\n");
-	for (i=0;i<nbv;i++){
+	for(i=0;i<nbv;i++){
 		vertices[i].r.x=x[i];
 		vertices[i].r.y=y[i];
@@ -5415,5 +4698,5 @@
 
 				/*Initialize variables*/
-				double Lstep=0,Lcurve=0;    // step between two points   (phase==1) 
+				double Lstep=0;             // step between two points   (phase==1) 
 				long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 
 
@@ -5548,5 +4831,4 @@
 						long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
 						Lstep = L/NbSegOnCurve; 
-						Lcurve = L;
 						NbCreatePointOnCurve = NbSegOnCurve-1;
 						NbOfNewEdge += NbSegOnCurve;
@@ -5679,11 +4961,11 @@
 				_error_("!v1 || !v2");
 			}
-			Icoor2 detss = 0,l=0,ks;
-			while ((ks=SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
+			Icoor2 detss = 0,l=0;
+			while ((SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
 			 if(l++ > 10000000) {
 				 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
 			 }
 			BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
-			if (( aa == &a ) && (bb == &b) ||  (bb ==  &a ) && (aa == &b)) {
+			if (((aa == &a ) && (bb == &b)) ||((bb ==  &a ) && (aa == &b))){
 				tc.SetLock();
 				a.Optim(1,0);
@@ -5827,10 +5109,10 @@
 			detsb=-detsb;
 
-			if (ToSwap) 
+			if(ToSwap){
 			 if (dets2 < 0) {// haut
 				 dets1 = (ToSwap ? dets1 : detsa) ;
 				 detsa = dets2; 
 				 tt1 =  Previous(tt2) ;}
-			 else if (dets2 > 0){// bas 
+			 else if(dets2 > 0){// bas 
 				 dets1 = (ToSwap ? dets1 : detsb) ;
 				 detsb =  dets2;
@@ -5840,4 +5122,5 @@
 				 tt1 = Next(tt2);
 				 ret =0;}
+			}
 
 		}
