Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2940) +++ /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2941) @@ -33,8 +33,8 @@ /*Bamg options*/ - int iso; + double anisomax; int maxnbv; double hmin,hmax; - double err,errg,coef; + double coef; int verbosity; int NbSmooth; @@ -46,5 +46,4 @@ int allquad=0; double costheta=2; - double anisomax = 1e6; int KeepBackVertices=1; double hminaniso=1e-100; @@ -56,7 +55,5 @@ /*Bamg options*/ - iso=bamgopts->iso; - err=bamgopts->err; - errg=bamgopts->errg; + anisomax=bamgopts->anisomax; NbSmooth=bamgopts->NbSmooth; omega=bamgopts->omega; @@ -71,6 +68,4 @@ 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,5 +154,5 @@ } - BTh.IntersectGeomMetric(errg,iso); + BTh.IntersectGeomMetric(bamgopts); if (verbosity>1) printf(" Smoothing Metric..."); if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation); Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2940) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2941) @@ -137,5 +137,5 @@ // ~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,5 +802,5 @@ void BuildMetric0(BamgOpts* bamgopts); void BuildMetric1(BamgOpts* bamgopts); - void IntersectGeomMetric(const Real8 err,const int iso); + void IntersectGeomMetric(BamgOpts* bamgopts); @@ -1078,6 +1078,5 @@ } - inline void Triangles::SetVertexFieldOn() - { + inline void Triangles::SetVertexFieldOn(){ for (Int4 i=0;imetric); - if (bamgopts->iso){ - bamgopts->metric=(double*)xmalloc(1*nbv*sizeof(double)); - for (i=0;imetric[i]=V.hmin(); - } - } - else { - bamgopts->metric=(double*)xmalloc(3*nbv*sizeof(double)); - for (i=0;imetric[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;imetric[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,5 +1059,5 @@ /*Compute Metric from Hessian*/ for ( iv=0;iverr[nusol],bamgopts); } @@ -1361,5 +1352,5 @@ /*Compute Metric from Hessian*/ for ( iv=0;iverr[nusol],bamgopts); } @@ -3242,10 +3233,6 @@ /*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,42 +3251,53 @@ /*}}}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;iR1tg(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;iR1tg(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 } Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 2940) +++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 2941) @@ -128,5 +128,5 @@ /*}}}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,5 +152,5 @@ */ if (Metrictype==0){ - ci= 2.0/9.0 * 1/(bamgopts->err*coef*coef); + ci= 2.0/9.0 * 1/(err*coef*coef); } @@ -163,5 +163,5 @@ */ 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,5 +173,5 @@ */ 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{ Index: /issm/trunk/src/c/objects/BamgOpts.h =================================================================== --- /issm/trunk/src/c/objects/BamgOpts.h (revision 2940) +++ /issm/trunk/src/c/objects/BamgOpts.h (revision 2941) @@ -17,5 +17,4 @@ int Metrictype; int nbjacobi; - int AbsError; double omega; double hmin; @@ -25,6 +24,5 @@ int splitcorners; int verbose; - double err; - double errg; + double* err; double coef; double* metric; Index: /issm/trunk/src/m/classes/public/bamg.m =================================================================== --- /issm/trunk/src/m/classes/public/bamg.m (revision 2940) +++ /issm/trunk/src/m/classes/public/bamg.m (revision 2941) @@ -76,11 +76,8 @@ % 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,4 +97,16 @@ %}}} +%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); Index: /issm/trunk/src/mex/Bamg/Bamg.cpp =================================================================== --- /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 2940) +++ /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 2941) @@ -32,10 +32,10 @@ /*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,5 +43,6 @@ double* metric=NULL; double* field=NULL; - int numfields=0; + double* err; + int numfields=0,numerr=0; /*Boot module: */ @@ -108,10 +109,4 @@ /*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,6 +119,4 @@ 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,4 +146,11 @@ 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: */