Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5530)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5531)
@@ -9,5 +9,4 @@
 
 	static const  Direction NoDirOfSearch=Direction();
-	long NbUnSwap =0;
 
 	/*Constructors/Destructors*/
@@ -4895,9 +4894,4 @@
 		int verbose=bamgopts->verbose;
 
-		//initialize Mesh
-		nbv=0;
-		NbVerticesOnGeomVertex=0;
-		NbVerticesOnGeomEdge=0;
-
 		//build background mesh flag (1 if background, else 0)
 		bool background=(&BTh != this);
@@ -4912,6 +4906,6 @@
 		VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];  
 		if(NbVerticesOnGeomVertex >= maxnbv) ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);
+		ISSMASSERT(nbv==0);
 		//Build VerticesOnGeomVertex
-		nbv=0;
 		for (i=0;i<Gh.nbv;i++){
 			/* Add vertex only if required*/
@@ -4941,7 +4935,7 @@
 		 * to the one present in the geometry provided.
 		 * We proceed in 2 steps
-		 *  -step 1: we count all the edges
-		 *           we allocate the number of edges at the end of step 1
-		 *  -step 2: the edges are created */
+		 *  -step 0: we count all the edges
+		 *           we allocate the number of edges at the end of step 0
+		 *  -step 1: the edges are created */
 		for (int step=0;step<2;step++){
 
@@ -4962,5 +4956,5 @@
 				for(int j=0;j<2;j++) {
 
-					/*Take only required vertices (corner)*/
+					/*Take only required vertices (corner->beginning of a new curve)*/
 					if (!ei.Mark() && ei[j].Required()){ 
 
@@ -4970,21 +4964,23 @@
 
 						/*If Edge is required (do that only once for the 2 vertices)*/
-						if(ei.Required() && j==0){
-							//do not create internal points if required (take it as is)
-							if(step==0) nbe++;
-							else{ 
-								e=&ei;
-								a=ei(0);
-								b=ei(1);
-
-								//check that edges has been allocated
-								ISSMASSERT(edges);
-								edges[nbe].v[0]=a->MeshVertexHook;
-								edges[nbe].v[1]=b->MeshVertexHook;;
-								edges[nbe].ReferenceNumber = e->ReferenceNumber;
-								edges[nbe].GeometricalEdgeHook = e;
-								edges[nbe].adj[0] = 0;
-								edges[nbe].adj[1] = 0;
-								nbe++;
+						if(ei.Required()){
+							if (j==0){
+								//do not create internal points if required (take it as is)
+								if(step==0) nbe++;
+								else{ 
+									e=&ei;
+									a=ei(0);
+									b=ei(1);
+
+									//check that edges has been allocated
+									ISSMASSERT(edges);
+									edges[nbe].v[0]=a->MeshVertexHook;
+									edges[nbe].v[1]=b->MeshVertexHook;;
+									edges[nbe].ReferenceNumber = e->ReferenceNumber;
+									edges[nbe].GeometricalEdgeHook = e;
+									edges[nbe].adj[0] = 0;
+									edges[nbe].adj[1] = 0;
+									nbe++;
+								}
 							}
 						}
@@ -4993,7 +4989,6 @@
 						else {
 							for (int kstep=0;kstep<=step;kstep++){
-								//kstep=0, do nothing
-								//kstep=1, compute the length of the curve
-								//kstep=2  create the points and edge
+								//kstep=0, compute number of edges (discretize curve)
+								//kstep=1  create the points and edge
 								PreviousNewEdge=0;
 								NbNewPoints=0;
@@ -5001,5 +4996,5 @@
 								if (nbvend>=maxnbv) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase maxnbv");
 								lcurve =0;
-								s = lstep;
+								s = lstep; //-1 initially, then length of each sub edge
 
 								/*reminder: i = edge number, j=[0;1] vertex index in edge*/
@@ -5010,5 +5005,5 @@
 								e->SetMark();   // Mark edge
 
-								/*If we have a Background mesh, we can use it to discretize the curve*/
+								/*Loop until we reach the end of the curve*/
 								for(;;){ 
 									k = 1-k;            // other vertx index of the curve
@@ -5016,6 +5011,6 @@
 									AB= b->r - a->r;   // AB = vector of the current edge
 									Metric MA = background ? BTh.MetricAt(a->r) :a->m ;  //Get metric associated to A
-									Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
-									double ledge = (MA(AB) + MB(AB))/2;                 //Get edge length in metric
+									Metric MB = background ? BTh.MetricAt(b->r) :b->m ;  //Get metric associated to B
+									double ledge = (MA(AB) + MB(AB))/2;                  //Get edge length in metric
 
 									/* We are now creating the mesh edges from the geometrical edge selected above.
@@ -5036,7 +5031,11 @@
 									//else, divide the edge
 									else {
-										//compute number of subedges (division of the edge)
+										//compute number of subedges (division of the edge), Maximum is 10
 										NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
-										//A and B are the position of points on the edge
+										/*Now, we are going to divide the edge according to the metric.
+										 * Get segment by sement along the edge.
+										 * Build lSubEdge, which holds the distance between the first vertex
+										 * of the edge and the next point on the edge according to the 
+										 * discretization (each SubEdge is AB)*/
 										R2 A,B;
 										A=a->r;
@@ -5049,5 +5048,5 @@
 											MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
 											AB = A-B;
-											lSubEdge[kk]= (ledge+=(MAs(AB)+MBs(AB))/2);
+											lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2);
 										}
 									}
@@ -5056,5 +5055,5 @@
 
 									/*Now, create corresponding points*/
-									while (lcurve<=s && s<=lcurveb && nbv<nbvend){
+									while(s>=lcurve && s<=lcurveb && nbv<nbvend){
 
 										double ss = s-lcurve;
@@ -5095,8 +5094,12 @@
 										va = vb;
 									}
+									
+									/*We just added one edge to the curve: Go to the next one*/
 									lcurve = lcurveb;
 									e->SetMark();
 									a=b;
-									if (b->Required() ) break;
+
+									/*If b is required, we are on a new curve->break*/
+									if (b->Required()) break;
 									int kprev=k;
 									k = e->AdjVertexIndex[kprev];// next vertices
@@ -5105,5 +5108,8 @@
 								}// for(;;)
 								vb = b->MeshVertexHook;
+
+								/*Number of edges in the last disretized curve*/
 								NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
+								/*Number of internal vertices in the last disretized curve*/
 								NbNewPoints = NbEdgeCurve-1;
 								if(!kstep){
@@ -5112,5 +5118,5 @@
 								}
 								nbvend=nbv+NbNewPoints; 
-								lstep = lcurve / NbEdgeCurve;
+								lstep = lcurve / NbEdgeCurve; //approximately one
 							}// end of curve --
 							if (edges) { // last edges of the curves