Changeset 5397


Ignore:
Timestamp:
08/19/10 09:10:32 (15 years ago)
Author:
Mathieu Morlighem
Message:

As usual

Location:
issm/trunk/src/c/objects/Bamg
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Bamg/Curve.cpp

    r5149 r5397  
    1111        /*Constructors/Destructors*/
    1212        /*FUNCTION Curve::Curve(){{{1*/
    13         Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0),master(true){
     13        Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0){
    1414
    1515        }
  • issm/trunk/src/c/objects/Bamg/Curve.h

    r5130 r5397  
    1616                        int kb,ke;  //  begin vetex and end vertex
    1717                        Curve* next; // next curve equi to this
    18                         bool master; // true => of equi curve point on this curve 
    1918
    2019                        //Methods
  • issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp

    r5151 r5397  
    152152        /*FUNCTION GeometricalEdge::SetCracked{{{1*/
    153153        void   GeometricalEdge::SetCracked()     {
    154                 type |= 1;
     154                type |= 1;/*=>1st digit to 1*/
    155155        }/*}}}*/
    156156        /*FUNCTION GeometricalEdge::SetTgA{{{1*/
    157157        void   GeometricalEdge::SetTgA()         {
    158                 type |=4;
     158                type |=4; /*=>2d digit to 1*/
    159159        }/*}}}*/
    160160        /*FUNCTION GeometricalEdge::SetTgB{{{1*/
    161161        void   GeometricalEdge::SetTgB()         {
    162                 type |=8;
     162                type |=8; /*=> 3d digit to 1*/
    163163        }/*}}}*/
    164164        /*FUNCTION GeometricalEdge::SetMark{{{1*/
    165165        void   GeometricalEdge::SetMark()        {
    166                 type |=16;
     166                type |=16;/*=> 4th digiy to 1*/
    167167        }/*}}}*/
    168168        /*FUNCTION GeometricalEdge::SetUnMark{{{1*/
    169169        void   GeometricalEdge::SetUnMark()      {
    170                 type &= 1007 /* 1023-16*/;
     170                type &= 1007 /* 1023-16 = 000111110111 => 4th digit to 0*/;
    171171        }/*}}}*/
    172172        /*FUNCTION GeometricalEdge::SetRequired{{{1*/
    173173        void   GeometricalEdge::SetRequired()    {
    174                 type |= 64;
     174                type |= 64;/*=>=>6th digit to 1*/
    175175        }/*}}}*/
    176176          /*FUNCTION GeometricalEdge::Tg{{{1*/
  • issm/trunk/src/c/objects/Bamg/Geometry.cpp

    r5340 r5397  
    2222                Init();
    2323                ReadGeometry(bamggeom,bamgopts);
    24                 AfterRead();
     24                PostRead();
    2525        }
    2626        /*}}}*/
     
    165165                        }
    166166
    167                         // definition  the default of the given mesh size
     167                        // definition the default of the given mesh size
    168168                        for (i=0;i<nbv;i++) {
    169169                                if (vertices[i].color > 0)
     
    400400
    401401        /*Methods*/
    402         /*FUNCTION Geometry::AfterRead{{{1*/
    403         void Geometry::AfterRead(){
     402        /*FUNCTION Geometry::PostRead{{{1*/
     403        void Geometry::PostRead(){
    404404                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
    405405
     
    606606                }
    607607
    608                 // generation of  all the tangents
     608                /* generation of  all the tangents*/
    609609                for (i=0;i<nbe;i++) {
    610610                        R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
     
    646646                }
    647647
     648                /* generation of  all curves (from corner to corner)*/
     649                /*We proceed in 2 steps: first allocate, second build*/
    648650                for (int step=0;step<2;step++){
    649651
     652                        //unmark all edges
    650653                        for (i=0;i<nbe;i++) edges[i].SetUnMark();
    651 
     654                        long  nb_marked_edges=0;
     655
     656                        //initialize number of curves
    652657                        nbcurves = 0;
    653                         long  nbgem=0;
    654 
    655                         for (int level=0;level < 2 && nbgem != nbe;level++){
     658
     659                        for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
    656660                                for (i=0;i<nbe;i++){
    657                                         GeometricalEdge & ei = edges[i];   
     661
     662                                        GeometricalEdge & ei=edges[i];   
    658663                                        for(j=0;j<2;j++){
     664                                                /*If current edge ei is unmarked and (level=1 or vertex i is required (corner)):
     665                                                 * we do have the first edge of a new curve*/
    659666                                                if (!ei.Mark() && (level || ei[j].Required())) {
    660                                                         // warning ei.Mark() can be change in loop for(j=0;j<2;j++)
    661667                                                        int k0=j,k1;
    662                                                         GeometricalEdge *e = & ei;
     668                                                        GeometricalEdge   *e=&ei;
    663669                                                        GeometricalVertex *a=(*e)(k0); // begin
    664670                                                        if(curves){
     
    667673                                                        }
    668674                                                        int nee=0;
    669                                                         for(;;) {
     675                                                        for(;;){
    670676                                                                nee++;
    671677                                                                k1 = 1-k0; // next vertex of the edge
    672678                                                                e->SetMark();
    673                                                                 nbgem++;
     679                                                                nb_marked_edges++;
    674680                                                                e->CurveNumber=nbcurves;
    675                                                                 if(curves) {
     681                                                                if(curves){
    676682                                                                        curves[nbcurves].ee=e;
    677683                                                                        curves[nbcurves].ke=k1;
     
    679685
    680686                                                                GeometricalVertex *b=(*e)(k1);
    681                                                                 if (a == b ||  b->Required() ) break;
    682                                                                 k0 = e->AdjVertexNumber[k1];//  vertex in next edge
    683                                                                 e = e->Adj[k1]; // next edge
     687
     688                                                                //break if we have reached the other end of the curve
     689                                                                if (a==b || b->Required()){
     690                                                                        break;
     691                                                                }
     692                                                                //else: go to next edge (adjacent)
     693                                                                else{
     694                                                                        k0 = e->AdjVertexNumber[k1];//  vertex in next edge
     695                                                                        e  = e->Adj[k1]; // next edge
     696                                                                }
    684697                                                        }
    685698                                                        nbcurves++;
     
    689702                                }
    690703                        }
    691                         ISSMASSERT(nbgem && nbe);
    692                         if(step==0){
    693                                 curves = new Curve[nbcurves];
    694                         }
     704                        ISSMASSERT(nb_marked_edges && nbe);
     705                        //allocate if first step
     706                        if(step==0) curves=new Curve[nbcurves];
    695707                }
    696708                for(int i=0;i<nbcurves ;i++){
    697709                        GeometricalEdge * be=curves[i].be, *eqbe=be;
    698710                        //GeometricalEdge * ee=curves[i].ee, *eqee=be;
    699                         curves[i].master=true;
    700711                }
    701712
     
    707718        }
    708719        /*}}}1*/
    709         /*FUNCTION Geometry::Containing{{{1*/
    710         GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
    711                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
    712 
    713                 GeometricalEdge* on =start,* pon=0;
    714                 // walk with the cos on geometry
    715                 int counter=0;
    716                 while(pon != on){ 
    717                         counter++;
    718                         ISSMASSERT(counter<100);
    719                         pon = on;
    720                         R2 A= (*on)[0];
    721                         R2 B= (*on)[1];
    722                         R2 AB = B-A;
    723                         R2 AP = P-A;
    724                         R2 BP = P-B;
    725                         if ( (AB,AP) < 0)
    726                          on = on->Adj[0];
    727                         else if ( (AB,BP)  > 0)
    728                          on = on->Adj[1];
    729                         else
    730                          return on;
    731                 }
    732                 return on;
    733         }
    734         /*}}}1*/
    735720        /*FUNCTION Geometry::Echo {{{1*/
    736 
    737721        void Geometry::Echo(void){
    738722
     
    774758        /*FUNCTION Geometry::MinimalHmin{{{1*/
    775759        double Geometry::MinimalHmin() {
     760                /* coeffIcoor = (2^30-1)/D
     761                 * We cannot go beyond hmin = D/2^30 because of the quadtree
     762                 * hmin is therefore approximately 2/coeffIcoor */
    776763                return 2.0/coefIcoor;
    777764        }/*}}}*/
     
    800787                return c - curves;
    801788        }/*}}}*/
     789        /*FUNCTION Geometry::Containing{{{1*/
     790        GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
     791                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
     792
     793                GeometricalEdge* on =start,* pon=0;
     794                // walk with the cos on geometry
     795                int counter=0;
     796                while(pon != on){ 
     797                        counter++;
     798                        ISSMASSERT(counter<100);
     799                        pon = on;
     800                        R2 A= (*on)[0];
     801                        R2 B= (*on)[1];
     802                        R2 AB = B-A;
     803                        R2 AP = P-A;
     804                        R2 BP = P-B;
     805                        if ( (AB,AP) < 0)
     806                         on = on->Adj[0];
     807                        else if ( (AB,BP)  > 0)
     808                         on = on->Adj[1];
     809                        else
     810                         return on;
     811                }
     812                return on;
     813        }
     814        /*}}}1*/
    802815        /*FUNCTION Geometry::ProjectOnCurve {{{1*/
    803816        GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
  • issm/trunk/src/c/objects/Bamg/Geometry.h

    r5180 r5397  
    5353                        void             ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
    5454                        void             Init(void);
    55                         void             AfterRead();
     55                        void             PostRead();
    5656                        long             GetId(const GeometricalVertex &t) const;
    5757                        long             GetId(const GeometricalVertex *t) const;
  • issm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp

    r5151 r5397  
    1111        /*FUNCTION MatVVP2x2::MatVVP2x2(const Metric M){{{1*/
    1212        MatVVP2x2::MatVVP2x2(const Metric M){
    13                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MatVVP2x2)*/
     13                /*From a metric (a11,a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
    1414
     15                /*Intermediaries*/
    1516                double a11=M.a11,a21=M.a21,a22=M.a22;
    16                 const double eps = 1.e-5;
    17                 double c11 = a11*a11, c22 = a22*a22, c21= a21*a21;
    18                 double b=-a11-a22,c=-c21+a11*a22;
    19                 double delta = b*b - 4 * c ;
    20                 double n2=(c11+c22+c21);
     17                double norm1,norm2,normM;
     18                double delta,b;
    2119
    22                 //if norm(M)<10^30 -> M=zeros(2,2)
    23                 if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0;
     20                /*To get the eigen values, we must solve the following equation:
     21                 *     | a11 - lambda    a21        |
     22                 * det |                            | = 0
     23                 *     | a21             a22-lambda |
     24                 *
     25                 * We have to solve the following polynom:
     26                 *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
    2427
    25                 //if matrix is already diagonal and has a 0 eigen value
    26                 else if (delta < eps*n2){
    27                         lambda1=lambda2=-b/2, v.x=1,v.y=0;
     28                /*Compute polynom determinant*/
     29                b=-a11-a22;
     30                delta=b*b - 4*(a11*a22-a21*a21);
     31
     32
     33                /*Compute norm of M to avoid round off errors*/
     34                normM=a11*a11 + a22*a22 + a21*a21;
     35
     36                /*1: normM too small: eigen values = 0*/
     37                if(normM<1.e-30){
     38                        lambda1=0;
     39                        lambda2=0;
     40                        v.x=1;
     41                        v.y=0;
     42                }
     43                /*2: delta is small -> double root*/
     44                else if (delta < 1.e-5*normM){
     45                        lambda1=-b/2;
     46                        lambda2=-b/2;
     47                        v.x=1;
     48                        v.y=0;
     49                }
     50                /*3: general case -> two roots*/
     51                else{
     52                        delta     = sqrt(delta);
     53                        lambda1   = (-b-delta)/2.0;
     54                        lambda2   = (-b+delta)/2.0;
     55
     56                        /*Now, one must find the eigen vectors. For that we use the following property of the inner product
     57                         *    <Ax,y> = <x,tAy>
     58                         * Here, M'(M-lambda*Id) is symmetrical, which gives:
     59                         *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
     60                         * And we have the following:
     61                         *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
     62                         * We have shown that
     63                         *    Im(M') ⊥ Ker(M')
     64                         *
     65                         * To find the eigen vectors of M, we only have to find two vectors
     66                         * of the image of M' and take their perpendicular as long as they are
     67                         * not 0.
     68                         * To do that, we take the images (1,0) and (0,1):
     69                         *  x1 = (a11 - lambda)      x2 = a21
     70                         *  y1 = a21                 y2 = (a22-lambda)
     71                         *
     72                         * We take the vector that has the larger norm and take its perpendicular.*/
     73
     74                        double norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21;
     75                        double norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
     76
     77                        if (norm2<norm1){
     78                                norm1=sqrt(norm1);
     79                                v.x = - a21/norm1;
     80                                v.y = (a11-lambda1)/norm1;
     81                        }
     82                        else{
     83                                norm2=sqrt(norm2);
     84                                v.x = - (a22-lambda1)/norm2;
     85                                v.y = a21/norm2;
     86                        }
    2887                }
    2988
    30                 //general case: construction of 2 eigen vectors
    31                 else {
    32                         delta     = sqrt(delta);
    33                         lambda1   = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;
    34                         double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;
    35                         double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
    36                         if(s1 < s0)
    37                          s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;
    38                         else
    39                          s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;
    40                 };
    4189        }
    4290        /*}}}1*/
  • issm/trunk/src/c/objects/Bamg/Mesh.cpp

    r5340 r5397  
    2121                if(bamggeom->Edges) {
    2222                        Gh.ReadGeometry(bamggeom,bamgopts);
    23                         Gh.AfterRead();
     23                        Gh.PostRead();
    2424                }
    2525
     
    3232                        printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
    3333                        BuildGeometryFromMesh(bamgopts);
    34                         Gh.AfterRead();
     34                        Gh.PostRead();
    3535                }
    3636
     
    141141                  //double cutoffradian = 10.0/180.0*Pi;
    142142                  BuildGeometryFromMesh(bamgopts);
    143                   Gh.AfterRead();
     143                  Gh.PostRead();
    144144                  SetIntCoor();
    145145                  ReconstructExistingMesh();
     
    300300                BuildGeometryFromMesh();
    301301                if (verbose) printf("Completing geometry\n");
    302                 Gh.AfterRead();
     302                Gh.PostRead();
    303303        }
    304304        /*}}}1*/
     
    53195319                                int jedge=bcurve[icurve]%2;
    53205320
    5321                                 /*Skip if we are on a equi curve (duplicate)*/
    5322                                 if(!Gh.curves[icurve].master) continue;
    5323 
    53245321                                /*Get edge of Bth with index iedge*/
    53255322                                Edge &ei = BTh.edges[iedge];
  • issm/trunk/src/c/objects/Bamg/Metric.cpp

    r5151 r5397  
    269269                        for(i=0;i<2;i++){
    270270                                /*Now, one must find the eigen vectors. For that we use the
    271                                  * following property of the scalare product
     271                                 * following property of the inner product
    272272                                 *    (Ax,b) = transp(b) Ax = transp(x) transp(A) b
    273273                                 *           = (transp(A) b ,x)
Note: See TracChangeset for help on using the changeset viewer.