Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.cpp	(revision 2781)
+++ /issm/trunk/src/c/Bamgx/Mesh2.cpp	(revision 2782)
@@ -4128,347 +4128,339 @@
 		}
 
-		void Triangles::FillHoleInMesh(){
-
+		void Triangles::FillHoleInMesh() {
 			Triangles * OldCurrentTh =CurrentTh;
-			long int verbosity=200;
-			
 			CurrentTh=this;
 			//  Int4 NbTold = nbt;
 			// generation of the integer coor
-
-			//  double coef = coefIcoor;
-			// recherche des extrema des vertices pmin,pmax
-			Int4 i;
-			if(verbosity>2)
-				cout << "  -- FillHoleInMesh: Nb of vertices =" << nbv 
+			  {
+
+				//  double coef = coefIcoor;
+				// recherche des extrema des vertices pmin,pmax
+				Int4 i;
+				if(verbosity>2)
+				 cout << "  -- FillHoleInMesh: Nb of vertices =" << nbv 
 					<< " Pmin = "<< pmin << " Pmax = "<< pmax << endl;
 
-			assert(ordre);
-			for (i=0;i<nbv;i++) 
-				ordre[i]= 0 ;
-
-
-			NbSubDomains =0;
-
-			// generation of the adjacence of the triangles
-			SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
-			Int4 * st = new Int4[nbt*3];
-			for (i=0;i<nbt*3;i++)
-				st[i]=-1;
-			Int4 kk =0;
-			for (i=0;i<nbe;i++)
-				kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
-			if (kk != nbe)
-			{ 
-				cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl;
-				MeshError(1002,this);
-			}
-			for (i=0;i<nbt;i++)
-				for (int j=0;j<3;j++)
-				{
-					// Int4 i0,i1;
-					Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
-							Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
-					Int4 invisible = triangles[i].Hidden(j);
-					if(st[k]==-1)
-						st[k]=3*i+j;
-					else if(st[k]>=0) {
-						assert( ! 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));
-						if (invisible)  triangles[i].SetHidden(j);
-						if (k<nbe) {
-							triangles[i].SetLocked(j);
-						}
-						st[k]=-2-st[k]; }
-					else {
-						cerr << " The edge (" 
+				assert(ordre);
+				for (i=0;i<nbv;i++) 
+				 ordre[i]= 0 ;
+
+
+				NbSubDomains =0;
+
+				// generation of the adjacence of the triangles
+				SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
+				Int4 * st = new Int4[nbt*3];
+				for (i=0;i<nbt*3;i++)
+				 st[i]=-1;
+				Int4 kk =0;
+				for (i=0;i<nbe;i++)
+				 kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
+				if (kk != nbe)
+				  { 
+					cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl;
+					MeshError(1002,this);
+				  }
+				for (i=0;i<nbt;i++)
+				 for (int j=0;j<3;j++)
+					{
+					 // Int4 i0,i1;
+					 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
+								 Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
+					 Int4 invisible = triangles[i].Hidden(j);
+					 if(st[k]==-1)
+					  st[k]=3*i+j;
+					 else if(st[k]>=0) {
+						 assert( ! 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));
+						 if (invisible)  triangles[i].SetHidden(j);
+						 if (k<nbe) {
+							 triangles[i].SetLocked(j);
+						 }
+						 st[k]=-2-st[k]; }
+					 else {
+						 cerr << " The edge (" 
 							<< Number(triangles[i][VerticesOfTriangularEdge[j][0]])
 							<< " , " 
 							<< Number(triangles[i][VerticesOfTriangularEdge[j][1]])
 							<< " ) is in more than 2 triangles " <<k <<endl;
-						cerr << " Edge " << j << " Of Triangle " << i << endl;
-						cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3  << endl;
-						cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3))
+						 cerr << " Edge " << j << " Of Triangle " << i << endl;
+						 cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3  << endl;
+						 cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3))
 							<< " Of Triangle " <<  Number(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) << endl;
-						MeshError(9999,this);}	
-
-
-				}
-			
-			if(verbosity>5) {
-				cout << "    On Mesh " << name << endl;
-				cout << "    - The number of Vertices  = " << nbv << endl;
-				cout << "    - The number of Triangles = " << nbt << endl;
-				cout << "    - The number of given edge = " << nbe << endl;
-				cout << "    - The number of all edges = " << edge4->nb() << endl;
-				cout << "    - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl; 
-			}
-
-
-			// check the consistant of edge[].adj and the geometrical required  vertex
-			Int4 k=0;
-			for (i=0;i<edge4->nb();i++)
-				if (st[i] >=0) // edge alone 
-					if (i < nbe) 
-					{
-						Int4 i0=edge4->i(i);ordre[i0] = vertices+i0;
-						Int4 i1=edge4->j(i);ordre[i1] = vertices+i1;
+						 MeshError(9999,this);}	
+
+
 					}
-					else {
-						k++;
-						if (verbosity>20 && k <20) 
-						{
-							Int4 i0=edge4->i(i);
-							Int4 i1=edge4->j(i);
-							cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl;
-						}
+				if(verbosity>5) {
+					cout << "    On Mesh " << name << endl;
+					cout << "    - The number of Vertices  = " << nbv << endl;
+					cout << "    - The number of Triangles = " << nbt << endl;
+					cout << "    - The number of given edge = " << nbe << endl;
+					cout << "    - The number of all edges = " << edge4->nb() << endl;
+					cout << "    - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl; }
+
+
+					// check the consistant of edge[].adj and the geometrical required  vertex
+					Int4 k=0;
+					for (i=0;i<edge4->nb();i++)
+					 if (st[i] >=0) // edge alone 
+					  if (i < nbe) 
+						 {
+						  Int4 i0=edge4->i(i);ordre[i0] = vertices+i0;
+						  Int4 i1=edge4->j(i);ordre[i1] = vertices+i1;
+						 }
+					  else {
+						  k++;
+						  if (verbosity>20 && k <20) 
+							 {
+							  Int4 i0=edge4->i(i);
+							  Int4 i1=edge4->j(i);
+							  cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl;
+							 }
+					  }
+
+					if(k != 0) {
+						if (verbosity>20)
+						  {
+							cout << " The given edge are " << endl;
+							for (int i=0;i< nbe;i++)
+							 cout <<  " Edge " << i << " : " <<  Number(edges[i][0]) << " " <<  Number(edges[i][1]) 
+								<< " " << edges[i].ref << endl; 
+						  }
+						cerr << k << " boundary edges  are not defined as edges " << endl;
+						MeshError(9998,this);
 					}
-
-			if(k != 0) {
-				if (verbosity>20)
-				{
-					cout << " The given edge are " << endl;
-					for (int i=0;i< nbe;i++)
-						cout <<  " Edge " << i << " : " <<  Number(edges[i][0]) << " " <<  Number(edges[i][1]) 
-							<< " " << edges[i].ref << endl; 
-				}
-				cerr << k << " boundary edges  are not defined as edges " << endl;
-				MeshError(9998,this);
-			}
-			printf("salope1\n");
-
-			// generation of the mesh with boundary points   
-			Int4 nbvb = 0;
-			for (i=0;i<nbv;i++)
-			{ 
-				vertices[i].t=0;
-				vertices[i].vint=0;
-				if (ordre[i]) 
-					ordre[nbvb++] = ordre[i];
-			}
-
-			Triangle *savetriangles= triangles;
-			Int4 savenbt=nbt;
-			Int4 savenbtx=nbtx;
-			SubDomain * savesubdomains = subdomains;
-			subdomains = 0;
-
-			Int4  Nbtriafillhole = 2*nbvb;
-			Triangle * triafillhole =new Triangle[Nbtriafillhole];
-			if (verbosity>9)
-				cout << " Nbtriafillhole triafillhole*" << triafillhole << endl; 
-			triangles =  triafillhole;
-
-			nbt=2;
-			nbtx= Nbtriafillhole;
-
-			for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
-				if  ( ++i >= nbvb) {
-					cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;
-					MeshError(998,this); 
-				}
-			}
-			Exchange( ordre[2], ordre[i]);
-
-			Vertex *  v0=ordre[0], *v1=ordre[1];
-
-
-			triangles[0](0) = 0; // sommet pour infini 
-			triangles[0](1) = v0;
-			triangles[0](2) = v1;
-
-			triangles[1](0) = 0;// sommet pour infini 
-			triangles[1](2) = v0;
-			triangles[1](1) = v1;
-			const int e0 = OppositeEdge[0];
-			const int e1 = NextEdge[e0];
-			const int e2 = PreviousEdge[e0];
-			triangles[0].SetAdj2(e0, &triangles[1] ,e0);
-			triangles[0].SetAdj2(e1, &triangles[1] ,e2);
-			triangles[0].SetAdj2(e2, &triangles[1] ,e1);
-
-			triangles[0].det = -1;  // faux triangles
-			triangles[1].det = -1;  // faux triangles
-
-			triangles[0].SetTriangleContainingTheVertex();
-			triangles[1].SetTriangleContainingTheVertex();
-
-			triangles[0].link=&triangles[1];
-			triangles[1].link=&triangles[0];
-
-			printf("salope2\n");
-			//  nbtf = 2;
-			if (  !quadtree ) delete  quadtree; // ->ReInitialise();
-
-			quadtree = new QuadTree(this,0);
-			quadtree->Add(*v0);
-			quadtree->Add(*v1);
-
-			// on ajoute les sommets un a un 
-			Int4 NbSwap=0;
-			for (Int4 icount=2; icount<nbvb; icount++) {
-
-				Vertex *vi  = ordre[icount];
-				//	  cout << " Add vertex " <<  Number(vi) << endl;
-				Icoor2 dete[3];
-				Triangle *tcvi = FindTriangleContening(vi->i,dete);
-				quadtree->Add(*vi); 
-				Add(*vi,tcvi,dete);
-				NbSwap += vi->Optim(1,1);
-
-			}// end loop on  icount	
-			printf("salope3\n");
-
-			//Int4 nbtfillhole = nbt;
-			// inforce the boundary 
-			TriangleAdjacent ta(0,0);
-			Int4 nbloss = 0,knbe=0;
-			for ( i = 0; i < nbe; i++) 
-				if (st[i] >=0)  // edge alone => on border ...  FH oct 2009
-				{
-					Vertex & a=edges[i][0], & b =    edges[i][1];
-					if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
-					{
-						knbe++;
-						if (ForceEdge(a,b,ta)<0)
-							nbloss++;
-					}
-				}
-			if(nbloss)
-			{
-				cerr << " we loss some  " << nbloss << " "  << " edges other " << knbe << endl;
-				MeshError(1100,this);
-			}
-			FindSubDomain(1);
-			// remove all the hole 
-			// remove all the good sub domain
-			Int4 krm =0;
-			for (i=0;i<nbt;i++)
-				if (triangles[i].link) // remove triangles
-				{
-					krm++;
-					for (int j=0;j<3;j++)
-					{
-						TriangleAdjacent ta =  triangles[i].Adj(j);
-						Triangle & tta = * (Triangle *) ta;
-						if(! tta.link) // edge between remove and not remove 
-						{ // change the link of ta;
-							int ja = ta;
-							Vertex *v0= ta.EdgeVertex(0);
-							Vertex *v1= ta.EdgeVertex(1);
-							Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
-							assert(st[k] >=0); 
-							tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
-							ta.SetLock();
-							st[k]=-2-st[k]; 
-						}
-					}
-				}
-			Int4 NbTfillHoll =0;
-			for (i=0;i<nbt;i++)
-				if (triangles[i].link) {
-					triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
-					triangles[i].color=-1;
-				}
-				else
-				{
-					triangles[i].color= savenbt+ NbTfillHoll++;
-				}
-			// cout <<      savenbt+NbTfillHoll << " " <<  savenbtx  << endl;
-			assert(savenbt+NbTfillHoll <= savenbtx );
-			// copy of the outside triangles in saveTriangles 
-			for (i=0;i<nbt;i++)
-				if(triangles[i].color>=0) 
-				{
-					savetriangles[savenbt]=triangles[i];
-					savetriangles[savenbt].link=0;
-					savenbt++;
-				}
-			// gestion of the adj
-			k =0;
-			Triangle * tmax = triangles + nbt;
-			for (i=0;i<savenbt;i++)  
-			{ 
-				Triangle & ti = savetriangles[i];
-				for (int j=0;j<3;j++)
-				{
-					Triangle * ta = ti.TriangleAdj(j);
-					int aa = ti.NuEdgeTriangleAdj(j);
-					int lck = ti.Locked(j);
-					if (!ta) k++; // bug 
-					else if ( ta >= triangles && ta < tmax) 
-					{
-						ta= savetriangles + ta->color;
-						ti.SetAdj2(j,ta,aa);
-						if(lck) ti.SetLocked(j);
-					}
-				}
-			}
-			//	 OutSidesTriangles = triangles;
-			//	Int4 NbOutSidesTriangles = nbt;
-
-			// restore triangles;
-			nbt=savenbt;
-			nbtx=savenbtx;
-			delete [] triangles;
-			delete [] subdomains;
-			triangles = savetriangles;
-			subdomains = savesubdomains;
-			//	 cout <<  triangles << " <> " << OutSidesTriangles << endl; 
-			/*	 k=0;
-				 for (i=0;i<nbt;i++)
-				 for (int j=0;j<3;j++)
-				 if (!triangles[i].TriangleAdj(j))
-				 k++;
-
-*/
-
-			printf("salope4\n");
-			if (k) {
-				cerr << "Error Nb of triangles edge alone = " << k << endl;
-				MeshError(9997,this);
-			}
-			printf("salope4a\n");
-			FindSubDomain();
-			// cout << " NbTOld = " << NbTold << " ==  " << nbt - NbOutT << " " << nbt << endl;
-
-			// 
-
-			printf("salope4b\n");
-			delete edge4;
-			printf("salope4ba\n");
-			delete [] st;
-			printf("salope4bb %i\n",nbv);
-			for (i=0;i<nbv;i++){
-				printf("%i\n",i);
-				quadtree->Add(vertices[i]);
-			}
-			printf("salope4c\n");
-
-			SetVertexFieldOn();
-			printf("salope5a\n");
-
-			for (i=0;i<nbe;i++){
-				if(edges[i].on){
-					for(int j=0;j<2;j++){
-						if (!edges[i].adj[j]){
-							if(!edges[i][j].on->IsRequiredVertex()) {
-								cerr << " Erreur adj et sommet requis edges [" << i <<  "][ " << j << "]= "
-									<<  Number(edges[i][j]) << " : "  << " on = " << Gh.Number(edges[i].on) ;
-								if (edges[i][j].on->OnGeomVertex())
+					// generation of the mesh with boundary points   
+					Int4 nbvb = 0;
+					for (i=0;i<nbv;i++)
+					  { 
+						vertices[i].t=0;
+						vertices[i].vint=0;
+						if (ordre[i]) 
+						 ordre[nbvb++] = ordre[i];
+					  }
+
+					Triangle *savetriangles= triangles;
+					Int4 savenbt=nbt;
+					Int4 savenbtx=nbtx;
+					SubDomain * savesubdomains = subdomains;
+					subdomains = 0;
+
+					Int4  Nbtriafillhole = 2*nbvb;
+					Triangle * triafillhole =new Triangle[Nbtriafillhole];
+					if (verbosity>9)
+					 cout << " Nbtriafillhole triafillhole*" << triafillhole << endl; 
+					triangles =  triafillhole;
+
+					nbt=2;
+					nbtx= Nbtriafillhole;
+
+					for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 
+					 if  ( ++i >= nbvb) {
+						 cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;
+						 MeshError(998,this); }
+						 Exchange( ordre[2], ordre[i]);
+
+						 Vertex *  v0=ordre[0], *v1=ordre[1];
+
+
+						 triangles[0](0) = 0; // sommet pour infini 
+						 triangles[0](1) = v0;
+						 triangles[0](2) = v1;
+
+						 triangles[1](0) = 0;// sommet pour infini 
+						 triangles[1](2) = v0;
+						 triangles[1](1) = v1;
+						 const int e0 = OppositeEdge[0];
+						 const int e1 = NextEdge[e0];
+						 const int e2 = PreviousEdge[e0];
+						 triangles[0].SetAdj2(e0, &triangles[1] ,e0);
+						 triangles[0].SetAdj2(e1, &triangles[1] ,e2);
+						 triangles[0].SetAdj2(e2, &triangles[1] ,e1);
+
+						 triangles[0].det = -1;  // faux triangles
+						 triangles[1].det = -1;  // faux triangles
+
+						 triangles[0].SetTriangleContainingTheVertex();
+						 triangles[1].SetTriangleContainingTheVertex();
+
+						 triangles[0].link=&triangles[1];
+						 triangles[1].link=&triangles[0];
+
+#ifdef DEBUG 
+						 triangles[0].check();
+						 triangles[1].check();
+#endif  
+						 //  nbtf = 2;
+						 if (  !quadtree ) 
+						  delete  quadtree; // ->ReInitialise();
+
+						 quadtree = new QuadTree(this,0);
+						 quadtree->Add(*v0);
+						 quadtree->Add(*v1);
+
+						 // on ajoute les sommets un a un 
+						 Int4 NbSwap=0;
+						 for (Int4 icount=2; icount<nbvb; icount++) {
+
+							 Vertex *vi  = ordre[icount];
+							 //	  cout << " Add vertex " <<  Number(vi) << endl;
+							 Icoor2 dete[3];
+							 Triangle *tcvi = FindTriangleContening(vi->i,dete);
+							 quadtree->Add(*vi); 
+							 Add(*vi,tcvi,dete);
+							 NbSwap += vi->Optim(1,1);
+
+#ifdef DRAWING2
+							 cout << Number(vi) << " " <<  NbSwap <<  endl;
+							 reffecran();
+							 Draw();
+							 vi->Draw();
+							 inquire();
+#endif
+						 }// end loop on  icount	
+#ifdef DRAWING1
+						 inquire();
+#endif
+
+						 //Int4 nbtfillhole = nbt;
+						 // inforce the boundary 
+						 TriangleAdjacent ta(0,0);
+						 Int4 nbloss = 0,knbe=0;
+						 for ( i = 0; i < nbe; i++) 
+						  if (st[i] >=0)  // edge alone => on border ...  FH oct 2009
+							 {
+							  Vertex & a=edges[i][0], & b =    edges[i][1];
+							  if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
+								 {
+								  knbe++;
+								  if (ForceEdge(a,b,ta)<0)
+									nbloss++;
+								 }
+							 }
+						 if(nbloss)
+							{
+							 cerr << " we loss some  " << nbloss << " "  << " edges other " << knbe << endl;
+							 MeshError(1100,this);
+							}
+						 FindSubDomain(1);
+						 // remove all the hole 
+						 // remove all the good sub domain
+						 Int4 krm =0;
+						 for (i=0;i<nbt;i++)
+						  if (triangles[i].link) // remove triangles
+							 {
+							  krm++;
+							  for (int j=0;j<3;j++)
+								 {
+								  TriangleAdjacent ta =  triangles[i].Adj(j);
+								  Triangle & tta = * (Triangle *) ta;
+								  if(! tta.link) // edge between remove and not remove 
+									 { // change the link of ta;
+									  int ja = ta;
+									  Vertex *v0= ta.EdgeVertex(0);
+									  Vertex *v1= ta.EdgeVertex(1);
+									  Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
+									  assert(st[k] >=0); 
+									  tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
+									  ta.SetLock();
+									  st[k]=-2-st[k]; 
+									 }
+								 }
+							 }
+						 Int4 NbTfillHoll =0;
+						 for (i=0;i<nbt;i++)
+						  if (triangles[i].link) {
+							  triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
+							  triangles[i].color=-1;
+						  }
+						  else
+							 {
+							  triangles[i].color= savenbt+ NbTfillHoll++;
+#ifdef DEBUG 
+							  triangles[i].check();
+#endif
+							 }
+						 // cout <<      savenbt+NbTfillHoll << " " <<  savenbtx  << endl;
+						 assert(savenbt+NbTfillHoll <= savenbtx );
+						 // copy of the outside triangles in saveTriangles 
+						 for (i=0;i<nbt;i++)
+						  if(triangles[i].color>=0) 
+							 {
+							  savetriangles[savenbt]=triangles[i];
+							  savetriangles[savenbt].link=0;
+							  savenbt++;
+							 }
+						 // gestion of the adj
+						 k =0;
+						 Triangle * tmax = triangles + nbt;
+						 for (i=0;i<savenbt;i++)  
+							{ 
+							 Triangle & ti = savetriangles[i];
+							 for (int j=0;j<3;j++)
+								{
+								 Triangle * ta = ti.TriangleAdj(j);
+								 int aa = ti.NuEdgeTriangleAdj(j);
+								 int lck = ti.Locked(j);
+								 if (!ta) k++; // bug 
+								 else if ( ta >= triangles && ta < tmax) 
+									{
+									 ta= savetriangles + ta->color;
+									 ti.SetAdj2(j,ta,aa);
+									 if(lck) ti.SetLocked(j);
+									}
+								}
+							}
+						 //	 OutSidesTriangles = triangles;
+						 //	Int4 NbOutSidesTriangles = nbt;
+
+						 // restore triangles;
+						 nbt=savenbt;
+						 nbtx=savenbtx;
+						 delete [] triangles;
+						 delete [] subdomains;
+						 triangles = savetriangles;
+						 subdomains = savesubdomains;
+						 //	 cout <<  triangles << " <> " << OutSidesTriangles << endl; 
+						 /*	 k=0;
+								 for (i=0;i<nbt;i++)
+								 for (int j=0;j<3;j++)
+								 if (!triangles[i].TriangleAdj(j))
+								 k++;
+								 */
+						 if (k) {
+							 cerr << "Error Nb of triangles edge alone = " << k << endl;
+							 MeshError(9997,this);
+						 }
+						 FindSubDomain();
+						 // cout << " NbTOld = " << NbTold << " ==  " << nbt - NbOutT << " " << nbt << endl;
+
+						 // 
+
+						 delete edge4;
+						 delete [] st;
+						 for (i=0;i<nbv;i++)
+						  quadtree->Add(vertices[i]);
+
+						 SetVertexFieldOn();
+
+						 for (i=0;i<nbe;i++)
+						  if(edges[i].on) 
+							for(int j=0;j<2;j++)
+							 if (!edges[i].adj[j])
+							  if(!edges[i][j].on->IsRequiredVertex()) {
+								  cerr << " Erreur adj et sommet requis edges [" << i <<  "][ " << j << "]= "
+									 <<  Number(edges[i][j]) << " : "  << " on = " << Gh.Number(edges[i].on) ;
+								  if (edges[i][j].on->OnGeomVertex())
 									cerr << " vertex " << Gh.Number(edges[i][j].on->gv);
-								else if (edges[i][j].on->OnGeomEdge())
+								  else if (edges[i][j].on->OnGeomEdge())
 									cerr << "Edges " << Gh.Number(edges[i][j].on->ge);
-								else
+								  else
 									cerr << " = " << edges[i][j].on ;
-								cerr << endl;
-							}
-						}
-					}
-				}
-			}
-
-			printf("salope5\n");
+								  cerr << endl;
+							  }
+			  }
 			CurrentTh=OldCurrentTh;
 		}