Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2954)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2955)
@@ -16,10 +16,4 @@
 #include <limits.h>
 #include <time.h>
-
-#if  (defined(unix) || defined(__unix)) && !defined(__AIX)
-#define SYSTIMES
-#include <sys/times.h>
-#include <unistd.h>
-#endif
 
 #include "meshtype.h"
@@ -124,6 +118,6 @@
 				Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed 
 				VertexOnGeom * onGeometry;     // if vint 8; // set with Triangles::SetVertexFieldOn()
-				Vertex * onbv; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
-				VertexOnEdge * onbe;   // if vint == 32 on Background edge
+				Vertex * onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
+				VertexOnEdge * onBackgroundEdge;   // if vint == 32 on Background edge
 			};
 			Int1 vint;  // the vertex number in triangle; varies between 0 and 2 in t
@@ -245,16 +239,13 @@
 	class GeometricalEdge {
 		public:
-			GeometricalVertex * v[2];
+			GeometricalVertex* v[2];
 			Int4 ref;
-			Int4  CurveNumber;
-			R2 tg[2]; // the 2 tangente 
-			//   if tg[0] =0 => no continuite 
-			GeometricalEdge * Adj [2]; 
+			Int4 CurveNumber;
+			R2   tg[2]; // the 2 tangentes
+			//   if tg[0] =0 => no continuity
+			GeometricalEdge* Adj[2]; 
 			int DirAdj[2];
-			//  private:
 			int flag ;
-		public: 
 			GeometricalEdge* link; // if   Cracked() or Equi()
-
 			// end of data 
 
@@ -578,5 +569,5 @@
 			VertexOnVertex() {};
 			inline void Set(const Triangles &,Int4,Triangles &);
-			void SetOnBTh(){v->onbv=bv;v->vint=IsVertexOnVertex;}
+			void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}
 	};
 	/*}}}1*/
@@ -590,5 +581,5 @@
 		VertexOnEdge(){}
 		inline void Set(const Triangles &,Int4,Triangles &);  
-		void SetOnBTh(){v->onbe=this;v->vint=IsVertexOnEdge;}  
+		void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;}  
 		Vertex & operator[](int i) const { return (*be)[i];}
 		operator Real8 () const { return abcisse;}
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2954)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 2955)
@@ -478,4 +478,6 @@
 
 		verbose=bamgopts->verbose;
+
+		//Build reft that holds the number the subdomain number of each triangle
 		Int4 *reft = new Int4[nbt];
 		Int4 nbInT = ConsRefTriangle(reft);
@@ -538,4 +540,5 @@
 			for (i=0;i<nbt;i++){
 				Triangle &t=triangles[i];
+				//reft[i]=-1 for outside triangle
 				if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
 					bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing
@@ -1770,19 +1773,36 @@
 		register Triangle *t0,*t;
 		register Int4 k=0, num;   
-		for (Int4 it=0;it<nbt;it++) reft[it]=-1; // outside triangle 
+
+		//initialize all triangles as -1 (outside)
+		for (Int4 it=0;it<nbt;it++) reft[it]=-1;
+
+		//loop over all subdomains
 		for (Int4 i=0;i<NbSubDomains;i++){ 
+
+			//first triangle of the subdomain i
 			t=t0=subdomains[i].head;
+
+			//check that the subdomain is not empty
 			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++;
+				throw ErrorException(__FUNCT__,exprintf("At least one subdomain is empty"));
+			}
+
+			//loop
+			do{
+				k++;
+
+				//get current triangle number
 				num = Number(t);
+
+				//check that num is in [0 nbt[
 				if (num<0 || num>=nbt){
 					throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt"));
 				}
+
+				//reft of this triangle is the subdomain number
 				reft[num]=i;
-			}
-			while (t0 != (t=t->link));
+			} while (t0 != (t=t->link));
+			//stop when all triangles of subdomains have been tagged
+
 			}
 			return k;   
@@ -3897,9 +3917,9 @@
 		Vertex * pvA=&vA, * pvB=&vB;
 		if (vA.vint == IsVertexOnVertex){
-			pA=vA.onbv;
+			pA=vA.onBackgroundVertex;
 		}
 		else if (vA.vint == IsVertexOnEdge){
-			pA=vA.onbe->be;
-			tA=vA.onbe->abcisse;
+			pA=vA.onBackgroundEdge->be;
+			tA=vA.onBackgroundEdge->abcisse;
 		}
 		else {
@@ -3908,9 +3928,9 @@
 
 		if (vB.vint == IsVertexOnVertex){
-			pB=vB.onbv;
+			pB=vB.onBackgroundVertex;
 		}
 		else if(vB.vint == IsVertexOnEdge){
-			pB=vB.onbe->be;
-			tB=vB.onbe->abcisse;
+			pB=vB.onBackgroundEdge->be;
+			tB=vB.onBackgroundEdge->abcisse;
 		}
 		else {
@@ -3934,5 +3954,5 @@
 		if( vA.vint == IsVertexOnEdge) 
 		  { // find the start edge 
-			e = vA.onbe->be;	 
+			e = vA.onBackgroundEdge->be;	 
 
 		  } 
@@ -3944,5 +3964,5 @@
 			Exchange(pvA,pvB);
 			Exchange(A,B);
-			e =  vB.onbe->be;
+			e =  vB.onBackgroundEdge->be;
 
 		  } 
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 2954)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 2955)
@@ -18,4 +18,6 @@
 bamg_geometry.NumEdges=0;
 bamg_geometry.Edges=zeros(0,3);
+bamg_geometry.NumCrackedEdges=0;
+bamg_geometry.CrackedEdges=zeros(0,1);
 bamg_geometry.hVertices=zeros(0,1);
 bamg_geometry.NumSubDomains=0;
@@ -36,6 +38,11 @@
 			bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  ones(nods,1)]];
 			if i>1,
-				clockwise=-1;
-				bamg_geometry.SubDomains=[2 count+1 clockwise 1];
+				%if closed : hole
+				if (domain(i).x(1)==domain(i).x(end) & domain(i).y(1)==domain(i).y(end)),
+					clockwise=-1;
+					bamg_geometry.SubDomains=[2 count+1 clockwise 1];
+				else
+					%rift
+				end
 			end
 			count=count+nods;
@@ -122,4 +129,5 @@
 
 %Fill in rest of fields:
+md.type='2d';
 md.numberofelements=length(md.elements);
 md.numberofgrids=length(md.x);
@@ -127,4 +135,5 @@
 md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
 md.gridonbed=ones(md.numberofgrids,1);
+md.gridonwater=zeros(md.numberofgrids,1);
 md.gridonsurface=ones(md.numberofgrids,1);
 md.elementonbed=ones(md.numberofelements,1);
@@ -138,13 +147,2 @@
 md.segments=findsegments(md);
 md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
-
-%Fill in rest of fields:
-md.z=zeros(md.numberofgrids,1);
-md.gridonbed=ones(md.numberofgrids,1);
-md.gridonsurface=ones(md.numberofgrids,1);
-md.elementonbed=ones(md.numberofelements,1);
-md.elementonsurface=ones(md.numberofelements,1);
-md.gridonwater=zeros(md.numberofgrids,1);
-
-%type of model
-md.type='2d';
Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 2954)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 2955)
@@ -171,7 +171,2 @@
 
 end
-
-
-
-
-
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2954)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 2955)
@@ -26,4 +26,6 @@
 	int     NumEdgesGeom;
 	double* EdgesGeom=NULL;
+	int     NumCrackedEdgesGeom;
+	double* CrackedEdgesGeom=NULL;
 	double* hVerticesGeom=NULL;
 	double  MaximalAngleOfCorner;
@@ -73,6 +75,8 @@
 	bamggeom.NumRequiredEdges=0;
 	bamggeom.RequiredEdges=NULL;
-	bamggeom.NumCrackedEdges=0;
-	bamggeom.CrackedEdges=NULL;
+	FetchData(&NumCrackedEdgesGeom,mxGetField(BAMGGEOMETRY,0,"NumCrackedEdges"));
+	bamggeom.NumCrackedEdges=NumCrackedEdgesGeom;
+	FetchData(&CrackedEdgesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"CrackedEdges"));
+	bamggeom.CrackedEdges=CrackedEdgesGeom;
 	FetchData(&NumSubDomainsGeom,mxGetField(BAMGGEOMETRY,0,"NumSubDomains"));
 	bamggeom.NumSubDomains=NumSubDomainsGeom;