Changeset 3277

Timestamp:

03/15/10 10:33:22 (15 years ago)

Author:

Mathieu Morlighem

Message:

Keep Bamg output in model for rifts

Location:

issm/trunk/src

Files:

• c/Bamgx/meshtype.h (modified) (1 diff)
• c/Bamgx/objects/Geometry.cpp (modified) (2 diffs)
• c/Bamgx/objects/Triangles.cpp (modified) (13 diffs)
• c/Bamgx/objects/Triangles.h (modified) (2 diffs)
• c/Bamgx/objects/VertexOnGeom.h (modified) (1 diff)
• c/objects/BamgMesh.cpp (modified) (1 diff)
• c/objects/BamgOpts.cpp (modified) (2 diffs)
• c/objects/BamgOpts.h (modified) (1 diff)
• m/classes/@bamggeom (added)
• m/classes/@bamggeom/bamggeom.m (added)
• m/classes/@bamggeom/display.m (added)
• m/classes/@bamggeom/subsasgn.m (added)
• m/classes/@bamggeom/subsref.m (added)
• m/classes/@bamgmesh (added)
• m/classes/@bamgmesh/bamgmesh.m (added)
• m/classes/@bamgmesh/display.m (added)
• m/classes/@bamgmesh/subsasgn.m (added)
• m/classes/@bamgmesh/subsref.m (added)
• m/classes/public/bamg.m (modified) (7 diffs)
• mex/Bamg/Bamg.cpp (modified) (2 diffs)

issm/trunk/src/c/Bamgx/meshtype.h

@@ -118 +118 @@
                 }
         }
+        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

issm/trunk/src/c/Bamgx/objects/Geometry.cpp

@@ -929 +929 @@
         /*}}}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 +969 @@
                                         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");

issm/trunk/src/c/Bamgx/objects/Triangles.cpp

@@ -401 +401 @@
                 }
 
-                //VerVerticesOnGeometricEdge
+                //VerticesOnGeometricEdge
                 if(bamgmesh->VerticesOnGeometricEdge){
                         if(verbose>5) printf("      processing VerticesOnGeometricEdge\n");
@@ -419 +419 @@
                 }
 
+                //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 +457 @@
                                 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 +473 @@
                                 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 +491 @@
                                         }
                                         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 +526 @@
                         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 +569 @@
                 /*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 +775 @@
                         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 +4316 @@
         }                  
         /*}}}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 +4390 @@
                 /*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 +4412 @@
                   }
 
+                printf("renumbering vertices on edge\n");
                 for (i=0;i< NbVerticesOnGeomEdge;i++)
                   {
@@ -4492 +4420 @@
                   }
 
+                printf("renumbering vertices on Bth vertex\n");
                 for (i=0;i< NbVertexOnBThVertex;i++)
                   {
@@ -4514 +4443 @@
                          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*/

issm/trunk/src/c/Bamgx/objects/Triangles.h

@@ -101 +101 @@
                         void Insert();
                         void ForceBoundary();
-                        void Renumber(BamgOpts* bamgopts);
                         void FindSubDomain(int OutSide=0);
                         long TriangleReferenceList(long *) const;
@@ -159 +158 @@
                           }
                         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();

issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h

@@ -34 +34 @@
 
                         //Operators
-                        operator Vertex * () const  {return mv;}
+                        operator Vertex*() const  {return mv;}
                         operator GeometricalVertex * () const  {return gv;}
                         operator GeometricalEdge * () const  {return ge;}

issm/trunk/src/c/objects/BamgMesh.cpp

@@ -146 +146 @@
 
         pfield=mxCreateDoubleMatrix(0,0,mxREAL);
-        mxSetM(pfield,2);
+        mxSetM(pfield,3);
         mxSetN(pfield,bamgmesh->NumVerticesOnGeometricEdge);
         mxSetPr(pfield,bamgmesh->VerticesOnGeometricEdge);

issm/trunk/src/c/objects/BamgOpts.cpp

@@ -11 +11 @@
         bamgopts->Metrictype=0;
         bamgopts->KeepVertices=0;
-        bamgopts->renumber=0;
         bamgopts->maxsubdiv=0;
         bamgopts->power=0;
@@ -42 +41 @@
         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"));

issm/trunk/src/c/objects/BamgOpts.h

@@ -14 +14 @@
         int     Metrictype;
         int     KeepVertices;
-        int     renumber;
         double  maxsubdiv;
         double  power;

issm/trunk/src/m/classes/public/bamg.m

@@ -31 +31 @@
 %   - 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 +45 @@
 %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 +175 @@
         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 +192 @@
 
 % 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 +237 @@
 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 +252 @@
 end
 
-
 %call Bamg
 [bamgmesh_out bamggeom_out]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
@@ -277 +257 @@
 % 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);

issm/trunk/src/mex/Bamg/Bamg.cpp

@@ -52 +52 @@
         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 +63 @@
         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"));