![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2984
- Timestamp:
- 02/08/10 13:33:59 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2983 r2984 32 32 //some functions 33 33 inline int BinaryRand(){ 34 34 #ifdef RAND_MAX 35 35 const long HalfRandMax = RAND_MAX/2; 36 36 return rand() < HalfRandMax; 37 37 #else 38 38 return rand() & 16384; //2^14 (for sun because RAND_MAX is not def in stdlib.h) 39 39 #endif 40 40 } 41 41 … … 45 45 typedef P2<Real8,Real8> R2; 46 46 typedef P2<Real4,Real8> R2xR2; 47 47 48 } 48 49 49 #include "Metric.h" 50 50 -
issm/trunk/src/c/Bamgx/Metric.h
r2981 r2984 6 6 namespace bamg { 7 7 8 typedef P2<double,double> D2;8 typedef P2<double,double> D2; 9 9 typedef P2xP2<double,double> D2xD2; 10 10 … … 13 13 14 14 15 class Metric{ public: 16 friend class MatVVP2x2; 17 Real8 a11,a21,a22; 18 Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} 19 Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){} 20 Metric() {}; // 21 Metric(const Real8 a[3],const Metric m0, 22 const Metric m1,const Metric m2 ); 23 R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} 24 Real8 det() const {return a11*a22-a21*a21;} 25 R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 26 R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 27 // D2 Orthogonal(const D2 x){return D2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 28 Metric( Real8 a,const Metric ma, 29 Real8 b,const Metric mb); 30 int IntersectWith(const Metric M2); 31 Metric operator*(Real8 c) const {Real8 c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} 32 Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} 33 operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} 15 class Metric{ 34 16 35 Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 36 // Real8 operator()(D2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 37 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;}; 38 inline void Box(Real4 &hx,Real4 &hy) const ; 39 Metric(const MatVVP2x2); 40 void Echo(); 17 public: 18 //fields 19 Real8 a11,a21,a22; 20 //friends 21 friend class MatVVP2x2; 22 //functions 23 Metric(){}; 24 Metric(const MatVVP2x2); 25 Metric(Real8 a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} 26 Metric(Real8 a,Real8 b,Real8 c) :a11(a),a21(b),a22(c){} 27 Metric( Real8 a,const Metric ma, Real8 b,const Metric mb); 28 Metric(const Real8 a[3],const Metric m0,const Metric m1,const Metric m2 ); 29 void Echo(); 30 R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} 31 Real8 det() const {return a11*a22-a21*a21;} 32 R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 33 R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 34 int IntersectWith(const Metric M2); 35 inline void Box(Real4 &hx,Real4 &hy) const ; 36 //operators 37 Metric operator*(Real8 c) const {Real8 c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} 38 Metric operator/(Real8 c) const {Real8 c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} 39 operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} 40 Real8 operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 41 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;}; 42 41 43 }; 42 44 43 45 class MatVVP2x2{ 44 friend class Metric;45 46 public: 46 double lambda1,lambda2; 47 D2 v; 48 49 void Echo(); 50 MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){} 51 52 void Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);} 53 void pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);} 54 void Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;} 55 void Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;} 56 57 void Minh(double h) {Min(1.0/(h*h));} 58 void Maxh(double h) {Max(1.0/(h*h));} 59 void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} 60 MatVVP2x2(const Metric ); 61 Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} 62 Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} 63 Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} 64 Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} 65 Real8 Aniso2() const { return lmax()/lmin();} 66 inline void BoundAniso2(const Real8 coef); 67 Real8 Aniso() const { return sqrt( Aniso2());} 68 void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));} 69 void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} 47 //fields 48 double lambda1,lambda2; 49 D2 v; 50 //friends 51 friend class Metric; 52 //functions 53 MatVVP2x2(const Metric ); 54 MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){} 55 void Echo(); 56 void Abs(){lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);} 57 void pow(double p){lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);} 58 void Min(double a) { lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;} 59 void Max(double a) { lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;} 60 void Minh(double h) {Min(1.0/(h*h));} 61 void Maxh(double h) {Max(1.0/(h*h));} 62 void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} 63 Real8 hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} 64 Real8 hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} 65 Real8 lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} 66 Real8 lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} 67 Real8 Aniso2() const { return lmax()/lmin();} 68 Real8 Aniso() const { return sqrt( Aniso2());} 69 void BoundAniso(const Real8 c){ BoundAniso2(1/(c*c));} 70 inline void BoundAniso2(const Real8 coef); 71 //operators 72 void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} 70 73 }; 71 74 75 class SaveMetricInterpole { 76 friend Real8 LengthInterpole(const Metric ,const Metric , R2 ); 77 friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8 ,int ); 78 int opt; 79 Real8 lab; 80 Real8 L[1024],S[1024]; 81 }; 82 83 extern SaveMetricInterpole LastMetricInterpole; // for optimization 84 //Functions 85 void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V); 86 Real8 LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); 87 Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8 s,int optim=0); 88 89 //inlines 72 90 inline void MatVVP2x2::BoundAniso2(const Real8 coef){ 73 91 if (coef<=1.00000000001){ … … 84 102 } 85 103 } 86 87 void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V);88 89 104 inline Metric::Metric(const MatVVP2x2 M) { 90 // recompose M in: M = V^t lambda V91 // V = ( v,v^\perp)92 // where v^\perp = (-v_1,v_0)93 105 double v00=M.v.x*M.v.x; 94 106 double v11=M.v.y*M.v.y; … … 97 109 a21=v01*(M.lambda1-M.lambda2); 98 110 a22=v00*M.lambda2+v11*M.lambda1; 99 } 100 111 } 101 112 inline void Metric::Box(Real4 &hx,Real4 &hy) const { 102 113 Real8 d= a11*a22-a21*a21; … … 104 115 hy = sqrt(a11/d); 105 116 } 106 107 class SaveMetricInterpole { 108 friend Real8 LengthInterpole(const Metric ,const Metric , R2 ); 109 friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8 ,int ); 110 int opt; 111 Real8 lab; 112 Real8 L[1024],S[1024]; 113 }; 114 115 extern SaveMetricInterpole LastMetricInterpole; // for optimization 116 117 Real8 LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); 118 Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8 s,int optim=0); 119 120 inline Real8 LengthInterpole(Real8 la,Real8 lb) 121 { return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb);} 122 123 inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s) 124 { 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);} 117 inline Real8 LengthInterpole(Real8 la,Real8 lb) { 118 return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb); 119 } 120 inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8 s){ 121 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); 122 } 125 123 126 124 } -
issm/trunk/src/c/Bamgx/QuadTree.h
r2981 r2984 5 5 namespace bamg { 6 6 7 typedef long IntQuad; 8 7 9 const int MaxDeep = 30; 8 typedef long IntQuad;9 //long int 1, bitwise operation: 8L = 00001000 << 2L -> 00100000 shifted left by 210 10 const IntQuad MaxISize = ( 1L << MaxDeep); 11 11 12 class Triangles; 12 13 class Vertex; 13 class QuadTree { 14 15 class QuadTree{ 14 16 private: 15 17 class QuadTreeBox { … … 39 41 40 42 public: 41 42 43 //fields 43 44 QuadTreeBox* root; 44 45 Triangles* th; 45 46 46 //functions 47 47 ~QuadTree(); -
issm/trunk/src/c/Bamgx/R2.h
r2981 r2984 4 4 template <class R,class RR> 5 5 class P2 { 6 7 6 public: 7 //fields 8 8 R x,y; 9 //functions 10 P2 () :x(0),y(0) {}; 11 P2 (R a,R b) :x(a),y(b) {} 12 P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {} 9 13 void Echo(){ 10 14 printf("Member of P2:\n"); … … 12 16 printf(" y: %g\n",y); 13 17 } 14 P2 () :x(0),y(0) {}; 15 P2 (R a,R b) :x(a),y(b) {} 16 P2 (P2 A,P2 B) : x(B.x-A.x),y(B.y-A.y) {} 18 //operators 17 19 RR operator,(const P2<R,RR> & cc) const {return (RR) x* (RR) cc.x+(RR) y* (RR) cc.y;} //scalar product 18 20 P2<R,RR> operator+(const P2<R,RR> & cc) const {return P2<R,RR>(x+cc.x,y+cc.y);} … … 29 31 template <class R,class RR> 30 32 class P2xP2 { 31 friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc){ 32 return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y); 33 } 33 private: 34 friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc){ 35 return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y); 36 } 34 37 public: 35 P2<R,RR> x,y; 36 37 //functions 38 void Echo(){ 39 printf("Member of P2xP2:\n"); 40 printf(" x.x: %g x.y: %g\n",x.x,x.y); 41 printf(" y.x: %g y.x: %g\n",y.x,y.y); 42 } 43 P2xP2 (): x(),y() {} 44 P2xP2 (P2<R,RR> a,P2<R,RR> b): x(a),y(b) {} 45 P2xP2 (P2<R,RR> a,P2<R,RR> b,P2<R,RR> c ): x(b-a),y(c-a) {} 46 P2xP2 (R xx,R xy,R yx,R yy) :x(xx,xy),y(yx,yy) {} 47 RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;} 48 P2xP2<R,RR> inv() const{ 49 RR d = (*this).det(); 50 return P2xP2<R,RR>((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) ); 51 }; 52 P2xP2<R,RR> t() {return P2xP2<R,RR>(x.x,y.x,x.y,y.y);} //transposer 53 P2<R,RR> tx() {return P2<R,RR>(x.x,y.x);} 54 P2<R,RR> ty() {return P2<R,RR>(x.y,y.y);} 55 56 //Operators 57 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);} 58 P2xP2<R,RR> operator*(P2xP2<R,RR> c) const{ 59 return P2xP2<R,RR>(x.x*c.x.x + x.y*c.y.x, 60 x.x*c.x.y + x.y*c.y.y, 61 y.x*c.x.x + y.y*c.y.x, 62 y.x*c.x.y + y.y*c.y.y); 63 } 38 //fields 39 P2<R,RR> x,y; 40 //functions 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 void Echo(){ 46 printf("Member of P2xP2:\n"); 47 printf(" x.x: %g x.y: %g\n",x.x,x.y); 48 printf(" y.x: %g y.x: %g\n",y.x,y.y); 49 } 50 RR det() const {return (RR) x.x* (RR) y.y - (RR) x.y * (RR) y.x;} 51 P2xP2<R,RR> inv() const{ 52 RR d = (*this).det(); 53 return P2xP2<R,RR>((R)( y.y /d) ,(R)(-x.y/d),(R)( -y.x/d) ,(R)( x.x/d) ); 54 }; 55 P2xP2<R,RR> t() {return P2xP2<R,RR>(x.x,y.x,x.y,y.y);} //transposer 56 P2<R,RR> tx() {return P2<R,RR>(x.x,y.x);} 57 P2<R,RR> ty() {return P2<R,RR>(x.y,y.y);} 58 //Operators 59 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);} 60 P2xP2<R,RR> operator*(P2xP2<R,RR> c) const{ 61 return P2xP2<R,RR>(x.x*c.x.x + x.y*c.y.x, 62 x.x*c.x.y + x.y*c.y.y, 63 y.x*c.x.x + y.y*c.y.x, 64 y.x*c.x.y + y.y*c.y.y); 65 } 64 66 }; 65 67 … … 69 71 return (RR) x.x * (RR) y.y - (RR) x.y * (RR) y.x ; 70 72 } 71 72 73 template <class R,class RR> 73 74 inline RR Area2 (const P2<R,RR> a,const P2<R,RR> b,const P2<R,RR> c) { 74 75 return Det(b-a,c-a) ; 75 76 } 76 77 77 template <class R,class RR> 78 78 inline R Norme1 (const P2<R,RR> x) { 79 79 return (Abs(x.x)+Abs(x.y)) ; 80 80 } 81 82 81 template <class R,class RR> 83 82 inline R NormeInfini (const P2<R,RR> x) { 84 83 return Max(Abs(x.x),Abs(x.y)) ; 85 84 } 86 87 85 template <class R,class RR> 88 86 inline RR Norme2_2 (const P2<R,RR> x) { 89 87 return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ; 90 88 } 91 92 89 template <class R,class RR> 93 90 inline RR Norme2 (const P2<R,RR> x) { 94 91 return sqrt((RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y) ; 95 92 } 96 97 93 template <class R,class RR> 98 94 inline P2<R,RR> Orthogonal (const P2<R,RR> x) { -
issm/trunk/src/c/Bamgx/SetOfE4.h
r2981 r2984 14 14 15 15 class SetOfEdges4 { 16 Int4 nx,nbax,NbOfEdges; 17 Int4* head; 18 Int4Edge* Edges; 16 17 private: 18 Int4 nx,nbax,NbOfEdges; 19 Int4* head; 20 Int4Edge* Edges; 19 21 20 22 public: 21 SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet22 ~SetOfEdges4() {delete [] head; delete [] Edges;}23 Int4 add (Int4 ii,Int4 jj);24 Int4 SortAndAdd (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 SortAndFind (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;}23 SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 24 ~SetOfEdges4() {delete [] head; delete [] Edges;} 25 Int4 add (Int4 ii,Int4 jj); 26 Int4 SortAndAdd (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;} 27 Int4 nb(){return NbOfEdges;} 28 Int4 find (Int4 ii,Int4 jj); 29 Int4 SortAndFind (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} 30 Int4 i(Int4 k){return Edges[k].i;} 31 Int4 j(Int4 k){return Edges[k].j;} 32 Int4 newarete(Int4 k){return NbOfEdges == k+1;} 31 33 32 //operators33 Int4Edge & operator[](Int4 k){return Edges[k];}34 //operators 35 Int4Edge & operator[](Int4 k){return Edges[k];} 34 36 }; 35 37 } -
issm/trunk/src/c/Bamgx/meshtype.h
r2885 r2984 3 3 #include <limits.h> 4 4 namespace bamg { 5 template<class T> inline T Square (const T &a) { return a*a;}6 template<class T> inline T Min (const T &a,const T &b){return a < b ? a : b;}7 template<class T> inline T Max (const T &a,const T & b){return a > b ? a : b;}8 template<class T> inline T Abs (const T &a){return a <0 ? -a : a;}9 template<class T> inline double Norme (const T &a){return sqrt(a*a);}10 template<class T> inline void Exchange (T& a,T& b) {T c=a;a=b;b=c;}11 // for pb on microsoft compiler12 template<class T> inline T Max3 (const T &a,const T & b,const T & c){return Max(Max(a,b),c);}13 template<class T> inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);}14 5 15 typedef float Real4; 16 typedef double Real8; 17 typedef short Int1; 18 typedef short Int2; 19 typedef long Int4; 6 //template functions 7 template<class T> inline T Square (const T &a) { return a*a;} 8 template<class T> inline T Min (const T &a,const T &b){return a < b ? a : b;} 9 template<class T> inline T Max (const T &a,const T & b){return a > b ? a : b;} 10 template<class T> inline T Abs (const T &a){return a <0 ? -a : a;} 11 template<class T> inline double Norme (const T &a){return sqrt(a*a);} 12 template<class T> inline void Exchange (T& a,T& b) {T c=a;a=b;b=c;} 13 14 // for pb on microsoft compiler 15 template<class T> inline T Max3 (const T &a,const T & b,const T & c){return Max(Max(a,b),c);} 16 template<class T> inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);} 17 18 typedef float Real4; 19 typedef double Real8; 20 typedef short Int1; 21 typedef short Int2; 22 typedef long Int4; 20 23 21 24 #if LONG_BIT > 63 22 // for alpha and silicon23 24 25 26 25 // for alpha and silicon 26 typedef int Icoor1; 27 typedef long Icoor2; 28 const Icoor1 MaxICoor = 1073741823; // 2^30-1 29 const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; 27 30 #else 28 29 30 31 31 typedef int Icoor1; 32 typedef double Icoor2; 33 const Icoor1 MaxICoor = 8388608; // 2^23 34 const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; 32 35 #endif 33 36 class Triangles; 34 37 } 35 38 #endif
Note:
See TracChangeset
for help on using the changeset viewer.
