Index: ../trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- ../trunk-jpl/src/c/bamg/Mesh.cpp	(revision 18492)
+++ ../trunk-jpl/src/c/bamg/Mesh.cpp	(revision 18493)
@@ -241,7 +241,7 @@
 			if (Gh.NbRef>0) Gh.NbRef--;
 			else if (Gh.NbRef==0) delete &Gh;
 		}
-		if(&BTh && (&BTh != this) && false){
+		if(&BTh && (&BTh != this)){
 			if (BTh.NbRef>0) BTh.NbRef--;
 			else if (BTh.NbRef==0) delete &BTh;
 		}
@@ -4504,295 +4504,308 @@
 		/*Get options*/
 		int verbose=bamgopts->verbose;
 
-		/*Intermediaries*/
-		int                i,k;
-		int                nbcurves    = 0;
-		int                NbNewPoints,NbEdgeCurve;
-		double             lcurve,lstep,s;
-		const int          MaxSubEdge  = 10;
+		Gh.NbRef++;// add a ref to Gh
 
-		R2          AB;
-		GeomVertex *a, *b;
-		BamgVertex *va,*vb;
-		GeomEdge   *e;
+		/************************************************************************* 
+		 * method in 2 steps
+		 * 1 - compute the number of new edges to allocate
+		 * 2 - construct the edges
+		 * remark:
+		 * in this part we suppose to have a background mesh with the same geometry 
+		 * 
+		 * To construct the discretization of the new mesh we have to 
+		 * rediscretize the boundary of background Mesh 
+		 * because we have only the pointeur from the background mesh to the geometry.
+		 * We need the abcisse of the background mesh vertices on geometry
+		 * so a vertex is 
+		 * 0 on GeomVertex ;
+		 * 1 on GeomEdge + abcisse
+		 * 2 internal 
+		 *************************************************************************/
 
-		// add a ref to GH to make sure that it is not destroyed by mistake
-		Gh.NbRef++;
+		//Check that background mesh and current mesh do have the same geometry
+		_assert_(&BTh.Gh==&Gh);
+		BTh.NbRef++; // add a ref to BackGround Mesh
 
-		//build background mesh flag (1 if background, else 0)
-		bool background=(&BTh != this);
+		//Initialize new mesh
+		BTh.SetVertexFieldOn();
+		int* bcurve = new int[Gh.nbcurves]; // 
 
-		/*Build VerticesOnGeomVertex*/
+		/* There are 2 ways to make the loop 
+		 * 1) on the geometry 
+		 * 2) on the background mesh
+		 *  if you do the loop on geometry, we don't have the pointeur on background,
+		 *  and if you do the loop in background we have the pointeur on geometry
+		 * so do the walk on  background */
 
-		//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++;
-		}
-		//allocate
-		VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];  
-		if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
-		_assert_(nbv==0);
-		//Build VerticesOnGeomVertex
-		for (i=0;i<Gh.nbv;i++){
-			/* Add vertex only if required*/
-			if (Gh[i].Required()) {//Gh  vertices Required
+		NbVerticesOnGeomVertex=0;
+		NbVerticesOnGeomEdge=0;
 
-				//Add the vertex
-				_assert_(nbv<maxnbv);
-				vertices[nbv]=Gh[i];
+		/*STEP 1 copy of Required vertices*/
 
-				//Add pointer from geometry (Gh) to vertex from mesh (Th)
-				Gh[i].MeshVertexHook=vertices+nbv;
+		int i; 
+		for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
+		printf("\n");
+		if(NbVerticesOnGeomVertex >= maxnbv){
+			_error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
+		}
 
-				//Build VerticesOnGeomVertex for current point
-				VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
+		VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
+		VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
 
-				//nbv increment
+		//At this point there is NO vertex but vertices should have been allocated by Init
+		_assert_(vertices);
+		for (i=0;i<Gh.nbv;i++){
+			if (Gh[i].Required()) {//Gh vertices Required
+				vertices[nbv]  =Gh[i];
+				vertices[nbv].i=I2(0,0);
+				Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th
+				VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
 				nbv++;
 			}
+			else Gh[i].MeshVertexHook=0;
+		} 
+		for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){ 
+			VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
+			if (vog.IsRequiredVertex()){
+				GeomVertex* gv=vog;
+				BamgVertex *bv = vog;
+				_assert_(gv->MeshVertexHook); // use of Geom -> Th
+				VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv);
+				gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh
+			}
 		}
+		_assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /*This might be due to MaxCornerAngle too small*/
 
-		/*Build VerticesOnGeomEdge*/
+		/*STEP 2: reseed boundary edges*/
 
-		//check that edges is still empty (Init)
-		_assert_(!edges);
+		//  find the begining of the curve in BTh
+		Gh.UnMarkEdges();	
+		int bfind=0;
+		for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1; 
 
-		/* Now we are going to create the first edges corresponding
-		 * to the one present in the geometry provided.
-		 * We proceed in 2 steps
-		 *  -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++){
+		/*Loop over the backgrounf mesh BTh edges*/
+		for (int iedge=0;iedge<BTh.nbe;iedge++){      
+			Edge &ei = BTh.edges[iedge];
 
-			//initialize number of edges and number of edges max
-			long nbex=0;
-			nbe=0;
-			long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
-			Gh.UnMarkEdges();	
-			nbcurves=0;
+			/*Loop over the 2 vertices of the current edge*/
+			for(int je=0;je<2;je++){
 
-			//go through the edges of the geometry
-			for (i=0;i<Gh.nbe;i++){
+				/* If one of the vertex is required we are in a new curve*/
+				if (ei[je].GeomEdgeHook->IsRequiredVertex()){ 
 
-				//ei = current Geometrical edge
-				GeomEdge &ei=Gh.edges[i];   
+					/*Get curve number*/
+					int nc=ei.GeomEdgeHook->CurveNumber;
 
-				//loop over the two vertices of the edge ei
-				for(int j=0;j<2;j++) {
+					//_printf_("Dealing with curve number " << nc << "\n");
+					//_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n");
+					//if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
+					//	_printf_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no\n");
+					//	_printf_("Do we have the right extremity? curve last  vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no\n");
+					//}
+					//BUG FIX from original bamg
+					/*Check that we are on the same edge and right vertex (0 or 1) */
+					if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge  && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
+						bcurve[nc]=iedge*2+je;
+						bfind++;	
+					}
+					else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge  && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
+						bcurve[nc]=iedge*2+je;
+						bfind++;	
+					}
+				}
+			}
+		} 
+		if (bfind!=Gh.nbcurves){
+			delete [] bcurve;
+			_error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)");
+		}
 
-					/*Take only required vertices (corner->beginning of a new curve)*/
-					if (!ei.Mark() && ei[j].Required()){ 
+		// method in 2 + 1 step 
+		//  0.0) compute the length and the number of vertex to do allocation
+		//  1.0) recompute the length
+		//  1.1) compute the  vertex 
 
-						long  nbvend=0;
-						Edge* PreviousNewEdge=NULL;
-						lstep = -1;
+		long nbex=0,NbVerticesOnGeomEdgex=0;
+		for (int step=0; step <2;step++){
 
-						/*If Edge is required (do that only once for the 2 vertices)*/
-						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);
+			long NbOfNewPoints=0;
+			long NbOfNewEdge=0;
+			long iedge;
+			Gh.UnMarkEdges();	
+			double L=0;
 
-									//check that edges has been allocated
-									_assert_(edges);
-									edges[nbe].v[0]=a->MeshVertexHook;
-									edges[nbe].v[1]=b->MeshVertexHook;;
-									edges[nbe].ReferenceNumber = e->ReferenceNumber;
-									edges[nbe].GeomEdgeHook = e;
-									edges[nbe].adj[0] = 0;
-									edges[nbe].adj[1] = 0;
-									nbe++;
-								}
-							}
-						}
+			/*Go through all geometrical curve*/
+			for (int icurve=0;icurve<Gh.nbcurves;icurve++){ 
 
-						/*If Edge is not required: we are on a curve*/
-						else {
-							for (int kstep=0;kstep<=step;kstep++){
-								//kstep=0, compute number of edges (discretize curve)
-								//kstep=1  create the points and edge
-								PreviousNewEdge=0;
-								NbNewPoints=0;
-								NbEdgeCurve=0;
-								if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
-								lcurve =0;
-								s = lstep; //-1 initially, then length of each sub edge
+				/*Get edge and vertex (index) of background mesh on this curve*/
+				iedge=bcurve[icurve]/2;
+				int jedge=bcurve[icurve]%2;
 
-								/*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->MeshVertexHook; // va = pointer toward mesh vertex associated
-								e->SetMark();   // Mark edge
+				/*Get edge of Bth with index iedge*/
+				Edge &ei = BTh.edges[iedge];
 
-								/*Loop until we reach the end of the curve*/
-								for(;;){ 
-									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
+				/*Initialize variables*/
+				double Lstep=0;             // step between two points   (phase==1) 
+				long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 
 
-									/* 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). */
+				/*Do phase 0 to step*/
+				for(int phase=0;phase<=step;phase++){
 
-									//By default, there is only one subedge that is the geometrical edge itself
-									int NbSubEdge = 1;
+					/*Current curve pointer*/
+					Curve *curve= Gh.curves+icurve;
 
-									//initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
-									double lSubEdge[MaxSubEdge];
+					/*Get index of current curve*/
+					int icurveequi= Gh.GetId(curve);
 
-									//Build Subedges according to the edge length
-									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 {
-										//compute number of subedges (division of the edge), Maximum is 10
-										NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
-										/*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;
-										Metric MAs=MA,MBs;
-										ledge=0; 
-										double x =0, xstep= 1./NbSubEdge;
-										for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
-											x += xstep;
-											B =  e->F(k ? x : 1-x);
-											MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
-											AB = A-B;
-											lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2);
-										}
-									}
+					/*For phase 0, check that we are at the begining of the curve only*/
+					if(phase==0 &&  icurveequi!=icurve)  continue;
 
-									double lcurveb = lcurve+ledge;
+					int   k0=jedge,k1;
+					Edge* pe=  BTh.edges+iedge;
+					int   iedgeequi=bcurve[icurveequi]/2;
+					int   jedgeequi=bcurve[icurveequi]%2;
 
-									/*Now, create corresponding points*/
-									while(s>=lcurve && s<=lcurveb && nbv<nbvend){
+					int k0equi=jedgeequi,k1equi;		  
+					Edge * peequi= BTh.edges+iedgeequi;
+					GeomEdge *ongequi = peequi->GeomEdgeHook;
 
-										/*Schematic of current curve
-										 *
-										 *  a                   vb                  b          // vertex
-										 *  0              ll0     ll1              ledge      // length from a
-										 *  + --- + - ... - + --S-- + --- + - ... - +          // where is S
-										 *  0              kk0     kk1              NbSubEdge  // Sub edge index
-										 *
-										 */
+					double sNew=Lstep;// abscisse of the new points (phase==1) 
+					L=0;// length of the curve
+					long i=0;// index of new points on the curve
+					GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook;
+					BamgVertex *A0;
+					A0 = GA0->MeshVertexHook;  // the vertex in new mesh
+					BamgVertex *A1;
+					VertexOnGeom *GA1;
+					Edge* PreviousNewEdge = 0;
 
-										double ss = s-lcurve;
+					// New Curve phase 
+					_assert_(A0-vertices>=0 && A0-vertices<nbv);
+					if(ongequi->Required()){
+						GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook;
+						A1 = GA1->MeshVertexHook;  //
+					}       
+					else {
+						for(;;){
+							Edge &ee=*pe; 
+							Edge &eeequi=*peequi; 
+							k1 = 1-k0; // next vertex of the edge 
+							k1equi= 1 - k0equi;
+							_assert_(pe && ee.GeomEdgeHook);
+							ee.GeomEdgeHook->SetMark();
+							BamgVertex & v0=ee[0], & v1=ee[1];
+							R2 AB=(R2)v1-(R2)v0;
+							double L0=L,LAB;
+							LAB=LengthInterpole(v0.m,v1.m,AB);
+							L+= LAB;
 
-										/*Find the SubEdge containing ss using Dichotomy*/
-										int kk0=-1,kk1=NbSubEdge-1,kkk;
-										double ll0=0,ll1=ledge,llk;
-										while (kk1-kk0>1){
-											if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
-											 kk1=kkk,ll1=llk;
-											else
-											 kk0=kkk,ll0=llk;
-										}
-										_assert_(kk1!=kk0);
+							if (phase){
+								// computation of the new points for the given curve
+								while ((i!=NbCreatePointOnCurve) && sNew<=L) { 
 
-										/*Curvilinear coordinate in [0 1] of ss in current edge*/
-										// WARNING: This is what we would do
-										// ssa = (ss-ll0)/(ll1-ll0);
-										// aa = (kk0+ssa)/NbSubEdge
-										// This is what Bamg does:
-										double sbb = (ss-ll0)/(ll1-ll0);
-										/*Curvilinear coordinate in [0 1] of ss in current curve*/
-										double bb = (kk1+sbb)/NbSubEdge;
-										double aa = 1-bb;
+									//some checks
+									_assert_(sNew>=L0);
+									_assert_(LAB);
+									_assert_(vertices && nbv<maxnbv);
+									_assert_(edges && nbe<nbex);
+									_assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
 
-										// new vertex on edge
-										vb = &vertices[nbv++];
-										vb->m = Metric(aa,a->m,bb,b->m);
-										vb->ReferenceNumber = e->ReferenceNumber;
-										vb->DirOfSearch =NoDirOfSearch;
-										double abcisse = k ? bb : aa;
-										vb->r =  e->F(abcisse);
-										VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);        
-
-										// to take into account the direction of the edge
-										s += lstep;
-										edges[nbe].v[0]=va;
-										edges[nbe].v[1]=vb;
-										edges[nbe].ReferenceNumber =e->ReferenceNumber;
-										edges[nbe].GeomEdgeHook = e;
-										edges[nbe].adj[0] = PreviousNewEdge;
-										if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
-										PreviousNewEdge=edges+nbe;
-										nbe++;
-										va = vb;
+									// new vertex on edge
+									A1=vertices+nbv++;
+									GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
+									Edge* e = edges + nbe++;
+									double se= (sNew-L0)/LAB;
+									if (se<0 || se>=1.000000001){
+										_error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
 									}
+									se = abscisseInterpole(v0.m,v1.m,AB,se,1);
+									if (se<0 || se>1){
+										_error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
+									}
+									se = k1         ? se : 1. - se;
+									se = k1==k1equi ? se : 1. - se;
+									VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 
+									ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 
+									A1->ReferenceNumber = eeequi.ReferenceNumber;
+									A1->DirOfSearch =NoDirOfSearch;
+									e->GeomEdgeHook = ongequi;
+									e->v[0]=A0;
+									e->v[1]=A1;
+									e->ReferenceNumber = eeequi.ReferenceNumber;
+									e->adj[0]=PreviousNewEdge;
 
-									/*We just added one edge to the curve: Go to the next one*/
-									lcurve = lcurveb;
-									e->SetMark();
-									a=b;
-
-									/*If b is required, we are on a new curve->break*/
-									if (b->Required()) break;
-									int kprev=k;
-									k = e->AdjVertexIndex[kprev];// next vertices
-									e = e->Adj[kprev];
-									_assert_(e);
-								}// 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){
-									NbVerticesOnGeomEdge0 += NbNewPoints;
-									nbcurves++;
+									if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
+									PreviousNewEdge=e;
+									A0=A1;
+									sNew += Lstep;
+									if (++i== NbCreatePointOnCurve) break;
 								}
-								nbvend=nbv+NbNewPoints; 
-								lstep = lcurve / NbEdgeCurve; //approximately one
-							}// end of curve --
-							if (edges) { // last edges of the curves 
-								edges[nbe].v[0]=va;
-								edges[nbe].v[1]=vb;
-								edges[nbe].ReferenceNumber = e->ReferenceNumber;
-								edges[nbe].GeomEdgeHook = e;
-								edges[nbe].adj[0] = PreviousNewEdge;
-								edges[nbe].adj[1] = 0;
-								if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];
-								nbe++;
 							}
-							else nbe += NbEdgeCurve;
-						} // end on  curve ---
+
+							//some checks
+							_assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
+							if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
+								_assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
+								GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
+								A1=GA1->MeshVertexHook;// the vertex in new mesh
+								_assert_(A1-vertices>=0 && A1-vertices<nbv);
+								break;
+							}
+							if (!ee.adj[k1]) {
+								_error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices);
+							}
+							pe = ee.adj[k1]; // next edge
+							k0 = pe->Intersection(ee); 
+							peequi= eeequi.adj[k1equi];  // next edge
+							k0equi=peequi->Intersection(eeequi);            
+						}// for(;;) end of the curve
 					}
+
+					if (phase){ // construction of the last edge
+						Edge* e=edges + nbe++;
+						e->GeomEdgeHook  = ongequi;
+						e->v[0]=A0;
+						e->v[1]=A1;
+						e->ReferenceNumber = peequi->ReferenceNumber;
+						e->adj[0]=PreviousNewEdge;
+						e->adj[1]=0;
+						if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
+						PreviousNewEdge = e;
+
+						_assert_(i==NbCreatePointOnCurve);
+					}
+
+					if (!phase)  { // 
+						long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
+						Lstep = L/NbSegOnCurve; 
+						NbCreatePointOnCurve = NbSegOnCurve-1;
+						NbOfNewEdge += NbSegOnCurve;
+						NbOfNewPoints += NbCreatePointOnCurve;
+					}
 				}
-			} // for (i=0;i<nbe;i++)
-			if(!step) {
-				_assert_(!edges);
-				_assert_(!VerticesOnGeomEdge);
+			}//  end of curve loop 
 
-				edges = new Edge[nbex=nbe];
-				if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
-
-				// do the vertex on a geometrical vertex
-				_assert_(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0);
-				NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;       
+			//Allocate memory
+			if(step==0){
+				if(nbv+NbOfNewPoints > maxnbv) {
+					_error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv);
+				}
+				edges = new Edge[NbOfNewEdge];
+				nbex = NbOfNewEdge;
+				if(NbOfNewPoints) {
+					VerticesOnGeomEdge    = new VertexOnGeom[NbOfNewPoints];
+					NbVertexOnBThEdge     = NbOfNewPoints;
+					VertexOnBThEdge       = new  VertexOnEdge[NbOfNewPoints];
+					NbVerticesOnGeomEdgex = NbOfNewPoints;
+				}
+				NbOfNewPoints =0;
+				NbOfNewEdge = 0;
 			}
-			else{
-				_assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
-			}
 		}
+		_assert_(nbe!=0);
+		delete [] bcurve;
 
-
 		//Insert points inside existing triangles
 		if (verbose>4) _printf_("      -- current number of vertices = " << nbv << "\n");
 		if (verbose>3) _printf_("      Creating initial Constrained Delaunay Triangulation...\n");