Changeset 2941

Timestamp:

02/01/10 14:57:22 (15 years ago)

Author:

Mathieu Morlighem

Message:

now different errors must be used to compute metric

Location:

issm/trunk/src

Files:

• c/Bamgx/Bamgx.cpp (modified) (5 diffs)
• c/Bamgx/Mesh2.h (modified) (4 diffs)
• c/Bamgx/objects/Triangles.cpp (modified) (5 diffs)
• c/Bamgx/objects/Vertex.cpp (modified) (4 diffs)
• c/objects/BamgOpts.h (modified) (2 diffs)
• m/classes/public/bamg.m (modified) (2 diffs)
• mex/Bamg/Bamg.cpp (modified) (5 diffs)

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

@@ -33 +33 @@
 
         /*Bamg options*/
-        int    iso;
+        double anisomax;
         int    maxnbv;
         double hmin,hmax;
-        double err,errg,coef;
+        double coef;
         int    verbosity;
         int    NbSmooth;
@@ -46 +46 @@
         int allquad=0;
         double costheta=2;
-        double anisomax = 1e6;
         int KeepBackVertices=1;
         double hminaniso=1e-100; 
@@ -56 +55 @@
 
         /*Bamg options*/
-        iso=bamgopts->iso;
-        err=bamgopts->err;
-        errg=bamgopts->errg;
+        anisomax=bamgopts->anisomax;
         NbSmooth=bamgopts->NbSmooth;
         omega=bamgopts->omega;
@@ -71 +68 @@
                 throw ErrorException(__FUNCT__,exprintf("maxsubdiv %g is should be in ]1,%g]",maxsubdiv,boundmaxsubdiv));
         }
-        // if iso, then anisomax = 1;
-        if (iso) anisomax=1;
         // no metric -> no smoothing 
         if (bamgopts->metric==NULL){
@@ -159 +154 @@
                 }
 
-                BTh.IntersectGeomMetric(errg,iso);
+                BTh.IntersectGeomMetric(bamgopts);
                 if (verbosity>1) printf("   Smoothing Metric...");
                 if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation);

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

@@ -137 +137 @@
                 //  ~Vertex(){}
                 Real8  Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,Real8 =1);
-                void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,BamgOpts* bamgopts);
+                void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
                 void Echo();
                 int ref() const { return ReferenceNumber;}
@@ -802 +802 @@
                                                  void BuildMetric0(BamgOpts* bamgopts);
                                                  void BuildMetric1(BamgOpts* bamgopts);
-                                                 void IntersectGeomMetric(const Real8 err,const int iso);
+                                                 void IntersectGeomMetric(BamgOpts* bamgopts);
 
 
@@ -1078 +1078 @@
           }
 
-        inline  void   Triangles::SetVertexFieldOn()
-          {
+        inline  void   Triangles::SetVertexFieldOn(){
                 for (Int4 i=0;i<nbv;i++) 
                  vertices[i].on=0;
@@ -1087 +1086 @@
                  VerticesOnGeomEdge[k].SetOn();
           }            
-        inline  void   Triangles::SetVertexFieldOnBTh()
-          {
+        inline  void   Triangles::SetVertexFieldOnBTh(){
                 for (Int4 i=0;i<nbv;i++) 
                  vertices[i].on=0;

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

@@ -712 +712 @@
                 int i;
                 xfree((void**)&bamgopts->metric);
-                if (bamgopts->iso){
-                        bamgopts->metric=(double*)xmalloc(1*nbv*sizeof(double));
-                        for (i=0;i<nbv;i++){
-                                MatVVP2x2 V=vertices[i].m;
-                                bamgopts->metric[i]=V.hmin();
-                        }
-                }
-                else {
-                        bamgopts->metric=(double*)xmalloc(3*nbv*sizeof(double));
-                        for (i=0;i<nbv;i++){
-                                bamgopts->metric[i*3+0]=vertices[i].m.a11;
-                                bamgopts->metric[i*3+1]=vertices[i].m.a21;
-                                bamgopts->metric[i*3+2]=vertices[i].m.a22;
-                        }
+                bamgopts->metric=(double*)xmalloc(3*nbv*sizeof(double));
+                for (i=0;i<nbv;i++){
+                        bamgopts->metric[i*3+0]=vertices[i].m.a11;
+                        bamgopts->metric[i*3+1]=vertices[i].m.a21;
+                        bamgopts->metric[i*3+2]=vertices[i].m.a22;
                 }
         }
@@ -1068 +1059 @@
                         /*Compute Metric from Hessian*/
                         for ( iv=0;iv<nbv;iv++){
-                                vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts);
+                                vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
                         }
 
@@ -1361 +1352 @@
                                 /*Compute Metric from Hessian*/
                                 for ( iv=0;iv<nbv;iv++){
-                                        vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts);
+                                        vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
                                 }
 
@@ -3242 +3233 @@
 /*FUNCTION Triangles::IntersectConsMetric{{{1*/
 void Triangles::IntersectConsMetric(BamgOpts* bamgopts){
-        //  the array of solution s is store    
-        // sol0,sol1,...,soln    on vertex 0
-        //  sol0,sol1,...,soln   on vertex 1
-        //  etc.
-        //  Hessiantype = 0 =>  H is computed with green formula
-        //  Hessiantype = 1 =>  H is computed from P2 on 4T 
+        //  Hessiantype = 0 =>  H is computed using double P2 projection
+        //  Hessiantype = 1 =>  H is computed with green formula
 
         /*Options*/
@@ -3264 +3251 @@
 /*}}}1*/
 /*FUNCTION Triangles::IntersectGeomMetric{{{1*/
-void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0){
-        long int verbosity=0;
+void Triangles::IntersectGeomMetric(BamgOpts* bamgopts){
+
+        /*Get options*/
+        int    verbosity=bamgopts->verbose;
+        double anisomax =bamgopts->anisomax;
+        double err      =bamgopts->err[0]; //to be changed...
+
         Real8 ss[2]={0.00001,0.99999};
         Real8 errC = 2*sqrt(2*err);
         Real8 hmax = Gh.MaximalHmax();
         Real8 hmin = Gh.MinimalHmin();
-        Real8 maxaniso = 1e6;
+
+        //check that hmax is positive
         if (hmax<=0){
                 throw ErrorException(__FUNCT__,exprintf("hmax<=0"));
         }
+
+        //errC cannot be higher than 1
+        if (errC > 1) errC = 1;
+
+        //Set all vertices to "on"
         SetVertexFieldOn();
-        if (errC > 1) errC = 1;
-        for (Int4  i=0;i<nbe;i++)
-         for (int j=0;j<2;j++)
-                {
-
-                 Vertex V;
-                 VertexOnGeom GV;
-                 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
-                        {
-                         GeometricalEdge * eg = GV;
-                         Real8 s = GV;
-                         R2 tg;
-                         Real8  R1= eg->R1tg(s,tg);
-                         Real8 ht = hmax;
-                         if (R1>1.0e-20) {  // err relative to the length of the edge
-                                 ht = Min(Max(errC/R1,hmin),hmax);
-                         }
-                         Real8 hn = iso? ht : Min(hmax,ht*maxaniso);
-                         if (ht<=0 || hn<=0){
-                                 throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
-                         }
-                         MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
-                         Metric MVp(Vp);
-                         edges[i][j].m.IntersectWith(MVp);
-                        }
-
-                }
+
+        //loop over all the vertices on edges
+        for (Int4  i=0;i<nbe;i++){
+                for (int j=0;j<2;j++){
+
+                        Vertex V;
+                        VertexOnGeom GV;
+                        Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
+
+                        GeometricalEdge * eg = GV;
+                        Real8 s = GV;
+                        R2 tg;
+                        Real8  R1= eg->R1tg(s,tg);
+                        Real8  ht=hmax;
+                        // err relative to the length of the edge
+                        if (R1>1.0e-20) {  
+                                ht = Min(Max(errC/R1,hmin),hmax);
+                        }
+                        Real8 hn=Min(hmax,ht*anisomax);
+                        if (ht<=0 || hn<=0){
+                                throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
+                        }
+                        MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
+                        Metric MVp(Vp);
+                        edges[i][j].m.IntersectWith(MVp);
+                }
+        }
         // the problem is for the vertex on vertex 
 }

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

@@ -128 +128 @@
         /*}}}1*/
         /*FUNCTION Vertex::MetricFromHessian{{{1*/
-        void Vertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,BamgOpts* bamgopts){
+        void Vertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){
                 /*Compute Metric from Hessian*/
 
@@ -152 +152 @@
                  */
                 if (Metrictype==0){
-                        ci= 2.0/9.0 * 1/(bamgopts->err*coef*coef);
+                        ci= 2.0/9.0 * 1/(err*coef*coef);
                 }
 
@@ -163 +163 @@
                  */
                 else if (Metrictype==1){
-                        ci= 2.0/9.0 * 1/(bamgopts->err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
+                        ci= 2.0/9.0 * 1/(err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
                 }
 
@@ -173 +173 @@
                  */
                 else if (Metrictype==2){
-                        ci= 2.0/9.0 * 1/(bamgopts->err*coef*coef) * 1/(smax-smin);
+                        ci= 2.0/9.0 * 1/(err*coef*coef) * 1/(smax-smin);
                 }
                 else{

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

@@ -17 +17 @@
         int     Metrictype;
         int     nbjacobi;
-        int     AbsError;
         double  omega;
         double  hmin;
@@ -25 +24 @@
         int     splitcorners;
         int     verbose;
-        double  err;
-        double  errg;
+        double*  err;
         double  coef;
         double* metric;

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

@@ -76 +76 @@
 
 % Bamg Options {{{1
-bamg_options.iso=getfieldvalue(options,'iso',0);
 bamg_options.err=getfieldvalue(options,'err',0.01);
-bamg_options.errg=getfieldvalue(options,'errg',0.1);
 bamg_options.coef=getfieldvalue(options,'coef',1);
 bamg_options.power=getfieldvalue(options,'power',1);
 bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1);
-bamg_options.AbsError=getfieldvalue(options,'AbsError',0);
 bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
 bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
@@ -100 +97 @@
 %}}}
 
+%check bamg options
+fieldsize=size(bamg_options.field);
+if fieldsize(1),
+        if fieldsize(1)~=md.numberofgrids,
+                error(['bamg error message, ''field'' should have ' num2str(md.numberofgrids) ' lines (numberofgrids)']);
+        end
+        if (fieldsize(2)~=0 & length(bamg_options.err)~=fieldsize(2)),
+                error(['bamg error message, ''err'' should be of length ' num2str(fieldsize(2)) ' (As many as fields)']);
+        end
+end
+
+
 %call Bamg
 [triangles vertices metric]=Bamg(bamg_mesh,bamg_geometry,bamg_options);

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

@@ -32 +32 @@
 
         /*Options inputs*/
-        int    iso,maxnbv,verbose,splitcorners;
+        int    maxnbv,verbose,splitcorners;
         double hmin,hmax,anisomax;
-        double err,errg,coef;
+        double coef;
         double power;
         int    Hessiantype,Metrictype,NbSmooth;
-        int    nbjacobi,AbsError;
+        int    nbjacobi;
         double omega;
         double gradation;
@@ -43 +43 @@
         double* metric=NULL;
         double* field=NULL;
-        int     numfields=0;
+        double* err;
+        int     numfields=0,numerr=0;
 
         /*Boot module: */
@@ -108 +109 @@
 
         /*create bamg options input*/
-        FetchData(&iso,mxGetField(BAMGOPTIONS,0,"iso"));
-        bamgopts.iso=iso;
-        FetchData(&err,mxGetField(BAMGOPTIONS,0,"err"));
-        bamgopts.err=err;
-        FetchData(&errg,mxGetField(BAMGOPTIONS,0,"errg"));
-        bamgopts.errg=errg;
         FetchData(&coef,mxGetField(BAMGOPTIONS,0,"coef"));
         bamgopts.coef=coef;
@@ -124 +119 @@
         FetchData(&nbjacobi,mxGetField(BAMGOPTIONS,0,"nbjacobi"));
         bamgopts.nbjacobi=nbjacobi;
-        FetchData(&AbsError,mxGetField(BAMGOPTIONS,0,"AbsError"));
-        bamgopts.AbsError=AbsError;
         FetchData(&NbSmooth,mxGetField(BAMGOPTIONS,0,"NbSmooth"));
         bamgopts.NbSmooth=NbSmooth;
@@ -153 +146 @@
         bamgopts.numfields=numfields;
         bamgopts.field=field;
+        FetchData(&err,NULL,&numerr,mxGetField(BAMGOPTIONS,0,"err"));
+        bamgopts.err=err;
+
+        /*Some checks*/
+        if (numfields!=0 && numerr!=numfields){
+                throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));
+        }
 
         /*!Generate internal degree of freedom numbers: */