Changeset 5130


Ignore:
Timestamp:
08/10/10 15:08:05 (15 years ago)
Author:
Mathieu Morlighem
Message:

Renamed some functions and moved routines from h to cpp

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

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Bamg/BamgVertex.h

    r5120 r5130  
    1919                public:
    2020
     21                        /*Fields*/
    2122                        I2        i;                 // integer coordinates
    2223                        R2        r;                 // real coordinates
    … …  
    3536                        };
    3637
    37                         //Operators
     38                        /*Operators*/
    3839                        operator I2() const {return i;}             // Cast operator
    3940                        operator const R2 & () const {return r;}    // Cast operator
    … …  
    4142                        double operator()(R2 x) const { return m(x);} // Get x in the metric m
    4243
    43                         //methods (No constructor and no destructors...)
     44                        /*methods (No constructor and no destructors...)*/
    4445                        double Smoothing(Mesh & ,const Mesh & ,Triangle  * & ,double =1);
    4546                        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);

  • issm/trunk/src/c/objects/Bamg/Curve.cpp

    r3913 r5130  
    1010
    1111        /*Constructors/Destructors*/
     12        /*FUNCTION Curve::Curve(){{{1*/
     13        Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0),master(true){
     14
     15        }
     16        /*}}}*/
    1217
    1318        /*Methods*/
     19        /*FUNCTION Curve::Reverse {{{1*/
     20        void Curve::Reverse() {
     21                /*reverse the direction of the curve */
     22                Exchange(be,ee);
     23                Exchange(kb,ke);
     24        }
     25        /*}}}*/
    1426        /*FUNCTION Curve::Set {{{1*/
    1527        void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew){

  • issm/trunk/src/c/objects/Bamg/Curve.h

    r4997 r5130  
    1919
    2020                        //Methods
    21                         Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {}
    22                         void Reverse() { Exchange(be,ee); Exchange(kb,ke);} //  revese the sens of the curse
     21                        Curve();
     22                        void Reverse(void);
    2323                        void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
    2424        };

  • issm/trunk/src/c/objects/Bamg/Direction.cpp

    r3913 r5130  
    99
    1010        /*Constructors/Destructors*/
     11        /*FUNCTION Direction() {{{1*/
     12        Direction::Direction():
     13                dir(MaxICoor){
     14
     15        }/*}}}*/
    1116        /*FUNCTION Direction(Icoor1 i,Icoor1 j) {{{1*/
    1217        Direction::Direction(Icoor1 i,Icoor1 j) {

  • issm/trunk/src/c/objects/Bamg/Direction.h

    r3913 r5130  
    1313                public:
    1414                        //Methods
    15                         Direction(): dir(MaxICoor){}; //  no direction set
     15                        Direction();
    1616                        Direction(Icoor1 i,Icoor1 j);
    1717                        int sens(Icoor1 i,Icoor1 j);

  • issm/trunk/src/c/objects/Bamg/Edge.cpp

    r5095 r5130  
    3333        }
    3434        /*}}}*/
     35        /*FUNCTION Edge::Renumbering{{{1*/
     36        void Edge::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
     37
     38                if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
     39                if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
     40
     41        }
     42        /*}}}*/
     43        /*FUNCTION Edge::Intersection{{{1*/
     44        int Edge::Intersection(const  Edge & e){
     45
     46                /*some shecks*/
     47                if (!(adj[0]==&e || adj[1]==&e)){ ISSMERROR("Intersection bug"); }
     48                ISSMASSERT(adj[0]==&e || adj[1]==&e);
     49
     50                return adj[0]==&e?0:1;
     51        }
     52        /*}}}*/
    3553
    3654}

  • issm/trunk/src/c/objects/Bamg/Edge.h

    r5124 r5130  
    2828
    2929                        //Methods
    30                         void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
    31                                 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
    32                                 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
    33                         }
    34                         int Intersection(const  Edge & e) const {
    35                                 if (!(adj[0]==&e || adj[1]==&e)){
    36                                         ISSMERROR("Intersection bug");
    37                                 }
    38                                 if (adj[0]!=&e && adj[1]!=&e){
    39                                         ISSMERROR("adj[0]!=&e && adj[1]!=&e");
    40                                 }
    41                                 return adj[0]==&e ? 0 : 1;
    42                         }
     30                        void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu);
     31                        int  Intersection(const  Edge & e);
    4332                        void Set(const Mesh &,long,Mesh &);
    4433                        void Echo(void);

  • issm/trunk/src/c/objects/Bamg/GeometricalEdge.h

    r3913 r5130  
    3030                        double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
    3131                        int    Tg(int i) const{return i==0 ? TgA() : TgB();}
    32                         int    Cracked() const    {return flag & 1;  }
    33                         int    TgA()const         {return flag & 4; }
    34                         int    TgB()const         {return flag & 8;  }
    35                         int    Mark()const        {return flag & 16; }
    36                         int    Required()         {return flag & 64; }
    37                         void   SetCracked()     { flag |= 1;  }
    38                         void   SetTgA()         { flag |=4;   }
    39                         void   SetTgB()         { flag |=8;   }
    40                         void   SetMark()        { flag |=16;  }
     32                        int    Cracked() const{return flag & 1;}
     33                        int    TgA()const{return flag & 4; }
     34                        int    TgB()const {return flag & 8;}
     35                        int    Mark()const {return flag & 16;}
     36                        int    Required() {return flag & 64;}
     37                        void   SetCracked()     { flag |= 1;}
     38                        void   SetTgA()         { flag |=4; }
     39                        void   SetTgB()         { flag |=8; }
     40                        void   SetMark()        { flag |=16;}
    4141                        void   SetUnMark()      { flag &= 1007 /* 1023-16*/;}
    4242                        void   SetRequired()    { flag |= 64; }

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

    r5124 r5130  
    1313
    1414        /*Constructors/Destructors*/
     15        /*FUNCTION Geometry::Geometry(){{{1*/
    1516        Geometry::Geometry(){
    1617                Init();
    1718        }
     19        /*}}}*/
     20        /*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{1*/
    1821        Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
    1922                Init();
    … …  
    2124                AfterRead();
    2225        }
    23         /*FUNCTION  Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/
     26        /*}}}*/
     27        /*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/
    2428        Geometry::Geometry(const Geometry & Gh) {
    2529                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/
    … …  
    405409
    406410        /*Methods*/
    407         /*FUNCTION  Geometry::AfterRead(){{{1*/
     411        /*FUNCTION Geometry::AfterRead(){{{1*/
    408412        void Geometry::AfterRead(){
    409413                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
    … …  
    719723        }
    720724        /*}}}1*/
    721         /*FUNCTION  Geometry::Containing{{{1*/
     725        /*FUNCTION Geometry::Containing{{{1*/
    722726        GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
    723727                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
    … …  
    745749        }
    746750        /*}}}1*/
    747         /*FUNCTION  Geometry::Echo {{{1*/
     751        /*FUNCTION Geometry::Echo {{{1*/
    748752
    749753        void Geometry::Echo(void){
    … …  
    768772        }
    769773        /*}}}*/
    770         /*FUNCTION  Geometry::Init{{{1*/
     774        /*FUNCTION Geometry::Init{{{1*/
    771775        void Geometry::Init(void){
    772776                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
    … …  
    786790        }
    787791        /*}}}1*/
    788         /*FUNCTION  Geometry::ProjectOnCurve {{{1*/
     792        /*FUNCTION Geometry::MinimalHmin{{{1*/
     793        double Geometry::MinimalHmin() {
     794                return 2.0/coefIcoor;
     795        }/*}}}*/
     796        /*FUNCTION Geometry::MaximalHmax{{{1*/
     797        double Geometry::MaximalHmax() {
     798                return Max(pmax.x-pmin.x,pmax.y-pmin.y);
     799        }/*}}}*/
     800        /*FUNCTION Geometry::Number(const GeometricalVertex &t){{{1*/
     801        long Geometry::Number(const GeometricalVertex & t) const  {
     802                return &t - vertices;
     803        }/*}}}*/
     804        /*FUNCTION Geometry::Number(const GeometricalVertex * t){{{1*/
     805        long Geometry::Number(const GeometricalVertex * t) const  {
     806                return t - vertices;
     807        }/*}}}*/
     808        /*FUNCTION Geometry::Number(const GeometricalEdge & t){{{1*/
     809        long Geometry::Number(const GeometricalEdge & t) const  {
     810                return &t - edges;
     811        }/*}}}*/
     812        /*FUNCTION Geometry::Number(const GeometricalEdge * t){{{1*/
     813        long Geometry::Number(const GeometricalEdge * t) const  {
     814                return t - edges;
     815        }/*}}}*/
     816        /*FUNCTION Geometry::Number(const Curve * c){{{1*/
     817        long Geometry::Number(const Curve * c) const  {
     818                return c - curves;
     819        }/*}}}*/
     820        /*FUNCTION Geometry::ProjectOnCurve {{{1*/
    789821        GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
    790822                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
    … …  
    844876                                        printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n");
    845877                                        ISSMERROR("see above");
    846                                   }
     878                                }
    847879                                NbTry++;
    848880                                goto retry;
    … …  
    920952        }
    921953        /*}}}1*/
    922 
     954        /*FUNCTION Geometry::toI2{{{1*/
     955        I2 Geometry::toI2(const R2 & P) const {
     956                return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
     957                                        ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
     958        }/*}}}*/
     959        /*FUNCTION Geometry::UnMarkEdges{{{1*/
     960        void Geometry::UnMarkEdges() {
     961                for (int i=0;i<nbe;i++) edges[i].SetUnMark();
     962        }/*}}}*/
    923963}

  • issm/trunk/src/c/objects/Bamg/Geometry.h

    r5120 r5130  
    4242
    4343                        //Operators
    44                         const GeometricalVertex & operator[]  (long i) const { return vertices[i];};
    45                         GeometricalVertex & operator[](long i) { return vertices[i];};
    46                         const  GeometricalEdge & operator()  (long i) const { return edges[i];};
    47                         GeometricalEdge & operator()(long i) { return edges[i];};
     44                        const GeometricalVertex &operator[](long i) const { return vertices[i]; };
     45                        GeometricalVertex       &operator[](long i) { return vertices[i];       };
     46                        const GeometricalEdge   &operator()(long i) const { return edges[i];    };
     47                        GeometricalEdge         &operator()(long  i) { return edges[i];                };
    4848
    4949                        //Methods
    50                         void Echo();
    51                         I2 toI2(const R2 & P) const {
    52                                 return  I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
    53                                                         ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
    54                         }
    55                         double MinimalHmin() {return 2.0/coefIcoor;}
    56                         double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
    57                         void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
    58                         void Init(void);
    59                         void AfterRead();
    60                         long Number(const GeometricalVertex & t) const  { return &t - vertices;}
    61                         long Number(const GeometricalVertex * t) const  { return t - vertices;}
    62                         long Number(const GeometricalEdge & t) const  { return &t - edges;}
    63                         long Number(const GeometricalEdge * t) const  { return t - edges;}
    64                         long Number(const Curve * c) const  { return c - curves;}
    65                         void UnMarkEdges() {for (int i=0;i<nbe;i++) edges[i].SetUnMark();}
    66                         GeometricalEdge *  ProjectOnCurve(const Edge & ,double,BamgVertex &,VertexOnGeom &) const ;
    67                         GeometricalEdge *  Containing(const R2 P,  GeometricalEdge * start) const;
    68                         void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
     50                        void             Echo();
     51                        I2               toI2(const R2 &P) const;
     52                        double           MinimalHmin();
     53                        double           MaximalHmax();
     54                        void             ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
     55                        void             Init(void);
     56                        void             AfterRead();
     57                        long             Number(const GeometricalVertex &t) const;
     58                        long             Number(const GeometricalVertex *t) const;
     59                        long             Number(const GeometricalEdge &t) const;
     60                        long             Number(const GeometricalEdge *t) const;
     61                        long             Number(const Curve *c) const;
     62                        void             UnMarkEdges();
     63                        GeometricalEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const;
     64                        GeometricalEdge *Containing(const R2 P, GeometricalEdge *start) const;
     65                        void             WriteGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
    6966        };
    7067       

  • issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp

    r5124 r5130  
    88namespace bamg {
    99
     10        /*Constructors Destructors*/
     11        /*FUNCTION ListofIntersectionTriangles::ListofIntersectionTriangles{{{1*/
     12        ListofIntersectionTriangles::ListofIntersectionTriangles(int n,int m)
     13          : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
     14          NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
     15          }
     16        /*}}}*/
     17        /*FUNCTION ListofIntersectionTriangles::~ListofIntersectionTriangles{{{1*/
     18        ListofIntersectionTriangles::~ListofIntersectionTriangles(){
     19                if (lIntTria) delete [] lIntTria,lIntTria=0;
     20                if (lSegsI) delete [] lSegsI,lSegsI=0;
     21        }
     22        /*}}}*/
     23
    1024        /*Methods*/
     25        /*FUNCTION ListofIntersectionTriangles::Init{{{1*/
     26        void ListofIntersectionTriangles::Init(void){
     27                state=0;
     28                len=0;
     29                Size=0;
     30        }
     31        /*}}}*/
     32        /*FUNCTION ListofIntersectionTriangles::Length{{{1*/
     33        double  ListofIntersectionTriangles::Length(){
     34                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Length)*/
     35
     36                // computation of the length
     37
     38                // check Size
     39                if (Size<=0){
     40                        ISSMERROR("Size<=0");
     41                }
     42
     43                Metric Mx,My;
     44                int ii,jj;
     45                R2 x,y,xy;
     46
     47                SegInterpolation* SegI=lSegsI;
     48                SegI=lSegsI;
     49                lSegsI[NbSeg].last=Size;
     50                int EndSeg=Size;
     51
     52                y = lIntTria[0].x;
     53                double sxy, s = 0;
     54                lIntTria[0].s =0;
     55                SegI->lBegin=s;
     56
     57                for (jj=0,ii=1;ii<Size;jj=ii++) { 
     58                        // seg jj,ii
     59                        x  = y;
     60                        y  = lIntTria[ii].x;
     61                        xy = y-x;
     62                        Mx = lIntTria[ii].m;
     63                        My = lIntTria[jj].m;
     64                        sxy= LengthInterpole(Mx,My,xy);
     65                        s += sxy;
     66                        lIntTria[ii].s = s;
     67                        if (ii == EndSeg){
     68                                SegI->lEnd=s;
     69                                SegI++;
     70                                EndSeg=SegI->last;
     71                                SegI->lBegin=s;
     72                        }
     73                }
     74                len = s;
     75                SegI->lEnd=s;
     76
     77                return s;
     78        }
     79        /*}}}1*/
     80        /*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {{{1*/
     81        int  ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {
     82                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
     83
     84                register int n;
     85                R2 x(0,0);
     86                if ( d0) x =      (*tt)[0].r * d0;
     87                if ( d1) x = x +  (*tt)[1].r * d1;
     88                if ( d2) x = x +  (*tt)[2].r * d2;
     89                // newer add same point
     90                if(!Size ||  Norme2_2(lIntTria[Size-1].x-x)) {
     91                        if (Size==MaxSize) ReShape();
     92                        lIntTria[Size].t=tt;
     93                        lIntTria[Size].bary[0]=d0;
     94                        lIntTria[Size].bary[1]=d1;
     95                        lIntTria[Size].bary[2]=d2;
     96                        lIntTria[Size].x = x;
     97                        Metric m0,m1,m2;
     98                        register BamgVertex * v;
     99                        if ((v=(*tt)(0))) m0    = v->m;
     100                        if ((v=(*tt)(1))) m1    = v->m;
     101                        if ((v=(*tt)(2))) m2    = v->m;
     102                        lIntTria[Size].m =  Metric(lIntTria[Size].bary,m0,m1,m2);
     103                        n=Size++;}
     104                else n=Size-1;
     105                return n;
     106        }
     107        /*}}}1*/
     108        /*FUNCTION ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm){{{1*/
     109        int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
     110                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
     111
     112                register int n;
     113                if(!Size ||  Norme2_2(lIntTria[Size-1].x-A)) {
     114                        if (Size==MaxSize) ReShape();
     115                        lIntTria[Size].t=0;
     116                        lIntTria[Size].x=A;
     117                        lIntTria[Size].m=mm;
     118                        n=Size++;
     119                }
     120                else  n=Size-1;
     121                return  n;
     122        }
     123        /*}}}1*/
     124        /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
     125        long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){
     126                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
     127
     128                //If length<1.5, do nothing
     129                double s=Length();
     130                if (s<1.5) return 0;
     131
     132                const long nbvold=nbv;
     133                int ii = 1 ;
     134                R2 y,x;
     135                Metric My,Mx ;
     136                double sx =0,sy;
     137                int nbi=Max(2,(int) (s+0.5));
     138                double sint=s/nbi;
     139                double si  =sint;
     140
     141                int EndSeg=Size;
     142                SegInterpolation* SegI=NULL;
     143                if (NbSeg) SegI=lSegsI,EndSeg=SegI->last;
     144
     145                for (int k=1;k<nbi;k++){
     146                        while ((ii < Size) && ( lIntTria[ii].s <= si )){
     147                                if (ii++ == EndSeg){
     148                                        SegI++;
     149                                        EndSeg=SegI->last;
     150                                }
     151                        }
     152
     153                        int ii1=ii-1;
     154                        x  =lIntTria[ii1].x;
     155                        sx =lIntTria[ii1].s;
     156                        Metric Mx=lIntTria[ii1].m;
     157                        y  =lIntTria[ii].x;
     158                        sy =lIntTria[ii].s;
     159                        Metric My=lIntTria[ii].m;
     160                        double lxy = sy-sx;
     161                        double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
     162
     163                        R2 C;
     164                        double cx = 1-cy;
     165                        C = SegI ? SegI->F(si): x * cx + y *cy;
     166                        //C.Echo();
     167                        //x.Echo();
     168                        //y.Echo();
     169                        //printf("cx = %g, cy=%g\n",cx,cy);
     170
     171                        si += sint;
     172                        if ( nbv<maxnbv) {
     173                                vertices[nbv].r = C;
     174                                vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
     175                        }
     176                        else return nbv-nbvold;
     177                  }
     178                return nbv-nbvold;
     179        }
     180        /*}}}1*/
     181        /*FUNCTION ListofIntersectionTriangles::NewSubSeg{{{1*/
     182        void  ListofIntersectionTriangles::NewSubSeg(GeometricalEdge *e,double s0,double s1){
     183                long int verbosity=0;
     184                if (NbSeg>=MaxNbSeg) {
     185                        int mneo= MaxNbSeg;
     186                        MaxNbSeg *= 2;
     187                        if (verbosity>3){
     188                                printf("   reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
     189                        }
     190                        SegInterpolation * lEn =  new SegInterpolation[MaxNbSeg];
     191                        if (!lSegsI || NbSeg>=MaxNbSeg){
     192                                ISSMERROR("!lSegsI || NbSeg>=MaxNbSeg");
     193                        }
     194                        for (int i=0;i< NbSeg;i++)
     195                         lEn[i] = lSegsI[MaxNbSeg]; // copy old to new           
     196                        delete []  lSegsI; // remove old
     197                        lSegsI = lEn;       
     198                }
     199                if (NbSeg) lSegsI[NbSeg-1].last=Size;
     200                lSegsI[NbSeg].e=e;
     201                lSegsI[NbSeg].sBegin=s0;
     202                lSegsI[NbSeg].sEnd=s1;     
     203                NbSeg++;           
     204        }
     205        /*}}}*/
     206        /*FUNCTION ListofIntersectionTriangles::ReShape{{{1*/
     207        void ListofIntersectionTriangles::ReShape(){
     208
     209                register int newsize = MaxSize*2;
     210                IntersectionTriangles* nw = new IntersectionTriangles[newsize];
     211                ISSMASSERT(nw);
     212
     213                // recopy
     214                for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];       
     215                long int verbosity=0;
     216                if(verbosity>3) printf("   ListofIntersectionTriangles  ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg);
     217                MaxSize = newsize;
     218                delete [] lIntTria;// remove old
     219                lIntTria = nw; // copy pointer
     220        }
     221        /*}}}*/
    11222        /*FUNCTION ListofIntersectionTriangles::SplitEdge{{{1*/
    12223        void ListofIntersectionTriangles::SplitEdge(const Mesh & Bh, const R2 &A,const R2  &B,int nbegin) {
    … …  
    33244                         ifirst = ilast;} 
    34245                else {// not optimisation
    35                         init();
     246                        Init();
    36247                        t=tbegin = Bh.TriangleFindFromCoord(a,deta);
    37248                        if( t->det>=0)
    … …  
    180391                                }
    181392
    182                                                 k = OppositeVertex[ocut];
    183 
    184                                                 Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
    185 
    186 
    187                                                 if (detbij >= 0) { //we find the triangle contening b
    188                                                         dt[0]=bamg::det((*t)[1],(*t)[2],b);
    189                                                         dt[1]=bamg::det((*t)[2],(*t)[0],b);
    190                                                         dt[2]=bamg::det((*t)[0],(*t)[1],b);
    191                                                         double dd = t->det;
    192                                                         NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);     
    193                                                         return ;}
    194                                                 else { // next triangle by  adjacent by edge ocut
    195                                                         deti = dt[i];
    196                                                         detj = dt[j];
    197                                                         double dij = detj-deti;
    198                                                         ba[i] =  detj/dij;
    199                                                         ba[j] = -deti/dij;
    200                                                         ba[3-i-j ] = 0;
    201                                                         ilast=NewItem(t, ba[0],ba[1],ba[2]);     
    202 
    203                                                         TriangleAdjacent ta =t->Adj(ocut);
    204                                                         t = ta;
    205                                                         iedge= ta;
    206                                                         if (t->det <= 0)  {
    207                                                                 double ba,bb;
    208                                                                 TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
    209                                                                 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
    210                                                                 return;
    211                                                         }
    212                                                 }// we  go outside of omega
     393                                k = OppositeVertex[ocut];
     394
     395                                Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
     396
     397
     398                                if (detbij >= 0) { //we find the triangle contening b
     399                                        dt[0]=bamg::det((*t)[1],(*t)[2],b);
     400                                        dt[1]=bamg::det((*t)[2],(*t)[0],b);
     401                                        dt[2]=bamg::det((*t)[0],(*t)[1],b);
     402                                        double dd = t->det;
     403                                        NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);     
     404                                        return ;}
     405                                else { // next triangle by  adjacent by edge ocut
     406                                        deti = dt[i];
     407                                        detj = dt[j];
     408                                        double dij = detj-deti;
     409                                        ba[i] =  detj/dij;
     410                                        ba[j] = -deti/dij;
     411                                        ba[3-i-j ] = 0;
     412                                        ilast=NewItem(t, ba[0],ba[1],ba[2]);     
     413
     414                                        TriangleAdjacent ta =t->Adj(ocut);
     415                                        t = ta;
     416                                        iedge= ta;
     417                                        if (t->det <= 0)  {
     418                                                double ba,bb;
     419                                                TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
     420                                                NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
     421                                                return;
     422                                        }
     423                                }// we  go outside of omega
    213424                } // for(;;)
    214425        }
    215426        /*}}}1*/
    216         /*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {{{1*/
    217         int  ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {
    218                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
    219 
    220                 register int n;
    221                 R2 x(0,0);
    222                 if ( d0) x =      (*tt)[0].r * d0;
    223                 if ( d1) x = x +  (*tt)[1].r * d1;
    224                 if ( d2) x = x +  (*tt)[2].r * d2;
    225                 // newer add same point
    226                 if(!Size ||  Norme2_2(lIntTria[Size-1].x-x)) {
    227                         if (Size==MaxSize) ReShape();
    228                         lIntTria[Size].t=tt;
    229                         lIntTria[Size].bary[0]=d0;
    230                         lIntTria[Size].bary[1]=d1;
    231                         lIntTria[Size].bary[2]=d2;
    232                         lIntTria[Size].x = x;
    233                         Metric m0,m1,m2;
    234                         register BamgVertex * v;
    235                         if ((v=(*tt)(0))) m0    = v->m;
    236                         if ((v=(*tt)(1))) m1    = v->m;
    237                         if ((v=(*tt)(2))) m2    = v->m;
    238                         lIntTria[Size].m =  Metric(lIntTria[Size].bary,m0,m1,m2);
    239                         n=Size++;}
    240                 else n=Size-1;
    241                 return n;
    242         }
    243         /*}}}1*/
    244         /*FUNCTION ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm){{{1*/
    245         int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
    246                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
    247 
    248                 register int n;
    249                 if(!Size ||  Norme2_2(lIntTria[Size-1].x-A)) {
    250                         if (Size==MaxSize) ReShape();
    251                         lIntTria[Size].t=0;
    252                         lIntTria[Size].x=A;
    253                         lIntTria[Size].m=mm;
    254                         n=Size++;
    255                 }
    256                 else  n=Size-1;
    257                 return  n;
    258         }
    259         /*}}}1*/
    260         /*FUNCTION ListofIntersectionTriangles::Length{{{1*/
    261         double  ListofIntersectionTriangles::Length(){
    262                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Length)*/
    263 
    264                 // computation of the length
    265 
    266                 // check Size
    267                 if (Size<=0){
    268                         ISSMERROR("Size<=0");
    269                 }
    270 
    271                 Metric Mx,My;
    272                 int ii,jj;
    273                 R2 x,y,xy;
    274 
    275                 SegInterpolation* SegI=lSegsI;
    276                 SegI=lSegsI;
    277                 lSegsI[NbSeg].last=Size;
    278                 int EndSeg=Size;
    279 
    280                 y = lIntTria[0].x;
    281                 double sxy, s = 0;
    282                 lIntTria[0].s =0;
    283                 SegI->lBegin=s;
    284 
    285                 for (jj=0,ii=1;ii<Size;jj=ii++) { 
    286                         // seg jj,ii
    287                         x  = y;
    288                         y  = lIntTria[ii].x;
    289                         xy = y-x;
    290                         Mx = lIntTria[ii].m;
    291                         My = lIntTria[jj].m;
    292                         sxy= LengthInterpole(Mx,My,xy);
    293                         s += sxy;
    294                         lIntTria[ii].s = s;
    295                         if (ii == EndSeg){
    296                                 SegI->lEnd=s;
    297                                 SegI++;
    298                                 EndSeg=SegI->last;
    299                                 SegI->lBegin=s;
    300                         }
    301                 }
    302                 len = s;
    303                 SegI->lEnd=s;
    304 
    305                 return s;
    306         }
    307         /*}}}1*/
    308         /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
    309         long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){
    310                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
    311 
    312                 //If length<1.5, do nothing
    313                 double s=Length();
    314                 if (s<1.5) return 0;
    315 
    316                 const long nbvold=nbv;
    317                 int ii = 1 ;
    318                 R2 y,x;
    319                 Metric My,Mx ;
    320                 double sx =0,sy;
    321                 int nbi=Max(2,(int) (s+0.5));
    322                 double sint=s/nbi;
    323                 double si  =sint;
    324 
    325                 int EndSeg=Size;
    326                 SegInterpolation* SegI=NULL;
    327                 if (NbSeg) SegI=lSegsI,EndSeg=SegI->last;
    328 
    329                 for (int k=1;k<nbi;k++){
    330                         while ((ii < Size) && ( lIntTria[ii].s <= si )){
    331                                 if (ii++ == EndSeg){
    332                                         SegI++;
    333                                         EndSeg=SegI->last;
    334                                 }
    335                         }
    336 
    337                         int ii1=ii-1;
    338                         x  =lIntTria[ii1].x;
    339                         sx =lIntTria[ii1].s;
    340                         Metric Mx=lIntTria[ii1].m;
    341                         y  =lIntTria[ii].x;
    342                         sy =lIntTria[ii].s;
    343                         Metric My=lIntTria[ii].m;
    344                         double lxy = sy-sx;
    345                         double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
    346 
    347                         R2 C;
    348                         double cx = 1-cy;
    349                         C = SegI ? SegI->F(si): x * cx + y *cy;
    350                         //C.Echo();
    351                         //x.Echo();
    352                         //y.Echo();
    353                         //printf("cx = %g, cy=%g\n",cx,cy);
    354 
    355                         si += sint;
    356                         if ( nbv<maxnbv) {
    357                                 vertices[nbv].r = C;
    358                                 vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
    359                         }
    360                         else return nbv-nbvold;
    361                   }
    362                 return nbv-nbvold;
    363         }
    364         /*}}}1*/
    365427
    366428}

  • issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h

    r5120 r5130  
    1111
    1212                class IntersectionTriangles {
     13
    1314                        public:
    14                                 Triangle* t;
    15                                 double  bary[3];  // use if t != 0
    16                                 R2 x;
    17                                 Metric m;
    18                                 double s; // curvilinear coordinate
    19                                 double sp;// length of the previous segment in m
    20                                 double sn;// length of the next segment in m
     15                                Triangle *t;
     16                                double    bary[3];   // use if t != 0
     17                                R2        x;
     18                                Metric    m;
     19                                double    s;         // curvilinear coordinate
     20                                double    sp;        // length of the previous segment in m
     21                                double    sn;        // length of the next segment in m
    2122                };
    2223
    2324                class SegInterpolation {
     25
    2426                        public:
    25                                 GeometricalEdge * e;
    26                                 double sBegin,sEnd; // abscisse of the seg on edge parameter
    27                                 double lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
    28                                 int last;// last index in ListofIntersectionTriangles for this Sub seg of edge
     27                                GeometricalEdge *e;
     28                                double           sBegin  ,sEnd; // abscisse of the seg on edge parameter
     29                                double           lBegin  ,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
     30                                int              last;          // last index in ListofIntersectionTriangles for this Sub seg of edge
    2931
    3032                                //Methods
    … …  
    4042                public:
    4143
    42                         int MaxSize;
    43                         int Size;
    44                         double len;
    45                         int state;
    46                         IntersectionTriangles * lIntTria;
    47                         int NbSeg;
    48                         int MaxNbSeg;
    49                         SegInterpolation * lSegsI;
     44                        int                    MaxSize;
     45                        int                    Size;
     46                        double                 len;
     47                        int                    state;
     48                        IntersectionTriangles *lIntTria;
     49                        int                    NbSeg;
     50                        int                    MaxNbSeg;
     51                        SegInterpolation      *lSegsI;
    5052
    5153                        //Constructors/Destructors
    52                         ListofIntersectionTriangles(int n=256,int m=16)
    53                           : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
    54                           NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
    55                                   long int verbosity=0;
    56                                   if (verbosity>9) printf("   construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
    57                         }
    58                         ~ListofIntersectionTriangles(){
    59                                 if (lIntTria) delete [] lIntTria,lIntTria=0;
    60                                 if (lSegsI) delete [] lSegsI,lSegsI=0;
    61                         }
     54                        ListofIntersectionTriangles(int n=256,int m=16);
     55                        ~ListofIntersectionTriangles();
    6256
    6357                        //Operators
    … …  
    6660
    6761                        //Methods
    68                         void  init(){state=0;len=0;Size=0;}
    69                         int   NewItem(Triangle * tt,double d0,double d1,double d2);
    70                         int   NewItem(R2,const Metric & );
    71                         void  SplitEdge(const Mesh & ,const R2 &,const R2  &,int nbegin=0);
    72                         double Length();
    73                         long  NewPoints(BamgVertex *,long & nbv,long maxnbv);
    74                         void  NewSubSeg(GeometricalEdge *e,double s0,double s1){
    75                                 long int verbosity=0;
    76                                 if (NbSeg>=MaxNbSeg) {
    77                                         int mneo= MaxNbSeg;
    78                                         MaxNbSeg *= 2;
    79                                         if (verbosity>3){
    80                                                 printf("   reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
    81                                         }
    82                                         SegInterpolation * lEn =  new SegInterpolation[MaxNbSeg];
    83                                         if (!lSegsI || NbSeg>=MaxNbSeg){
    84                                                 ISSMERROR("!lSegsI || NbSeg>=MaxNbSeg");
    85                                         }
    86                                         for (int i=0;i< NbSeg;i++)
    87                                          lEn[i] = lSegsI[MaxNbSeg]; // copy old to new           
    88                                         delete []  lSegsI; // remove old
    89                                         lSegsI = lEn;       
    90                                 }
    91                                 if (NbSeg) lSegsI[NbSeg-1].last=Size;
    92                                 lSegsI[NbSeg].e=e;
    93                                 lSegsI[NbSeg].sBegin=s0;
    94                                 lSegsI[NbSeg].sEnd=s1;     
    95                                 NbSeg++;           
    96                         }
    97                         void ReShape(){
    98                                 register int newsize = MaxSize*2;
    99                                 IntersectionTriangles* nw = new IntersectionTriangles[newsize];
    100                                 if (!nw) ISSMERROR("!nw");
    101                                 // recopy
    102                                 for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];       
    103                                 long int verbosity=0;
    104                                 if(verbosity>3) printf("   ListofIntersectionTriangles  ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg);
    105                                 MaxSize = newsize;
    106                                 delete [] lIntTria;// remove old
    107                                 lIntTria = nw; // copy pointer
    108                         }
     62                        void   Init();
     63                        int    NewItem(Triangle *tt,double d0,double d1,double d2);
     64                        int    NewItem(R2 ,const Metric &);
     65                        void   SplitEdge(const Mesh &,const R2 &,const R2 &,int nbegin=0);
     66                        double Length();
     67                        long   NewPoints(BamgVertex *,long &nbv,long maxnbv);
     68                        void   NewSubSeg(GeometricalEdge *e,double s0,double s1);
     69                        void   ReShape();
    10970        };
    11071

  • issm/trunk/src/c/objects/Bamg/Mesh.cpp

    r5124 r5130  
    219219
    220220          }
     221        /*}}}1*/
     222        /*FUNCTION Mesh::Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices){{{1*/
     223        Mesh::Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices) :Gh(BT.Gh),BTh(BT) {
     224                TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);
     225        }
     226        /*}}}1*/
     227        /*FUNCTION Mesh::Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts){{{1*/
     228        Mesh::Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
     229                TriangulateFromGeom0(maxnbv,bamgopts);
     230        }
    221231        /*}}}1*/
    222232        /*FUNCTION Mesh::~Mesh(){{{1*/

  • issm/trunk/src/c/objects/Bamg/Mesh.h

    r5124 r5130  
    6262                        Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator
    6363                        Mesh(const Mesh &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
    64                         Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
    65                                 try {TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);}
    66                                 catch(...) { this->~Mesh(); throw; }
    67                         }
    68                         Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){
    69                                 try { TriangulateFromGeom0(maxnbv,bamgopts);}
    70                                 catch(...) { this->~Mesh(); throw; }
    71                         }
     64                        Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1);
     65                        Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts);
    7266                        ~Mesh();
    7367
    7468                        //Operators
    75                         const BamgVertex & operator[]  (long i) const { return vertices[i];};
    76                         BamgVertex & operator[](long i) { return vertices[i];};
    77                         const Triangle & operator()  (long i) const { return triangles[i];};
    78                         Triangle & operator()(long i) { return triangles[i];};
     69                        const BamgVertex &operator[](long i) const { return vertices[i];  };
     70                        BamgVertex       &operator[](long i) { return vertices[i];        };
     71                        const Triangle   &operator()(long i) const { return triangles[i]; };
     72                        Triangle         &operator()(long  i) { return triangles[i];             };
    7973
    8074                        //Methods
    … …  
    159153                        void TriangulateFromGeom0(long maxnbv,BamgOpts* bamgopts);// the real constructor mesh generator
    160154                        void Init(long);
    161 
    162155        };
    163156

  • issm/trunk/src/c/objects/Bamg/include.h

    r3913 r5130  
    66#define  _INCLUDE2_H_
    77
    8 
    98#include "./macros.h"
    109#include "./typedefs.h"
    1110
    12 
    1311#endif //ifndef _INCLUDE2_H_
    14 

  • issm/trunk/src/c/objects/Bamg/macros.h

    r3913 r5130  
    1919        static const short NextVertex[3] = {1,2,0};
    2020        static const short PreviousVertex[3] = {2,0,1};
    21 #if LONG_BIT > 63
     21        #if LONG_BIT > 63
    2222        const  Icoor1 MaxICoor   = 1073741823; // 2^30-1 =111...111 (29 times)
    23 #else
     23        #else
    2424        const  Icoor1 MaxICoor   = 8388608;    // 2^23
    25 #endif
     25        #endif
    2626        const  Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ;
    2727}
Note: See TracChangeset for help on using the changeset viewer.