Changeset 3022


Ignore:
Timestamp:
02/11/10 15:28:10 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added some doc and some checks

Location:
issm/trunk/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r2981 r3022  
    9292                //generate mesh
    9393                if (verbosity>1) printf("   Generating Mesh...\n");
    94                 Triangles Th(maxnbv,Gh);
     94                Triangles Th(maxnbv,Gh,bamgopts);
    9595                if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
    9696                Th.ReNumberingTheTriangleBySubDomain();
    … …  
    103103                if (verbosity>1) printf("   Write Geometry...\n");
    104104                Gh.WriteGeometry(bamggeom,bamgopts);
     105
     106                //clean up
     107        //      delete &Th;
     108        //      delete &Gh;
    105109                /*}}}*/
    106110        }
    … …  
    183187
    184188                /*clean up*/
    185                 delete &Th; //TEST crash
     189                delete &Th;
     190                //delete &BTh;
    186191                /*}}}*/
    187192        }

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

    r2992 r3022  
    655655
    656656                        //Genetare mesh from geometry
    657                         Triangles(Int4 nbvx,Geometry & G) :Gh(G),BTh(*this){
    658                                 try { GeomToTriangles0(nbvx);}
     657                        Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
     658                                try { GeomToTriangles0(nbvx,bamgopts);}
    659659                                catch(...) { this->~Triangles(); throw; }
    660660                        }
    … …  
    680680                        void Insert();
    681681                        void ForceBoundary();
    682                         void Heap();
     682                        //void Heap();
     683                        void Renumber(BamgOpts* bamgopts);
    683684                        void FindSubDomain(int OutSide=0);
    684685                        Int4 ConsRefTriangle(Int4 *) const;
    … …  
    742743                private:
    743744                        void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption
    744                         void GeomToTriangles0(Int4 nbvx);                   // the real constructor mesh generator
     745                        void GeomToTriangles0(Int4 nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
    745746                        void PreInit(Int4);
    746747
    … …  
    755756                        Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux,
    756757                        Int4 NbSubDomains; //
    757                         Int4 NbEquiEdges;
    758                         Int4 NbCrackedEdges;
    759758                        //  Int4 nbtf;//  de triangle frontiere
    760759                        Int4 NbOfCurves;
    … …  
    810809        // to sort in descending order
    811810        template<class T>  inline void  HeapSort(T *c,long n){
    812                 long l,j,r,i;
    813                 T crit;
    814                 c--; // on decale de 1 pour que le tableau commence a 1
    815                 if( n <= 1) return;
    816                 l = n/2 + 1;
    817                 r = n;
    818                 while (1) { // label 2
    819                         if(l <= 1 ) { // label 20
    820                                 crit = c[r];
    821                                 c[r--] = c[1];
    822                                 if ( r == 1 ) { c[1]=crit; return;}
    823                         } else  crit = c[--l];
     811                /*Intermediary*/
     812                int l,j,r,i;
     813                T   crit;
     814                c--;                    //the array must starts at 1 and not 0
     815                if(n<=1) return;        //return if size <=1
     816                l=n/2+1;                //initialize l and r
     817                r=n;
     818                for(;;){
     819                        if(l<=1){
     820                                crit  =c[r];
     821                                c[r--]=c[1];
     822                                if (r==1){c[1]=crit; return;}
     823                        }
     824                        else  crit = c[--l];
    824825                        j=l;
    825                         while (1) {// label 4
     826                        for(;;){
    826827                                i=j;
    827828                                j=2*j;
    828                                 if  (j>r) {c[i]=crit;break;} // L8 -> G2
    829                                 if ((j<r) && (c[j] < c[j+1])) j++; // L5
    830                                 if (crit < c[j]) c[i]=c[j]; // L6+1 G4
    831                                 else {c[i]=crit;break;} //L8 -> G2
     829                                if  (j>r) {c[i]=crit;break;}
     830                                if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
     831                                if (crit < c[j]) c[i]=c[j];       //c[j]  > crit -> stock this large value in i(<j)
     832                                else{c[i]=crit;break;}            //c[j]  < crit -> stock crit in i (<j)
     833                        }
     834                }
     835        }
     836        // to sort in descending order and return ordering
     837        template<class T>  inline void  HeapSort(int** porder,T* c,int n){
     838                /*Intermediary*/
     839                int  l,j,r,i;
     840                T    crit;
     841                int  pos;
     842                int* order=NULL;
     843                order=(int*)xmalloc(n*sizeof(int));
     844                for(i=0;i<n;i++) order[i]=i+1;
     845                c--;                    //the array must starts at 1 and not 0
     846                order--;
     847                if(n<=1) return;        //return if size <=1
     848                l=n/2+1;                //initialize l and r
     849                r=n;
     850                for(;;){
     851                        if(l<=1){
     852                                crit  =c[r]; pos=order[r];
     853                                c[r--]=c[1]; order[r+1]=order[1];
     854                                if (r==1){
     855                                        c[1]=crit; order[1]=pos;
     856                                        order++;
     857                                        *porder=order;
     858                                        return;
     859                                }
     860                        }
     861                        else  {crit=c[--l]; pos=order[l];}
     862                        j=l;
     863                        for(;;){
     864                                i=j;
     865                                j=2*j;
     866                                if  (j>r) {c[i]=crit;order[i]=pos;break;}
     867                                if ((j<r) && (c[j] < c[j+1]))j++;
     868                                if (crit < c[j]) {c[i]=c[j];order[i]=order[j];}
     869                                else{c[i]=crit;order[i]=pos;break;}
    832870                        }
    833871                }
    … …  
    958996        /*INLINE functions of CLASS Triangles{{{1*/
    959997        inline Triangles::Triangles(Int4 i) :Gh(*new Geometry()),BTh(*this){PreInit(i);}
    960         inline  void  Triangles::ReMakeTriangleContainingTheVertex()
    961           {
     998        inline  void  Triangles::ReMakeTriangleContainingTheVertex(){
    962999                register Int4 i;
    963                 for ( i=0;i<nbv;i++)
    964                   {
    965                         vertices[i].vint = 0;
     1000                for ( i=0;i<nbv;i++){
     1001                        vertices[i].vint=0;
    9661002                        vertices[i].t=0;
    9671003                  }
    968                 for ( i=0;i<nbt;i++)
    969                  triangles[i].SetTriangleContainingTheVertex();
     1004                for ( i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
    9701005          }
    9711006

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

    r2966 r3022  
    893893                printf("   nbv  (number of initial vertices) : %i\n",nbiv);
    894894                printf("   NbSubDomains: %i\n",NbSubDomains);
    895                 printf("   NbEquiEdges: %i\n",NbEquiEdges);
    896                 printf("   NbCrackedEdges: %i\n",NbCrackedEdges);
    897895                printf("   NbOfCurves: %i\n",NbOfCurves);
    898896                printf("   vertices: %p\n",vertices);

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

    r2981 r3022  
    479479
    480480                int i,j,k,num;
    481                 int verbose;
    482 
    483                 verbose=bamgopts->verbose;
     481                int verbose=bamgopts->verbose;
     482
     483                //renumber
     484                if (bamgopts->renumber){
     485                        if(verbose>3) printf("      Renumbering...");
     486                        Renumber(bamgopts);
     487                        if(verbose>3) printf(" done\n");
     488                }
    484489
    485490                //Build reft that holds the number the subdomain number of each triangle
    486                 Int4 *reft = new Int4[nbt];
     491                Int4* reft = new Int4[nbt];
    487492                Int4 nbInT = ConsRefTriangle(reft);
    488493
    … …  
    524529
    525530                //Segments
     531                bamgmesh->NumSegments=0;
    526532                if(verbose>3) printf("      writing Segments\n");
     533
    527534                //chaining algorithm
    528                 int head_v[nbv];
    529                 int next_p[3*nbt];
     535                int* head_v=NULL;
     536                head_v=(int*)xmalloc(nbv*sizeof(int));
     537                int* next_p=NULL;
     538                next_p=(int*)xmalloc(3*nbt*sizeof(int));
     539
    530540                for (i=0;i<nbv;i++) head_v[i]=-1;
    531541                k=0;
    … …  
    535545                                for (j=0;j<3;j++){
    536546                                        int v=Number(triangles[i][j]); //jth vertex of the ith triangle
     547                                        if (k>3*nbt-1 || k<0) throw ErrorException(__FUNCT__,exprintf("k = %i, nbt = %i",k,nbt));
    537548                                        next_p[k]=head_v[v];
     549                                        if (v>nbv-1 || v<0) throw ErrorException(__FUNCT__,exprintf("v = %i, nbv = %i",v,nbv));
    538550                                        head_v[v]=k++;
    539551                                }
    … …  
    582594                        }
    583595                }
     596                //clean up
     597                xfree((void**)&head_v);
     598                xfree((void**)&next_p);
    584599
    585600                //CrackedEdges
    … …  
    615630                                        bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
    616631                                        num=num+1;
    617                                         /*if (reft[i]==-1){
    618                                           if (!t(0)) printf("%i %i    %i\n",Number(t[1])+1,Number(t[2]),i);
    619                                           if (!t(1)) printf("%i %i    %i\n",Number(t[2])+1,Number(t[0]),i);
    620                                           if (!t(2)) printf("%i %i    %i\n",Number(t[0])+1,Number(t[1]),i);
    621                                           }*/
    622632                                }
    623633                        }
    … …  
    746756                        bamgmesh->EdgesOnGeometricEdge=NULL;
    747757                }
     758
     759
     760                //clean up
     761                delete [] reft;
    748762        }
    749763        /*}}}1*/
    … …  
    14621476
    14631477                //first, build the chains that will be used for the Hessian computation, as weel as the area of each element
    1464                 int head_s[nbv];
    1465                 int next_p[3*nbt];
     1478                int* head_s=NULL;
     1479                head_s=(int*)xmalloc(nbv*sizeof(int));
     1480                int* next_p=NULL;
     1481                next_p=(int*)xmalloc(3*nbt*sizeof(int));
    14661482                int  p=0;
    14671483                //initialization
    … …  
    16021618
    16031619                //clean up
     1620                xfree((void**)&head_s);
     1621                xfree((void**)&next_p);
    16041622                delete [] detT;
    16051623                delete [] alpha;
    … …  
    27692787        /*}}}1*/
    27702788        /*FUNCTION Triangles::GeomToTriangles0{{{1*/
    2771         void Triangles::GeomToTriangles0(Int4 inbvx) {
     2789        void Triangles::GeomToTriangles0(Int4 inbvx,BamgOpts* bamgopts) {
    27722790                /*Generate mesh from geometry*/
    27732791
    … …  
    27842802                GeometricalEdge * e;
    27852803
     2804                /*Get options*/
     2805                int verbose=bamgopts->verbose;
     2806
    27862807                //initialize this
     2808                if (verbose>3) printf("      Generating vertices\n");
    27872809                PreInit(inbvx);
    27882810                nbv=0;
    … …  
    30643086
    30653087                //Insert points inside existing triangles
     3088                if (verbose>3) printf("      Inserting vertices in mesh\n");
    30663089                Insert();
    30673090
    30683091                //Force the boundary
     3092                if (verbose>3) printf("      Forcing boundaries\n");
    30693093                ForceBoundary();
    30703094
    30713095                //Extract SubDomains
     3096                if (verbose>3) printf("      Extracting subdomains\n");
    30723097                FindSubDomain();
    30733098
    30743099                // NewPointsOld(*this) ;
     3100                if (verbose>3) printf("      Generating mesh properties\n");
    30753101                NewPoints(*this,0) ;
    30763102        }
    … …  
    33813407                                        }
    33823408                                  }//for(phase;;)
    3383                                 /*  modif FH add Curve class             
    3384                                          }}//for (iedge=0;iedge<BTh.nbe;iedge++)
    3385                                          */
    3386                                 // new code Curve class         
    33873409                  } //  end of curve loop
    33883410                // end new code     
    … …  
    38863908        /*FUNCTION Triangles::NewPoints{{{1*/
    38873909        void  Triangles::NewPoints(Triangles & Bh,int KeepVertices) {
     3910
    38883911                long int verbosity=2;
    38893912                Int4 nbtold(nbt),nbvold(nbv);
    3890                 if (verbosity>2)  printf("   Triangles::NewPoints\n");
    38913913                Int4 i,k;
    38923914                int j;
    38933915                Int4 *first_np_or_next_t = new Int4[nbtx];
    38943916                Int4 NbTSwap =0;
    3895                 //    insert old point
     3917
     3918                //display info
     3919                if (verbosity>2)  printf("   Triangles::NewPoints\n");
     3920
     3921                /*First, insert old points*/
     3922
    38963923                nbtold = nbt;
    38973924
    38983925                if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){
    38993926                        //   Bh.SetVertexFieldOn();
    3900                         for (i=0;i<Bh.nbv;i++)
    3901                           {
    3902                                 Vertex & bv = Bh[i];
    3903                                 if (!bv.onGeometry) {
     3927                        for (i=0;i<Bh.nbv;i++){
     3928                                Vertex &bv=Bh[i];
     3929                                if (!bv.onGeometry){
    39043930                                        vertices[nbv].r = bv.r;
    3905                                         vertices[nbv++].m = bv.m;}
    3906                           }
     3931                                        vertices[nbv++].m = bv.m;
     3932                                }
     3933                        }
    39073934                        int nbv1=nbv;
    39083935                        Bh.ReMakeTriangleContainingTheVertex();     
    39093936                        InsertNewPoints(nbvold,NbTSwap)   ;           
    39103937                } 
    3911                 else
    3912                  Bh.ReMakeTriangleContainingTheVertex();     
     3938                else Bh.ReMakeTriangleContainingTheVertex();     
    39133939
    39143940                Triangle *t;
    … …  
    39223948                // the next traingle on i is  -first_np_or_next_t[i]
    39233949                int iter=0;
    3924                 // Big loop
     3950
     3951                // Big loop (most time consuming)
    39253952                do {
    39263953                        iter++;
    … …  
    39293956
    39303957                        // default size of  IntersectionTriangle
    3931 
    39323958                        i=Headt;
    39333959                        next_t=-first_np_or_next_t[i];
    … …  
    39994025                Int4 NbSwapf =0,NbSwp;
    40004026
    4001                 // bofbof
    4002 
    4003 
    40044027                NbSwp = NbSwapf;
    40054028                for (i=0;i<nbv;i++)
    40064029                 NbSwapf += vertices[i].Optim(0);
    4007                 /*
    4008                         for (i=0;i<nbv;i++)
    4009                         NbSwapf += vertices[i].Optim(0);
    4010                         for (i=0;i<nbv;i++)
    4011                         NbSwapf += vertices[i].Optim(0);
    4012                         for (i=0;i<nbv;i++)
    4013                         NbSwapf += vertices[i].Optim(0);
    4014                         for (i=0;i<nbv;i++)
    4015                         NbSwapf += vertices[i].Optim(0);
    4016                         */
    40174030                NbTSwap +=  NbSwapf ;
    40184031        }
    … …  
    42124225        }                 
    42134226        /*}}}1*/
     4227/*FUNCTION Triangles::Renumber {{{1*/
     4228void Triangles::Renumber(BamgOpts* bamgopts){
     4229
     4230        /*Intermediary*/
     4231        int i,j,k;
     4232        int*    r=NULL;
     4233        int*    rprime=NULL;
     4234        rprime  =(int*)   xmalloc(nbv*sizeof(int));
     4235
     4236        //renumbering with respect to lower left corner
     4237        if (bamgopts->renumber==1){
     4238
     4239                //compute distance to pmin
     4240                double* distance;
     4241                distance=(double*)xmalloc(nbv*sizeof(double));
     4242                for (i=0;i<nbv;i++) distance[i]=pow(vertices[i].r.x - pmin.x,2)+pow(vertices[i].r.y - pmin.y,2);
     4243
     4244                //get ordering
     4245                HeapSort(&r,distance,nbv);
     4246                xfree((void**)&distance);
     4247        }
     4248        //renumbering randomly
     4249        else if (bamgopts->renumber==2){
     4250
     4251                //allocate memory to r
     4252                r=(int*)xmalloc(nbv*sizeof(int));
     4253                int PrimeNumber= AGoodNumberPrimeWith(nbv) ;
     4254                int k0=rand()%nbv;
     4255                for (int i=0; i<nbv; i++){
     4256                        r[i]=(k0=(k0+PrimeNumber)%nbv)+1;
     4257                }
     4258        }
     4259        //else: error
     4260        else{
     4261                throw ErrorException(__FUNCT__,exprintf("renumbering option %i not supported yet",bamgopts->renumber));
     4262        }
     4263
     4264        //Keep copy of old vertices
     4265        Vertex* oldvertices=NULL;
     4266        oldvertices=(Vertex*)xmalloc(nbv*sizeof(Vertex));
     4267        for (i=0;i<nbv;i++){
     4268                oldvertices[i]=vertices[i];
     4269        }
     4270
     4271        //renumber vertices
     4272        for (i=0;i<nbv;i++){
     4273                //check r[i]
     4274                if (r[i]>nbv){
     4275                        throw ErrorException(__FUNCT__,exprintf("r[i]>nbv (r[i]=%i and nbv=%i)",r[i],nbv));
     4276                }
     4277                if (r[i]<=0){
     4278                        throw ErrorException(__FUNCT__,exprintf("r[i]<=0 (r[i]=%i)",r[i]));
     4279                }
     4280                vertices[i]=oldvertices[r[i]-1];
     4281                rprime[r[i]-1]=i+1;
     4282        }
     4283        //clean up
     4284        xfree((void**)&oldvertices);
     4285
     4286        //renumber triangles
     4287        for (i=0; i<nt;i++){
     4288                for (j=0;j<3;j++){
     4289                        if (triangles[i](j)){
     4290                                triangles[i](j)=&vertices[rprime[Number(triangles[i](j))] -1 ];
     4291                        }
     4292                }
     4293        }
     4294
     4295        //renumber edges
     4296        for (i=0; i<nbe;i++){
     4297                for (j=0;j<2;j++){
     4298                        edges[i].v[j]=&vertices[rprime[Number(edges[i].v[j])] -1 ];
     4299                }
     4300        }
     4301
     4302        //clean up
     4303        xfree((void**)&r);
     4304        xfree((void**)&rprime);
     4305}
     4306/*}}}1*/
    42144307        /*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/
    42154308        void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) {

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

    r2942 r3022  
    1414        int     Metrictype;
    1515        int     KeepVertices;
     16        int     renumber;
    1617        double  maxsubdiv;
    1718        double  power;

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

    r2977 r3022  
    22%BAMG - mesh generation
    33%
    4 %   Usage:
    5 %      md=bamg(md,varargin);
     4%   Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):
     5%
     6%   - domain : followed by an ARGUS file that prescribes the domain outline (rifts and holes)
     7%   - hmin   : minimum edge length (default is 10^-100)
     8%   - hmax   : maximum esge length (default is 10^100)
     9%   - hVertices : imposed edge length for each vertex (geometry or mesh)
     10%
     11%   - anisomax    : maximum ration between the smallest and largest edges (default is 10^30)
     12%   - coef        : coefficient applied to the metric (2-> twice as many elements, default is 1)
     13%   - cutoff      : scalar used to compute the metric when metric type 2 or 3 are applied
     14%   - err         : error used to generate the metric from a field
     15%   - errg        : geometrical error (default is 0.1)
     16%   - field       : field of the model that will be used to compute the metric
     17%                   to apply several fields, use one column per field
     18%   - gradation   : maximum ration between two adjacent edges
     19%   - Hessiantype : 0 -> use double P2 projection (default)
     20%                   1 -> use Green formula
     21%   - KeepVertices: try to keep initial vertices when adaptation is done on an existing mesh (default 1)
     22%   - MaximalAngleOfCorner: maximal angle of corners in degree (default is 10)
     23%   - maxnbv      : maximum number of vertices used to allocate memory (default is 10^6)
     24%   - maxsubdiv   : maximum subdivision of exisiting elements (default is 10)
     25%   - metric      : matrix (numberofgrids x 3) used as a metric
     26%   - Metrictype  : 1 -> absolute error          c/(err coeff^2) * Abs(H)        (default)
     27%                   2 -> relative error          c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s))
     28%                   3 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin)
     29%   - nbjacoby    : correction used by Hessiantype=1 (default is 1)
     30%   - NbSmooth    : number of metric smoothing procedure (default is 3)
     31%   - omega       : relaxation parameter of the smoothing procedure (default is 1.8)
     32%   - power       : power applied to the metric (default is 1)
     33%   - renumber    : 0 -> no renumbering
     34%                   1 -> renumber by distance to lower left corner (default)
     35%                   2 -> random renumbering
     36%   - splitcorner : split triangles whuch have 3 vertices on the outline (default is 1)
     37%   - verbose     : level of verbosity (default is 1)
     38%
     39%   Examples:
     40%      md=bamg(md,'domain','DomainOutline.exp','hmax',3000);
     41%      md=bamg(md,'field',[md.vel_obs md.thickness],'hmax',20000,'hmin',1000);
     42%      md=bamg(md,'metric',A,'hmin',1000,'hmax',20000,'gradation',3,'anisomax',1);
    643
    744%process options
    … …  
    104141
    105142% Bamg Options {{{1
     143bamg_options.anisomax=getfieldvalue(options,'anisomax',10^30);
     144bamg_options.coef=getfieldvalue(options,'coef',1);
     145bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
    106146bamg_options.err=getfieldvalue(options,'err',0.01);
    107147bamg_options.errg=getfieldvalue(options,'errg',0.1);
    108 bamg_options.coef=getfieldvalue(options,'coef',1);
     148bamg_options.field=getfieldvalue(options,'field',[]);
     149bamg_options.gradation=getfieldvalue(options,'gradation',1.5);
     150bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
     151bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
     152bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
     153bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1);
     154bamg_options.MaximalAngleOfCorner=getfieldvalue(options,'MaximalAngleOfCorner',10);
     155bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6);
    109156bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10);
    110 bamg_options.power=getfieldvalue(options,'power',1);
     157bamg_options.metric=getfieldvalue(options,'metric',[]);
     158bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
    111159bamg_options.nbjacobi=getfieldvalue(options,'nbjacobi',1);
    112 bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
    113 bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
    114 bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1);
    115160bamg_options.NbSmooth=getfieldvalue(options,'NbSmooth',3);
    116161bamg_options.omega=getfieldvalue(options,'omega',1.8);
    117 bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6);
    118 bamg_options.MaximalAngleOfCorner=getfieldvalue(options,'MaximalAngleOfCorner',10);
    119 bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
    120 bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
    121 bamg_options.anisomax=getfieldvalue(options,'anisomax',10^30);
    122 bamg_options.gradation=getfieldvalue(options,'gradation',1.5);
    123 bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
     162bamg_options.power=getfieldvalue(options,'power',1);
     163bamg_options.renumber=getfieldvalue(options,'renumber',1);
     164bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
    124165bamg_options.verbose=getfieldvalue(options,'verbose',1);
    125 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
    126 bamg_options.metric=getfieldvalue(options,'metric',[]);
    127 bamg_options.field=getfieldvalue(options,'field',[]);
    128166%}}}
    129167

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

    r2977 r3022  
    1010        int   i;
    1111        int   nods=0; //to be removed
     12        int   lines,cols;
    1213        BamgOpts bamgopts;
    1314        BamgMesh bamgmesh;
    … …  
    4142        double power;
    4243        int    Hessiantype,Metrictype,NbSmooth;
    43         int    nbjacobi,KeepVertices;
     44        int    nbjacobi,KeepVertices,renumber;
    4445        double omega;
    4546        double gradation,errg;
    … …  
    6566        FetchData(&EdgesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges"));
    6667        bamggeom.Edges=EdgesGeom;
    67         FetchData(&hVerticesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"hVertices"));
     68        FetchData(&hVerticesGeom,&lines,&cols,mxGetField(BAMGGEOMETRY,0,"hVertices"));
     69        if (hVerticesGeom && (cols!=1 || lines!=NumVerticesGeom)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",NumVerticesGeom,1));}
    6870        bamggeom.hVertices=hVerticesGeom;
    6971        bamggeom.MetricVertices=NULL;
    … …  
    9597        FetchData(&TrianglesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles"));
    9698        bamgmesh.Triangles=TrianglesMesh;
    97         FetchData(&hVerticesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"hVertices"));
     99        FetchData(&hVerticesMesh,&lines,&cols,mxGetField(BAMGMESH,0,"hVertices"));
     100        if (hVerticesMesh && (cols!=1 || lines!=NumVerticesMesh)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",NumVerticesMesh,1));}
    98101        bamgmesh.hVertices=hVerticesMesh;
    99102        bamgmesh.NumQuadrilaterals=0;
    … …  
    121124        /*create bamg options input*/
    122125        FetchData(&coef,mxGetField(BAMGOPTIONS,0,"coef"));
     126        if (coef==0) throw ErrorException(__FUNCT__,exprintf("'coef' should be positive"));
    123127        bamgopts.coef=coef;
    124128        FetchData(&maxsubdiv,mxGetField(BAMGOPTIONS,0,"maxsubdiv"));
     129        if (maxsubdiv<1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >=1"));
    125130        bamgopts.maxsubdiv=maxsubdiv;
    126131        FetchData(&Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype"));
     132        if (Hessiantype!=0 && Hessiantype!=1) throw ErrorException(__FUNCT__,exprintf("'Hessiantype' supported options are 0 and 1"));
    127133        bamgopts.Hessiantype=Hessiantype;
    128134        FetchData(&Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype"));
     135        if (Metrictype!=0 && Metrictype!=1 && Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2"));
    129136        bamgopts.Metrictype=Metrictype;
    130137        FetchData(&KeepVertices,mxGetField(BAMGOPTIONS,0,"KeepVertices"));
     138        if (KeepVertices!=0 && KeepVertices!=1) throw ErrorException(__FUNCT__,exprintf("'KeepVertices' supported options are 0 and 1"));
    131139        bamgopts.KeepVertices=KeepVertices;
     140        FetchData(&renumber,mxGetField(BAMGOPTIONS,0,"renumber"));
     141        if (renumber!=0 && renumber!=1 && renumber!=2) throw ErrorException(__FUNCT__,exprintf("'renumber' supported options are 0, 1 and 2"));
     142        bamgopts.renumber=renumber;
    132143        FetchData(&power,mxGetField(BAMGOPTIONS,0,"power"));
    133144        bamgopts.power=power;
    134145        FetchData(&errg,mxGetField(BAMGOPTIONS,0,"errg"));
     146        if (errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0"));
    135147        bamgopts.errg=errg;
    136148        FetchData(&nbjacobi,mxGetField(BAMGOPTIONS,0,"nbjacobi"));
     149        if (nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0"));
    137150        bamgopts.nbjacobi=nbjacobi;
    138151        FetchData(&NbSmooth,mxGetField(BAMGOPTIONS,0,"NbSmooth"));
     152        if (NbSmooth<=0) throw ErrorException(__FUNCT__,exprintf("'NbSmooth' option should be >0"));
    139153        bamgopts.NbSmooth=NbSmooth;
    140154        FetchData(&omega,mxGetField(BAMGOPTIONS,0,"omega"));
    141155        bamgopts.omega=omega;
    142156        FetchData(&maxnbv,mxGetField(BAMGOPTIONS,0,"maxnbv"));
     157        if (maxnbv<3) throw ErrorException(__FUNCT__,exprintf("'maxnbv' option should be >3"));
    143158        bamgopts.maxnbv=maxnbv;
    144159        FetchData(&hmin,mxGetField(BAMGOPTIONS,0,"hmin"));
     160        if (hmin<=0) throw ErrorException(__FUNCT__,exprintf("'hmin' option should be >0"));
    145161        bamgopts.hmin=hmin;
    146162        FetchData(&hmax,mxGetField(BAMGOPTIONS,0,"hmax"));
     163        if (hmax<=0 || hmax<hmin) throw ErrorException(__FUNCT__,exprintf("'hmax' option should be between 0 and hmin=%g",hmin));
    147164        bamgopts.hmax=hmax;
    148165        FetchData(&anisomax,mxGetField(BAMGOPTIONS,0,"anisomax"));
     166        if (anisomax<1) throw ErrorException(__FUNCT__,exprintf("'anisomax' option should be >=1"));
    149167        bamgopts.anisomax=anisomax;
    150168        FetchData(&gradation,mxGetField(BAMGOPTIONS,0,"gradation"));
     169        if (anisomax<1) throw ErrorException(__FUNCT__,exprintf("'anisomax' option should be >=1"));
    151170        bamgopts.gradation=gradation;
    152171        FetchData(&cutoff,mxGetField(BAMGOPTIONS,0,"cutoff"));
    … …  
    158177        FetchData(&MaximalAngleOfCorner,mxGetField(BAMGOPTIONS,0,"MaximalAngleOfCorner"));
    159178        bamgopts.MaximalAngleOfCorner=MaximalAngleOfCorner;
    160         FetchData(&metric,NULL,NULL,mxGetField(BAMGOPTIONS,0,"metric"));
     179        FetchData(&metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric"));
     180        if (metric && (cols!=3 || lines!=NumVerticesMesh)){throw ErrorException(__FUNCT__,exprintf("the size of 'metric' should be [%i %i]",NumVerticesMesh,3));}
    161181        bamgopts.metric=metric;
    162         FetchData(&field,NULL,&numfields,mxGetField(BAMGOPTIONS,0,"field"));
     182        FetchData(&field,&lines,&numfields,mxGetField(BAMGOPTIONS,0,"field"));
     183        if (field && lines!=NumVerticesMesh){throw ErrorException(__FUNCT__,exprintf("the size of 'field' should be [%i %i]",NumVerticesMesh,numfields));}
    163184        bamgopts.numfields=numfields;
    164185        bamgopts.field=field;
    165         FetchData(&err,NULL,&numerr,mxGetField(BAMGOPTIONS,0,"err"));
     186        FetchData(&err,NULL,&cols,mxGetField(BAMGOPTIONS,0,"err"));
     187        if (numfields!=0 && cols!=numfields){throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));}
     188        for (i=0;i<numfields;i++) {if (err[i]<=0) throw ErrorException(__FUNCT__,exprintf("'err' option should be >0"));};
    166189        bamgopts.err=err;
    167 
    168         /*Some checks*/
    169         if (numfields!=0 && numerr!=numfields){
    170                 throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));
    171         }
    172190
    173191        /*!Generate internal degree of freedom numbers: */
    … …  
    192210{
    193211        _printf_("\n");
    194         _printf_("   usage: [elements,x,y]=%s(bamgmesh,bamggeom,bamgoptions,nbv);\n",__FUNCT__);
     212        _printf_("   usage: [elements,vertices,segments,segmentmarkers,metric]=%s(bamgmesh,bamggeom,bamgoptions);\n",__FUNCT__);
    195213        _printf_("\n");
    196214}
Note: See TracChangeset for help on using the changeset viewer.