Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2968)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2969)
@@ -314,5 +314,5 @@
 		inline void Set(const Triangle &,const Triangles &,Triangles &);
 		inline int In(Vertex *v) const { return ns[0]==v || ns[1]==v || ns[2]==v ;}
-		TriangleAdjacent FindBoundaryEdge(int ) const;
+		TriangleAdjacent FindBoundaryEdge(int i) const;
 
 		void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 2968)
+++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp	(revision 2969)
@@ -20,26 +20,46 @@
 	/*FUNCTION Triangle::FindBoundaryEdge{{{1*/
 	TriangleAdjacent Triangle::FindBoundaryEdge(int i) const{
-		// turn around  the vertex ns[i] also call  s
+
+		/*Intermediary*/
+		Triangle* ttc=NULL;
+		int k,j,jc;
+
+		// call current triangle t
 		Triangle* t = (Triangle*)this;
-		Triangle* ttc=NULL;
 
-		int k=0,j = EdgesVertexTriangle[i][0],jc;
-		int exterieur  = !link  ;
+		//is the current triangle inside or outside?
+		int outside=!link  ;
 
+		// EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}};
+		// initialize j as the first vertex of the ith edge
+		j=EdgesVertexTriangle[i][0];
+
+		//Loop over the adjacent triangle of t
+		k=0;
 		do{
-			int exterieurp = exterieur;
-			k++; 
+			//keep track of outside
+			int outsidep = outside;
+			//increment k
+			k++;
+			//Get ttc, adjacent triangle of t with respect to vertex j
 			ttc =  t->at[j];
-			exterieur = !ttc->link;
-			if (exterieur+exterieurp == 1) 
-			 return TriangleAdjacent(t,j);
+			//is the current triangle inside or outside?
+			outside = !ttc->link;
+			//if both previous triangle are outside, return
+			if (outside+outsidep == 1) return TriangleAdjacent(t,j);
+
+			//update t and j
+			t = ttc;
+			//NextEdge[3] = {1,2,0};
 			jc = NextEdge[t->aa[j]&3];
-			t = ttc;
 			j = NextEdge[jc];
+
+			//check number of iterations
 			if (k>=2000){
 				throw ErrorException(__FUNCT__,exprintf("too many iteration in Triangle::FindBoundaryEdge (k>=2000)"));
 			}
-		} while (this!= t); 
-		return TriangleAdjacent(0,0);
+		} while (this!= t);
+		//not found, return empty triangle
+		return TriangleAdjacent(NULL,0);
 	}
 	/*}}}1*/
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2968)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2969)
@@ -513,4 +513,7 @@
 				bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing
 				bamgmesh->Edges[i*3+2]=edges[i].ref;
+				/*if(edges[i].onGeometry){
+					printf("ref = %i [%i %i]\n",edges[i].ref,Number(edges[i][0])+1,Number(edges[i][1])+1);
+				}*/
 			}
 		}
@@ -551,4 +554,9 @@
 					bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
 					num=num+1;
+					/*if (reft[i]==-1){
+					  if (!t(0)) printf("%i %i    %i\n",Number(t[1])+1,Number(t[2]),i);
+					  if (!t(1)) printf("%i %i    %i\n",Number(t[2])+1,Number(t[0]),i);
+					  if (!t(2)) printf("%i %i    %i\n",Number(t[0])+1,Number(t[1]),i);
+					  }*/
 				}
 			}
@@ -3489,8 +3497,8 @@
 
 		nbt = 2;
-		triangles[0](0) = 0; //infinite vertex
+		triangles[0](0) = NULL; //infinite vertex
 		triangles[0](1) = v0;
 		triangles[0](2) = v1;
-		triangles[1](0) = 0;//infinite vertex
+		triangles[1](0) = NULL;//infinite vertex
 		triangles[1](2) = v0;
 		triangles[1](1) = v1;
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 2968)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 2969)
@@ -39,5 +39,5 @@
 				nods=domain(i).nods-1; %the domain are closed 1=end;
 				bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
-				bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  ones(nods,1)]];
+				bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  1*ones(nods,1)]];
 				if i>1,
 					%if closed : hole
@@ -51,5 +51,5 @@
 				nods=domain(i).nods-1;
 				bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(:) domain(i).y(:) ones(nods+1,1)]];
-				bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods+1])  ones(nods,1)]];
+				bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods+1])  2*ones(nods,1)]];
 
 			else
@@ -65,6 +65,5 @@
 		end
 		bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
-		bamg_geometry.NumEdges=size(bamg_geometry.Edges,1)
-		bamg_geometry.Edges
+		bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
 		bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1);
 	end