Index: /issm/trunk/src/c/Bamgx/meshtype.h
===================================================================
--- /issm/trunk/src/c/Bamgx/meshtype.h	(revision 3276)
+++ /issm/trunk/src/c/Bamgx/meshtype.h	(revision 3277)
@@ -118,4 +118,39 @@
 		}
 	}
+	template<class T> inline void  HeapSort(long** porder,T* c,int n){
+		long l,j,r,i;
+		T    crit;
+		long pos;
+		long* order=NULL;
+		order=(long*)xmalloc(n*sizeof(int));
+		for(i=0;i<n;i++) order[i]=i+1;
+		c--;                    //the array must starts at 1 and not 0 
+		order--;
+		if(n<=1) return;        //return if size <=1
+		l=n/2+1;                //initialize l and r
+		r=n;
+		for(;;){
+			if(l<=1){
+				crit  =c[r]; pos=order[r];
+				c[r--]=c[1]; order[r+1]=order[1];
+				if (r==1){
+					c[1]=crit; order[1]=pos;
+					order++;
+					*porder=order;
+					return;
+				}
+			}
+			else  {crit=c[--l]; pos=order[l];}
+			j=l;
+			for(;;){
+				i=j;
+				j=2*j;
+				if  (j>r) {c[i]=crit;order[i]=pos;break;}
+				if ((j<r) && (c[j] < c[j+1]))j++;
+				if (crit < c[j]) {c[i]=c[j];order[i]=order[j];} 
+				else{c[i]=crit;order[i]=pos;break;}
+			}
+		}
+	}
 
 	//Inline functions
Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3276)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 3277)
@@ -929,17 +929,17 @@
 	/*}}}1*/
 	/*FUNCTION  Geometry::ProjectOnCurve {{{1*/
-	GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,double s,Vertex &V,VertexOnGeom &GV ) const {
+	GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,Vertex &V,VertexOnGeom &GV) const {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
 
 		double save_s=s;
 		int NbTry=0;
-retry:    
+retry:
 		s=save_s;
-		GeometricalEdge* on = e.onGeometry;
+		GeometricalEdge* on=e.onGeometry;
 		if (!on){
-			throw ErrorException(__FUNCT__,exprintf("!on"));
+			throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: edge provided shonld be on geometry"));
 		}
 		if (!e[0].onGeometry ||  !e[1].onGeometry){
-			throw ErrorException(__FUNCT__,exprintf("!e[0].on ||  !e[1].on"));
+			throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry"));
 		}
 		const Vertex &v0=e[0],&v1=e[1];
@@ -969,5 +969,5 @@
 					printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n");
 					printf("That bug might come from:\n");
-					printf(" 1)  a mesh edge  contening more than %i geometrical edges\n",mxe/2);
+					printf(" 1)  a mesh edge  containing more than %i geometrical edges\n",mxe/2);
 					printf(" 2)  code bug : be sure that we call   Triangles::SetVertexFieldOn() before\n");
 					printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n");
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3276)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3277)
@@ -401,5 +401,5 @@
 		}
 
-		//VerVerticesOnGeometricEdge
+		//VerticesOnGeometricEdge
 		if(bamgmesh->VerticesOnGeometricEdge){
 			if(verbose>5) printf("      processing VerticesOnGeometricEdge\n");
@@ -419,16 +419,31 @@
 		}
 
+		//VerticesOnGeometricVertex
+		if(bamgmesh->NumVerticesOnGeometricVertex){
+			if(verbose>5) printf("      processing VerticesOnGeometricVertex\n");
+			NbVerticesOnGeomVertex=bamgmesh->NumVerticesOnGeometricVertex;
+			VerticesOnGeomVertex  = new  VertexOnGeom[NbVerticesOnGeomVertex] ;
+			for (i=0;i<NbVerticesOnGeomVertex;i++){
+				long  i1,i2;
+				i1=(long)bamgmesh->VerticesOnGeometricVertex[i*2+0]-1; //for C indexing
+				i2=(long)bamgmesh->VerticesOnGeometricVertex[i*2+1]-1; //for C indexing
+				VerticesOnGeomVertex[i]=VertexOnGeom(vertices[i1],Gh.vertices[i2]);
+			}
+		}
+		else{
+			if(verbose>5) printf("      no VertexOnGeometricVertex found\n");
+		}
+
 		//Edges
 		if (bamgmesh->Edges){
 			int i1,i2;
-			R2 zero2(0,0);
-			double *len =0;
+			double* len=NULL;
 
 			if(verbose>5) printf("      processing Edges\n");
 			nbe=bamgmesh->NumEdges;
-			edges = new Edge[nbe];
+			edges= new Edge[nbe];
 
 			if (!hvertices) {
-				len = new double[nbv];
+				len= new double[nbv];
 				for(i=0;i<nbv;i++)
 				 len[i]=0;
@@ -442,12 +457,12 @@
 				edges[i].adj[0]=NULL;
 				edges[i].adj[1]=NULL;
-				R2 x12 = vertices[i2].r-vertices[i1].r;
-				double l12=sqrt( (x12,x12));        
-
-				if (!hvertices) {
+				R2 x12=vertices[i2].r-vertices[i1].r;
+				double l12=sqrt((x12,x12));
+
+				if (!hvertices){
 					vertices[i1].color++;
 					vertices[i2].color++;
-					len[i1]+= l12;
-					len[i2] += l12;
+					len[i1]+=l12;
+					len[i2]+=l12;
 				}
 				Hmin = Min(Hmin,l12);
@@ -458,7 +473,7 @@
 				for (i=0;i<nbv;i++) 
 				 if (vertices[i].color > 0) 
-				  vertices[i].m=  Metric(len[i] /(double) vertices[i].color);
+				  vertices[i].m=Metric(len[i] /(double) vertices[i].color);
 				 else 
-				  vertices[i].m=  Metric(Hmin);
+				  vertices[i].m=Metric(Hmin);
 				delete [] len;
 			}
@@ -476,11 +491,9 @@
 					}
 					else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
-						j0 =  i0%2;
-						i0 =  i0/2;
-						if (v!=edges[i0 ].v[j0]){
-							throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]"));
-						}
-						edges[i ].adj[ j ] =edges +i0;
-						edges[i0].adj[ j0] =edges +i ;
+						j0 = i0%2;
+						i0 = i0/2;
+						assert(v==edges[i0 ].v[j0]);
+						edges[i ].adj[j ] =edges +i0;
+						edges[i0].adj[j0] =edges +i ;
 						v->color = -3;
 					}
@@ -513,10 +526,11 @@
 			i2=bamgmesh->NumEdgesOnGeometricEdge;
 			for (i1=0;i1<i2;i1++) {
-				i=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+0];
-				j=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+1];
-				if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe)) {
-					throw ErrorException(__FUNCT__,exprintf("ReadMesh error: EdgesOnGeometricEdge edge provided (line %i: [%i %i]) is incorrect (must be positive, [0<i<nbe=%i 0<j<Gh.nbe=%i]",i1+1,i,j,nbe,Gh.nbe));
-				}
-				edges[i-1].onGeometry = Gh.edges + j-1;
+				i=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+0]-1; //C indexing
+				j=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+1]-1; //C indexing
+				//Check value
+				if(!(i>=0 && j>=0 && i<nbe && j<Gh.nbe)) {
+					throw ErrorException(__FUNCT__,exprintf("ReadMesh error: EdgesOnGeometricEdge edge provided (line %i: [%i %i]) is incorrect (must be positive, [0<i<nbe=%i 0<j<Gh.nbe=%i]",i1+1,i+1,j+1,nbe,Gh.nbe));
+				}
+				edges[i].onGeometry=Gh.edges+j;
 			}
 		}
@@ -555,12 +569,4 @@
 		/*Initialize output*/
 		BamgMeshInit(bamgmesh);
-
-
-		//renumber
-		if (bamgopts->renumber){
-			if(verbose>5) printf("      Renumbering...");
-			Renumber(bamgopts);
-			if(verbose>5) printf(" done\n"); 
-		}
 
 		//Build reft that holds the number the subdomain number of each triangle
@@ -769,9 +775,7 @@
 			for (i=0;i<NbVerticesOnGeomVertex;i++){
 				VertexOnGeom &v=VerticesOnGeomVertex[i];
-				if (!v.OnGeomVertex()){
-					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not"));
-				}
+				assert(v.OnGeomVertex());
 				bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing
-				bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number((GeometricalVertex*)v)+1; //back to Matlab indexing
 			}
 		}
@@ -4312,84 +4316,4 @@
 	}                  
 	/*}}}1*/
-/*FUNCTION Triangles::Renumber {{{1*/
-void Triangles::Renumber(BamgOpts* bamgopts){
-
-	/*Intermediary*/
-	int i,j,k;
-	int*    r=NULL;
-	int*    rprime=NULL;
-	rprime  =(int*)   xmalloc(nbv*sizeof(int));
-
-	//renumbering with respect to lower left corner
-	if (bamgopts->renumber==1){
-
-		//compute distance to pmin
-		double* distance;
-		distance=(double*)xmalloc(nbv*sizeof(double));
-		for (i=0;i<nbv;i++) distance[i]=pow(vertices[i].r.x - pmin.x,2)+pow(vertices[i].r.y - pmin.y,2);
-
-		//get ordering
-		HeapSort(&r,distance,nbv);
-		xfree((void**)&distance);
-	}
-	//renumbering randomly
-	else if (bamgopts->renumber==2){
-
-		//allocate memory to r
-		r=(int*)xmalloc(nbv*sizeof(int));
-		int PrimeNumber= BigPrimeNumber(nbv) ;
-		int k0=rand()%nbv; 
-		for (int i=0; i<nbv; i++){
-			r[i]=(k0=(k0+PrimeNumber)%nbv)+1;
-		}
-	}
-	//else: error
-	else{
-		throw ErrorException(__FUNCT__,exprintf("renumbering option %i not supported yet",bamgopts->renumber)); 
-	}
-
-	//Keep copy of old vertices
-	Vertex* oldvertices=NULL;
-	oldvertices=(Vertex*)xmalloc(nbv*sizeof(Vertex));
-	for (i=0;i<nbv;i++){
-		oldvertices[i]=vertices[i];
-	}
-
-	//renumber vertices
-	for (i=0;i<nbv;i++){
-		//check r[i]
-		if (r[i]>nbv){
-			throw ErrorException(__FUNCT__,exprintf("r[i]>nbv (r[i]=%i and nbv=%i)",r[i],nbv));
-		}
-		if (r[i]<=0){
-			throw ErrorException(__FUNCT__,exprintf("r[i]<=0 (r[i]=%i)",r[i]));
-		}
-		vertices[i]=oldvertices[r[i]-1];
-		rprime[r[i]-1]=i+1;
-	}
-	//clean up
-	xfree((void**)&oldvertices);
-
-	//renumber triangles
-	for (i=0; i<nt;i++){
-		for (j=0;j<3;j++){
-			if (triangles[i](j)){
-				triangles[i](j)=&vertices[rprime[Number(triangles[i](j))] -1 ];
-			}
-		}
-	}
-
-	//renumber edges
-	for (i=0; i<nbe;i++){
-		for (j=0;j<2;j++){
-			edges[i].v[j]=&vertices[rprime[Number(edges[i].v[j])] -1 ];
-		}
-	}
-
-	//clean up
-	xfree((void**)&r);
-	xfree((void**)&rprime);
-}
-/*}}}1*/
 	/*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/
 	void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) {
@@ -4466,16 +4390,19 @@
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
 
-		// warning be carfull because pointeur
+		// warning be carfull because pointer
 		// from on mesh to over mesh 
-		//  --  so do ReNumbering a the beginning
+		//  --  so do ReNumbering at the beginning
 		Vertex * ve = vertices+nbv;
 		long it,ie,i;
 
+		printf("renumbering triangles\n");
 		for ( it=0;it<nbt;it++) 
 		 triangles[it].ReNumbering(vertices,ve,renu);
 
+		printf("renumbering edges\n");
 		for ( ie=0;ie<nbe;ie++) 
 		 edges[ie].ReNumbering(vertices,ve,renu);
 
+		printf("renumbering vertices on geom\n");
 		for (i=0;i< NbVerticesOnGeomVertex;i++)
 		  {
@@ -4485,4 +4412,5 @@
 		  }
 
+		printf("renumbering vertices on edge\n");
 		for (i=0;i< NbVerticesOnGeomEdge;i++)
 		  {
@@ -4492,4 +4420,5 @@
 		  }
 
+		printf("renumbering vertices on Bth vertex\n");
 		for (i=0;i< NbVertexOnBThVertex;i++)
 		  {
@@ -4514,20 +4443,19 @@
 			 i=it;
 			 Vertex ti=vertices[i],tj;
-			 while ( (j=renu[i]) >= 0) 
-				{ // i is old, and j is new 
+			 while ( (j=renu[i]) >= 0){
+				 // i is old, and j is new 
 				 renu[i] = -1-renu[i]; // mark 
-				 tj = vertices[j]; // save new
-				 vertices[j]= ti; // new <- old
+				 tj = vertices[j];     // save new
+				 vertices[j]= ti;      // new <- old
 				 i=j;     // next 
 				 ti = tj;
 				}  
 			}
-		if (quadtree) 
-		  {  delete quadtree;
+		if (quadtree) {delete quadtree;
 			quadtree = new QuadTree(this);
-		  }
-		for ( it=0;it<nbv;it++)
-		 renu[i]= -renu[i]-1;
-
+		}
+
+		for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
+		printf("end renumbervertex\n");
 	}
 	/*}}}1*/
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3276)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3277)
@@ -101,5 +101,4 @@
 			void Insert();
 			void ForceBoundary();
-			void Renumber(BamgOpts* bamgopts);
 			void FindSubDomain(int OutSide=0);
 			long TriangleReferenceList(long *) const;
@@ -159,5 +158,5 @@
 			  }
 			inline  void  SetVertexFieldOn(){
-				for (int i=0;i<nbv;i++)                    vertices[i].onGeometry=0;
+				for (int i=0;i<nbv;i++)                    vertices[i].onGeometry=NULL;
 				for (int j=0;j<NbVerticesOnGeomVertex;j++) VerticesOnGeomVertex[j].SetOn();
 				for (int k=0;k<NbVerticesOnGeomEdge;k++ )  VerticesOnGeomEdge[k].SetOn();
Index: /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h	(revision 3276)
+++ /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h	(revision 3277)
@@ -34,5 +34,5 @@
 
 			//Operators
-			operator Vertex * () const  {return mv;}
+			operator Vertex*() const  {return mv;}
 			operator GeometricalVertex * () const  {return gv;}
 			operator GeometricalEdge * () const  {return ge;}
Index: /issm/trunk/src/c/objects/BamgMesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/BamgMesh.cpp	(revision 3276)
+++ /issm/trunk/src/c/objects/BamgMesh.cpp	(revision 3277)
@@ -146,5 +146,5 @@
 
 	pfield=mxCreateDoubleMatrix(0,0,mxREAL);
-	mxSetM(pfield,2);
+	mxSetM(pfield,3);
 	mxSetN(pfield,bamgmesh->NumVerticesOnGeometricEdge);
 	mxSetPr(pfield,bamgmesh->VerticesOnGeometricEdge);
Index: /issm/trunk/src/c/objects/BamgOpts.cpp
===================================================================
--- /issm/trunk/src/c/objects/BamgOpts.cpp	(revision 3276)
+++ /issm/trunk/src/c/objects/BamgOpts.cpp	(revision 3277)
@@ -11,5 +11,4 @@
 	bamgopts->Metrictype=0;
 	bamgopts->KeepVertices=0;
-	bamgopts->renumber=0;
 	bamgopts->maxsubdiv=0;
 	bamgopts->power=0;
@@ -42,5 +41,4 @@
 	if (bamgopts->Metrictype!=0   && bamgopts->Metrictype!=1 && bamgopts->Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2"));
 	if (bamgopts->KeepVertices!=0 && bamgopts->KeepVertices!=1) throw ErrorException(__FUNCT__,exprintf("'KeepVertices' supported options are 0 and 1"));
-	if (bamgopts->renumber!=0     && bamgopts->renumber!=1   && bamgopts->renumber!=2)   throw ErrorException(__FUNCT__,exprintf("'renumber' supported options are 0, 1 and 2"));
 	if (bamgopts->errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0"));
 	if (bamgopts->nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0"));
Index: /issm/trunk/src/c/objects/BamgOpts.h
===================================================================
--- /issm/trunk/src/c/objects/BamgOpts.h	(revision 3276)
+++ /issm/trunk/src/c/objects/BamgOpts.h	(revision 3277)
@@ -14,5 +14,4 @@
 	int     Metrictype;
 	int     KeepVertices;
-	int     renumber;
 	double  maxsubdiv;
 	double  power;
Index: /issm/trunk/src/m/classes/@bamggeom/bamggeom.m
===================================================================
--- /issm/trunk/src/m/classes/@bamggeom/bamggeom.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamggeom/bamggeom.m	(revision 3277)
@@ -0,0 +1,54 @@
+function bg = bamggeom(varargin)
+%ICEFRONT - constructor for bamggeom object
+%
+%   Usage:
+%      bamggeom = bamggeom(varargin)
+
+switch nargin
+case 0
+	% if no input arguments, create a default object
+	bg.NumVertices=0;
+	bg.Vertices=[];
+
+	bg.NumEdges=0;
+	bg.Edges=zeros(0,3);
+
+	bg.NumTangentAtEdges=0;
+	bg.TangentAtEdges=[];
+
+	bg.NumCorners=0;
+	bg.Corners=[];
+
+	bg.NumRequiredVertices=0;
+	bg.RequiredVertices=[];
+
+	bg.NumRequiredEdges=0;
+	bg.RequiredEdges=[];
+
+	bg.NumCrackedEdges=0;
+	bg.CrackedEdges=zeros(0,3);
+
+	bg.hVertices=[];
+	bg.MetricVertices=zeros(0,3);
+	bg.h1h2VpVertices=[];
+
+	bg.NumSubDomains=0;
+	bg.SubDomains=zeros(0,4);
+
+	bg=class(bg,'bamggeom');
+
+case 1
+
+	bg=bamggeom;
+	object=varargin{1};
+	fields=fieldnames(object);
+	for i=1:length(fields)
+		field=fields{i};
+		if isfield(struct(bg),field),
+			bg.(field)=object.(field);
+		end
+	end
+
+otherwise
+	error('bamggeom constructor error message: unknown type of constructor call');
+end
Index: /issm/trunk/src/m/classes/@bamggeom/display.m
===================================================================
--- /issm/trunk/src/m/classes/@bamggeom/display.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamggeom/display.m	(revision 3277)
@@ -0,0 +1,7 @@
+function display(object)
+%DISPLAY - displays the fields of an bamggeom 
+%
+%   echo function for 'bamggeom' class
+
+disp(sprintf('\n%s = \n',inputname(1)));
+disp(struct(object))
Index: /issm/trunk/src/m/classes/@bamggeom/subsasgn.m
===================================================================
--- /issm/trunk/src/m/classes/@bamggeom/subsasgn.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamggeom/subsasgn.m	(revision 3277)
@@ -0,0 +1,12 @@
+function bg = subsasgn(bg,index,val)
+%SUBSASGN - handles indexed assignments to objects
+%
+%   the first argument is the object, the second argument a structure array
+%   the third one is the value
+%
+%   Usage:
+%      bg = subsasgn(bg,index,val)
+%
+%   See also SUBSREF
+
+bg=builtin('subsasgn',bg,index,val);
Index: /issm/trunk/src/m/classes/@bamggeom/subsref.m
===================================================================
--- /issm/trunk/src/m/classes/@bamggeom/subsref.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamggeom/subsref.m	(revision 3277)
@@ -0,0 +1,10 @@
+function bg = subsref(bg,index)
+%SUBSREF - handles indexed references to objects
+%
+%   the first argument is the object and the second one a structure array
+%   Usage:
+%      bg = subsref(bg,index)
+% 
+%   See also SUBSASGN
+
+bg=builtin('subsref',bg,index);
Index: /issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m
===================================================================
--- /issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m	(revision 3277)
@@ -0,0 +1,62 @@
+function bm = bamgmesh(varargin)
+%ICEFRONT - constructor for bamgmesh object
+%
+%   Usage:
+%      bamgmesh = bamgmesh(varargin)
+
+switch nargin
+case 0
+	% if no input arguments, create a default object
+	bm.NumVertices=0;
+	bm.Vertices=[];
+
+	bm.NumEdges=0;
+	bm.Edges=zeros(0,3);
+
+	bm.NumTriangles=0;
+	bm.Triangles=[];
+
+	bm.NumQuadrilaterals=0;
+	bm.Quadrilaterals=[];
+
+	bm.NumSegments=0;
+	bm.Segments=zeros(0,3);
+	bm.SegmentsMarkers=[];
+
+	bm.NumVerticesOnGeometricVertex=0;
+	bm.VerticesOnGeometricVertex=zeros(0,2);
+
+	bm.NumVerticesOnGeometricEdge=0;
+	bm.VerticesOnGeometricEdge=zeros(0,2);
+
+	bm.NumEdgesOnGeometricEdge=0;
+	bm.EdgesOnGeometricEdge=zeros(0,2);
+
+	bm.NumCrackedEdges=0;
+	bm.CrackedEdges=zeros(0,2);
+
+	bm.NumSubDomains=0;
+	bm.SubDomains=zeros(0,4);
+
+	bm.NumSubDomainsFromGeom=0;
+	bm.SubDomainsFromGeom=zeros(0,4);
+
+	bm.hVertices=[];
+
+	bm=class(bm,'bamgmesh');
+
+case 1
+
+	bm=bamgmesh;
+	object=varargin{1};
+	fields=fieldnames(object);
+	for i=1:length(fields)
+		field=fields{i};
+		if isfield(struct(bm),field),
+			bm.(field)=object.(field);
+		end
+	end
+
+otherwise
+	error('bamgmesh constructor error message: unknown type of constructor call');
+end
Index: /issm/trunk/src/m/classes/@bamgmesh/display.m
===================================================================
--- /issm/trunk/src/m/classes/@bamgmesh/display.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamgmesh/display.m	(revision 3277)
@@ -0,0 +1,7 @@
+function display(object)
+%DISPLAY - displays the fields of an bamgmesh 
+%
+%   echo function for 'bamgmesh' class
+
+disp(sprintf('\n%s = \n',inputname(1)));
+disp(struct(object))
Index: /issm/trunk/src/m/classes/@bamgmesh/subsasgn.m
===================================================================
--- /issm/trunk/src/m/classes/@bamgmesh/subsasgn.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamgmesh/subsasgn.m	(revision 3277)
@@ -0,0 +1,12 @@
+function icefront = subsasgn(icefront,index,val)
+%SUBSASGN - handles indexed assignments to objects
+%
+%   the first argument is the object, the second argument a structure array
+%   the third one is the value
+%
+%   Usage:
+%      icefront = subsasgn(icefront,index,val)
+%
+%   See also SUBSREF
+
+icefront=builtin('subsasgn',icefront,index,val);
Index: /issm/trunk/src/m/classes/@bamgmesh/subsref.m
===================================================================
--- /issm/trunk/src/m/classes/@bamgmesh/subsref.m	(revision 3277)
+++ /issm/trunk/src/m/classes/@bamgmesh/subsref.m	(revision 3277)
@@ -0,0 +1,10 @@
+function icefront = subsref(icefront,index)
+%SUBSREF - handles indexed references to objects
+%
+%   the first argument is the object and the second one a structure array
+%   Usage:
+%      icefront = subsref(icefront,index)
+% 
+%   See also SUBSASGN
+
+icefront=builtin('subsref',icefront,index);
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 3276)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 3277)
@@ -31,7 +31,4 @@
 %   - omega       : relaxation parameter of the smoothing procedure (default is 1.8)
 %   - power       : power applied to the metric (default is 1)
-%   - renumber    : 0 -> no renumbering
-%                   1 -> renumber by distance to lower left corner (default)
-%                   2 -> random renumbering
 %   - splitcorner : split triangles whuch have 3 vertices on the outline (default is 1)
 %   - verbose     : level of verbosity (default is 1)
@@ -48,17 +45,8 @@
 %initialize the structures required as input of Bamg
 bamg_options=struct();
-bamg_geometry=struct();
-bamg_mesh=struct();
+bamg_geometry=bamggeom;
+bamg_mesh=bamgmesh;
 
 % Bamg Geometry parameters {{{1
-bamg_geometry.NumVertices=0;
-bamg_geometry.Vertices=zeros(0,3);
-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;
-bamg_geometry.SubDomains=zeros(0,3);
 if exist(options,'domain'),
 
@@ -187,10 +175,11 @@
 	end
 
+elseif isstruct(md.bamg),
+
+	bamg_geometry=bamggeom(md.bamg.geometry); %TEST
+
 else
 
-	%Use segments as the geometry
-	%bamg_geometry.Vertices=[md.x md.y ones(md.numberofgrids,1)];
-	%bamg_geometry.Edges=[md.segments(:,1:2) md.segmentmarkers];
-	%bamg_geometry.Edges
+	%do nothing...
 
 end
@@ -203,22 +192,15 @@
 
 % Bamg Mesh parameters {{{1
-bamg_mesh.NumVertices=0;
-bamg_mesh.Vertices=[];
-bamg_mesh.NumTriangles=0;
-bamg_mesh.Triangles=[];
-bamg_mesh.NumEdges=0;
-bamg_mesh.Edges=zeros(0,3);
-bamg_mesh.NumEdgesOnGeometricEdge=0;
-bamg_mesh.EdgesOnGeometricEdge=zeros(0,2);
-bamg_mesh.NumVerticesOnGeometricEdge=0;
-bamg_mesh.VerticesOnGeometricEdge=zeros(0,2);
-bamg_mesh.hVertices=[];
-bamg_mesh.NumCrackedEdges=0;
-bamg_mesh.CrackedEdges=zeros(0,2);
 if (~exist(options,'domain') & md.numberofgrids~=0 & strcmpi(md.type,'2d')),
-	bamg_mesh.NumVertices=md.numberofgrids;
-	bamg_mesh.Vertices=[md.x md.y ones(md.numberofgrids,1)];
-	bamg_mesh.NumTriangles=md.numberofelements;
-	bamg_mesh.Triangles=[md.elements ones(md.numberofelements,1)];
+
+	if isstruct(md.bamg),%TEST
+		bamg_mesh=bamgmesh(md.bamg.mesh);
+	else
+		bamg_mesh.NumVertices=md.numberofgrids;
+		bamg_mesh.Vertices=[md.x md.y ones(md.numberofgrids,1)];
+		bamg_mesh.NumTriangles=md.numberofelements;
+		bamg_mesh.Triangles=[md.elements ones(md.numberofelements,1)];
+	end
+
 	if exist(options,'hVertices'),
 		bamg_mesh.hVertices=getfieldvalueerr(options,'hVertices');
@@ -255,5 +237,4 @@
 bamg_options.omega=getfieldvalue(options,'omega',1.8);
 bamg_options.power=getfieldvalue(options,'power',1);
-bamg_options.renumber=getfieldvalue(options,'renumber',1);
 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
 bamg_options.verbose=getfieldvalue(options,'verbose',1);
@@ -271,5 +252,4 @@
 end
 
-
 %call Bamg
 [bamgmesh_out bamggeom_out]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
@@ -277,6 +257,6 @@
 % plug results onto model
 md.bamg=struct();
-md.bamg.mesh=bamgmesh_out;
-md.bamg.geometry=bamggeom_out;
+md.bamg.mesh=bamgmesh(bamgmesh_out);
+md.bamg.geometry=bamggeom(bamggeom_out);
 md.x=bamgmesh_out.Vertices(:,1);
 md.y=bamgmesh_out.Vertices(:,2);
Index: /issm/trunk/src/mex/Bamg/Bamg.cpp
===================================================================
--- /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3276)
+++ /issm/trunk/src/mex/Bamg/Bamg.cpp	(revision 3277)
@@ -52,4 +52,6 @@
 	FetchData(&bamgmesh_in.NumVerticesOnGeometricEdge,mxGetField(BAMGMESH,0,"NumVerticesOnGeometricEdge"));
 	FetchData(&bamgmesh_in.VerticesOnGeometricEdge,NULL,NULL,mxGetField(BAMGMESH,0,"VerticesOnGeometricEdge"));
+	FetchData(&bamgmesh_in.NumVerticesOnGeometricVertex,mxGetField(BAMGMESH,0,"NumVerticesOnGeometricVertex"));
+	FetchData(&bamgmesh_in.VerticesOnGeometricVertex,NULL,NULL,mxGetField(BAMGMESH,0,"VerticesOnGeometricVertex"));
 	if (bamgmesh_in.hVertices && (cols!=1 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",bamgmesh_in.NumVertices,1));}
 
@@ -61,5 +63,4 @@
 	FetchData(&bamgopts.Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype"));
 	FetchData(&bamgopts.KeepVertices,mxGetField(BAMGOPTIONS,0,"KeepVertices"));
-	FetchData(&bamgopts.renumber,mxGetField(BAMGOPTIONS,0,"renumber"));
 	FetchData(&bamgopts.power,mxGetField(BAMGOPTIONS,0,"power"));
 	FetchData(&bamgopts.errg,mxGetField(BAMGOPTIONS,0,"errg"));