Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2786)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2787)
@@ -173,8 +173,9 @@
 		Triangles BTh(bamgmesh,bamgargs); // read the background mesh 
 
+		throw ErrorException(__FUNCT__,exprintf("TESTTESTTESTTEST"));
+		printf("ok1\n"); 
 		hmin = Max(hmin,BTh.MinimalHmin());
+		printf("ok2\n"); 
 		hmax = Min(hmax,BTh.MaximalHmax());
-
-		throw ErrorException(__FUNCT__," done");
 
 		BTh.MakeQuadTree();
Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.cpp	(revision 2786)
+++ /issm/trunk/src/c/Bamgx/Mesh2.cpp	(revision 2787)
@@ -26,13 +26,7 @@
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
    */
-#ifdef __MWERKS__
-#ifdef __INTEL__
-//#pragma global_optimizer off
-//#pragma inline_depth(0)
-//#pragma optimization_level 2
-#endif
-//#pragma inline_depth 0
-#endif
 extern bool withrgraphique;
+
+#include "../shared/shared.h"
 #include <stdio.h>
 #include <string.h>
@@ -40,4 +34,5 @@
 #include <time.h>
 #include <iostream>
+
 using namespace std; 
 
@@ -4129,270 +4124,235 @@
 
 		void Triangles::FillHoleInMesh() {
-			Triangles * OldCurrentTh =CurrentTh;
+
+			Triangles* OldCurrentTh =CurrentTh;
 			CurrentTh=this;
-			int verbosity=1;
-			//  Int4 NbTold = nbt;
-			// generation of the integer coor
+
+			int verbosity=10;
+
+			// generation of the integer coordinate
 			  {
 
-				//  double coef = coefIcoor;
-				// recherche des extrema des vertices pmin,pmax
+				// find extrema coordinates of vertices pmin,pmax
 				Int4 i;
-				if(verbosity>2)
-				 cout << "  -- FillHoleInMesh: Nb of vertices =" << nbv 
-					<< " Pmin = "<< pmin << " Pmax = "<< pmax << endl;
-
+				if(verbosity>2) printf("      Filling holes in mesh of %i vertices\n",nbv); 
+
+				//initialize ordre
 				assert(ordre);
-				for (i=0;i<nbv;i++) 
-				 ordre[i]= 0 ;
-
+				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;
+				/* generation of the adjacence of the triangles*/
+
+				SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
+
+				//initialize st
+				Int4* st = new Int4[nbt*3];
+				for (i=0;i<nbt*3;i++) st[i]=-1;
+
+				//check number of edges
 				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);
+				for (i=0;i<nbe;i++){
+					kk=kk+(i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
+				}
+				if (kk != nbe) { 
+					throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe));
+				}
+
+				//
+				for (i=0;i<nbt;i++){
+					for (int j=0;j<3;j++) {
+						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 {
+							throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",
+											Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));
+						}
+					}
+				}
+				if(verbosity>5) {
+					printf("         info on Mesh %s:\n",name);
+					printf("            - number of vertices    = %i \n",nbv); 
+					printf("            - number of triangles   = %i \n",nbt); 
+					printf("            - number of given edges = %i \n",nbe); 
+					printf("            - number of all edges   = %i \n"  ,edge4->nb()); 
+					printf("            - Euler number 1 - nb of holes = %i \n"  ,nbt-edge4->nb()+nbv); 
+				}
+
+				// 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=k+1;
+							if (k <20) {
+								throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));
+							}
+						}
+					}
+				}
+				if(k != 0) {
+					throw ErrorException(__FUNCT__,exprintf("%i boundary edges  are not defined as edges",k));
+				}
+
+				/* mesh generation 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];
+				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];
+
+					 if (!quadtree ) 
+					  delete  quadtree; // ->ReInitialise();
+
+					 quadtree = new QuadTree(this,0);
+					 quadtree->Add(*v0);
+					 quadtree->Add(*v1);
+
+					 // We add the vertices one by one
+					 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);
+					 }
+
+					 // 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++;
+								}
 						 }
-						 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))
-							<< " 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;
+					 }
+					 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]; 
+									}
+								}
 						 }
-					  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;
+					 }
+						 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;
 							 }
-					  }
-
-					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);
-					}
-					// 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
-							 }
+							 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++;
+						 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;
@@ -4439,9 +4399,6 @@
 						 // cout << " NbTOld = " << NbTold << " ==  " << nbt - NbOutT << " " << nbt << endl;
 
-						 // 
-
-						 //delete edge4; //TESTESTESTSTSTSTSTST
-						 //delete [] st; //TESTESTESTSTSTSTSTST
-						 printf("ok1\n"); 
+						 delete edge4;
+						 delete [] st;
 						 for (i=0;i<nbv;i++)
 						  quadtree->Add(vertices[i]);
@@ -4465,5 +4422,4 @@
 							  }
 			  }
-			printf("ok2\n"); 
 			CurrentTh=OldCurrentTh;
 		}
Index: /issm/trunk/src/c/Bamgx/MeshRead.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/MeshRead.cpp	(revision 2786)
+++ /issm/trunk/src/c/Bamgx/MeshRead.cpp	(revision 2787)
@@ -70,5 +70,5 @@
 			vertices[i].DirOfSearch =NoDirOfSearch;
 			vertices[i].m=M1;
-			vertices[i].color=bamgmesh->Vertices[i*3+2];
+			vertices[i].color=(Int4)bamgmesh->Vertices[i*3+2];
 		}
 		nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals 
@@ -81,5 +81,6 @@
 	if(bamgmesh->Triangles){
 		if(verbose>3) printf("      processing Triangles\n");
-		triangles =new Triangle[nbt];
+		triangles =new Triangle[nbtx]; //we cannot allocate only nbt triangles since 
+												 //other triangles will be added for each edge
 		for (i=0;i<nbt;i++){
 			Triangle & t = triangles[i];
@@ -92,5 +93,5 @@
 	}
 	else{
-		if(verbose>3) throw ErrorException(__FUNCT__,exprintf("no Vertices found in the initial mesh"));
+		if(verbose>3) throw ErrorException(__FUNCT__,exprintf("no Triangles found in the initial mesh"));
 	}
 
@@ -252,5 +253,5 @@
 			if (i3!=23) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 3"));
 			if (head<0 || head>=nbt) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: head should in [1 %i] (triangle number)",nbt));
-			//subdomains[i].head=i3;   //TESTTESTETETETETETETETEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
+			subdomains[i].head = triangles+head;
 		}
 	}
@@ -263,5 +264,5 @@
 		printf("Warning: mesh present but no geometry found. Reconstructing...\n");
 		/*Recreate geometry: */
-		ConsGeometry(-1);	//MaximalAngleOfCorner is actually in BamgGeom...
+		ConsGeometry(-1);	//TESTTESTTESTTEST MaximalAngleOfCorner is actually in BamgGeom...
 		Gh.AfterRead();
 	}
@@ -1148,5 +1149,5 @@
   SetIntCoor();
   FillHoleInMesh();
-  throw ErrorException(__FUNCT__,exprintf("TEST"));
+  printf("ok0\n");
 }