Changeset 2899

Timestamp:

01/22/10 17:14:29 (15 years ago)

Author:

Mathieu Morlighem

Message:

some improvements in Bamg

Location:

issm/trunk/src

Files:

• c/Bamgx/Bamgx.cpp (modified) (5 diffs)
• c/Bamgx/Mesh2.h (modified) (1 diff)
• c/Bamgx/Metric.h (modified) (3 diffs)
• c/Bamgx/objects/Geometry.cpp (modified) (3 diffs)
• c/Bamgx/objects/MatVVP2x2.cpp (modified) (1 diff)
• c/Bamgx/objects/MetricAnIso.cpp (modified) (1 diff)
• c/Bamgx/objects/Triangles.cpp (modified) (16 diffs)
• c/objects/BamgOpts.h (modified) (1 diff)
• m/classes/public/bamg.m (modified) (6 diffs)
• mex/Bamg/Bamg.cpp (modified) (6 diffs)
• mex/Bamg/Bamg.h (modified) (1 diff)

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

@@ -52 +52 @@
         double omega=1.8;
         int NbSmooth=3;
+        int ChoiseHessien = 0;
+        double power=1;
+        int Rescaling=1;
+
         Triangles* Thr=NULL;
         Triangles* Thb=NULL;
@@ -79 +83 @@
                 if (verbosity>1) printf("   Processing geometry...\n");
                 Geometry Gh(bamggeom,bamgopts);
+                if (verbosity>10){
+                        Gh.Echo();
+                }
+
+                //get hmin and hmax from geometry to generate the metric
                 hmin = Max(hmin,Gh.MinimalHmin());
                 hmax = Min(hmax,Gh.MaximalHmax());
@@ -118 +127 @@
                 hmax = Min(hmax,BTh.MaximalHmax());
 
-                //?????TEST
+                //Make Quadtree from background mesh
                 BTh.MakeQuadTree();
 
                 //build metric if not given in input
-                if (verbosity>1) printf("   Processing Metric...\n");
                 if (bamgopts->metric){
+                        if (verbosity>1) printf("   Processing Metric...\n");
                         BTh.ReadMetric(bamgopts,hmin,hmax,coef);
+                }
+                else if (bamgopts->field){
+                        if (verbosity>1) printf("   Generating Metric from solution field...\n");
+                        BTh.IntersectConsMetric(bamgopts->field,1,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:cutoff,nbjacoby,Rescaling,power,ChoiseHessien);
                 }
                 else { // init with hmax 
                         Metric Mhmax(hmax);
-                        for (Int4 iv=0;iv<BTh.nbv;iv++)
-                         BTh[iv].m = Mhmax;
+                        for (Int4 iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;
                 }
 
@@ -138 +150 @@
 
                 //Build new mesh Thr
+                if (verbosity>1) printf("   Generating Mesh...\n");
                 Thr=&BTh,Thb=0;
                 Triangles & Th( *(0 ?  new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) :  new Triangles(maxnbv,BTh,KeepBackVertices)));
@@ -163 +176 @@
                 if (verbosity>1) printf("   Write Geometry...\n");
                 Th.Gh.WriteGeometry(bamggeom,bamgopts);
+                if (verbosity>1) printf("   Write Metric...\n");
+                BTh.WriteMetric(bamgopts);
 
                 /*clean up*/

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

@@ -845 +845 @@
                         Real8 MaximalAngleOfCorner;
 
-                        //  end of data 
-
+                        void Echo();
 
                         I2 toI2(const R2 & P) const {

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

@@ -73 +73 @@
                 inline void  Box(Real4 &hx,Real4 &hy) const ;  
                 MetricAnIso(const MatVVP2x2);
+                void  Echo();
         };
 
@@ -82 +83 @@
                 D2 v;
 
+                void  Echo();
                 MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
 
@@ -103 +105 @@
                 void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
                 void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
-          };
+        };
 
         inline void  MatVVP2x2::BoundAniso2(const Real8 coef){

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

@@ -580 +580 @@
                         }
                 }
-
-                // bulle sort on the angle of edge  
+                // in our case, for 4 points and 4 edges (in parenthesis from now on):
+                // hv =  7 2 4 6 
+                // ev =  -1 -1 1 -1 3 -1 5 0 
+
+                //compute angle of edges
                 for (i=0;i<nbv;i++) {
+
                         int exch=1, ord=0;      
                         while (exch){
-                                exch  =0;
-                                Int4* p=hv+i;
-                                Int4* po=p;
-                                Int4  n=*p;
+                                Int4 *p=hv+i;                   // pointer to hv[vertex i]  (p=hv)
+                                Int4 *po=p;                     // copy of pointer p        (po=p)
+                                Int4  n=*p;                     // value hv[vertex i]       (n=7)
                                 register float angleold=-1000 ; // angle = - infinity
-                                ord = 0;
-                                while (n >=0){
+                                ord=0; exch=0;
+                                while (n >=0){// loop as long as ...
                                         ord++;
-                                        register Int4 i1=n/2;
-                                        register Int4 j1=n%2;
-                                        register Int4 *pn=ev+n;
+                                        register Int4  i1=n/2;       // i1 = floor (n/2)          (i1 = 7/2 = 3)
+                                        register Int4  j1=n%2;       // j1 = 1 if n is odd        (j1 = 1)
+                                        register Int4* pn=ev+n;      // pointer to ev[n]          (pn = &ev[7])
+                                        //compute angle from eangle:
+                                        // if n odd: angle = eangle - Pi [2*Pi]
+                                        // else    : angle = eangle
                                         float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1];
-                                        n = *pn;
-                                        if (angleold > angle) // exch to have : po -> pn -> p 
-                                         exch=1,*pn = *po,*po=*p,*p=n,po = pn;
-                                        else //  to have : po -> p -> pn 
-                                         angleold =  angle, po = p,p  = pn;
-                                }
-                        }
-
-                        // angulare test to find a corner 
+                                        n=*pn;                       //  n = ev[n]                (n = ev[7] = 0)
+                                        if (angleold > angle){       // exch to have : po -> pn -> p 
+                                                exch=1, *pn=*po, *po=*p, *p=n, po=pn;
+                                        }
+                                        else{                        //  to have : po -> p -> pn
+                                                angleold=angle, po=p, p=pn;
+                                        }
+                                }
+                        }
+                        //first iteration: i=0
+                        //hv 0 2 4 6 
+                        //ev 7 -1 1 -1 3 -1 5 -1 
+
+                        // angular test on current vertex to guess whether it is a corner
                         if(ord == 2) { 
-                                Int4 n1 = hv[i];
-                                Int4 n2 = ev[n1];
-                                Int4 i1 = n1 /2, i2 = n2/2; // edge number
-                                Int4  j1 = n1 %2, j2 = n2%2; // vertex in the edge 
+                                Int4  n1 = hv[i];
+                                Int4  n2 = ev[n1];
+                                Int4  i1 = n1/2, i2 = n2/2; // edge number
+                                Int4  j1 = n1%2, j2 = n2%2; // vertex in the edge 
+                                /*//first iteration i=0
+                                //n1=hv[i]=0
+                                //n2=ev[i]=7
+                                //i1 = 0  i2 = 3
+                                //j1 = 0  j2 = 1
+                                printf("i=%i\n",i);
+                                printf("n1=hv[i]=%i\n",n1);
+                                printf("n2=ev[i]=%i\n",n2);
+                                printf("i1 = %i  i2 = %i\n",i1,i2);
+                                printf("j1 = %i  j2 = %i\n",j1,j2);*/
                                 float angle1=  j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
                                 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
@@ -638 +659 @@
                 k =0;
                 for (i=0;i<nbe;i++){
-                        for (jj=0;jj<2;jj++){
+                        for (j=0;j<2;j++){
                                 Int4 n1 = ev[k++]; 
                                 Int4 i1 = n1/2 ,j1=n1%2;
-                                if( edges[i1].v[j1] != edges[i].v[jj]) {
+                                if( edges[i1].v[j1] != edges[i].v[j]) {
                                         throw ErrorException(__FUNCT__,exprintf("Bug Adj edge"));
                                 }
                                 edges[i1].Adj[j1] = edges + i;
-                                edges[i1].DirAdj[j1] = jj;
+                                edges[i1].DirAdj[j1] = j;
                         }
                 }
@@ -776 +797 @@
         }
         /*}}}1*/
+        /*FUNCTION  Geometry::Echo {{{1*/
+
+        void Geometry::Echo(void){
+
+                printf("Geometry:\n");
+                printf("   name: %s\n",name);
+                printf("   nbvx (maximum number of vertices) : %i\n",nbvx);
+                printf("   nbtx (maximum number of triangles): %i\n",nbtx);
+                printf("   nbv  (number of vertices) : %i\n",nbv);
+                printf("   nbt  (number of triangles): %i\n",nbt);
+                printf("   nbe  (number of edges)    : %i\n",nbe);
+                printf("   nbv  (number of initial vertices) : %i\n",nbiv);
+                printf("   NbSubDomains: %i\n",NbSubDomains);
+                printf("   NbEquiEdges: %i\n",NbEquiEdges);
+                printf("   NbCrackedEdges: %i\n",NbCrackedEdges);
+                printf("   NbOfCurves: %i\n",NbOfCurves);
+                printf("   vertices: %p\n",vertices);
+                printf("   triangles: %p\n",triangles);
+                printf("   edges: %p\n",edges);
+                printf("   quadtree: %p\n",quadtree);
+                printf("   subdomains: %p\n",subdomains);
+                printf("   curves: %p\n",curves);
+                printf("   pmin (x,y): (%g %g)\n",pmin.x,pmin.y);
+                printf("   pmax (x,y): (%g %g)\n",pmax.x,pmax.y);
+                printf("   coefIcoor: %g\n",coefIcoor);
+                printf("   MaximalAngleOfCorner: %g\n",MaximalAngleOfCorner);
+
+                return;
+        }
+        /*}}}*/
         /*FUNCTION  Geometry::EmptyGeometry(){{{1*/
         void Geometry::EmptyGeometry() {

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

@@ -47 +47 @@
         /*}}}1*/
 
+        /*Methods*/
+        /*FUNCTION MatVVP2x2::Echo {{{1*/
+
+        void MatVVP2x2::Echo(void){
+
+                printf("MatVVP2x2:\n");
+                printf("   lambda1: %g\n",lambda1);
+                printf("   lambda2: %g\n",lambda2);
+                printf("   v.x: %g\n",v.x);
+                printf("   v.y: %g\n",v.y);
+
+                return;
+        }
+        /*}}}*/
+
 } 

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

@@ -59 +59 @@
 
         /*Methods*/
+        /*FUNCTION MetricAnIso::Echo {{{1*/
+
+        void MetricAnIso::Echo(void){
+
+                printf("MetricAnIso:\n");
+                printf("   [a11 a21 a22]: [%g %g %g]\n",a11,a21,a22);
+
+                return;
+        }
+        /*}}}*/
         /*FUNCTION MetricAnIso::IntersectWith{{{1*/
         int MetricAnIso::IntersectWith(const MetricAnIso M2) {

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

@@ -325 +325 @@
                         hvertices=1;
                         for (i=0;i< nbv;i++){
-                                vertices[i].m=Metric((Real4)bamgmesh->hVertices[i]);
+                                if (!isnan(bamgmesh->hVertices[i])){
+                                        vertices[i].m=Metric((Real4)bamgmesh->hVertices[i]);
+                                }
                         }
                 }
@@ -689 +691 @@
                         }
                         else if (j==3){
-                                Real8 a,b,c;         
-                                a=bamgopts->metric[i*3+0];
-                                b=bamgopts->metric[i*3+1];
-                                c=bamgopts->metric[i*3+2];
-                                MetricAnIso M(a,b,c);
-                                MatVVP2x2 Vp(M/coef);
-
-                                Vp.Maxh(hmin);
-                                Vp.Minh(hmax);
-                                vertices[i].m = Vp;
+                                //do not erase metric computed by hVertices
+                                if (vertices[i].m.a11==1 && vertices[i].m.a21==0 && vertices[i].m.a22==1){
+                                        Real8 a,b,c;         
+                                        a=bamgopts->metric[i*3+0];
+                                        b=bamgopts->metric[i*3+1];
+                                        c=bamgopts->metric[i*3+2];
+                                        MetricAnIso M(a,b,c);
+                                        MatVVP2x2 Vp(M/coef);
+
+                                        Vp.Maxh(hmin);
+                                        Vp.Minh(hmax);
+                                        vertices[i].m = Vp;
+                                }
                         }
                 }
@@ -2723 +2728 @@
         /*}}}1*/
 /*FUNCTION Triangles::IntersectConsMetric{{{1*/
-void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols,
+void Triangles::IntersectConsMetric(const double* s,const Int4 nbsol,const int * typsols,
                         const  Real8 hmin1,const Real8 hmax1,const Real8 coef,
                         const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,
-                        const int DoNormalisation,const double power,const int choice)
-{ //  the array of solution s is store    
+                        const int DoNormalisation,const double power,const int choice){
+        //  the array of solution s is store    
         // sol0,sol1,...,soln    on vertex 0
         //  sol0,sol1,...,soln   on vertex 1
         //  etc.
-        //  choise = 0 =>  H is computed with green formule
-        //   otherwise  => H is computed from P2 on 4T 
+        //  choise = 0 =>  H is computed with green formula
+        //  otherwise  =>  H is computed from P2 on 4T 
+
         const int dim = 2;
-
         long int verbosity=0;
 
@@ -2741 +2746 @@
         // computation of the nb of field 
         Int4 ntmp = 0;
-        if (typsols)
-          {
-                for (Int4 i=0;i<nbsol;i++)
-                 ntmp += sizeoftype[typsols[i]];
+        if (typsols){
+                for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]];
           }
-        else
-         ntmp = nbsol;
+        else ntmp = nbsol;
 
         // n is the total number of fields
@@ -2768 +2770 @@
                 else printf("      Absolute error\n");
         }
-        double *ss=(double*)s;
+        double* ss=(double*)s;
 
         double sA,sB,sC;
@@ -2781 +2783 @@
         Real8 *workV= new Real8[nbv];
         int *OnBoundary = new int[nbv];
-        for (iv=0;iv<nbv;iv++)
-          {
+        for (iv=0;iv<nbv;iv++){
                 Mmass[iv]=0;
                 OnBoundary[iv]=0;
                 Mmassxx[iv]=0;
-          }
+        }
 
         for (i=0;i<nbt;i++) 
-         if(triangles[i].link) // the real triangles 
-                {
+         if(triangles[i].link){ // the real triangles 
                  const Triangle &t=triangles[i];
                  // coor of 3 vertices 
@@ -2823 +2823 @@
                  Mmass[iC] += dett;
 
-                 if((nbb==0)|| !choice)
-                        {
+                 if((nbb==0)|| !choice){
                          Mmassxx[iA] += dett;
                          Mmassxx[iB] += dett;
@@ -2830 +2829 @@
                         }
                 }
-         else
-          workT[i]=-1;
-
-        for (Int4 nusol=0;nusol<nbsol;nusol++)
-          { //for all Solution  
+         else workT[i]=-1;
+
+        for (Int4 nusol=0;nusol<nbsol;nusol++) {
+                //for all Solution  
 
                 Real8 smin=ss[0],smax=ss[0];
@@ -2843 +2841 @@
                 int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 
                 if (nbfield == 1) 
-                 for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-                        {
+                 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){
                          dxdx[iv]=dxdy[iv]=dydy[iv]=0;
                          smin=Min(smin,ss[k]);
                          smax=Max(smax,ss[k]);
-                        }
-                else
-                  {
+                 }
+                else{
                         //  cas vectoriel 
-                        for ( iv=0,k=0; iv<nbv; iv++,k+=n )
-                          {     
+                        for ( iv=0,k=0; iv<nbv; iv++,k+=n ){    
                                 double v=0;                  
                                 for (int i=0;i<nbfield;i++) 
@@ -2860 +2855 @@
                                 smin=Min(smin,v);
                                 smax=Max(smax,v);
-                          }
-                  }
+                        }
+                }
                 Real8 sdelta = smax-smin;
                 Real8 absmax=Max(Abs(smin),Abs(smax));
@@ -2874 +2869 @@
 
                 double *sf  = ss; 
-                for (Int4 nufield=0;nufield<nbfield;nufield++,ss++) 
-                  {
+                for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){
                         for ( iv=0,k=0; iv<nbv; iv++,k+=n )
                          dxdx[iv]=dxdy[iv]=dydy[iv]=0;
@@ -3035 +3029 @@
                          else kk++;
 
-
-                        // correction of second derivate
+                        // correction of second derivative
                         // by a laplacien
-
                         Real8 *d2[3] = { dxdx, dxdy, dydy};
                         Real8 *dd;
-                        for (int xy = 0;xy<3;xy++)
-                          {
+                        for (int xy = 0;xy<3;xy++) {
                                 dd = d2[xy];
                                 // do leat 2 iteration for boundary problem
-                                for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++)
-                                  {
+                                for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
                                         for (i=0;i<nbt;i++) 
                                          if(triangles[i].link) // the real triangles 
@@ -3078 +3068 @@
                                          if( ijacobi<NbJacobi || OnBoundary[iv])
                                           dd[iv] = workV[iv]/(Mmass[iv]*6);
-
-
-                                  }
-
-
-                          }
+                                }
+                        }
 
                         // constuction  of the metrix from the Hessian dxdx. dxdy,dydy
-
                         Real8 rCutOff=CutOff*absmax;// relative cut off 
 
@@ -3107 +3092 @@
 
                                 Vp.Abs();
-                                if(power!=1.0) 
-                                 Vp.pow(power);
-
-
+                                if(power!=1.0) Vp.pow(power);
 
                                 h1=Min(h1,Vp.lmin());
@@ -3128 +3110 @@
                                 Metric MVp(Vp);
                                 vertices[iv].m.IntersectWith(MVp);
+                                //vertices[iv].m=Vp; TEST: work like that
                           }// for all vertices 
-                        if (verbosity>2) { 
+                        //vertices[nbv-1].m.Echo();
+                        if (verbosity>-1) { 
                                 printf("      Field %i of solution %i\n",nufield,nusol);
                                 printf("         before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5));
@@ -3135 +3119 @@
                         }
                   } //  end of for all field
-          }// end for all solution 
+        }// end for all solution 
 
         delete [] detT;

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

@@ -18 +18 @@
         int     verbose;
         double* metric;
+        double* field;
 
 };

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

@@ -35 +35 @@
                         bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
                         bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  ones(nods,1)]];
-                        bamg_geometry.hVertices=[];%[bamg_geometry.hVertices; resolution*ones(nods,1)];
                         if i>1,
                                 clockwise=-1;
@@ -41 +40 @@
                         end
                         count=count+nods;
+                end
+                if exist(options,'hVertices'),
+                        bamg_geometry.hVertices=getfieldvalueerr(options,'hVertices');
+                        if length(bamg_geometry.hVertices)~=nods,
+                                error(['hVertices option should be of size [' num2str(nods) ',1]']);
+                        end
                 end
                 bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
@@ -55 +60 @@
 bamg_mesh.Triangles=[];
 bamg_mesh.index=zeros(0,3);
+bamg_mesh.hVertices=[];
 if (~exist(options,'domain') & md.numberofgrids~=0 & strcmpi(md.type,'2d')),
         bamg_mesh.NumVertices=md.numberofgrids;
@@ -60 +66 @@
         bamg_mesh.NumTriangles=md.numberofelements;
         bamg_mesh.Triangles=[md.elements ones(md.numberofelements,1)];
+        if exist(options,'hVertices'),
+                bamg_mesh.hVertices=getfieldvalueerr(options,'hVertices');
+                if length(bamg_mesh.hVertices)~=md.numberofgrids,
+                        error(['hVertices option should be of size [' num2str(md.numberofgrids) ',1]']);
+                end
+        end
 end
 %}}}
@@ -74 +86 @@
 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
 bamg_options.metric=getfieldvalue(options,'metric',zeros(0,3));
+bamg_options.field=getfieldvalue(options,'field',zeros(0,1));
 %}}}
 
 %call Bamg
-[triangles vertices]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
+[triangles vertices metric]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
 
 % plug results onto model
@@ -83 +96 @@
 md.y=vertices(:,2);
 md.elements=triangles(:,1:3);
+md.dummy=metric;
 
 %Fill in rest of fields:

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

@@ -9 +9 @@
         int   noerr=1;
         int   i;
+        int   nods=0; //to be removed
         BamgOpts bamgopts;
         BamgMesh bamgmesh;
@@ -18 +19 @@
         int     NumTrianglesMesh;
         double* TrianglesMesh=NULL;
+        double* hVerticesMesh=NULL;
 
         /*Geom inputs: */
@@ -34 +36 @@
         double gradation;
         double cutoff;
-        double* metric;
+        double* metric=NULL;
+        double* field=NULL;
 
         /*Boot module: */
@@ -79 +82 @@
         FetchData(&TrianglesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles"));
         bamgmesh.Triangles=TrianglesMesh;
-        bamgmesh.hVertices=NULL;
+        FetchData(&hVerticesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"hVertices"));
+        bamgmesh.hVertices=hVerticesMesh;
         bamgmesh.NumQuadrilaterals=0;
         bamgmesh.Quadrilaterals=NULL;
@@ -118 +122 @@
         FetchData(&metric,NULL,NULL,mxGetField(BAMGOPTIONS,0,"metric"));
         bamgopts.metric=metric;
+        FetchData(&field,NULL,NULL,mxGetField(BAMGOPTIONS,0,"field"));
+        bamgopts.field=field;
 
         /*!Generate internal degree of freedom numbers: */
+        nods=bamgmesh.NumVertices;
         Bamgx(&bamgmesh,&bamggeom,&bamgopts);
 
@@ -125 +132 @@
         WriteData(TRIANGLESOUT,bamgmesh.Triangles,bamgmesh.NumTriangles,4);
         WriteData(VERTICESOUT,bamgmesh.Vertices,bamgmesh.NumVertices,3);
+        WriteData(METRICOUT,bamgopts.metric,nods,3);
 
         /*Free ressources: */

issm/trunk/src/mex/Bamg/Bamg.h

@@ -25 +25 @@
 #define TRIANGLESOUT (mxArray**)&plhs[0]
 #define VERTICESOUT  (mxArray**)&plhs[1]
+#define METRICOUT    (mxArray**)&plhs[2]
 
 /* serial arg counts: */
 #undef NLHS
-#define NLHS  2
+#define NLHS  3
 #undef NRHS
 #define NRHS  3