Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5447)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5448)
@@ -4875,21 +4875,23 @@
 		/*Generate mesh from geometry*/
 
-		Gh.NbRef++;// add a ref to GH
-
-		int i,j,k;
-		int nbcurves=0,NbNewPoints,NbEdgeCurve;
-		double lcurve,lstep,s;
-		const int MaxSubEdge = 10;
-
-		R2 AB;
-		GeometricalVertex *a,*b;
-		BamgVertex *va,*vb;
-		GeometricalEdge *e;
+		/*Intermediaries*/
+		int                i,j,k;
+		int                nbcurves    = 0;
+		int                NbNewPoints,NbEdgeCurve;
+		double             lcurve,lstep,s;
+		const int          MaxSubEdge  = 10;
+
+		R2                 AB;
+		GeometricalVertex *a, *b;
+		BamgVertex        *va, *vb;
+		GeometricalEdge   *e;
+
+		// add a ref to GH to make sure that it is not destroyed by mistake
+		Gh.NbRef++;
 
 		/*Get options*/
 		int verbose=bamgopts->verbose;
 
-		//initialize this
-		if (verbose>3) printf("      Generating Boundary vertices\n");
+		//initialize Mesh
 		Init(imaxnbv);
 		nbv=0;
@@ -4898,18 +4900,16 @@
 
 		//build background mesh flag (1 if background, else 0)
-		int  background=(&BTh != this);
-
-		//Compute number of vertices on geometrical vertex
+		bool background=(&BTh != this);
+
+		/*Build VerticesOnGeomVertex*/
+
+		//Compute the number of geometrical vertices that we are going to use to mesh
 		for (i=0;i<Gh.nbv;i++){
 			if (Gh[i].Required()) NbVerticesOnGeomVertex++;
 		}
-
-		//initialize VerticesOnGeomVertex
+		//allocate
 		VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];  
-		if( NbVerticesOnGeomVertex >= maxnbv) {
-			ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);
-		}
-
-		//Add all the geometrical vertices to the mesh
+		if(NbVerticesOnGeomVertex >= maxnbv) ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);
+		//Build VerticesOnGeomVertex
 		nbv=0;
 		for (i=0;i<Gh.nbv;i++){
@@ -4917,11 +4917,7 @@
 			if (Gh[i].Required()) {//Gh  vertices Required
 
-				//Add the vertex (provided that nbv<maxnbv)
-				if (nbv<maxnbv){
-					vertices[nbv]=Gh[i];
-				}
-				else{
-					ISSMERROR("Maximum number of vertices (maxnbv = %i) too small",maxnbv);
-				}
+				//Add the vertex
+				ISSMASSERT(nbv<maxnbv);
+				vertices[nbv]=Gh[i];
 				
 				//Add pointer from geometry (Gh) to vertex from mesh (Th)
@@ -4929,5 +4925,5 @@
 
 				//Build VerticesOnGeomVertex for current point
-				VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]);
+				VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].to,Gh[i]);
 
 				//nbv increment
@@ -4936,8 +4932,8 @@
 		}
 
-		//check that edges has been allocated
-		if (edges){
-			ISSMERROR("edges is empty");
-		}
+		/*Build VerticesOnGeomEdge*/
+
+		//check that edges is still empty (Init)
+		ISSMASSERT(!edges);
 
 		/* Now we are going to create the first edges corresponding
@@ -4957,48 +4953,46 @@
 
 			//go through the edges of the geometry
-			for (i=0;i<Gh.nbe;i++) {
+			for (i=0;i<Gh.nbe;i++){
 
 				//ei = current Geometrical edge
 				GeometricalEdge &ei=Gh.edges[i];   
 
+				//loop over the two vertices of the edge ei
 				for(int j=0;j<2;j++) {
 
-					/*The first time, no edge is marked but this might change during the loop*/
+					/*Take only required vertices (corner)*/
 					if (!ei.Mark() && ei[j].Required()){ 
 
 						long  nbvend=0;
 						Edge* PreviousNewEdge=NULL;
-
 						lstep = -1;
 
-						/*If Edge is required*/
-						if(ei.Required()){
+						/*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 (j==0){
-								if(step==0) nbe++;
-								else{ 
-									e=&ei;
-									a=ei(0);
-									b=ei(1);
-
-									//check that edges has been allocated
-									if (!edges) ISSMERROR("edges has not been allocated...");
-									edges[nbe].v[0]=a->to;
-									edges[nbe].v[1]=b->to;;
-									edges[nbe].ReferenceNumber = e->ReferenceNumber;
-									edges[nbe].GeometricalEdgeHook = e;
-									edges[nbe].adj[0] = 0;
-									edges[nbe].adj[1] = 0;
-									nbe++;
-								}
+							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->to;
+								edges[nbe].v[1]=b->to;;
+								edges[nbe].ReferenceNumber = e->ReferenceNumber;
+								edges[nbe].GeometricalEdgeHook = e;
+								edges[nbe].adj[0] = 0;
+								edges[nbe].adj[1] = 0;
+								nbe++;
 							}
 						}
 
-						/*If Edge is not required: on a curve*/
+						/*If Edge is not required: we are on a curve*/
 						else {
 							for (int kstep=0;kstep<=step;kstep++){
-								//step=0, do nothing
-								//step=1, compute the length of the curve
-								//step=2  create the points and edge
+								//kstep=0, do nothing
+								//kstep=1, compute the length of the curve
+								//kstep=2  create the points and edge
 								PreviousNewEdge=0;
 								NbNewPoints=0;
@@ -5008,27 +5002,23 @@
 								s = lstep;
 
-								// i = edge number, j=[0;1] vertex number in edge
-
-								k=j;            // k = vertex number in edge (0 or 1)
+								/*reminder: i = edge number, j=[0;1] vertex index in edge*/
+								k=j;            // k = vertex index in edge (0 or 1)
 								e=&ei;          // e = reference of current edge
 								a=ei(k);        // a = pointer toward the kth vertex of the current edge
-								va = a->to;     // va = pointer toward vertex associated
+								va = a->to;     // va = pointer toward mesh vertex associated
 								e->SetMark();   // Mark edge
 
-								//if SameGeo We have go to the background geometry 
-								//to find the discretisation of the curve
+								/*If we have a Background mesh, we can use it to discretize the curve*/
 								for(;;){ 
-									k = 1-k;
-									b = (*e)(k);// b = pointer toward the other vertex of the current edge
+									k = 1-k;            // other vertx index of the curve
+									b = (*e)(k);        // b = pointer toward the other vertex of the current edge
 									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
-
-									/* We are now creating the edges of the mesh from the
-									 * geometrical edge selected above.
-									 * The edge will be divided according to the metric
-									 * previously computed and cannot be divided more
-									 * than 10 times (MaxSubEdge). */
+									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.
+									 * The edge will be divided according to the metric previously computed and cannot
+									 * be divided more than 10 times (MaxSubEdge). */
 
 									//By default, there is only one subedge that is the geometrical edge itself
@@ -5039,6 +5029,8 @@
 
 									//Build Subedges according to the edge length
-									//if ledge < 1.5 (between one and 2), take the edge as is
-									if (ledge < 1.5) lSubEdge[0] = ledge;
+									if (ledge < 1.5){
+										//if ledge < 1.5 (between one and 2), take the edge as is
+										lSubEdge[0] = ledge;
+									}
 									//else, divide the edge
 									else {
@@ -5054,5 +5046,5 @@
 											x += xstep;
 											B =  e->F(k ? x : 1-x);
-											MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB);
+											MBs= background ? BTh.MetricAt(B) : Metric(1-x, MA, x ,MB);
 											AB = A-B;
 											lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);
@@ -5138,21 +5130,16 @@
 			} // for (i=0;i<nbe;i++)
 			if(!step) {
-				if (edges){
-					ISSMERROR("edges");
-				}
-				if (VerticesOnGeomEdge){
-					ISSMERROR("VerticesOnGeomEdge");
-				}
+				ISSMASSERT(!edges);
+				ISSMASSERT(!VerticesOnGeomEdge);
+
 				edges = new Edge[nbex=nbe];
-				if(NbVerticesOnGeomEdge0)
-				 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
-				if (!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0){
-					ISSMERROR("!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0");
-				}
+				if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
+
 				// do the vertex on a geometrical vertex
+				ISSMASSERT(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0);
 				NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;       
 			}
-			else if (NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0){
-				ISSMERROR("NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0");
+			else{
+				ISSMASSERT(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
 			}
 		}