Changeset 2862


Ignore:
Timestamp:
01/15/10 14:50:32 (15 years ago)
Author:
Mathieu Morlighem
Message:

removed useless code

Location:
issm/trunk/src/c/Bamgx
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r2861 r2862  
    66namespace bamg {
    77
    8 typedef P2<double,double> D2;
    9 typedef P2xP2<double,double> D2xD2;
     8        typedef P2<double,double> D2;
     9        typedef P2xP2<double,double> D2xD2;
    1010
    11 class  MetricAnIso;
    12 class MatVVP2x2;
    13 class MetricIso;
     11        class  MetricAnIso;
     12        class MatVVP2x2;
     13        class MetricIso;
    1414
    15 typedef TYPEMETRIX Metric;
     15        typedef TYPEMETRIX Metric;
    1616
    1717
    18 class MetricIso{
    19   friend class MatVVP2x2;
    20    Real4 h;
    21 public:
    22   MetricIso(Real4 a): h(a){}
    23   MetricIso(const MetricAnIso M);// {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
    24   MetricIso(Real8 a11,Real8 a21,Real8 a22);// {*this=MetricAnIso(a11,a21,a22));}
    25   MetricIso(): h(1) {}; //
    26   MetricIso(const Real8  a[3],const  MetricIso m0,
    27            const  MetricIso m1,const  MetricIso m2 )
    28     : h(hinterpole
    29         ? (a[0]*m0.h+a[1]*m1.h+a[2]*m2.h)
    30         : 1/sqrt(a[0]/(m0.h*m0.h)+a[1]/(m1.h*m1.h)+a[2]/(m2.h*m2.h))) {}
    31     MetricIso(const Real8  a,const  MetricIso ma,
    32                   const Real8  b,const  MetricIso mb)
    33     : h(hinterpole
    34         ? a*ma.h+b*mb.h
    35         :1/sqrt(a/(ma.h*ma.h)+b/(mb.h*mb.h))) {}
    36   R2 Orthogonal(const R2 A)const {return R2(-h*A.y,h*A.x);}
    37   R2 Orthogonal(const I2 A)const {return R2(-h*A.y,h*A.x);}
    38 //  D2 Orthogonal(const D2 A)const {return D2(-h*A.y,h*A.x);}
    39   Real8 operator()(R2 x) const { return sqrt((x,x))/h;};
    40   Real8 operator()(R2 x,R2 y) const { return ((x,y))/(h*h);};
    41   int  IntersectWith(MetricIso M) {int r=0;if (M.h<h) r=1,h=M.h;return r;}
    42   MetricIso operator*(Real8 c) const {return  MetricIso(h/c);}
    43   MetricIso operator/(Real8 c) const {return  MetricIso(h*c);}
    44   Real8 det() const {return 1./(h*h*h*h);}   
    45   operator D2xD2(){ return D2xD2(1/(h*h),0.,0.,1/(h*h));}
    46   void     Box(Real4 & hx,Real4 & hy) { hx=h,hy=h;}
    47 };
     18        class MetricIso{
     19                friend class MatVVP2x2;
     20                Real4 h;
     21                public:
     22                MetricIso(Real4 a): h(a){}
     23                MetricIso(const MetricAnIso M);// {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
     24                MetricIso(Real8 a11,Real8 a21,Real8 a22);// {*this=MetricAnIso(a11,a21,a22));}
     25                MetricIso(): h(1) {}; //
     26                MetricIso(const Real8  a[3],const  MetricIso m0,
     27                                        const  MetricIso m1,const  MetricIso m2 )
     28                  : h(hinterpole
     29                                          ? (a[0]*m0.h+a[1]*m1.h+a[2]*m2.h)
     30                                          : 1/sqrt(a[0]/(m0.h*m0.h)+a[1]/(m1.h*m1.h)+a[2]/(m2.h*m2.h))) {}
     31                MetricIso(const Real8  a,const  MetricIso ma,
     32                                        const Real8  b,const  MetricIso mb)
     33                  : h(hinterpole
     34                                          ? a*ma.h+b*mb.h
     35                                          :1/sqrt(a/(ma.h*ma.h)+b/(mb.h*mb.h))) {}
     36                R2 Orthogonal(const R2 A)const {return R2(-h*A.y,h*A.x);}
     37                R2 Orthogonal(const I2 A)const {return R2(-h*A.y,h*A.x);}
     38                //  D2 Orthogonal(const D2 A)const {return D2(-h*A.y,h*A.x);}
     39                Real8 operator()(R2 x) const { return sqrt((x,x))/h;};
     40                Real8 operator()(R2 x,R2 y) const { return ((x,y))/(h*h);};
     41                int  IntersectWith(MetricIso M) {int r=0;if (M.h<h) r=1,h=M.h;return r;}
     42                MetricIso operator*(Real8 c) const {return  MetricIso(h/c);}
     43                MetricIso operator/(Real8 c) const {return  MetricIso(h*c);}
     44                Real8 det() const {return 1./(h*h*h*h);}   
     45                operator D2xD2(){ return D2xD2(1/(h*h),0.,0.,1/(h*h));}
     46                void     Box(Real4 & hx,Real4 & hy) { hx=h,hy=h;}
     47        };
    4848
    4949
    50 class MetricAnIso{ public:
    51   friend class MatVVP2x2;
    52   Real8 a11,a21,a22;
    53   MetricAnIso(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}
    54   MetricAnIso(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){}
    55   MetricAnIso()  {}; //
    56   MetricAnIso(const Real8  a[3],const  MetricAnIso m0,
    57               const  MetricAnIso m1,const  MetricAnIso m2 );
    58   R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);}
    59   Real8 det() const {return a11*a22-a21*a21;} 
    60   R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
    61   R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
    62 //  D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
    63   MetricAnIso( Real8  a,const  MetricAnIso ma,
    64                Real8  b,const  MetricAnIso mb);
    65   int  IntersectWith(const MetricAnIso M);
    66   MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return  MetricAnIso(a11*c2,a21*c2,a22*c2);}
    67   MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return  MetricAnIso(a11*c2,a21*c2,a22*c2);}
    68   operator D2xD2(){ return D2xD2(a11,a21,a21,a22);}
     50        class MetricAnIso{ public:
     51                friend class MatVVP2x2;
     52                Real8 a11,a21,a22;
     53                MetricAnIso(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}
     54                MetricAnIso(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){}
     55                MetricAnIso()  {}; //
     56                MetricAnIso(const Real8  a[3],const  MetricAnIso m0,
     57                                        const  MetricAnIso m1,const  MetricAnIso m2 );
     58                R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);}
     59                Real8 det() const {return a11*a22-a21*a21;} 
     60                R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
     61                R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
     62                //  D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);}
     63                MetricAnIso( Real8  a,const  MetricAnIso ma,
     64                                        Real8  b,const  MetricAnIso mb);
     65                int  IntersectWith(const MetricAnIso M);
     66                MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return  MetricAnIso(a11*c2,a21*c2,a22*c2);}
     67                MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return  MetricAnIso(a11*c2,a21*c2,a22*c2);}
     68                operator D2xD2(){ return D2xD2(a11,a21,a21,a22);}
    6969
    70   Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
    71 //  Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
    72   Real8 operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};
    73   inline void  Box(Real4 &hx,Real4 &hy) const ; 
    74   MetricAnIso(const MatVVP2x2);
    75 };
     70                Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
     71                //  Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
     72                Real8 operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};
     73                inline void  Box(Real4 &hx,Real4 &hy) const ; 
     74                MetricAnIso(const MatVVP2x2);
     75        };
    7676
     77        class MatVVP2x2{
     78                friend  class MetricAnIso;
     79                friend  class MetricIso;
     80                public:
     81                double lambda1,lambda2;
     82                D2 v;
    7783
    78 class MatVVP2x2
    79 {
    80   friend  class MetricAnIso;
    81   friend  class MetricIso;
    82 public:
    83   double lambda1,lambda2;
    84   D2 v;
     84                MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
    8585
     86                void  Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);}
     87                void  pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);}
     88                void  Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;}
     89                void  Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;}
    8690
    87   MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){}
    88  
    89   void  Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);}
    90   void  pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);}
    91   void  Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;}
    92   void  Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;}
     91                void Minh(double h) {Max(1.0/(h*h));}
     92                void Maxh(double h) {Min(1.0/(h*h));}
     93                void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
     94                MatVVP2x2(const MetricAnIso );
     95                MatVVP2x2(const MetricIso M) :  lambda1(1/(M.h*M.h)),lambda2(1/(M.h*M.h)),v(1,0) {}
     96                Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));}
     97                Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));}
     98                Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);}
     99                Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);}
     100                Real8 Aniso2() const  { return lmax()/lmin();}
     101                inline void BoundAniso2(const Real8 coef);
     102                Real8 Aniso() const  { return sqrt( Aniso2());}
     103                void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
     104                void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
     105          };
    93106
    94   void Minh(double h) {Max(1.0/(h*h));}
    95   void Maxh(double h) {Min(1.0/(h*h));}
    96   void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);}
    97   MatVVP2x2(const MetricAnIso );
    98   MatVVP2x2(const MetricIso M) :  lambda1(1/(M.h*M.h)),lambda2(1/(M.h*M.h)),v(1,0) {}
    99   Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));}
    100   Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));}
    101   Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);}
    102   Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);}
    103   Real8 Aniso2() const  { return lmax()/lmin();}
    104   inline void BoundAniso2(const Real8 coef);
    105   Real8 Aniso() const  { return sqrt( Aniso2());}
    106   void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));}
    107   void operator *=(double coef){ lambda1*=coef;lambda2*=coef;}
    108 };
     107        inline void  MatVVP2x2::BoundAniso2(const Real8 coef){
     108                if (coef<=1.00000000001)
     109                 if (lambda1 < lambda2)
     110                  lambda1 = bamg::Max(lambda1,lambda2*coef);
     111                 else
     112                  lambda2 = bamg::Max(lambda2,lambda1*coef);
     113                else  // a verifier
     114                 if (lambda1 > lambda2)
     115                  lambda1 = bamg::Min(lambda1,lambda2*coef);
     116                 else
     117                  lambda2 = bamg::Min(lambda2,lambda1*coef);
     118          }
    109119
    110 inline void  MatVVP2x2::BoundAniso2(const Real8 coef)
    111 {
    112   if (coef<=1.00000000001)
    113     if (lambda1 < lambda2)
    114       lambda1 = bamg::Max(lambda1,lambda2*coef);
    115     else
    116       lambda2 = bamg::Max(lambda2,lambda1*coef);
    117   else  // a verifier
    118     if (lambda1 > lambda2)
    119       lambda1 = bamg::Min(lambda1,lambda2*coef);
    120     else
    121       lambda2 = bamg::Min(lambda2,lambda1*coef);
    122 }
     120        void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
     121        MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
     122
     123        inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) {
     124                //     recompose M in: M = V^t lambda V
     125                //     V = ( v,v^\perp)
     126                //  where v^\perp = (-v_1,v_0)
     127                double v00=M.v.x*M.v.x;
     128                double v11=M.v.y*M.v.y;
     129                double v01=M.v.x*M.v.y;
     130                a11=v00*M.lambda1+v11*M.lambda2;
     131                a21=v01*(M.lambda1-M.lambda2);
     132                a22=v00*M.lambda2+v11*M.lambda1;
     133          }
     134
     135        inline   void  MetricAnIso::Box(Real4 &hx,Real4 &hy) const {
     136                Real8 d=  a11*a22-a21*a21;
     137                hx = sqrt(a22/d);
     138                hy = sqrt(a11/d);
     139        }
     140
     141        inline MetricIso::MetricIso(const MetricAnIso M)
     142          {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
     143
     144        inline  MetricIso::MetricIso(Real8 a11,Real8 a21,Real8 a22)
     145          {MatVVP2x2 vp(MetricAnIso(a11,a21,a22));h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
    123146
    124147
    125148
     149        class SaveMetricInterpole {
     150                friend  Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
     151                friend Real8 abscisseInterpole(const MetricAnIso ,const  MetricAnIso , R2 ,Real8 ,int );
     152                int opt;
     153                Real8 lab;
     154                Real8 L[1024],S[1024];
     155        };
    126156
    127 void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
    128 MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
     157        extern SaveMetricInterpole  LastMetricInterpole; // for optimization
    129158
    130 inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) 
    131 {
    132  //     recompose M in: M = V^t lambda V
    133   //     V = ( v,v^\perp)
    134   //  where v^\perp = (-v_1,v_0)
    135   double v00=M.v.x*M.v.x;
    136   double v11=M.v.y*M.v.y;
    137   double v01=M.v.x*M.v.y;
    138   a11=v00*M.lambda1+v11*M.lambda2;
    139   a21=v01*(M.lambda1-M.lambda2);
    140   a22=v00*M.lambda2+v11*M.lambda1;
    141 }
     159        Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB);
     160        Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim=0);
    142161
     162        inline Real8 LengthInterpole(Real8 la,Real8 lb)
     163          {   return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ?  (la+lb)/2  : la*lb*log(la/lb)/(la-lb);}
    143164
    144 inline   void  MetricAnIso::Box(Real4 &hx,Real4 &hy) const {
    145   Real8 d=  a11*a22-a21*a21;
    146   hx = sqrt(a22/d);
    147   hy = sqrt(a11/d);
    148 }
    149 
    150 
    151 inline MetricIso::MetricIso(const MetricAnIso M)
    152   {MatVVP2x2 vp(M);h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
    153  
    154 inline  MetricIso::MetricIso(Real8 a11,Real8 a21,Real8 a22)
    155   {MatVVP2x2 vp(MetricAnIso(a11,a21,a22));h=1/sqrt(Max(vp.lambda1,vp.lambda2));}
    156 
    157 
    158 
    159 class SaveMetricInterpole {
    160   friend  Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
    161   friend Real8 abscisseInterpole(const MetricAnIso ,const  MetricAnIso , R2 ,Real8 ,int );
    162   int opt;
    163   Real8 lab;
    164   Real8 L[1024],S[1024];
    165 };
    166 
    167 extern SaveMetricInterpole  LastMetricInterpole; // for optimization
    168 
    169 
    170  Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB);
    171  Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim=0);
    172 
    173 
    174 
    175 inline Real8 LengthInterpole(Real8 la,Real8 lb)
    176 {   return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ?  (la+lb)/2  : la*lb*log(la/lb)/(la-lb);}
    177 
    178 inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s)
    179 { return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20))  ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb);}
     165        inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s)
     166          { return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20))  ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb);}
    180167
    181168}

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

    r2861 r2862  
    55namespace bamg {
    66
    7 const int MaxDeep = 30;
    8 typedef  long  IntQuad;
    9 const IntQuad MaxISize = ( 1L << MaxDeep);
     7        const int MaxDeep = 30;
     8        typedef  long  IntQuad;
     9        const IntQuad MaxISize = ( 1L << MaxDeep);
    1010
    1111
    12 class Triangles;
    13 class Vertex;
     12        class Triangles;
     13        class Vertex;
    1414
    15 class QuadTree {
    16  public:
     15        class QuadTree {
     16                public:
    1717
    18   class QuadTreeBox {
    19   public:
     18                        class QuadTreeBox {
     19                                public:
    2020
    21     long n; // if n < 4 => Vertex else =>  QuadTreeBox;
    22     union {
    23       QuadTreeBox *b[4];
    24       Vertex * v[4];
    25     };
    26    
    27 
    28   }; // end class QuadTreeBox  /////////////////
    29 
    30   class StorageQuadTreeBox {
    31   public:
    32     QuadTreeBox *b,*bc,*be;
    33     long len;
    34     StorageQuadTreeBox *n; // next StorageQuadTreeBox
    35     StorageQuadTreeBox(long ,StorageQuadTreeBox * =0);
    36     ~StorageQuadTreeBox() {
    37       if(n) delete n;
    38       delete [] b;
    39     }
    40     long  SizeOf() const {
    41       return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);
    42     }
    43   }; // end class  StorageQuadTreeBox
    44  
    45   StorageQuadTreeBox * sb;
    46  
    47  
    48   long  lenStorageQuadTreeBox;
    49 
    50 public:
    51   QuadTreeBox * root;
    52   Triangles *th;
    53   long NbQuadTreeBox,NbVertices;
    54   long NbQuadTreeBoxSearch,NbVerticesSearch;
    55   Vertex * NearestVertex(Icoor1 i,Icoor1 j);
    56   Vertex *  NearestVertexWithNormal(Icoor1 i,Icoor1 j); // new version 
    57   Vertex * ToClose(Vertex & ,Real8 ,Icoor1,Icoor1);
    58   long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();}
     21                                        long n; // if n < 4 => Vertex else =>  QuadTreeBox;
     22                                        union {
     23                                                QuadTreeBox *b[4];
     24                                                Vertex * v[4];
     25                                        };
    5926
    6027
    61   void  Add( Vertex & w);
     28                        }; // end class QuadTreeBox  /////////////////
    6229
    63   QuadTreeBox* NewQuadTreeBox(){
    64     if(! (sb->bc<sb->be))
    65           sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
     30                        class StorageQuadTreeBox {
     31                                public:
     32                                        QuadTreeBox *b,*bc,*be;
     33                                        long len;
     34                                        StorageQuadTreeBox *n; // next StorageQuadTreeBox
     35                                        StorageQuadTreeBox(long ,StorageQuadTreeBox * =0);
     36                                        ~StorageQuadTreeBox() {
     37                                                if(n) delete n;
     38                                                delete [] b;
     39                                        }
     40                                        long  SizeOf() const {
     41                                                return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);
     42                                        }
     43                        }; // end class  StorageQuadTreeBox
    6644
    67          if (!sb || (sb->bc->n != 0)){
    68                  throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));
    69          }
    70     NbQuadTreeBox++;
    71     return sb->bc++;
    72   }
    73   ~QuadTree();
    74   QuadTree(Triangles * t,long nbv=-1);
    75   QuadTree();
    76        
    77 };
     45                        StorageQuadTreeBox * sb;
     46
     47
     48                        long  lenStorageQuadTreeBox;
     49
     50                public:
     51                        QuadTreeBox * root;
     52                        Triangles *th;
     53                        long NbQuadTreeBox,NbVertices;
     54                        long NbQuadTreeBoxSearch,NbVerticesSearch;
     55                        Vertex * NearestVertex(Icoor1 i,Icoor1 j);
     56                        Vertex *  NearestVertexWithNormal(Icoor1 i,Icoor1 j); // new version 
     57                        Vertex * ToClose(Vertex & ,Real8 ,Icoor1,Icoor1);
     58                        long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();}
     59
     60
     61                        void  Add( Vertex & w);
     62
     63                        QuadTreeBox* NewQuadTreeBox(){
     64                                if(! (sb->bc<sb->be))
     65                                 sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
     66
     67                                if (!sb || (sb->bc->n != 0)){
     68                                        throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));
     69                                }
     70                                NbQuadTreeBox++;
     71                                return sb->bc++;
     72                        }
     73                        ~QuadTree();
     74                        QuadTree(Triangles * t,long nbv=-1);
     75                        QuadTree();
     76        };
    7877}

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

    r2861 r2862  
    22
    33namespace bamg {
    4 template <class R,class RR> class P2xP2;
     4        template <class R,class RR> class P2xP2;
    55
    6 template <class R,class RR>
    7 class P2 {
     6        template <class R,class RR>
     7          class P2 {
    88
    9 public: 
    10   R x,y;
    11   P2 () :x(0),y(0) {};
    12   P2 (R a,R b)  :x(a),y(b)  {}
    13   P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {}
    14   P2<R,RR>   operator+(const P2<R,RR> & cc) const {return P2<R,RR>(x+cc.x,y+cc.y);}
    15   P2<R,RR>   operator-(const P2<R,RR> & cc) const {return P2<R,RR>(x-cc.x,y-cc.y);}
    16   P2<R,RR>   operator-()  const{return P2<R,RR>(-x,-y);}
    17 //  RR   operator*(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} // produit scalaire
    18   RR   operator,(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} // produit scalaire
    19   P2<R,RR>   operator*(R  cc) const {return P2<R,RR>(x*cc,y*cc);}
    20  // P2<R,RR>   operator*(RR  cc) const {return P2<R,RR>((R)(x*cc),(R)(y*cc));}
    21   P2<R,RR>   operator/(R  cc) const {return P2<R,RR>(x/cc,y/cc);}
    22   P2<R,RR>  operator+=(const  P2<R,RR> & cc) {x += cc.x;y += cc.y;return *this;}
    23   P2<R,RR>  operator/=(const  R r) {x /= r;y /= r;return *this;}
    24   P2<R,RR>  operator*=(const  R r) {x *= r;y *= r;return *this;}
    25   P2<R,RR>  operator-=(const  P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;}
    26 //  P2<R,RR> Orthogonal(const   P2<R,RR> r) {return P2<R,RR>(-r.y,r.x);}
    27  };
     9                  public: 
     10                          R x,y;
     11                          P2 () :x(0),y(0) {};
     12                          P2 (R a,R b)  :x(a),y(b)  {}
     13                          P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {}
     14                          P2<R,RR>   operator+(const P2<R,RR> & cc) const {return P2<R,RR>(x+cc.x,y+cc.y);}
     15                          P2<R,RR>   operator-(const P2<R,RR> & cc) const {return P2<R,RR>(x-cc.x,y-cc.y);}
     16                          P2<R,RR>   operator-()  const{return P2<R,RR>(-x,-y);}
     17                          RR   operator,(const P2<R,RR> & cc) const {return  (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} // produit scalaire
     18                          P2<R,RR>   operator*(R  cc) const {return P2<R,RR>(x*cc,y*cc);}
     19                          P2<R,RR>   operator/(R  cc) const {return P2<R,RR>(x/cc,y/cc);}
     20                          P2<R,RR>  operator+=(const  P2<R,RR> & cc) {x += cc.x;y += cc.y;return *this;}
     21                          P2<R,RR>  operator/=(const  R r) {x /= r;y /= r;return *this;}
     22                          P2<R,RR>  operator*=(const  R r) {x *= r;y *= r;return *this;}
     23                          P2<R,RR>  operator-=(const  P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;}
     24          };
    2825
    29 template <class R,class RR>
    30 class P2xP2 { // x ligne 1 y ligne2
     26        template <class R,class RR>
     27          class P2xP2 { // x ligne 1 y ligne2
    3128
    32   friend std::ostream& operator <<(std::ostream& f, const P2xP2<R,RR> & c)
    33     { f << '[' << c.x << ',' << c.y << ']' <<std::flush ; return f; }
    34      
    35   friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc)
    36     {return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y);}
     29                  friend std::ostream& operator <<(std::ostream& f, const P2xP2<R,RR> & c)
     30                        { f << '[' << c.x << ',' << c.y << ']' <<std::flush ; return f; }
     31
     32                  friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc)
     33                        {return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y);}
    3734
    3835
    39  public:
    40   P2<R,RR> x,y;
    41   P2xP2 (): x(),y()  {}
    42   P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {}
    43   P2xP2 (P2<R,RR> a,P2<R,RR> b,P2<R,RR> c ): x(b-a),y(c-a) {}
    44   P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {}
    45   P2<R,RR> operator*(const P2<R,RR> c) const {return P2<R,RR>(x.x*c.x + x.y*c.y, y.x*c.x + y.y*c.y);}
    46   P2xP2<R,RR>  operator*(P2xP2<R,RR> c) const
    47     { return  P2xP2<R,RR>(x.x*c.x.x + x.y*c.y.x,
    48                           x.x*c.x.y + x.y*c.y.y,
    49                           y.x*c.x.x + y.y*c.y.x,
    50                           y.x*c.x.y + y.y*c.y.y);}
    51   RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;}
    52   P2xP2<R,RR> inv()  const
    53     { RR d = (*this).det();
    54        return P2xP2<R,RR>((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) );
    55     };
    56    P2xP2<R,RR> t() {return P2xP2<R,RR>(x.x,y.x,x.y,y.y);} //transposer
    57    P2<R,RR>tx() {return P2<R,RR>(x.x,y.x);}
    58    P2<R,RR>ty() {return P2<R,RR>(x.y,y.y);}
     36                  public:
     37                  P2<R,RR> x,y;
     38                  P2xP2 (): x(),y()  {}
     39                  P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {}
     40                  P2xP2 (P2<R,RR> a,P2<R,RR> b,P2<R,RR> c ): x(b-a),y(c-a) {}
     41                  P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {}
     42                  P2<R,RR> operator*(const P2<R,RR> c) const {return P2<R,RR>(x.x*c.x + x.y*c.y, y.x*c.x + y.y*c.y);}
     43                  P2xP2<R,RR>  operator*(P2xP2<R,RR> c) const
     44                        { return  P2xP2<R,RR>(x.x*c.x.x + x.y*c.y.x,
     45                                                x.x*c.x.y + x.y*c.y.y,
     46                                                y.x*c.x.x + y.y*c.y.x,
     47                                                y.x*c.x.y + y.y*c.y.y);}
     48                  RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;}
     49                  P2xP2<R,RR> inv()  const
     50                        { RR d = (*this).det();
     51                          return P2xP2<R,RR>((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) );
     52                        };
     53                  P2xP2<R,RR> t() {return P2xP2<R,RR>(x.x,y.x,x.y,y.y);} //transposer
     54                  P2<R,RR>tx() {return P2<R,RR>(x.x,y.x);}
     55                  P2<R,RR>ty() {return P2<R,RR>(x.y,y.y);}
    5956
    60 }; 
     57          }; 
    6158
    62 template  <class R,class RR> 
    63 inline RR Det(const P2<R,RR> x,const P2<R,RR> y) {
    64   return (RR) x.x * (RR) y.y - (RR) x.y * (RR) y.x ;}
     59        template  <class R,class RR> 
     60          inline RR Det(const P2<R,RR> x,const P2<R,RR> y) {
     61                  return (RR) x.x * (RR) y.y - (RR) x.y * (RR) y.x ;}
    6562
    66 template  <class R,class RR> 
    67 inline RR Area2 (const P2<R,RR> a,const P2<R,RR> b,const P2<R,RR> c) {
    68   return Det(b-a,c-a) ;}
     63        template  <class R,class RR> 
     64          inline RR Area2 (const P2<R,RR> a,const P2<R,RR> b,const P2<R,RR> c) {
     65                  return Det(b-a,c-a) ;}
    6966
    70 template  <class R,class RR> 
    71 inline R Norme1 (const P2<R,RR> x) {
    72   return (Abs(x.x)+Abs(x.y)) ;}
     67        template  <class R,class RR> 
     68          inline R Norme1 (const P2<R,RR> x) {
     69                  return (Abs(x.x)+Abs(x.y)) ;}
    7370
    74 template  <class R,class RR> 
    75 inline R NormeInfini (const P2<R,RR> x) {
    76   return Max(Abs(x.x),Abs(x.y)) ;}
     71        template  <class R,class RR> 
     72          inline R NormeInfini (const P2<R,RR> x) {
     73                  return Max(Abs(x.x),Abs(x.y)) ;}
    7774
    78 template  <class R,class RR> 
    79 inline RR Norme2_2 (const P2<R,RR> x) {
    80   return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ;}
     75        template  <class R,class RR> 
     76          inline RR Norme2_2 (const P2<R,RR> x) {
     77                  return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ;}
    8178
    82 template  <class R,class RR> 
    83 inline RR Norme2 (const P2<R,RR> x) {
    84   return sqrt((RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y) ;}
     79        template  <class R,class RR> 
     80          inline RR Norme2 (const P2<R,RR> x) {
     81                  return sqrt((RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y) ;}
    8582
    86 template  <class R,class RR> 
    87 inline P2<R,RR> Orthogonal (const P2<R,RR> x) {
    88   return  P2<R,RR>(-x.y,x.x);}
     83        template  <class R,class RR> 
     84          inline P2<R,RR> Orthogonal (const P2<R,RR> x) {
     85                  return  P2<R,RR>(-x.y,x.x);}
    8986
    90 template <class R,class RR>
    91 inline  std::ostream& operator <<(std::ostream& f, const P2<R,RR> & c)
    92   { f << '[' << c.x << ',' << c.y <<']' <<std::flush ; return f; }
     87        template <class R,class RR>
     88          inline  std::ostream& operator <<(std::ostream& f, const P2<R,RR> & c)
     89                { f << '[' << c.x << ',' << c.y <<']' <<std::flush ; return f; }
    9390}

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

    r2856 r2862  
    44namespace bamg {
    55
    6 class SetOfEdges4 ;
    7 class Int4Edge{
    8 friend class SetOfEdges4;
    9 public:
    10   Int4 i,j;
    11   Int4 next;
    12 };
     6        class SetOfEdges4 ;
     7        class Int4Edge{
     8                friend class SetOfEdges4;
     9                public:
     10                Int4 i,j;
     11                Int4 next;
     12        };
    1313
    14 class SetOfEdges4 {
    15   Int4 nx,nbax,NbOfEdges;
    16   Int4 * tete;
    17   Int4Edge * Edges;
     14        class SetOfEdges4 {
     15                Int4 nx,nbax,NbOfEdges;
     16                Int4 * tete;
     17                Int4Edge * Edges;
    1818
    19 public:
    20   SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet
    21   ~SetOfEdges4() {
    22   delete [] tete; delete [] Edges;}
    23    Int4 add (Int4 ii,Int4 jj);
    24   Int4 addtrie (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
    25   Int4  nb(){return NbOfEdges;}
    26   Int4 find (Int4 ii,Int4 jj);
    27   Int4 findtrie (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj)  : find (jj,ii) ;}
    28   Int4 i(Int4 k){return Edges[k].i;}
    29   Int4 j(Int4 k){return Edges[k].j;}
    30   Int4 newarete(Int4 k){return NbOfEdges == k+1;}
    31   Int4Edge & operator[](Int4 k){return  Edges[k];}
    32 };
     19                public:
     20                SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet
     21                ~SetOfEdges4() {delete [] tete; delete [] Edges;}
     22                Int4 add (Int4 ii,Int4 jj);
     23                Int4 addtrie (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
     24                Int4  nb(){return NbOfEdges;}
     25                Int4 find (Int4 ii,Int4 jj);
     26                Int4 findtrie (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj)  : find (jj,ii) ;}
     27                Int4 i(Int4 k){return Edges[k].i;}
     28                Int4 j(Int4 k){return Edges[k].j;}
     29                Int4 newarete(Int4 k){return NbOfEdges == k+1;}
     30                Int4Edge & operator[](Int4 k){return  Edges[k];}
     31        };
    3332}
    3433
Note: See TracChangeset for help on using the changeset viewer.